Symbolic eigenvectors returned by eig are incorrect
6 views (last 30 days)
Show older comments
Kamil Dylewicz
on 25 Apr 2023
Commented: Christine Tobler
on 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 Comment
Christine Tobler
on 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.
Accepted Answer
Steven Lord
on 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 Comments
Christine Tobler
on 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
More Answers (0)
See Also
Categories
Find more on Linear Algebra in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!