反復法によるAx=bを解く

行列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 $$

1 準備

行列を作成します。

解をMATLABの命令で計算する。

2. 3つの方法で解を求める

2.1 ヤコビ方法

この方法では、$M=D, N=E+F$という$A$の分解を使います。

反復計算

a) 行列の計算で考える反復計算式


x = M\(b-Nx);

b) ベクトルの成分毎の計算で考える反復計算式

注意:ここに書いているアルゴリズムは「疑似コード」です。実際の計算では、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ノルムを計算します。

演習1

初期値 $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 | ) } $$

2.2 Gauss-Seidel法

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

演習2

初期値 $x=0$からスタートして、Gauss-Seidel法の反復計算を20回行い、計算の結果を確認しなさい。

2.3 SOR法

SOR法では以下の行列の分解を使っています。

M=D/w+E; N=(w-1)/w*D+F; (0<w<2)

演習3

初期値 $x=0$からスタートして、SOR法の反復計算を20回行い、計算の結果を確認しなさい。 ただし、wはw=0.5, w=1.1, w=1.8の値を使います。

注意:行列を用いる反復計算、またはベクトルの成分毎の反復計算を利用してください。

演習4:3つの方法の比較

ヤコビ方法、ガウス・ザイデル法、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を使用してください。