Errors in using null command due to truncation error

3 ビュー (過去 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