反復法による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 $$
  • D:対角成分
  • E:対角以下の部分
  • F:対角以上の部分

1 準備

行列を作成します。

In [17]:
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の命令で計算する。

In [ ]:
sol_x=A\b
In [16]:
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

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

2.1 ヤコビ方法

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

反復計算

  • 初期値ベクトル $x^{0}$
  • 反復計算式 $$ x^{(k)} = M^{-1} (b - Nx^{(k-1)}), ~~ (k=1,2, \cdots) $$

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ノルムを計算します。

In [7]:
% Jacobi method
M=D; N=A-M;
In [8]:
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 
In [9]:
hold off; %現在のグラウをクリアする
plot(1:k,log(errlist), '-o');
grid on
legend('Jacobi method');
Gnuplot Produced by GNUPLOT 5.0 patchlevel 3 -5 -4 -3 -2 -1 0 0 5 10 15 20 Jacobi method Jacobi method
In [15]:
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

演習1

初期値 $x=0$からスタートして、ベクトルの成分毎の計算で考える反復計算式を行い、計算の結果を確認しなさい。

  • 毎回の反復計算で得られる近似解の誤差をerrlistに格納して、計算後に、k - log(errlist(k)) の関係図を描いてください。
  • 収束のオーダーと収束の定数を計算しなさい。

復習

$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 | ) } $$
In [45]:

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回行い、計算の結果を確認しなさい。

  • 毎回の反復計算で得られる近似解の誤差をerrlistに格納して、計算後に、k - log(errlist(k)) の関係図を描いてください。
  • 収束のオーダーと収束の定数を計算しなさい。

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の値を使います。

  • 毎回の反復計算で得られる近似解の誤差をerrlistに格納して、計算後に、k - log(errlist(k)) の関係図を描いてください。
  • 収束のオーダーと収束の定数を計算しなさい。
  • w=0.5, w=1.1, w=1.8の値により、もっと収束の速いwを試してみてください。

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

In [ ]:

演習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を使用してください。

In [ ]: