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;
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;
Author And Source
この問題について(Sakurai-Sugiura法をMATLABで実装してみた), 我々は、より多くの情報をここで見つけました https://qiita.com/HSND/items/a907446807f815382cd6著者帰属：元の著者の情報は、元のURLに含まれています。著作権は原作者に属する。
Content is automatically searched and collected through network algorithms . If there is a violation . Please contact us . We will adjust (correct author information ,or delete content ) as soon as possible .