Errors in using null command due to truncation error

2023 8월 17
1 답변
업데이트 시간: 2023 8월 18
조회 수: 2 (30일)
Ran in:

0 개 추천

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일

2 개 추천

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개 표시 이전 댓글 4개 숨기기

Walter Roberson
Walter Roberson 2023년 8월 17일
Ran in:
Ran in:
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 = (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
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 NaN + 0i 0 + 0i 1 + 0i NaN + 0i 1.00659554466799 + 0i 0 + 0i NaN + 0i 1 + 0i 0 + 0i NaN + 0i 0 + 0i 0 + 0i NaN + 0i 0 + 0i 0 + 0i NaN + 0i 0 + 0i 0 + 0i NaN + 0i 0 + 0i 0 + 0i NaN + 0i 0 + 0i 0 + 0i NaN + 0i Columns 7 through 9 NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 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.996416379214725 + 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 NaN + 0i 0 + 0i 1 + 0i NaN + 0i 1.0132786693931 + 0i 0 + 0i NaN + 0i 1 + 0i 0 + 0i NaN + 0i 0 + 0i 0 + 0i NaN + 0i 0 + 0i 0 + 0i NaN + 0i 0 + 0i 0 + 0i NaN + 0i 0 + 0i 0 + 0i NaN + 0i 0 + 0i 0 + 0i NaN + 0i Columns 7 through 9 NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 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 NaN + 0i 0 + 0i 1 + 0i NaN + 0i 1.02005113025445 + 0i 0 + 0i NaN + 0i 1 + 0i 0 + 0i NaN + 0i 0 + 0i 0 + 0i NaN + 0i 0 + 0i 0 + 0i NaN + 0i 0 + 0i 0 + 0i NaN + 0i 0 + 0i 0 + 0i NaN + 0i 0 + 0i 0 + 0i NaN + 0i Columns 7 through 9 NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i NaN + 0i
Walter Roberson
Walter Roberson 2023년 8월 17일
Note that in this above code, the null() calculation is working on the symbolic solutions, so there is no truncation error going on.
Seung Hyeop Hyun
Seung Hyeop Hyun 2023년 8월 18일
Walter,
Thank you for leaving your comments. The issue here is that there must be eigenvectors corresponding to the eigenvalues, which can be obtained from the 'null', and also the dimension of the eigenvectors is definitely 9x1. The 'null' cannot find the eigenvectors and provides 0X1 empty vectors. I believe that it's because of the approximation due to using 'vpa'. The 'a' is already replaced by one of eigenvalues, and therefore, the null is not working on the symbolic solutions. Alternatively, I can do 'eig' for M after one of eigenvalues is entered.
Matt J
Matt J 2023년 8월 18일
편집: Matt J 2023년 8월 18일
There doesn't appear to be any reason for you to be using symbolic operations for this problem. Eigenvalue problems generally do not have closed form solutions. Also, all of your problem input parameters (e.g. R,r,s) have a priori known values.
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:
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

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

카테고리

도움말 센터 및 File Exchange에서 Linear Algebra에 대해 자세히 알아보기

제품

릴리스

R2017b

태그

질문:

Seung Hyeop Hyun
Seung Hyeop Hyun
2023년 8월 17일

댓글:

Walter Roberson
Walter Roberson
2023년 8월 18일

Community Treasure Hunt

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

Start Hunting!

Translated by