eig() gives wrong results for complex eigenvalues in comparison to eigs()
4 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
lennart
le 18 Mai 2015
Commenté : Nalini Vishnoi
le 20 Mai 2015
Hi,
I'm working on a eigenvalue problem of the following two matrices A and B
A = [M zeros(N); zeros(N) -K];
B = [zeros(N) M; M D];
M, D and K are submatrices of size NxN. Let's say:
N = 1;
K = rand(N); D = rand(N); M = rand(N);
The eigenvalue problem can be formulated as:
A*x = lambda*B*x
or
y.' *A = lambda * y.'*B
with the right and left eigenvectors or modal matrices x, y and X, Y.
Obviously, there are 6 ways to calculate the eigenvalue problem:
[X1, L1, Y1] = eig(A, B);
[X2, L21] = eigs(A, B);
[Y2, L22] = eigs(A.', B.');
%
[X3, L3, Y3] = eig(B \ A);
Y3 = (Y3.' / B).';
[X4, L41] = eigs(B \ A);
[Y4, L42] = eigs((B \ A).');
Y4 = (Y4.' / B).';
%
[X5, L5, Y5] = eig(A / B);
X5 = B \ X5;
[X6, L61] = eigs(A / B);
[Y6, L62] = eigs((A / B).');
X6 = B \ X6;
The corrections terms of the modal matrices Y3, Y4, X5 and X6 are required in case of using left and right multiplication of the matrix B.
We will find the eigenvalues lambda1-lambda6 as diagonal entries of the matrices L1-L6 to be equal.
But, the following equation should give us a zero solution vector
(A - lambda_i * B) * x_i = 0
y_j.' * (A - lambda_j * B) = 0
Unfortunately, this is not the case for (X1, Y1), (X3, Y3) and (X5, Y5) for complex eigenvalues. However, for completely real eigenvalues, they give the correct solution.
for qq = 1:length(l1)
%
right1(:, qq) = (A - l1(qq) * B) * X1(:, qq);
left1(qq, :) = Y1(:, qq).' * (A - l1(qq) * B);
%
right2(:, qq) = (A - l2(qq) * B) * X2(:, qq);
left2(qq, :) = Y2(:, qq).' * (A - l2(qq) * B);
%
right3(:, qq) = (A - l3(qq) * B) * X3(:, qq);
left3(qq, :) = Y3(:, qq).' * (A - l3(qq) * B);
%
right4(:, qq) = (A - l4(qq) * B) * X4(:, qq);
left4(qq, :) = Y4(:, qq).' * (A - l4(qq) * B);
%
right5(:, qq) = (A - l5(qq) * B) * X5(:, qq);
left5(qq, :) = Y5(:, qq).' * (A - l5(qq) * B);
%
right6(:, qq) = (A - l6(qq) * B) * X6(:, qq);
left6(qq, :) = Y6(:, qq).' * (A - l6(qq) * B);
end
[left1 right1;
left2 right2;
left3 right3;
left4 right4;
left5 right5;
left6 right6]
Can you please tell me, why "eig()" gives "wrong" eigenvectors in this case?
Thank you very much for your help, regards
Lennart
0 commentaires
Réponse acceptée
Nalini Vishnoi
le 19 Mai 2015
Hi Lennart,
It seems there are a few things that are not correct with the last part of the code:
1. Instead of li(qq), it should be li(qq,qq).
2. You should not use (.'), instead (') should be used. (.') computes the nonconjugate transpose of a matrix whereas (.) computes complex conjugate transpose. Since the eigenvalues and eigenvectors could be complex in nature even for a real matrix, the latter is the correct form to be used.
3. The code below gives you the desired results:
for qq = 1:length(L1)
%
right1(:, qq) = (A - L1(qq,qq) * B) * X1(:, qq);
left1(qq, :) = Y1(:, qq)' * (A - L1(qq,qq) * B);
%
right2(:, qq) = (A - L21(qq,qq) * B) * X2(:, qq);
left2(qq, :) = Y2(:, qq)' * (A - L22(qq,qq)' * B);
%
right3(:, qq) = (A - L3(qq,qq) * B) * X3(:, qq);
left3(qq, :) = Y3(:, qq)' * (A - L3(qq,qq) * B);
%
right4(:, qq) = (A - L41(qq,qq) * B) * X4(:, qq);
left4(qq, :) = Y4(:, qq)' * (A - L42(qq,qq)' * B);
%
right5(:, qq) = (A - L5(qq,qq) * B) * X5(:, qq);
left5(qq, :) = Y5(:, qq)' * (A - L5(qq,qq) * B);
%
right6(:, qq) = (A - L61(qq,qq) * B) * X6(:, qq);
left6(qq, :) = Y6(:, qq)' * (A - L62(qq,qq)' * B);
end
[left1 right1;
left2 right2;
left3 right3;
left4 right4;
left5 right5;
left6 right6]
Please note that I have changed Li2(qq,qq) to Li2(qq,qq)' for the computation of left2, left4 and left6 in the 2nd, 4th and 6th case. Here is why: As per the documentation of eig,
[Y2, L22] = eigs(A', B');
ensures that the solution satisfies the following:
A'*Y2 = B'*Y2*L22 % which gives
Y2'*A = L22'*Y2'*B % (taking transpose of both sides)
Y2'*A - L22'*Y2'*B = 0
Y2'*(A-L22'*B) = 0
(since values of L22 are taken one by one in your calculations, so it can behave as a scalar and moving its position does not change the result). However, taking its (complex conjugate) transpose is crucial since it is a complex number. Instead of computing the residuals inside a loop, it can be written as one single statement (treating L22 as a matrix):
left2 = Y2'*A - L22'*Y2'*B;
The same reasoning applies for cases 4 and 6. Please note this reasoning does not apply to cases 1, 3 and 5 and the original form of equation stays valid.
I hope the above helps.
Nalini
2 commentaires
Nalini Vishnoi
le 20 Mai 2015
Hi Lennart,
You are right, if L22 is a matrix, this is not true. However, I did mention in my previous answer that for this specific case I am treating L22 as a scalar (since you are calculating values for each eigenvalue one by one so shifting a scalar value would not change the answer). If L22 is treated as a matrix, the above won't hold true.
Plus de réponses (0)
Voir également
Catégories
En savoir plus sur Linear Algebra dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!