Symbolic eigenvectors returned by eig are incorrect
6 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Kamil Dylewicz
am 25 Apr. 2023
Kommentiert: Christine Tobler
am 25 Apr. 2023
Hi Community,
I am trying to work with symbolic eigenvalues and eigenvectors and a small, 5-by-5, matrix.
I am able to recover correct eigenvalues, that satisfy the definition, but when I compute the eigenvectors they do not satisfy the definition of eigenvectors.
According to the documentation (https://uk.mathworks.com/help/symbolic/eig.html), given [vecR, lambda] = eig(A), if vecR is the same size as A, then the matrix A has a full set of linearly independent eigenvectors that satisfy A*vecR = vecR*lambda.
Firstly, I check that all eigenvectors are linearly independent and follow that with check of eigenvalues according to the definition. Yet, 'Check 1' in eigenvectors fail. Do you have any ideas why?
Full code is attached below.
Thanks in advance.
clc; clear;
% Define variables
syms U R T E Gamma Ma
% Define Jacobian
A = [0, 0, R, 0, 0; 0, 0, R*U, 0, 0; T/(Gamma*Ma^2), 0, 0, 0, R/(Gamma*Ma^2); ...
0, 0, 0, 0, 0; 0, 0, R*(E + T/(Gamma*Ma^2)), 0, 0];
%% Find eigenvalues and right eigenvectors
[vecR, lambda, p] = eig(A);
if (length(A) == length(p))
fprintf('All eigenvectors are linearly independent. \n')
end
%% Check eigenvalues
% Check 1
fprintf('Checking eigenvalues ... ')
msg = 'Eigenvalues do not satisfy the characteristic polynomial!';
assert(~any(det(A - lambda)), msg)
fprintf('[PASS] \n')
%% Check eigenvectors
% Check 1
fprintf('Checking right eigenvectors ... ')
dummy1 = A*vecR;
dummy2 = vecR*lambda;
cond = isequaln(dummy1, dummy2);
msg = 'Right eigenvectors do not satisfy the eigenvalue problem!';
assert(cond, msg)
fprintf('[PASS] \n')
1 Kommentar
Christine Tobler
am 25 Apr. 2023
Calling
>> simplify(dummy1 - dummy2)
ans =
[0, 0, 0, 0, 0]
[0, 0, 0, 0, 0]
[0, 0, 0, 0, 0]
[0, 0, 0, 0, 0]
[0, 0, 0, 0, 0]
seems to show that the eigenvectors are correct, the problem would be that isequaln doesn't recognize that the two symbolic expressions in dummy1 and dummy2 are equal.
Akzeptierte Antwort
Steven Lord
am 25 Apr. 2023
clc; clear;
% Define variables
syms U R T E Gamma Ma
% Define Jacobian
A = [0, 0, R, 0, 0; 0, 0, R*U, 0, 0; T/(Gamma*Ma^2), 0, 0, 0, R/(Gamma*Ma^2); ...
0, 0, 0, 0, 0; 0, 0, R*(E + T/(Gamma*Ma^2)), 0, 0];
%% Find eigenvalues and right eigenvectors
[vecR, lambda, p] = eig(A);
if (length(A) == length(p))
fprintf('All eigenvectors are linearly independent. \n')
end
%% Check eigenvalues
% Check 1
fprintf('Checking eigenvalues ... ')
msg = 'Eigenvalues do not satisfy the characteristic polynomial!';
assert(~any(det(A - lambda)), msg)
fprintf('[PASS] \n')
%% Check eigenvectors
% Check 1
fprintf('Checking right eigenvectors ... ')
dummy1 = A*vecR;
dummy2 = vecR*lambda;
Checking with isequaln is not the right tool to use. Check that the difference between the two is all zero.
cond = all(simplify(dummy1-dummy2) == 0, 'all');
msg = 'Right eigenvectors do not satisfy the eigenvalue problem!';
assert(cond, msg)
fprintf('[PASS] \n')
2 Kommentare
Christine Tobler
am 25 Apr. 2023
The following works for me:
[vecL, lam_ch] = eig(A.');
dummy1 = vecL.'*A;
dummy2 = lam_ch*vecL.';
cond = all(simplify(dummy1 - dummy2) == 0, 'all')
Note .' corresponds to calling transpose
Weitere Antworten (0)
Siehe auch
Kategorien
Mehr zu Linear Algebra finden Sie in Help Center und File Exchange
Produkte
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!