计算科学 - 八度音阶的 eigs 例程 - 吾爱随笔录

我正在使用八度音阶并观察到非对称矩阵的 eigs 例程存在问题。使用 GNU octave 版本 3.8.1，下面的代码给出了特征值的显着差异，尽管执行了相同的命令。通过减去eig- 和 -特征值来计算特征值的差异eigs。

% solves the Eigenvalue Problem 
%   -y''= lambda^2 y on the interval (a,b)
%
% with mixed boundary conditions
%   -- dirichlet boundary condition y(a)=0
%   -- robin boundary condition   y'(b)=y(b)
clear all;
close all;
p=6;
%% number of points
n=2^p;
%% grid space
h=1/2^p;
%% numbers of eigenvalues to find
kmax=8;
%% assemble symmetric 2nd finite difference matrix
Lh=-gallery('tridiag',n,1,-2,1)./(h^2);
%% apply robin boundary in last line
Lh(n,:)=0.0;
Lh(n,n-1)=-2/h^2;
Lh(n,n)=-1/h^2*(2*h-2);
%% options for eigs
% opts.tol=1e-10;
opts.p=max(floor(n/2),2*kmax);
comp=0;
while comp<10 % do the same thing 10 times
    %% use eig to get all eigenvalues
    [~,LambdaEig]=eig(Lh);
    if isreal(diag(LambdaEig))
      LambdaEig=sort(diag(LambdaEig));
      lambdaEig=LambdaEig(1:kmax);
    else
      warning('eigenvalues of eig are not real')
    end
    %% use eigs routine to get 1st kmax eigenvalues of smallest magnitude
    [~,LambdaEigs,fl]=eigs(Lh,kmax,'sr');
    if isreal(diag(LambdaEigs))
      LambdaEig=sort(diag(LambdaEig));
      lambdaEigs=sort(diag(LambdaEigs));
    else
      warning('eigenvalues of eigs are not real')
    end
    errmax=max(abs(lambdaEig-lambdaEigs));
    disp(['eigs flag: ' num2str(fl) ', maximum error of eigenvalues: '...
            num2str(errmax)])
    comp=comp+1;
end

你对此有什么解释吗？GNU octave 的手册报告eigs基于 ARPACK 包。我猜 ARPACK 被广泛用于大特征值问题并且经过良好测试。但是 GNU octave 的eigs-routine 似乎对于这个问题并不可靠。

我想在 GNU 八度音程中使用eigs的较大一般特征值问题由于矩阵大小，Krylow 子空间方法对我来说是合理的，而且据我所知是修改后的 Arnoldi 迭代。 $A x=\lambda B x$ $A\in\mathbb{R}^{2000\times2000}$ eigs

您对如何继续有任何建议？