How to find the best solution to make eigs function converge?

Hello everyone,

I wrote this code, it is a simple linear eigenvalue problem:

Lx = 50; Nx = 501;                                   
Ly = 50; Ny = 501;
dx = Lx / (Nx - 1);                                     
dy = Ly / (Ny - 1);
x = linspace(-Lx/2, Lx/2, Nx)';                   
y = linspace(-Ly/2, Ly/2, Ny)';
[X, Y] = meshgrid(x, y);

Data = load('mode(049).dat');
Data = reshape(Data(:,1), Ny, Nx);
Pot = load('Potential (r=1.60).dat');
Pot = reshape(Pot(:,1), Ny, Nx);
R = Pot(:);
u = Data(:);
b = 0.8;

% Construct DX2 matrix
diagonals_DX2 = [(1) * ones(Ny * Nx, 1), (-2) * ones(Ny * Nx, 1), (1) * ones(Ny * Nx, 1)];
positions_DX2 = [-Ny, 0, Ny];
DX2 = spdiags(diagonals_DX2, positions_DX2, Ny * Nx, Ny * Nx);

% Construct DY2 matrix
e1 = ones(Ny * Nx, 1);
e2 = ones(Ny * Nx, 1);
e1(Ny:Ny:end) = 0;
e2(Ny+1:Ny:end) = 0;
diagonals_DY2 = [(1) * e1, (-2) * ones(Ny * Nx, 1), (1) * e2];
positions_DY2 = -1:1;
DY2 = spdiags(diagonals_DY2, positions_DY2, Ny * Nx, Ny * Nx);

H1 = 1/2 * 1i * ( DX2/dx^2 + DY2/dy^2 ) + 1i * spdiags(R,0,Ny*Nx,Ny*Nx) - ...
     1i * b * speye(Ny * Nx) + 2 * 1i * spdiags(abs(u).^2, 0, Ny * Nx, Ny * Nx);
H2 =  + 1i * spdiags(u.^2, 0, Ny * Nx, Ny * Nx);
H3 = - 1/2 * 1i * ( DX2/dx^2 + DY2/dy^2 ) - 1i * spdiags(R,0,Ny*Nx,Ny*Nx) + ...
     1i * b * speye(Ny * Nx) - 2 * 1i * spdiags(abs(u).^2, 0, Ny * Nx, Ny * Nx);
H4 = - 1i * spdiags(conj(u).^2, 0, Ny * Nx, Ny * Nx);

H = [H1 , H2;...
     H4 , H3];
H = sparse(H);
opts = struct('maxit', 4000);
Delta = eigs(H,1,'lr', opts);

The goal is to find the Largest real of Delta, however every time the eigenvalue does not converge.

I increased the maximum iteration but that didn't help.

I think the problem is because H is too large.

Any suggestions to handle such a problem?

Thank you!

[William Rose]

@Khalil Sabour, please attach file mode(049).data, so others can run your code. Thank you.

[Torsten]

... and Potential (r=1.60).dat

[Khalil Sabour]

https://drive.google.com/file/d/1DGlv6H4UpLSEKSTwqLjHs0LmHz82iqeF/view?usp=sharing

Answer 1

[Christine Tobler]

It seems that the eigenvalues of H are pure imaginary (I'm seeing real parts of magnitude about 1e-14, and the maximum imaginary part is between -400 and 400).

eigs will converge better the more isolated the eigenvalue being searched is from other eigenvalues. With this 'largestreal' search, it's basically not able to decide which eigenvalue to converge to, since all their real parts are identical.

Two first steps I recommend when searching for eigenvalues of a new matrix:

1) Find a few largest eigenvalues by magnitude: For this H, I turned down the tolerance, since each iteration is quite expensive at this matrix size. Here, I got values

e = eigs(H, 3, 'lm', Display=true, Tolerance=1e-4)
e =
    1.0e+02 *
    0.0000 - 4.0076i
    0.0000 + 4.0076i
   -0.0000 + 4.0074i

I'm using display here to get an understand of how close to convergence the algorithm is. That made me realize I should lower the tolerance.

2) Get a rough estimate of the spectrum:

U = randn(size(H, 1), 500);
[U, ~] = qr(U, "econ");
plot(eig(U'*H*U), 'o')

Note none of the values here are likely to match an eigenvalue of H, but they should lie in the same area of the complex plane as the actual eigenvalues. All of these have real values of about 1e-14, which makes it very likely that all eigenvalues of H are pure imaginary.