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!