I can not say wether Eigs is giving expected results or not

56 次查看（过去 30 天）
Hugo
Hugo 2025-9-4，17:19
评论： Hugo 2025-9-5，8:35
Hello,
I am building a little mode solver for optical fibers based on this article : "Jesper Riishede et al 2003 J. Opt. A: Pure Appl. Opt. 5 534". As explained in the article I need to get the eigenvalues and eigenvectors of my matrix Phi. Since it is a large sized one I turned Phi into a sparse matrix and used eigs on it. But when I compute a simple "Phi*eigenvectors-eigenvectors*eigenvalues" I do not get small values <<1. The flag says eigs converged toward the solution without issue? I know based on what I expect to get from the simulation that the result is completely false but I do not know yet what to blame in my code.
So here is my question : Is there a scaling factor in the Phi*eigenvectors-eigenvectors*eigenvalues based on how big the initial values were, and eigs indeed converged or did it simply failed to reach an accurate solution? Or is there another issue?
Here is a sample of my code (the loop does not make sense but I had to break it to look for where I might have gone wrong a little bit faster since it fails on every iteration).
%%Function topology for DCF design
format long
%%Define the space for the study-----------------------------------
sizeX=25e-6;
sizeY=25e-6;
l=100;
m=100;
d=sizeX/l;
x=-sizeX/2:d:sizeX/2;
y=-sizeY/2:d:sizeY/2;
[xx,yy]=meshgrid(x,y);
n=zeros(size(xx));
Phi=zeros(m*l);
dPhi=sizeX/(l*m);
xX=-sizeX/2:dPhi:sizeX/2;
yY=-sizeY/2:dPhi:sizeY/2;
%%%%%%%%%%%%%%%%%%%%%%%%Iterate over 3 wl%%%%%%%%%%%%%%%%%%%%%%%%%%
%Basic physical parameters-----------------------------------------
A=readmatrix("C:\Users\Tessonneau Hugo\Desktop\Sellmeier Du goat2.txt");
lambdac=1030e-9;
dlambda=0.5e-9;
rcoreIni=5e-6;
rclad=125e-6;
nbg=1;
neff=[];
lambdac0=1535e-9;
lambdacmax=1575e-9;
for lambda0=lambdac0:5e-9:lambdac0%lambdac-dlambda:dlambda:lambdac+dlambda
n_SiO2=interp1(A(:,1)*1e-9,A(:,9),lambda0,'spline');
k0=2*pi/lambda0;
lambda0
%Initial design guess----------------------------------------------
nmin=n_SiO2;
nmax=sqrt(0.14^2+n_SiO2^2);
gammaicore=1;
gammaiclad=0.25;
% nicore=nmin+gammaicore*(nmax-nmin);
% niclad=nmin+gammaiclad*(nmax-nmin);
nicore=n_SiO2+0.38;
niclad=n_SiO2;
n=n*0+nbg;
n(xx.^2+yy.^2<rclad^2)=niclad;
n(xx.^2+yy.^2<rcoreIni^2)=nicore;
% figure(1);
% surf(xx,yy,n,'EdgeColor','none')
%Define Phi based on Helmoltz equation---------------------------------
% N=zeros(size(Phi));
for ii=1:1:l*m
for jj=1:1:l*m
% N(ii,jj)=n(floor(ii/l)+1,floor(jj/m)+1);
if(ii==jj)
Phi(ii,jj)=-4/dPhi^2+n(floor(ii/l)+1,floor(jj/m)+1)^2*k0^2;
elseif(ii==jj+1 || ii==jj-1 || jj==ii-1 || jj==ii+1)
Phi(ii,jj)=1/dPhi^2;
end
end
end
%%Eigenvalues calc--------------------------------------------------
Atest=sparse(rand(l*m));
Phi=sparse(Phi);
[eigenvector,eigenvalues,flag]=eigs(Phi,100,'largestabs','MaxIterations',1000000);
beta2=diag(eigenvalues);
Psi=eigenvector;
neff=[neff sqrt(beta2)./k0];
end
Thanks!
  3 个评论
显示 1更早的评论
Torsten
Torsten 2025-9-4，19:51
The matrix A is missing.
Where do you test
Phi*eigenvectors-eigenvectors*eigenvalues = 0
in your code ?
Hugo
Hugo 2025-9-4，22:05
Sorry, I forgot to add it, I do not have it with me right now but since it is a list of evolving refractive index based on the wavelength and my issue can be studied under only one wavelength you can just get rid of this issue using :
n_SiO2 = 1.45;
That matches my data. As for the following line :
Phi*eigenvectors-eigenvectors*eigenvalues = 0
I tested outside of the main code, directly in matlab command window, sorry for that.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

John D'Errico
John D'Errico 2025-9-4，20:03
编辑：John D'Errico 2025-9-4，20:14
Ran in:
It is literally trivial to get large differences there. We need not even use eigs. Eig will suffice.
A0 = randn(6); A0 = A0 + A0'; % I'll be nice and make the eigenvalues at least real.
[V0,D0] = eig(A0);
I'm sure you would agree that this problem is easy to solve using eig.
A0*V0-V0*D0
ans = 6×6
1.0e-14 * -0.0111 0.0222 -0.0517 0.0139 0.0444 0.0389 0.1776 0.0666 0.1034 -0.0333 -0.1110 -0.0222 -0.0222 0 -0.0645 0 0.0860 -0.0555 0.0444 0.0222 0.1017 -0.0222 -0.2665 0.1665 0.0083 0.1193 0.0534 0.0416 0.0555 0.2665 0.1776 0.0042 -0.0066 0.0222 -0.1554 0
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
And so we see tiny numbers for a difference.
The probem is, you have numbers with wide ranges of magnitudes. And that suggests it is very likely that your problem is not terribly well-scaled as is my A0. As such, I'll modify A0 by just a very simple scale factor.
A1 = A0*1e15;
[V1,D1] = eig(A1);
A1*V1-V1*D1
ans = 6×6
-1.0000 -1.4219 -0.6611 0.0312 0.5000 -0.7500 0 -0.8750 -0.7734 -0.6562 2.1250 -1.4375 -0.2500 -0.4375 0.0039 -0.6250 -0.5938 0.8125 0.7500 -0.7500 -0.0938 0.1562 0.8750 -1.8750 0.8906 -1.6875 0.8730 -0.3438 -1.6250 0.5000 0.5000 -0.2812 -0.4336 0 0.2500 -0.3125
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
Note that now the differences are on the order of 1.
So look carefully at your matrix. Look at the eigenvalues it found.
diag(D1)
ans = 6×1
1.0e+15 * -4.0766 -1.6740 -0.0536 1.2641 3.0689 3.2348
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
You might expect to see differences in there on the order of roughly
eps(sum(abs(diag(D1))))
ans = 2
For the original array, we had
eps(sum(abs(diag(D0))))
ans = 1.7764e-15
As such, you CANNOT just compare that difference to 1, and hope the result will be much smaller than 1. You need to consider the magnitude of what is in the array to choose a "tolerance". (I could arguably have used other similar measures, that arguaby would have also been reasonable. But in any case, 1 is not a reasonable value to compare to.)
In the case of a sparse matrix, you cannot compare to the eigenvalues themselves, but we might as well have used something like this:
tol = eps(sum(abs(A1),'all'))
tol = 4
And that would have given a reasonable comparison value, Probably a little conservative, but not by a lot.
  3 个评论
显示 1更早的评论
John D'Errico
John D'Errico 2025-9-4，23:19
But scaling Phi by some costant does not do anything meaningful. That scaling simply gets applied to the eigenvalues. The eigenvectors would stay the same. In a relative sense, it does nothing to improve anything about the problem.
It might suggest some questions. First, there may be a problem in the code, building the matrices incorrectly. Not like I've never made a mistake in code. :) So this is one thing you could check. Or possibly assumptions you made in the model may have been inadequate. I've not a clue as to what the code models.
Next, that you only got a limited subset of the eigensystem from eigs may be an issue. You got the 100 largest eigenvalues. Was that sufficient? Maybe, maybe not. So you might check to see if more eigenvalues from eigs would help.
Hugo
Hugo 2025-9-5，8:35
Thank you, there is definitely an issue with my array's definition and I will look further into that, what I really wanted to know here was : "what to blame", now I know.

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

产品


版本

R2024a

Community Treasure Hunt

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

Start Hunting!

Translated by