MATLAB Answers

eigs gives wrong eigenvalues

46 views (last 30 days)
Chiao-Ting Li
Chiao-Ting Li on 17 Jan 2020
Commented: Chiao-Ting Li on 3 Feb 2020 at 6:42
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');
eig_vs_eigs.png

  2 Comments

Steven Lord
Steven Lord on 17 Jan 2020
Which release of MATLAB are you using?
Chiao-Ting Li
Chiao-Ting Li on 20 Jan 2020 at 0:07
I have MATLAB R2019b.

Sign in to comment.

Answers (1)

Christine Tobler
Christine Tobler on 21 Jan 2020 at 12:16
Edited: Christine Tobler on 23 Jan 2020 at 18:43
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.

  2 Comments

Christine Tobler
Christine Tobler on 23 Jan 2020 at 18:42
After looking at this a bit more, I was going in the wrong direction in my answer above.
The problem here is that K is badly conditioned. Internally, a linear system with K has to be solved, and this is causing large round-off error in that computation. On a higher level, 'sm' is only safe to use if none of the eigenvalues of (K, M) are zero or extremely close to zero.
To address this problem, replace
[V, D, flag] = eigs(K, M, 523, 'sm');
with
[V, D, flag] = eigs(K, M, 523, 1e-2);
which instead of computing the eigenvalues of smallest magnitude, is computes those closest to 1e-2. The effect of this is that the linear system to be solved is K - 1e-2*M, which is a well conditioned matrix, because this is not within round-off (~ 1e-15) of an eigenvalue of (K, M).
Usually, eigs would give a warning in case the input matrix is badly conditioned, but for this particular matrix, it estimated the condition number as just too small to cause that warning to be given.
Chiao-Ting Li
Chiao-Ting Li on 3 Feb 2020 at 6:42
Hi Christine,
Thank you for looking into my problems. Your solution to address the ill-conditioning in my K matrix works!
However, I wish to follow up a few more things:
1) I computed the condition number of K and K - 1e-1*M, both of which have huge orders of magnitude and are larger then 1/eps. How do you tell that K is ill-conditioned, but K-1e-2*M is well-conditioned?
condest(K) % 6.267333282663795e+18
condest(K - 1e-1*M) % 1.277502807524298e+16
eps % Roundoff error: 2.220446049250313e-16
1/eps % Reciprocal of eps: 4.503599627370496e+15
2) I figure that my M & K matrices from ANSYS will mostly be ill-conditioned, because they are both stored in as space matrices and have 2~5% entries being non-zero and the non-zero entries often have the order of magnitudes ~1e5-~1e8, although M is relatively better conditioned than K. Can I always use your tweak to scale M by 1e-2 to re-condition “K-1e-2*M”? It is magical that the scaling of 1e-2 does not change the eigenvalues of the system, but why? Also, how do you figure out 1e-2 is a good scaling at the very beginning?
Thank you very much for your help!
Chiao-Ting

Sign in to comment.

Sign in to answer this question.

Products


Release

R2019b