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

5 次查看（过去 30 天）
Khalil Sabour
Khalil Sabour 2024-9-2
编辑： Walter Roberson 2024-9-3
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 个）

Christine Tobler
Christine Tobler 2024-9-3
编辑：Christine Tobler 2024-9-3
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.

类别

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

标签

产品


版本

R2020a

Community Treasure Hunt

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

Start Hunting!

Translated by