前回扱ったべき乗法のアルゴリズムを参考にして逆べき乗法のアルゴリズムを実現するコードを書き、
$$ A= \begin{pmatrix} 2 & -1 & 0 & 0 & 0\\ -1 & 2 & -1 & 0 & 0\\ 0 & -1 & 2 & -1 & 0\\ 0 & 0 & -1 & 2 & -1\\ 0 & 0 & 0 & -1 & 2 \end{pmatrix} $$の絶対値最小固有値および対応する固有ベクトルの近似値を求めてください。ただし、
$$ x^{(0)}= \begin{pmatrix} 1\\ 0\\ 0\\ 0\\ 0 \end{pmatrix} ,\quad EPS=10^{-2},\quad K=100 $$とします。逆行列 $A^{-1}$ の計算には、関数invを使用すること。
A = [2,-1,0,0,0; -1,2,-1,0,0; 0,-1,2,-1,0; 0,0,-1,2,-1; 0,0,0,-1,2];
EPS = 1E-2;
K = 100;
x = [1;0;0;0;0];
x = x/norm(x); %初期値ベクトルの正規化（しなくてもよい）
inv_A = inv(A); %Aの逆行列の計算
%inv_Aに対してべき乗法を実行する
r = x'*inv_A*x;
k = 0;
while norm(inv_A*x-r*x) >= EPS && k <= K
x = inv_A*x/norm(inv_A*x);
r = x'*inv_A*x;
k = k+1;
end
r = 1/r; %固有値の逆数を取る
if k > K
printf("%d回の反復では、収束しなかった。\n",K)
else
printf("%d回の反復で、絶対値最小固有値の近似値%.5f\n",k,r)
printf("および対応する固有ベクトルの近似値\n")
printf("%.5f\n",x)
printf("を求めることができた。")
end
5回の反復で、絶対値最小固有値の近似値0.26795 および対応する固有ベクトルの近似値 0.28993 0.50119 0.57730 0.49881 0.28752 を求めることができた。
ガウスの消去法の前進消去によるLU分解と前進代入・後退代入を利用した逆べき乗法によって、演習1と同じく
$$ A= \begin{pmatrix} 2 & -1 & 0 & 0 & 0\\ -1 & 2 & -1 & 0 & 0\\ 0 & -1 & 2 & -1 & 0\\ 0 & 0 & -1 & 2 & -1\\ 0 & 0 & 0 & -1 & 2 \end{pmatrix} $$の絶対値最小固有値および対応する固有ベクトルの近似値を求めてください。ただし、
$$ x^{(0)}= \begin{pmatrix} 1\\ 0\\ 0\\ 0\\ 0 \end{pmatrix} ,\quad EPS=10^{-2},\quad K=100 $$とします。
次のコードでは、理解の手助けとなるよう、ひとまず連立一次方程式を演算子\で解くことによる逆べき乗法を実装している。
A = [2,-1,0,0,0; -1,2,-1,0,0; 0,-1,2,-1,0; 0,0,-1,2,-1; 0,0,0,-1,2];
EPS = 1E-2;
K = 100;
x = [1;0;0;0;0];
x = x/norm(x); %初期値ベクトルの正規化（しなくてもよい）
y = A\x; %連立一次方程式Ay=xの解yを求める
r = x'*y;
k = 0;
while norm(y-r*x) >= EPS && k <= K
x = y/norm(y);
y = A\x; %連立一次方程式Ay=xの解yを求める
r = x'*y;
k = k+1;
end
r = 1/r; %固有値の逆数を取る
if k > K
printf("%d回の反復では、収束しなかった。\n",K)
else
printf("%d回の反復で、絶対値最小固有値の近似値%.5f\n",k,r)
printf("および対応する固有ベクトルの近似値\n")
printf("%.5f\n",x)
printf("を求めることができた。")
end
5回の反復で、絶対値最小固有値の近似値0.26795 および対応する固有ベクトルの近似値 0.28993 0.50119 0.57730 0.49881 0.28752 を求めることができた。
ガウスの消去法の前進消去による $A$ のLU分解を行っておき、y = A\xの代わりに前進代入と後退代入で連立一次方程式を解けばよいので、解答のコードは次のようになる。
function value = LU_decomposition(A) %与えられた行列のLU分解を行い、結果を一つの行列にして返す関数
n = size(A,1);
for k = 1:n-1
for i = k+1:n
alpha = -A(i,k)/A(k,k);
A(i,k+1:n) = A(i,k+1:n)+alpha*A(k,k+1:n);
A(i,k) = -alpha;
end
end
value = A;
end
function value = fb_substitution(A,b) %LU分解済みの行列とベクトルに対して、前進代入と後退代入により連立一次方程式の解を返す関数
n = size(A,1);
%前進代入（L（Aの対角より下）とbを用いて、Ly=bの解yを求める）
y = zeros(n,1);
y(1) = b(1);
for i = 2:n
y(i) = b(i)-A(i,1:i-1)*y(1:i-1);
end
%後退代入（U（Aの対角以上）とyを用いて、Ux=yの解xを求める）
x = zeros(n,1);
x(n) = y(n)/A(n,n);
for i = n-1:-1:1
x(i) = (y(i)-A(i,i+1:n)*x(i+1:n))/A(i,i);
end
value = x;
end
%---------- ここまでが関数の定義（変数はローカルであることに注意）----------
%---------- 実際に実行されるコードはここから ----------
A = [2,-1,0,0,0; -1,2,-1,0,0; 0,-1,2,-1,0; 0,0,-1,2,-1; 0,0,0,-1,2];
EPS = 1E-2;
K = 100;
x = [1;0;0;0;0];
x = x/norm(x); %初期値ベクトルの正規化（しなくてもよい）
A = LU_decomposition(A); %AのLU分解を一度だけ行う（結果はAにまとめて格納する）
y = fb_substitution(A,x); %前進代入と後退代入により、連立一次方程式Ay=xの解yを求める
r = x'*y;
k = 0;
while norm(y-r*x) >= EPS && k <= K
x = y/norm(y);
y = fb_substitution(A,x); %前進代入と後退代入により、連立一次方程式Ay=xの解yを求める
r = x'*y;
k = k+1;
end
r = 1/r; %固有値の逆数を取る
if k > K
printf("%d回の反復では、収束しなかった。\n",K)
else
printf("%d回の反復で、絶対値最小固有値の近似値%.5f\n",k,r)
printf("および対応する固有ベクトルの近似値\n")
printf("%.5f\n",x)
printf("を求めることができた。")
end
5回の反復で、絶対値最小固有値の近似値0.26795 および対応する固有ベクトルの近似値 0.28993 0.50119 0.57730 0.49881 0.28752 を求めることができた。
逆べき乗法を利用して、
$$ A= \begin{pmatrix} 2 & -1 & 0 & 0 & 0\\ -1 & 2 & -1 & 0 & 0\\ 0 & -1 & 2 & -1 & 0\\ 0 & 0 & -1 & 2 & -1\\ 0 & 0 & 0 & -1 & 2 \end{pmatrix} $$の $\mu=2.3$ に最も近い固有値および対応する固有ベクトルの近似値を求めてください。ただし、
$$ x^{(0)}= \begin{pmatrix} 1\\ 0\\ 0\\ 0\\ 0 \end{pmatrix} ,\quad EPS=10^{-2},\quad K=100 $$とします。関数invを使用して逆行列の計算を行ってもよいです。
A = [2,-1,0,0,0; -1,2,-1,0,0; 0,-1,2,-1,0; 0,0,-1,2,-1; 0,0,0,-1,2];
EPS = 1E-2;
K = 100;
x = [1;0;0;0;0];
mu = 2.3; %指定された値μ
x = x/norm(x); %初期値ベクトルの正規化（しなくてもよい）
n = size(A,1);
B = A-mu*eye(n); %B=A-μInとおく
%Bに対して逆べき乗法を実行する
inv_B = inv(B); %Bの逆行列の計算
%inv_Bに対してべき乗法を実行する
r = x'*inv_B*x;
k = 0;
while norm(inv_B*x-r*x) >= EPS && k <= K
x = inv_B*x/norm(inv_B*x);
r = x'*inv_B*x;
k = k+1;
end
r = 1/r+mu; %固有値の逆数を取り、μだけシフトする
%（元のrはinv_Bの絶対値最大固有値、1/rはBの絶対値最小固有値、1/r+muはAのμに最も近い固有値）
if k > K
printf("%d回の反復では、収束しなかった。\n",K)
else
printf("%d回の反復で、%.5fに最も近い固有値の近似値%.5f\n",k,mu,r)
printf("および対応する固有ベクトルの近似値\n")
printf("%.5f\n",x)
printf("を求めることができた。")
end
8回の反復で、2.30000に最も近い固有値の近似値2.00000 および対応する固有ベクトルの近似値 0.57785 -0.00049 -0.57735 0.00049 0.57685 を求めることができた。