Help us understand the problem. What is going on with this article?

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;
Why do not you register as a user and use Qiita more conveniently?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
Comments
Sign up for free and join this conversation.
If you already have a Qiita account
Why do not you register as a user and use Qiita more conveniently?
You need to log in to use this function. Qiita can be used more conveniently after logging in.
You seem to be reading articles frequently this month. Qiita can be used more conveniently after logging in.
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away