行列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) $$
不動点式を解くために、反復計算法はよく使用されています。今回の授業では、3つの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より)
// xから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
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
初期値 $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
ベクトルの成分毎の計算で考える反復計算式
// xの各成分を更新します。(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法の反復計算を20回行い、計算の結果を確認しなさい。
初期値 $x=0$からスタートして、SOR法の反復計算を20回行い、計算の結果を確認しなさい。 ただし、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を使用してください。