Errors in using null command due to truncation error

조회 수: 5 (최근 30일)
Seung Hyeop Hyun
Seung Hyeop Hyun 2023년 8월 17일
댓글: Walter Roberson 2023년 8월 18일
Ran in:
Hi,
I'm trying to find eigenvectors of a 9-by-9 square matrix, corresponding to its eigenvalues. The matrix consists of components with complex numbers and one symbolic 'a', so I found nine eigenvalues ('a'), via solving the determinant of the matrix. For some eigenvalues, I used the 'vpa' command since, without 'vpa', they are obtained as a form of 'root(eqn, z, integer)'. Here the issue seems to arise. Due to the truncation error, the 'null' shows empty eigenvectors corresponding to the eigenvalues. FYI, I don't know how to assign a specific variable in 'eig' and it takes forever to run. Is there a breakthrough other than Gauss elimination method with suppressing close-to-zero values?
clc;
close all;
clear;
syms a
R = 0.5234;
r = 0.0054;
s = 0.0084;
for p = 1:3
M = [a*R*r 1i*r 0 0 0 0 a*R*8.78+1i*(8.78-p*R) 0 0;
0 0 a*R*s 1i*s 0 0 0 a*R*79.88+1i*(79.88-p*R) 0;
0 0 0 0 a*R*0.4542 1i*0.4542 -(a*R*291.06+1i*(291.06+p*R)) -(a*R*291.06+1i*(291.06+p*R)) 0;
a*R*8.78+1i*(8.78-p*R) 0 0 0 0 0 0 0 a*R;
0 a*R*8.78+1i*(8.78-p*R) 0 0 0 0 0 0 1i;
0 0 a*R*79.88+1i*(79.88-p*R) 0 0 0 0 0 a*R;
0 0 0 a*R*79.88+1i*(79.88-p*R) 0 0 0 0 1i;
0 0 0 0 a*R*291.06+1i*(291.06+p*R) 0 0 0 a*R;
0 0 0 0 0 a*R*291.06+1i*(291.06+p*R) 0 0 1i];
Det = det(M);
DetEqn = Det == 0;
EigenVal1 = solve(DetEqn,a);
EigVal = vpa(EigenVal1);
for j=1:rank(M)
M_temp = subs(M,a,EigVal(j));
EigVec(:,j) = null(M_temp)
end
end
Unable to perform assignment because the indices on the left side are not compatible with the size of the right side.

Error in sym/privsubsasgn (line 1168)
L_tilde2 = builtin('subsasgn',L_tilde,struct('type','()','subs',{varargin}),R_tilde);

Error in indexing (line 999)
C = privsubsasgn(L,R,inds{:});

이 질문에 답변하려면 로그인하십시오.

답변 (1개)

Walter Roberson
Walter Roberson 2023년 8월 17일
Your code assumes that the null space is the same size each time, but most of the time the null space is empty. You cannot store an empty vector into a definite vector location.
You need to decide what you want to do when the null space is empty.
  댓글 수: 6
이전 댓글 4개 표시
Torsten
Torsten 2023년 8월 18일
@Seung Hyeop Hyun
I don't know why you talk about "eigenvalues", but I agree that if "a" gives det(M(a)) = 0, null(M(a)) should be at least 1-dimensional and not empty.
Walter Roberson
Walter Roberson 2023년 8월 18일
Ran in:
@Seung Hyeop Hyun
You did not take into account that you use floating point constants and that some of the calculations take place in floating point instead of as symbolic numbers.
When you use symbolic numbers consistently then the problem does not show up.
Q = @(v) sym(v);
syms a
R = Q(5234)/Q(10)^4;
r = Q(54)/Q(10)^4;
s = Q(84)/Q(10)^4;
n8_78 = Q(878)/Q(10)^2;
n79_88 = Q(7988)/Q(10)^2;
n_4542 = Q(4542)/Q(10)^4;
n291_06 = Q(29106)/Q(10)^2;
for p = 1:3
M = [a*R*r 1i*r 0 0 0 0 a*R*n8_78+1i*(n8_78-p*R) 0 0;
0 0 a*R*s 1i*s 0 0 0 a*R*n79_88+1i*(n79_88-p*R) 0;
0 0 0 0 a*R*n_4542 1i*n_4542 -(a*R*n291_06+1i*(n291_06+p*R)) -(a*R*n291_06+1i*(n291_06+p*R)) 0;
a*R*n8_78+1i*(n8_78-p*R) 0 0 0 0 0 0 0 a*R;
0 a*R*n8_78+1i*(n8_78-p*R) 0 0 0 0 0 0 1i;
0 0 a*R*n79_88+1i*(n79_88-p*R) 0 0 0 0 0 a*R;
0 0 0 a*R*n79_88+1i*(n79_88-p*R) 0 0 0 0 1i;
0 0 0 0 a*R*n291_06+1i*(n291_06+p*R) 0 0 0 a*R;
0 0 0 0 0 a*R*n291_06+1i*(n291_06+p*R) 0 0 1i];
Det = det(M);
DetEqn = Det == 0;
EigenVal1 = solve(DetEqn,a);
EigVal = EigenVal1;
for j=1:rank(M)
M_temp = subs(M,a,EigVal(j));
EV = null(M_temp);
if isempty(EV)
EigVec(:,j,p) = sym(NaN(size(EV,1),1));
else
EigVec(:,j,p) = EV;
end
end
end
EigVec
EigVec(:,:,1) = 
EigVec(:,:,2) = 
EigVec(:,:,3) = 
format long g
EigVec = double(EigVec)
EigVec =
EigVec(:,:,1) = Columns 1 through 3 -0.0604434080666857 - 0.0568402099138265i -0.0604434080666857 + 0.0568402099138265i 0 + 0i -0.0568402099138265 + 0.0604434080666857i -0.0568402099138265 - 0.0604434080666857i 0 + 0i -0.00630053702136478 - 0.00625925383312013i -0.00630053702136478 + 0.00625925383312013i 0 + 0i -0.00625925383312013 + 0.00630053702136478i -0.00625925383312013 - 0.00630053702136478i 0 + 0i -0.00171477249055503 - 0.00171785608816912i -0.00171477249055503 + 0.00171785608816912i 0.998204973259795 + 0i -0.00171785608816912 + 0.00171477249055503i -0.00171785608816912 - 0.00171477249055503i 1 + 0i 0 + 0i 0 + 0i 0 + 0i 0 + 0i 0 + 0i 0 + 0i 1 + 0i 1 + 0i 0 + 0i Columns 4 through 6 0 + 0i 1.06339171087373 + 0i -2.07917517265388 - 0.12335357734933i 0 + 0i 1 + 0i 2.08986241502569 + 0.131173171659899i 1.00659554466799 + 0i 0 + 0i -1.17431482491886 - 3.23249177899391i 1 + 0i 0 + 0i 1.16945872447306 + 3.25381182291119i 0 + 0i 0 + 0i 0.389692214298109 - 0.187973926292369i 0 + 0i 0 + 0i -0.392422256869245 + 0.187636508068213i 0 + 0i 0 + 0i 0.0467535877540691 + 0.00573056182938062i 0 + 0i 0 + 0i -0.15363302189684 + 0.127699515578846i 0 + 0i 0 + 0i 1 + 0i Columns 7 through 9 -2.07917517265388 + 0.12335357734933i -1.91125259969198 - 0.520193365405037i -1.91125259969198 + 0.520193365405037i 2.08986241502569 - 0.131173171659899i 1.91129494286941 + 0.553169312823223i 1.91129494286941 - 0.553169312823223i -1.17431482491886 + 3.23249177899391i -0.105661537118353 - 0.770245234313409i -0.105661537118353 + 0.770245234313409i 1.16945872447306 - 3.25381182291119i 0.0937570861782642 + 0.775325421161631i 0.0937570861782642 - 0.775325421161631i 0.389692214298109 + 0.187973926292369i 0.0716575680173222 - 0.185960430138347i 0.0716575680173222 + 0.185960430138347i -0.392422256869245 - 0.187636508068213i -0.0749584912828433 + 0.185626626193629i -0.0749584912828433 - 0.185626626193629i 0.0467535877540691 - 0.00573056182938062i 0.0363384221349967 + 0.0221560988269227i 0.0363384221349967 - 0.0221560988269227i -0.15363302189684 - 0.127699515578846i -0.00986540081559587 + 0.00258850280209401i -0.00986540081559587 - 0.00258850280209401i 1 + 0i 1 + 0i 1 + 0i EigVec(:,:,2) = Columns 1 through 3 -0.0641387232677706 - 0.0564917511132487i -0.0641387232677706 + 0.0564917511132487i 0 + 0i -0.0564917511132487 + 0.0641387232677706i -0.0564917511132487 - 0.0641387232677706i 0 + 0i -0.00634195365255762 - 0.00625884452532305i -0.00634195365255762 + 0.00625884452532305i 0 + 0i -0.00625884452532305 + 0.00634195365255762i -0.00625884452532305 - 0.00634195365255762i 0 + 0i -0.00171169167540565 - 0.00171784779045346i -0.00171169167540565 + 0.00171784779045346i 0.996416379214726 + 0i -0.00171784779045346 + 0.00171169167540565i -0.00171784779045346 - 0.00171169167540565i 1 + 0i 0 + 0i 0 + 0i 0 + 0i 0 + 0i 0 + 0i 0 + 0i 1 + 0i 1 + 0i 0 + 0i Columns 4 through 6 0 + 0i 1.13536440283453 + 0i -1.03424396514104 - 0.0577669915193808i 0 + 0i 1 + 0i 1.04493120751284 + 0.0655865858299493i 1.0132786693931 + 0i 0 + 0i -0.589585462682334 - 1.60558586753832i 1 + 0i 0 + 0i 0.58472936223653 + 1.62690591145559i 0 + 0i 0 + 0i 0.193481085863487 - 0.0941556722582628i 0 + 0i 0 + 0i -0.196211128434622 + 0.0938182540341066i 0 + 0i 0 + 0i 0.0116310771845003 + 0.00138540923306503i 0 + 0i 0 + 0i -0.0380957587743605 + 0.0318852196758995i 0 + 0i 0 + 0i 1 + 0i Columns 7 through 9 -1.03424396514104 + 0.0577669915193808i -0.955605128257274 - 0.243608708993425i -0.955605128257274 + 0.243608708993425i 1.04493120751284 - 0.0655865858299493i 0.955647471434705 + 0.276584656411611i 0.955647471434705 - 0.276584656411611i -0.589585462682334 + 1.60558586753832i -0.058782994029221 - 0.382582523732593i -0.058782994029221 + 0.382582523732593i 0.58472936223653 - 1.62690591145559i 0.0468785430891321 + 0.387662710580816i 0.0468785430891321 - 0.387662710580816i 0.193481085863487 + 0.0941556722582628i 0.0341783223759006 - 0.0931471170415329i 0.0341783223759006 + 0.0931471170415329i -0.196211128434622 - 0.0938182540341066i -0.0374792456414217 + 0.0928133130968145i -0.0374792456414217 - 0.0928133130968145i 0.0116310771845003 - 0.00138540923306503i 0.00912923446797742 + 0.00536880051870816i 0.00912923446797742 - 0.00536880051870816i -0.0380957587743605 - 0.0318852196758995i -0.00244438930219064 + 0.000683128476076533i -0.00244438930219064 - 0.000683128476076533i 1 + 0i 1 + 0i 1 + 0i EigVec(:,:,3) = Columns 1 through 3 -0.0680253050205014 - 0.0558597772365388i -0.0680253050205014 + 0.0558597772365388i 0 + 0i -0.0558597772365388 + 0.0680253050205014i -0.0558597772365388 - 0.0680253050205014i 0 + 0i -0.00638363887049742 - 0.00625815577392186i -0.00638363887049742 + 0.00625815577392186i 0 + 0i -0.00625815577392186 + 0.00638363887049742i -0.00625815577392186 - 0.00638363887049742i 0 + 0i -0.0017086164151074 - 0.00171783399737567i -0.0017086164151074 + 0.00171783399737567i 0.99463418334813 + 0i -0.00171783399737567 + 0.0017086164151074i -0.00171783399737567 - 0.0017086164151074i 1 + 0i 0 + 0i 0 + 0i 0 + 0i 0 + 0i 0 + 0i 0 + 0i 1 + 0i 1 + 0i 0 + 0i Columns 4 through 6 0 + 0i 1.21778690116231 + 0i -0.685933562636758 - 0.0359047962427311i 0 + 0i 1 + 0i 0.696620805008562 + 0.0437243905532995i 1.02005113025445 + 0i 0 + 0i -0.394675675270157 - 1.06328389705312i 1 + 0i 0 + 0i 0.389819574824353 + 1.0846039409704i 0 + 0i 0 + 0i 0.12807737638528 - 0.0628829209135606i 0 + 0i 0 + 0i -0.130807418956415 + 0.0625455026894044i 0 + 0i 0 + 0i 0.00514395587663915 + 0.000594946349441334i 0 + 0i 0 + 0i -0.0167933653854523 + 0.0141531959048798i 0 + 0i 0 + 0i 1 + 0i Columns 7 through 9 -0.685933562636758 + 0.0359047962427311i -0.637055971112372 - 0.151413823522888i -0.637055971112372 + 0.151413823522888i 0.696620805008562 - 0.0437243905532995i 0.637098314289803 + 0.184389770941074i 0.637098314289803 - 0.184389770941074i -0.394675675270157 + 1.06328389705312i -0.0431568129995103 - 0.253361620205655i -0.0431568129995103 + 0.253361620205655i 0.389819574824353 - 1.0846039409704i 0.0312523620594214 + 0.258441807053877i 0.0312523620594214 - 0.258441807053877i 0.12807737638528 + 0.0628829209135606i 0.0216852404954267 - 0.0622093460092614i 0.0216852404954267 + 0.0622093460092614i -0.130807418956415 - 0.0625455026894044i -0.0249861637609478 + 0.061875542064543i -0.0249861637609478 - 0.061875542064543i 0.00514395587663915 - 0.000594946349441334i 0.0040759677296243 + 0.00231048172032896i 0.0040759677296243 - 0.00231048172032896i -0.0167933653854523 - 0.0141531959048798i -0.0010764184837051 + 0.000319388109428762i -0.0010764184837051 - 0.000319388109428762i 1 + 0i 1 + 0i 1 + 0i

댓글을 달려면 로그인하십시오.

카테고리

Help Center 및 File Exchange에서 Linear Algebra에 대해 자세히 알아보기

태그

제품


릴리스

R2017b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by