What is matlab doing under the hood when I solve this generalized eigenvalue problem?
조회 수: 7 (최근 30일)
이전 댓글 표시
I have the following code for a non-hermitian generalized eigenvalue problem that runs extremely fast on matlab. This seems like it should be a hard problem, given its size and lack of nice properties such as symmetry and hermiticity. Can someone illucidate to me what's going on here to make this work so nicely?
dx = 0.0025;
xMin = -100;
xMax = 20;
a = 0.9;
mu = 0.01;
n = 1;
rp = 1 + sqrt(1 - a^2);
rm = 1 - sqrt(1 - a^2);
omegac = a / (2 * rp);
omegan = mu * (1 - mu^2 / (2 * (n + 2)^2));
x = (xMin : dx : xMax)';
xLeft = x + dx;
xRight = x - dx;
r = rp * (exp(x) + 1);
rLeft = rp * (exp(xLeft ) + 1);
rRight = rp * (exp(xRight) + 1);
N = length(x);
e = ones(N,1);
HElements = ((r - rp)./(r - rm));
HElementsLeft = ((rLeft - rp)./(rLeft - rm));
HElementsRight = ((rRight - rp)./(rRight - rm));
Vt0Elements = (a^2)./((r - rm).^2) - HElements .*(mu^2 * r.^2 + 2);
Vt1Elements = -(4 * a * r)./((r - rm).^2);
Vt2Elements = ((r.^2 + a^2).^2)./((r - rm).^2) - HElements * a^2;
B0 = 1 + dx * omegac * 2 * 1i * (rp/(rp - rm));
B1 = - dx * 2 * 1i * (rp/(rp - rm));
D2Left = (1/(dx^2) - HElementsLeft /(2 * dx));
D2Middle = (Vt0Elements - 2 /(dx^2));
D2Right = (1/(dx^2) + HElementsRight/(2 * dx));
DLeft = D2Left;
DMiddle = D2Middle;
DRight = D2Right;
D2 = spdiags([DLeft DMiddle DRight] ,-1:1,N,N);
D2(1,1) = 0;
D2(1,2) = 0;
Vt1 = spdiags(Vt1Elements,0,N,N);
Vt1(1,1) = 0;
Vt2 = spdiags(Vt2Elements,0,N,N);
Vt2(1,1) = 0;
Id = speye(N,N);
Z = sparse(N,N);
S11 = Z;
S12 = Z;
S11(1,1) = 1;
S12(1,2) = 1;
A11 = S11 * B0 - S12 + D2;
A12 = S11 * B1 + Vt1 + omegan * Vt2;
A21 = Id * omegan;
A22 = Id * (-1);
A = [A11 A12;A21 A22];
B11 = Z;
B12 = Vt2 * (-1);
B21 = Id * (-1);
B22 = Z;
B = [B11 B12;B21 B22];
[V,D] = eigs(A,B,1,'smallestabs');
disp(D)
댓글 수: 0
채택된 답변
Christine Tobler
2022년 9월 7일
You can use the Display option to get some more information on what's going on inside of eigs.
load matrices.mat
[V,D] = eigs(A,B,1,"smallestabs", Display=true);
D
norm(A*V - B*V*D)
So for the case of a non-hermitian generalized eigenvalue problem A*x = lambda*B*x, eigs iteratively solves the simple eigenvalue problems A\(B*x) = mu*x.
This works well unless A is close to singular. Unfortunately, that's the case for your matrix A, and it looks like the shifted matrix A - sigma*B is also singular for any value of sigma, so there's not much to do about this unless you have some way to project out the singularity of this system.
However, based on the residual the found eigenvalue and eigenvector are quite accurate.
For cases where eigs gives a warning about the first input matrix being singular, I would recommended
a) trying to choose a shift sigma for which A is not singular (not possible here)
b) making sure to double-check the output to see the effect of this ill-conditioning on the computed result.
추가 답변 (0개)
참고 항목
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!