Errors in using null command due to truncation error

4 views (last 30 days)
Seung Hyeop Hyun
Seung Hyeop Hyun on 17 Aug 2023
Commented: Walter Roberson on 18 Aug 2023
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{:});

Sign in to answer this question.

Answers (1)

Walter Roberson
Walter Roberson on 17 Aug 2023
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 Comments
Show 4 older comments
Torsten
Torsten on 18 Aug 2023
@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 on 18 Aug 2023
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

Sign in to comment.

Categories

Find more on Linear Algebra in Help Center and File Exchange

Tags

Products


Release

R2017b

Community Treasure Hunt

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

Start Hunting!