行列Aの分解$A=M+N$ (M:正則行列)によって、$Ax=b$という方程式を不動点式に変形することができます。
$$ Ax=b \to (M+N)x=b \to Mx=b-Nx \to x = M^{-1}(b-Nx) $$関数$g(x):=M^{-1}(b-Nx)$を定義します。このとき、$Ax=b$の解がg(x)の不動点となります。即ち、
$$ x=g(x) $$不動点式を解くために、反復計算法はよく使用されています。今回の授業では、３つのAの分解方法によって、$Ax=b$の反復計算方法を検討します。
準備として、Aの成分を以下のように表します。
$$ A = D + E + F $$行列を作成します。
n=4;
A=[2, -1, 0, 0;
-1, 2, -1, 0;
0, -1, 2, -1;
0, 0, -1, 2];
b = ones(n,1); % Vector b can also be inputed by: b=[1,1,1,1]';
A
b
解をMATLABの命令で計算する。
sol_x=A\b
D=diag(diag(A));
E=[0, 0, 0, 0;
-1, 0, 0, 0;
0, -1, 0, 0;
0, 0, -1, 0];
F=A-D-E;
D
E
F
この方法では、$M=D, N=E+F$という$A$の分解を使います。
反復計算
x = M\(b-Nx);
注意：ここに書いているアルゴリズムは「疑似コード」です。実際の計算では、MATLAB（Octave）の言語に書き直すのは必要です。
擬似コード （ぎじコード、英: pseudocode）とは、アルゴリズムなどを、架空の非常に高水準なプログラミング言語（擬似言語）で記述したものである。 Pascal、Fortran、C言語などの既存のプログラミング言語の構文と、自然言語に近い表現を組み合わせて記述することが多い。（Wikipediaより）
// ｘからnew_xを作ります.
for i from 1 to n
//b - Nxを計算します。
tmp = b(i)
for j from 1 to (i-1)
tmp = tmp - A(i,j) * x(j)
next
for j from (i+1) to n
tmp = tmp - A(i,j) * x(j)
next
// 1/A(i,i)*(b-Nx)の計算
new_x(i) = 1/A(i,i) * tmp
next
// xをnew_xにします。
for i from 1 to n
x(i) = new_x(i)
next
以下は行列の計算で書かれる反復計算式を利用しています。
補足： MATLABのnorm関数は与えられるベクトルのL2ノルムを計算します。
% Jacobi method
M=D; N=A-M;
x=zeros(n,1); %初期値ベクトルx.
errlist=[];
k=0;
while k<20
x = M\(b-N*x);
solution_error = norm(x-sol_x)/norm(sol_x);
k = k+1;
errlist(k) = solution_error;
printf("%2d interations: error = %f \n",k, solution_error );
end
1 interations: error = 0.808608 2 interations: error = 0.654129 3 interations: error = 0.529196 4 interations: error = 0.428128 5 interations: error = 0.346362 6 interations: error = 0.280213 7 interations: error = 0.226697 8 interations: error = 0.183402 9 interations: error = 0.148375 10 interations: error = 0.120038 11 interations: error = 0.097113 12 interations: error = 0.078566 13 interations: error = 0.063561 14 interations: error = 0.051422 15 interations: error = 0.041601 16 interations: error = 0.033656 17 interations: error = 0.027228 18 interations: error = 0.022028 19 interations: error = 0.017821 20 interations: error = 0.014418
hold off; %現在のグラウをクリアする
plot(1:k,log(errlist), '-o');
grid on
legend('Jacobi method');
length(errlist)
for k=3:length(errlist)
alpha = log( errlist(k-1)/errlist(k))/log(errlist(k-2)/errlist(k-1));
lambda=errlist(k)/errlist(k-1);
printf("Order: %f, Lambda: %f \n", alpha, lambda);
end
ans = 20
初期値 $x=0$からスタートして、ベクトルの成分毎の計算で考える反復計算式を行い、計算の結果を確認しなさい。
$x$を$Ax=b$の解または精度の良い近似解とする。 反復回数$n$が十分大きい時，反復計算で得られる数列$\{x_i\}$が以下の式を満たすことは考えられる。
$$ \frac{| x_n - x|}{|x_{n-1} - x |^{\alpha}} \approx \frac{| x_{n-1} - x|}{|x_{n-2} - x |^{\alpha}} \approx \lambda $$ここで、$\lambda$は収束の定数である。
よって，$n\ge 2$のとき，次の式で収束オーダー $\alpha$ を計算できる。一般的に，$n$が大きい時，$\alpha$の計算結果が一定の値となる。
$$ \alpha \approx \frac{ \log (| x_{n-1} - x | / | x_n - x | )}{ \log (| x_{n-2} - x | / | x_{n-1} - x | ) } $$
Gauss-Seidel法では以下の行列の分解を使っています。
M=D+E; N=A-M
ベクトルの成分毎の計算で考える反復計算式
// ｘの各成分を更新します。（new_xという新しい変数は不要です。）
for i from 1 to n
//b - Nxを計算します。
tmp = b(i)
for j from 1 to (i-1)
tmp = tmp - A(i,j) * x(j)
next
for j from (i+1) to n
tmp = tmp - A(i,j) * x(j)
next
// 1/A(i,i)*(b-Nx)の計算
x(i) = 1/A(i,i) * tmp
next
初期値 $x=0$からスタートして、Gauss-Seidel法の反復計算を２０回行い、計算の結果を確認しなさい。
初期値 $x=0$からスタートして、SOR法の反復計算を２０回行い、計算の結果を確認しなさい。 ただし、wはw=0.5, w=1.1, w=1.8の値を使います。
注意：行列を用いる反復計算、またはベクトルの成分毎の反復計算を利用してください。
ヤコビ方法、ガウス・ザイデル法、SOR法について、それぞれの計算誤差は error_list_1, error_list_2, error_list_3に格納されています。 error_list_1, error_list_2, error_list_3を同じグラフで描いて比較してください。
特に、SOR法はw=0.5, w=1.1, w=1.8の中で収束の一番速いw、またはw=0.5, w=1.1, w=1.8よりもっと速いwを使用してください。