現在この質問をフォロー中です
- フォローしているコンテンツ フィードに更新が表示されます。
- コミュニケーション基本設定に応じて電子メールを受け取ることができます。
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.
ma Jack
2022 年 7 月 7 日
Sir thank you for your suggestion, but I think Mr. Matt J's answer is better.
回答 (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!エラーが発生しました
ページに変更が加えられたため、アクションを完了できません。ページを再度読み込み、更新された状態を確認してください。
Web サイトの選択
Web サイトを選択すると、翻訳されたコンテンツにアクセスし、地域のイベントやサービスを確認できます。現在の位置情報に基づき、次のサイトの選択を推奨します:
また、以下のリストから Web サイトを選択することもできます。
最適なサイトパフォーマンスの取得方法
中国のサイト (中国語または英語) を選択することで、最適なサイトパフォーマンスが得られます。その他の国の 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 (한국어)