Sakurai-Sugiura法をMATLABで実装してみた

Sakurai-Sugiura法

指定した領域内の固有値を求める解法であるSakurai-Sugiura法を実装してみました．
ムダな部分もあるかと思いますが，とりあえず動きます．
固有値の求解まではコーディングしてますが，固有ベクトルは求めていませんので注意してください．

引用文献

詳細なアルゴリズムや，原理は以下の原著論文をご参考ください．
1. Tetsuya Sakurai, Hiroshi Sugiura, A projection method for generalized eigenvalue problems using numerical integration, Journal of Computational and Applied Mathematics, Volume 159, Issue 1, 2003, Pages 119-128,
2. 日本語 -> https://repository.kulib.kyoto-u.ac.jp/dspace/bitstream/2433/43096/1/1320_24.pdf

MATLABサンプルコード

ss.m

% テスト行列の準備
A = (diag(99:-1:0) + circshift(eye(100),1,2))./100;
B = [zeros(80,80) zeros(80,20);zeros(20,80) eye(20)];
[mA nA] = size(A);

% 適当な値
N = 64;
j = 0:N-1;

% ジョルダン閉曲線の中心
r = 0.015;
% ジョルダン閉曲線の半径
p  = 0.02;

% 上記の領域内の固有値数
m = 4;

% ランダムベクトル
u = rand(mA,1);
v = rand(mA,1);

% (1) ジョルダン閉曲線（今回は円）
w = r + p .* exp(2 * pi * 1i * j ./ N);
plot(real(w),imag(w));
grid on;

% (2) 
for lp = 1:N
    y(:,lp) = v' / (w(lp) * B - A);
end

% (3)
f = u' * y;

% (4)
for lp2 = 1:2 * m
    mu(lp2) = 0;
    for lp = 1:N
        mu(lp2) = mu(lp2) + ( w(lp)-r ).^lp2 * f(lp);
    end
    mu(lp2) = mu(lp2) / N;
end

% (5) Hankel行列の計算と固有値の求解
H1 = hankel(mu(1:m),mu(m:end-1));
H2 = hankel(mu(2:m+1),mu(m+1:end));

t = eig(H2,H1);
lambda = r + t
figure;
scatter(real(lambda),imag(lambda));
hold on;
scatter(real(w),imag(w));
grid on;