eigs gives wrong eigenvalues

20 views (last 30 days)
Chiao-Ting Li
Chiao-Ting Li on 17 Jan 2020
Commented: Chiao-Ting Li on 7 Jun 2020
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
I have MATLAB R2019b.

Sign in to comment.

Answers (1)

Christine Tobler
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
Andrew Knyazev on 5 Jun 2020
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...
Chiao-Ting Li
Chiao-Ting Li on 7 Jun 2020
Thanks, Antrew, for suggesting the fileexchange post on locally optimal block pcg. I will look into it and update results when available to keep people with similar problems informed.
Two more thoughts to add here:
1) I totally agree that it is possible that my toy example in ANSYS is spitting out ill-conditioned matrices because I created the toy example with overly simplified boundary conditions for debugging purposes. Due to the simplicity in my BCs setups, the toy example has rigid-body motions, i.e. eigen-frequencies at 0 Hz, which is fundamentally hard to sovle numerically. With more realistic structure without rigid-body motions, I may obtain more well-conditioned matrices.
2) On a different note, I have worked out two ways around the problem at hand for now. The first is with Christine's suggestion to scale K and M matrices. The second is to make ANSYS splitting out only the first several eigen-frequencies and eigen-values, and in addition also something called the residual vectors, then Ifollow the method of calculating the flexibility matrix (the inverse of K matrix) in traditional FE literature to solve for modal displacement etc. The calculation with the residual vector allows me to avoid using MATLAB to calculate eigen-frequencies completly, but still preserve the low-frequency modes. Just FYI.
Regards,
Chiao-Ting

Sign in to comment.

Categories

Find more on Simulink in Help Center and File Exchange

Products


Release

R2019b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!