Symbolic eigenvectors returned by eig are incorrect
6 次查看(过去 30 天)
显示 更早的评论
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 个评论
Christine Tobler
2023-4-25
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.
采纳的回答
Steven Lord
2023-4-25
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 个评论
Christine Tobler
2023-4-25
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
更多回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Linear Algebra 的更多信息
产品
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!