Filter löschen
Filter löschen

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

4 Ansichten (letzte 30 Tage)
Khalil Sabour
Khalil Sabour am 2 Sep. 2024
Bearbeitet: Walter Roberson am 3 Sep. 2024
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 Kommentare
Torsten
Torsten am 2 Sep. 2024
... and Potential (r=1.60).dat

Melden Sie sich an, um zu kommentieren.

Antworten (1)

Christine Tobler
Christine Tobler am 3 Sep. 2024
Bearbeitet: Christine Tobler am 3 Sep. 2024
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.

Kategorien

Mehr zu Linear Algebra finden Sie in Help Center und File Exchange

Produkte


Version

R2020a

Community Treasure Hunt

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

Start Hunting!

Translated by