How to find the best solution to make eigs function converge?
4 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
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
Antworten (1)
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.
0 Kommentare
Siehe auch
Kategorien
Mehr zu Linear Algebra finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!