46 views (last 30 days)

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');

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.

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.

Sign in to comment.

Sign in to answer this question.

Opportunities for recent engineering grads.

Apply Today
## 2 Comments

## Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/500699-eigs-gives-wrong-eigenvalues#comment_787076

⋮## Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/500699-eigs-gives-wrong-eigenvalues#comment_787076

## Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/500699-eigs-gives-wrong-eigenvalues#comment_787721

⋮## Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/500699-eigs-gives-wrong-eigenvalues#comment_787721

Sign in to comment.