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

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!

[the cyclist]

Can you upload the input file Sellmeier Du goat2.txt so that we can run your code?

[Torsten]

The matrix A is missing.

Where do you test

Phi*eigenvectors-eigenvectors*eigenvalues = 0

in your code ?

[Hugo]

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.

Answer 1

[John D'Errico]

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.

[Hugo]

Thanks, that explanation helps a lot, so if I understood correctly the overall results when compared to the scale of my Phi Matrix confirm the fact that eigs converged and that my results not matching the physical theories is linked to an error in the design of Phi.

Still one last question, let's imagine "I can" play with the magnitude of Phi's coefficients making it so that I get better scaled values. Would the overall accuracy of the eigs algorithm increase despite a factor between the tolerance and the values in Phi being on the same order ?

[John D'Errico]

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]

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.