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

Question

Khalil Sabour about 23 hours 前

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2149554-how-to-find-the-best-solution-to-make-eigs-function-converge

编辑： Christine Tobler about 4 hours 前

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!

3 个评论
显示 1更早的评论隐藏 1更早的评论

Torsten about 22 hours 前

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

Khalil Sabour about 5 hours 前

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

请先登录，再进行评论。

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

Answer 1

Christine Tobler about 4 hours 前

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2149554-how-to-find-the-best-solution-to-make-eigs-function-converge#answer_1509709

编辑：Christine Tobler about 4 hours 前

在 MATLAB Online 中打开

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.