What is matlab doing under the hood when I solve this generalized eigenvalue problem?

8 次查看（过去 30 天）
David Cyncynates
David Cyncynates 2022-9-5
评论： David Cyncynates 2022-9-7
Ran in:
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');
Warning: First input matrix is close to singular or badly scaled (RCOND = 2.060539e-16) and results may be inaccurate. Consider specifying a small nonzero numeric sigma value instead of 'smallestabs' to improve the condition of the matrix.
disp(D)
-5.1680e-12 + 5.8500e-21i

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

采纳的回答

Christine Tobler
Christine Tobler 2022-9-7
Ran in:
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);
=== Generalized eigenvalue problem A*x = lambda*B*x === The eigenvalue problem is complex non-Hermitian. Matrix B is not Hermitian positive (semi-)definite. Computing 1 eigenvalues of type 'smallestabs'. Parameters passed to Krylov-Schur method: Maximum number of iterations: 300 Tolerance: 1e-14 Subspace Dimension: 20 Find eigenvalues of A\(B*x) = mu*x. Compute decomposition of A...
Warning: First input matrix is close to singular or badly scaled (RCOND = 2.060539e-16) and results may be inaccurate. Consider specifying a small nonzero numeric sigma value instead of 'smallestabs' to improve the condition of the matrix.
--- Start of Krylov-Schur method --- Iteration 1: 1 of 1 eigenvalues converged.
D
D = -5.1680e-12 + 5.8500e-21i
norm(A*V - B*V*D)
ans = 1.5547e-06
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 个）

类别

在 Help Center 和 File Exchange 中查找有关 Large Files and Big Data 的更多信息

标签

产品


版本

R2021b

Community Treasure Hunt

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

Start Hunting!

Translated by