이 질문을 팔로우합니다.
- 팔로우하는 게시물 피드에서 업데이트를 확인할 수 있습니다.
- 정보 수신 기본 설정에 따라 이메일을 받을 수 있습니다.
How to avoid using det, when looking for the complex root w with det(M(w)) = 0
조회 수: 3 (최근 30일)
이전 댓글 표시
ma Jack
2022년 7월 6일
Hi all,
I want to find a complex number w such that det(M(w)) converges to 0 (note that M is a matrix and it is a function of w), but the det function seems to have a large error, how can I avoid using it?M is an 8*8 matrix, so it would be very complicated to write out its determinant expression.
Thanks in advance.
댓글 수: 6
ma Jack
2022년 7월 6일
GvXY is the matrix M,Zero_seek_Main() is the main function
function []=Zero_seek_Main()
clear;clc;close all
%=====================
d=500*1e-9;
rx=d/4;
ry=rx;
num_Re=5;
num_Im=5;
%====================
%==========%
lsqoptions = optimoptions(@lsqnonlin,'Display','none');
%==========%
Real_w0=linspace(0.8,1.2,num_Re);
Imag_w0=linspace(-0.2,0.2,num_Im);
Result_stoXY=zeros(num_Re*num_Im,1);
count=0;
for uu=1:num_Re
for vv=1:num_Im
count=count+1;
mark=count/(num_Re*num_Im)
Re_w=Real_w0(uu);
Im_w=Imag_w0(vv);
w=Re_w+1j*Im_w;
FUNXY=@(x)det(GvXY(x,d,rx,ry));
%======================
[find_x,~,~,~]=lsqnonlin(FUNXY,w,[],[],lsqoptions);
Result_stoXY(count,1)=find_x;
%======================
end
end
figure
Result_index=[1:length(Result_stoXY)].';
plot(Result_index,Result_stoXY,'bo')
end
function [ModeXY]=GvXY(w,d,rx,ry)
c0=299792458;%light
wsp=3.742*1e15;%rad/s
eps0=8.8541*1e-12;
k=w*wsp/c0;
M1=RegSigGtXY(k,d);%size:2*2
M2=SigGtXY(k,rx,ry,d);%size:2*2
ModeXY=[M1,M2,M2,M2;...
M2,M1,M2,M2;...
M2,M2,M1,M2;...
M2,M2,M2,M1];%size:(2*4,2*4)
ModeXY=(k^2)/eps0.*ModeXY;
ModeXY(isnan(ModeXY))=0;%set NAN to 0
ModeXY(isinf(ModeXY))=0;%set INF to 0
%ModeXY=ModeXY(1:4,1:4);%If I enable this command. Then lsqnonlin will not
%report an error, if I increase the
%size of the matrix (for example, 5*5) lsqnonlin will report an error
end
function [o4]=RegSigGtXY(k,d)
N=1;
index_sto=zeros((2*N+1)^2-1,2);
cc=0;
for m=-N:N
for n=-N:N
if m==0&&n==0
cc=cc;
else
cc=cc+1;
index_sto(cc,1:2)=[m,n];
end
end
end
Sum_sto11=zeros((2*N+1)^2-1,1);%xx
Sum_sto22=zeros((2*N+1)^2-1,1);%yy
Sum_sto12=zeros((2*N+1)^2-1,1);%xy
Sum_sto21=zeros((2*N+1)^2-1,1);%yx
parfor ii=1:size(index_sto,1)
m=index_sto(ii,1);
n=index_sto(ii,2);
%==================
nR=abs(m*d+1j*n*d);
Rx=m*d;
Ry=n*d;
B=BB(k,nR);
A=AA(k,nR);
term=exp(1j*k*nR)/(4*pi*nR);
%==================
Sum_sto11(ii,1)=B/(nR^2)*(Rx^2)*term+A*term;
Sum_sto22(ii,1)=B/(nR^2)*(Ry^2)*term+A*term;
Sum_sto12(ii,1)=B/(nR^2)*Rx*Ry*term;
Sum_sto21(ii,1)=B/(nR^2)*Ry*Rx*term;
end
o4=zeros(2,2);
o4(1,1)=sum(Sum_sto11);
o4(2,2)=sum(Sum_sto22);
o4(1,2)=sum(Sum_sto12);
o4(2,1)=sum(Sum_sto21);
end
function [o3]=SigGtXY(k,rx,ry,d)
N=1;
index_sto=zeros((2*N+1)^2,2);
cc=0;
for m=-N:N
for n=-N:N
cc=cc+1;
index_sto(cc,1:2)=[m,n];
end
end
Sum_sto11=zeros((2*N+1)^2,1);%xx
Sum_sto22=zeros((2*N+1)^2,1);%yy
Sum_sto12=zeros((2*N+1)^2,1);%xy
Sum_sto21=zeros((2*N+1)^2,1);%yx
parfor ii=1:size(index_sto,1)
m=index_sto(ii,1);
n=index_sto(ii,2);
%=====================
nR=abs(m*d+n*d*1j+rx+ry*1j);
Rx=m*d+rx;
Ry=n*d+ry;
Term=exp(1j*k*nR)/(4*pi*nR);
A=AA(k,nR);
B=BB(k,nR);
%=====================
Sum_sto11(ii,1)=B/(nR^2)*(Rx^2)*Term+A*Term;
Sum_sto22(ii,1)=B/(nR^2)*(Ry^2)*Term+A*Term;
Sum_sto12(ii,1)=B/(nR^2)*Rx*Ry*Term;
Sum_sto21(ii,1)=B/(nR^2)*Ry*Rx*Term;
end
o3=zeros(2,2);
o3(1,1)=sum(Sum_sto11);
o3(2,2)=sum(Sum_sto22);
o3(1,2)=sum(Sum_sto12);
o3(2,1)=sum(Sum_sto21);
end
function [o1]=AA(k,R)
o1=1+(1j*k*R-1)/((k*R)^2);
end
function [o2]=BB(k,R)
o2=(3-3*1j*k*R-(k*R)^2)/((k*R)^2);
end
ma Jack
2022년 7월 6일
x is w, which is defined in the following code:
FUNXY=@(x)det(GvXY(x,d,rx,ry));
%======================
[find_x,~,~,~]=lsqnonlin(FUNXY,w,[],[],lsqoptions);
Walter Roberson
2022년 7월 6일
As discussed in your previous question, you have the difficulty that you are working with a matrix whose determinant is on the order of 10^250
On the first ModeXY matrix that is generated, there are only four unique values. If you construct a representative symbolic matrix in terms of the pattern of unique values, and take det() of it, and substitute in the unique values, then that is a lot faster than calculating det() of the original matrix symbolically. The idea of calculating it symbolically being to reduce the error in the calculation of det()
ModeXY = [M1,M2,M2,M2;...
M2,M1,M2,M2;...
M2,M2,M1,M2;...
M2,M2,M2,M1];%size:(2*4,2*4)
That promises that the pattern continues of there being only 4 unique elements in the matrix, so you can
pre-calculate the determinant as
(V1 + V2 - V3 - V4)^3*(V1 - V2 - V3 + V4)^3*(V1 + V2 + 3*V3 + 3*V4)*(V1 - V2 + 3*V3 - 3*V4)
which would be 0 if and only if any one of the sub-expressions is 0, which happens if
V1 + V2 = V3 + V4
V1 + V4 = V2 + V3
V1 + 3*V3 = V2 + 3*V4
V1 + V2 + 3*V3 + 3*V4 == 0
You might be able to take advantage of those to seek for a zero with a lower range.
답변 (1개)
Matt J
2022년 7월 6일
편집: Matt J
2022년 7월 7일
Because your matrix appears to be symmetric, I suggest minimizing instead norm(M(w)) which is the maximum absolute eigenvalue of M. This is the same as forcing M to be singular.
EDIT: rcond(M) is probably more appropriate than norm(M)
Additionally, I suggest using fminsearch instead of lsqnonlin, since you only have a small number of variables and a non-differentiable cost function. Be mindful, however, that you must express your objective function in terms of a vector z of real variables.
w=@(z) complex(z(1),z(2));
zopt=fminsearch(@(z) norm(M(w(z))) ,z0)
wopt=w(zopt)
댓글 수: 21
ma Jack
2022년 7월 6일
Sorry, I don't quite understand why norm(M(w)) is used instead of det(M(w)), is this some kind of mathematical knowledge?
Also, is there any advantage of fminsearch over lsqnonlin? Why should we use it
Matt J
2022년 7월 6일
편집: Matt J
2022년 7월 6일
I don't quite understand why norm(M(w)) is used instead of det(M(w)),
John demonstrated the problems with det() in your earlier post. Let's see how norm(M) does in similar examples:
M1=eye(1200)/2;
M2=0*M1;
det(M1), det(M2)
ans = 0
ans = 0
norm(M1), norm(M2)
ans = 0.5000
ans = 0
Clearly norm(M) is better able to distinguish between singular and non-singular matrices than det. rcond() would be even better because it is independent of the scaling of M1:
rcond(M1),rcond(1e-100*M1), rcond(M2)
ans = 1
ans = 1
ans = 0
is there any advantage of fminsearch over lsqnonlin?
lsqnonlin assumes that the cost function is differentiable whereas fminsearch does not. It's not clear to me that your M(w) is a differentiable function of w, so fminsearch might therefore be safer. But you can also try lsqnonlin if you wish.
Matt J
2022년 7월 6일
ma Jack's comment moved here:
Thank you very much for your reply sir, but I think I need to spend some time to understand and test them.
ma Jack
2022년 7월 6일
Sorry sir, I have one more doubt, if we find a complex w such that norm(M(w)) -> 0,will this also make det(M(w)) -> 0? Please excuse my tardiness.
Paul
2022년 7월 7일
For a matrix M, norm(M) is the maximum singular value of M. M can be singular even though norm(M) ~=0. Of course, norm(M) == 0 means M is singular, but that condition seems overly restrictive (I think norm(M) == 0 is only true for the zero matrix?) .
M = [1 0;0 0];
norm(M)
ans = 1
det(M)
ans = 0
Am I missing something?
Bruno Luong
2022년 7월 7일
편집: Bruno Luong
2022년 7월 7일
And of course det(M(w)) can converge to 0, "independent" of smalest/largest sinigular value
M1(w) = w*eye(n)
and
M2(w) = diag(1/w, w)
The question as it stands has very little interest mathematically, and if M is numercial approximated, one can never estimate accuracy smalleste eigen value below eps(1)*norm(M); which prevent ANY numerical algorithm to solve the problem OP asks.
In other word asking question on det is mostly a dead end numericaly in gerenal.
Matt J
2022년 7월 7일
편집: Matt J
2022년 7월 7일
I've edited my post to recommend rcond(M) rather than norm(M). As Bruno's examples suggest, rcond would also be problematic without an upper bound on the largest eigenvalue. However, from the OP's code, it appears that the M(i,j) in this case are all globally bounded functions of w. They are just sums of complex exponential terms whos phases are functions of w. That would be enough to globally bound the largest eigenvalue.
Bruno Luong
2022년 7월 7일
On a side note, numerically it is more reasonable to work with log(det(A)). It has some physical interpretation in some circumstance and equivalent to det(log(A)) such as A is definite positive matrix.
If the eigen values are purely complex phase, log(det) transforms to the sim of the phases.
ma Jack
2022년 7월 7일
Thanks for the discussion, I have learned a lot and I would like to add a fact here which might help to solve my problem better:
My goal is to find a complex w such that det(ModeXY)->0 (since this ensures that the eigenvector of an eigenequation like H\psi=lam\psi is not a 0-vector), but note that the w found is in the range where its real part is around 1 (roughly 0.6 to 1.5) and its imaginary part is negative (roughly -0.1 to 0), since my problem is relevant to the physical world, the correct w must be within this range, based on experimental observations.
Bruno Luong
2022년 7월 7일
편집: Bruno Luong
2022년 7월 7일
What is the unit and physical interpretation of the determinant of M for you?
If you can't answer my question accurately I would bet you stil not totally clear about the claim that you know how to translate the physics into maths.
Ifthe unit is (eV)^8 you might also have trouble and think about the physics again.
The only possibility of valid answer is det(M) has unitless unit. Is it the case?
ma Jack
2022년 7월 7일
Uh .... I'm sorry, I didn't think too hard about this, but these formulas should be fine, I've seen people looking for complex roots on these formulas, but I'm not sure what method they used to find them.
Matt J
2022년 7월 7일
편집: Matt J
2022년 7월 7일
My goal is to find a complex w such that det(ModeXY)->0 (since this ensures that the eigenvector of an eigenequation like H\psi=lam\psi is not a 0-vector)
Eigenvectors are never 0-vectors, regardless of the determinant of H.
If you meant that you want to make sure H has a trivial null-space, det(H)=0 would be the opposite of what you want.
ma Jack
2022년 7월 8일
Sorry I did't make it clear, not det(H), but det(H-eye*lam), in addition I use norm found that the effect is not det good (I scaled the problem so that det can work).
Matt J
2022년 7월 8일
편집: Matt J
2022년 7월 8일
Sorry I did't make it clear, not det(H), but det(H-eye*lam),
Still not clear. The condition det(H-eye*lam)=0 ensures that lam is an eigenvalue. Is that what you want?
in addition I use norm found that the effect is not det good
Of course not. Above, we said the idea of using norm(M) should be abandoned. What about rcond?
Bruno Luong
2022년 7월 8일
det(H-eye*lam) (= 0)
So your problem is find the eigen values/vectors or H. Sorry but why not simply tell us that at the first place instead of asking a question with innacurate description that is not mathematically equivalent?
ma Jack
2022년 7월 8일
Sorry please give me some time to redescribe the problem again, also the above code does not reflect my problem very accurately (due to simplification), I am now testing the difference between the various methods, it will take a while.
Bruno Luong
2022년 7월 8일
"due to simplification"
Please never do that when asking question. You would think it helps? Generally no, you should ask accurately question, since a simplified problem can have completily solution.
ma Jack
2022년 7월 12일
Hi everyone, interested in following my latest issue?I encountered a very, very strange problem when using the fminsearch+rcond function(Same code, same computer, same matlab version, different running times and different results).
Rosalinda
2024년 4월 9일
hello ma jack..i encounter the same problem for 8by 8 matrix..if you could please tell me how you resolve your issue
참고 항목
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!오류 발생
페이지가 변경되었기 때문에 동작을 완료할 수 없습니다. 업데이트된 상태를 보려면 페이지를 다시 불러오십시오.
웹사이트 선택
번역된 콘텐츠를 보고 지역별 이벤트와 혜택을 살펴보려면 웹사이트를 선택하십시오. 현재 계신 지역에 따라 다음 웹사이트를 권장합니다:
또한 다음 목록에서 웹사이트를 선택하실 수도 있습니다.
사이트 성능 최적화 방법
최고의 사이트 성능을 위해 중국 사이트(중국어 또는 영어)를 선택하십시오. 현재 계신 지역에서는 다른 국가의 MathWorks 사이트 방문이 최적화되지 않았습니다.
미주
- América Latina (Español)
- Canada (English)
- United States (English)
유럽
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
아시아 태평양
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)