eigs gives wrong eigenvalues
20 views (last 30 days)
Show older comments
Hi,
I am using MATLAB R2019b to compute eigenvalues of sparse matrices and discovered that eigs gives wrong results if I only ask for the first several smallest eigenvlues.
The sparse matrices come from the mass and stiff matrices of an ANSYS model I build for modal analysis. Here I provide the toy example of a mass connected to the ground with a translational spring, and I know for sure that this toy example will have eigenvalues with real parts and NO imaginary parts.
Here is my question:
If I compute all the 1047 eigenvalues in my toy example by eig and eigs, I can get matching and correct results. However, if I use eigs to compute only the first several smallest eigenvalues, it gives wrong results with imaginary parts! To be specific, if I ask eigs to give out eigenvalues for the first 523 smallest eigenvalues (the length being equal or shorter than half of 1047), the results are WRONG. eigs only computes correctly when I ask for more than half of the available eigenvalues. I have read many posts on this forum about why eigs is different than eig, due to the former uses the random start, but this does not solve my problem.
My real ANSYS model is way more complex than this toy example, and I cannot afford to use eig to find all eigenvalues, instead, I will use eigs to find the first 20~50 smallest eigenvalues. What can I do with eigs being wrong when I ask it to give eigenvalues shorter than half of the total available eigenvalues?
I have included my data and script blow.
Thank you,
Chiao-Ting
clear; close all; clc;
K = load('stiff_mtx_NoMass.txt'); % Symmetric sparse matrix in MMF format
K = spconvert(K(2:end,:));
K = K + K' - speye(size(K,1)).*diag(diag(K));
M = load('mass_mtx_NoMass.txt');
M = spconvert(M(2:end,:));
M = M + M' - speye(size(M,1)).*diag(diag(M));
[V, D, flag] = eigs(K, M, 523, 'sm'); % NG
% [V, D, flag] = eigs(K, M, 524, 'sm'); % OK
% [V, D, flag] = eigs(K, M, 1047, 'sm'); % OK
D = diag(D);
vecWn = sqrt(D)/2/pi;
K2 = full(K);
M2 = full(M);
[V2, D2] = eig(K2, M2);
D2 = diag(D2);
D2 = sort(D2);
vecWn2 = sqrt(D2)/2/pi;
figure(1); clf;
plot(real(vecWn2), 's-', 'markersize', 4); hold on;
% plot(imag(vecWn2), '+-', 'markersize', 4);
plot(real(vecWn), 'x-', 'markersize', 4);
% plot(imag(vecWn), 'o-', 'markersize', 4);
grid on;
xlabel('DOF');
ylabel('Eigenfrequency (Hz)');
legend('eig; Re(\omega)', 'eigs; Re(\omega)');
set(legend, 'location', 'northwest');
2 Comments
Answers (1)
Christine Tobler
on 21 Jan 2020
Edited: Christine Tobler
on 23 Jan 2020
Edit: Please see the comment below, the first answer I gave here was going in the wrong direction.
Thank you for the detailed description of the problem. It looks like something is going wrong in the symmetric branch of eigs. If I make a small modification to K, so that it's not exactly symmetric anymore, the residual error is as small as the one returned by eig:
clear; close all; clc;
K = load('stiff_mtx_NoMass.txt'); % Symmetric sparse matrix in MMF format
K = spconvert(K(2:end,:));
K = K + K' - speye(size(K,1)).*diag(diag(K));
M = load('mass_mtx_NoMass.txt');
M = spconvert(M(2:end,:));
M = M + M' - speye(size(M,1)).*diag(diag(M));
% Rescale the matrices (this will only change the scale of the eigenvalues,
% making it easier to tell if the residuals are correct)
K = K / norm(full(K));
M = M / norm(full(M));
[V2, D2] = eig(full(K), full(M));
resEig = norm(K*V2 - M*V2*D2) % 4.3904e-15
[V, D, flag] = eigs(K, M, 150, 'smallestabs');
resEigs = norm(K*V - M*V*D) % 0.3924
Kmod = K;
Kmod(1, 2) = Kmod(1, 2) + eps(Kmod(1, 2));
[V, D, flag] = eigs(Kmod, M, 150, 'smallestabs');
resEigsMod = norm(K*V - M*V*D) % 9.2051e-15
If I plot the eigenvalues, these also look the same between eig and eigs with the modified input K.
As an alternative to modifying K slightly, you can also pass in a function handle and specify that this should not be treated as symmetric in its name-value pair:
dK = decomposition(K, 'CheckCondition', false);
[V, D, flag] = eigs(@(x) dK\x, 315, M, 150, 'smallestabs', 'IsFunctionSymmetric', false);
resEigsDK = norm(K*V - M*V*D) % 1.4698e-14
I will look into what's happening inside of eigs in the symmetric case. For now, this workaround seems to fix the issue.
4 Comments
Andrew Knyazev
on 5 Jun 2020
You may want to try https://www.mathworks.com/matlabcentral/fileexchange/48-locally-optimal-block-preconditioned-conjugate-gradient
But generally speaking, no method can be expested to work well in double precision for such ill conditioned matrtices, unless you can find some way to scale them. If not, it is basically a bug in ANSYS of in your setup that it generates matrices unsuitable for numerical solution...
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!