Nonorthogonal eigenvectors for general eigenvalue problem with eig() and eigs()
29 views (last 30 days)
Show older comments
William
on 5 Jun 2024
Answered: Christine Tobler
on 5 Jun 2024
I am working with a finite element model and want to decompose my system using the eigenvectors of the system. The eigenvalue problem is of the form: K*V = M*V*D, where V is the eigenvector matrix and D is the eigenvalue matrix. Both M and K are block diagonal symmetric matrices. For the generalized eigenvalue problem it should be the case that V'*M*V = a matrix whose only nonzero terms are along the diagonal, but I am getting nonzero terms on either side of these. Using the code included below for something like N =5, I get (after rounding to nearest 3rd decimal to remove near zero terms)
abs( Orthgonal_Check_eig) = [ 78.8050 0 0 0 0
0 43.4540 0 0 0
0 0 0 34.9310 0
0 0 34.9310 0 0
0 0 0 0 29.6430]
Where I am taking the absolute value because this case has complex eigenvectors (though I believe I've seen it with real eigenvector matrices as well). If I rerun the code it does not always happen, but I am trying to understand why this can occur for both eig() and eigs() so that I can determine if it is the conditioning of the matrices or if it is inherent to the approximations inside of eig() and eigs(). Is there a better one of the two to use to avoid this issue?
N = 5;
M = zeros(N,N);
K = zeros(N,N);
for ii = 1:N
% Diagonal Terms
M(ii,ii) = rand*50;
K(ii,ii) = rand*500;
if ii ~=N
% Off Diagonal Terms
off_diagonal = rand*50;
M(ii,ii+1) = off_diagonal;
K(ii,ii+1) = off_diagonal*10;
M(ii+1,ii) = off_diagonal;
K(ii+1,ii) = off_diagonal*10;
end
end
M= round(M);
K= round(K);
[Veig,~] = eig(K,M);
[Veigs,~] = eigs(K,M,3);
Orthgonal_Check_eig = Veig'*M*Veig;
Orthgonal_Check_eigs =Veigs'*M*Veigs;
0 Comments
Accepted Answer
Christine Tobler
on 5 Jun 2024
For simple eigenvalue problem A*x = lambda*x, the eigenvalues are real and the eigenvectors can form an orthogonal basis only if A is symmetric. In the generalized case, this condition becomes that A must be symmetric and B must be symmetric positive definite.
So what you're expecting is true if K is symmetric (Hermitian if K is complex) and M is symmetric positive definite (Hermitian if M is complex). However, in this example M is not positive definite:
>> eig(M)
ans =
-32.1828
30.7556
37.5342
50.4526
110.4404
When M is not positive definite, there is no guarantee that the eigenvalues are all real, which you can see by computing them:
eig(K, M)
ans =
11.0502 + 0.0000i
11.1566 + 5.5626i
11.1566 - 5.5626i
6.8143 + 0.0000i
6.4945 + 0.0000i
This is also why your proof in the comment above doesn't work out: It's assuming that the eigenvectors and eigenvalues of (K, M) are real (and therefore a transpose instead of a conjugate transpose should be used). If you use conjugate transpose there, the terms don't cancel anymore there.
0 Comments
More Answers (1)
Torsten
on 5 Jun 2024
Edited: Torsten
on 5 Jun 2024
From the documentation:
[V,D] = eig(A,B) returns diagonal matrix D of generalized eigenvalues and full matrix V whose columns are the corresponding right eigenvectors, so that A*V = B*V*D.
I don't see where it is guaranteed that V'*B*V is a diagonal matrix. I also don't see that V should be a unitary matrix ( V' = inv(V) ) (maybe it is if B^(-1)*A is Hermitian).
If the inverse exists,
inv(B*V)*A*V
should be diagonal and equal D.
2 Comments
Torsten
on 5 Jun 2024
I already wrote the V' does not need to be inv(V).
From the equation
A*V = B*V*D
that MATLAB says to be true it follows that
inv(B*V)*A*V = D
if B*V is invertible - nothing less and nothing more.
See Also
Categories
Find more on Linear Algebra in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!