4 連立一次方程式(反復法)

実用的なプログラムでは，非常に大きな連立方程式を計算しなくてはならない．数百万元 に及ぶことも珍しくない．これを，ガウス・ジョルダン法で計算するの時間的にほとんど 不可能である．そこで，これよりは格段に計算の速い方法が用いられる．ここでは，その 一つとして反復法を簡単に説明する．

当然ここでも，連立方程式

$$\boldsymbol{A}\boldsymbol{x}=\boldsymbol{b}$$ (31)

を満たす $\boldsymbol{x}$を数値計算により，求めることになる．ここで，真の解 $\boldsymbol{x}$とする．

ここで，ある計算によりn回目で求められたものを $\boldsymbol{x}^(n)$とする．そして，計算回数を増やして，

$$\lim_{n\rightarrow\infty}\boldsymbol{x}^{(n)}=\boldsymbol{x}$$ (32)

になったとする．この様に計算回数を増やして，真の解に近づける方法を反復法という．

この様な方法は，行列 $\boldsymbol{A}$を $\boldsymbol{S}-\boldsymbol{T}$と分解するだけで，容易に作ることが できる．たとえば，

$$\boldsymbol{S}\boldsymbol{x}^{(k+1)}=\boldsymbol{T}\boldsymbol{x}^{(k)}+\boldsymbol{b}$$ (33)

とすればよい．ここで， $\boldsymbol{x}^{(k)}$が $\boldsymbol{\alpha}$に収束するとする．すると，式 (33)と式(31)を比べれば， $\boldsymbol{\alpha}$と $\boldsymbol{x}$は等しいことがわかる．すなわち，式(33)で元の方程式 (31)を表した場合， $\boldsymbol{x}^{(k)}$が収束すれば，必ず真の解 $\boldsymbol{x}$に収束するのである．別の解に収束することはなく，真の解に収束するか，発散 するかのいずれかである．

4.1 ヤコビ法

係数行列 $\boldsymbol{A}$の対角行列を反復計算の行列 $\boldsymbol{S}$としたものがヤコビ(Jacobi)法で ある．係数行列を

$$\left[ \begin{array}{@{\,}ccccc@{\,}} a_{11} & a_{12} & a_{13} & ... ... & \ddots & \vdots \\ a_{n1} & a_{n2} & a_{n3} & \ldots & 0 \end{array} \right]$$ (34)

と分解し，右辺第1項が行列 $\boldsymbol{S}$で第2項が $-\boldsymbol{T}$となる． $\boldsymbol{x}^{(k+1)}$の解の 計算に必要な $\boldsymbol{S}$の逆行列は，それが対角行列なので，

$$\boldsymbol{S}^{-1}= \left[ \begin{array}{@{\,}ccccc@{\,}} a_{11}... ...vdots & \ddots & \vdots \\ 0 & 0 & 0 & \ldots & a_{nn}^{-1} \end{array} \right]$$ (35)

と簡単に求めることができる．$k+1$番目の近似解は $\boldsymbol{x}^{(k+1)}=\boldsymbol{S}^{-1}\left(\boldsymbol{b}+\boldsymbol{T}\boldsymbol{x}^{(k)}\right)$となるなので，容易に 求めることができる．実際，$k$番目の解を

$$x_1^{(k)},\,x_2^{(k)},\,x_3^{(k)},\,\cdots,\,x_n^{(k)}$$

とすると，k+1番目の解は

$$\begin{align*}\begin{aligned}&x_1^{(k+1)}=a_{11}^{-1}\left\{b_1-\left( a_{12}x_2... ...^{(k)}+ \cdots+a_{nn-1}x_{n-1}^{(k)}\right)\right\} \\ \end{aligned}\end{align*}$$ (36)

と計算できる．これを繰り返して連立方程式の解を求める方法が，ヤコビ法である．

4.2 ガウス・ザイデル法

ヤコビ法の特徴では， $\boldsymbol{x}^{(k+1)}$の近似値は，すべてその前の値 $\boldsymbol{x}^{(k)}$を 使うことにある．大きな行列を扱う場合，全ての $\boldsymbol{x}^{(k+1)}$と $\boldsymbol{x}^{(k)}$を記 憶する必要があり，大きなメモリーが必要となり問題が生じる．今では，個人で大きなメ モリーを使い計算することは許されるが，ちょっと前まではできるだけメモリーを節約し たプログラムを書かなくてはならなかった．

そこで， $\boldsymbol{x}^{(k+1)}$の各成分の計算が終わると，それを直ちに使うことが考えば， メモリーは半分で済む．即ち， $x_i^{(k+1)}$を計算するときに，

$$x_i^{(k+1)} =a_{ii}^{-1}\left\{b_i-( a_{i1}x_1^{(k+1)}+a_{i2}x_2^{(k+1)}+a_{i3}x_3^{(k+1)}+\cdots+ a_{ii-1}x_{i-1}^{(k+1)}+\right.$$

$$\qquad\qquad\qquad\qquad\qquad\qquad \left. a_{ii+1}x_{i+1}^{(k)}+a_{ii+2}x_{i+2}^{(k)}+a_{ii+3}x_{i+3}^{(k)}+\cdots+ a_{in}x_n^{(k)})\right\}$$ (37)

とするのである．実際の計算では，$k+1$番目の解は

$$\begin{align*}\begin{aligned}&x_1^{(k+1)}=a_{11}^{-1}\left\{b_1-\left( a_{12}x_2... ...k+1)}+\cdots+a_{nn-1}x_{n-1}^{(k+1)}\right)\right\} \\ \end{aligned}\end{align*}$$ (38)

と計算できる．これを繰り返して連立方程式の解を求める方法が，ガウス・ザイデル法である．

4.3 プログラム例(ガウス・ザイデル法)

ガウス・ザイデル法のような反復法は大きな連立方程式の計算に敵している．しかし，こ こではその計算原理を分かり易くするため，次の連立方程式を計算する．

$$\begin{bmatrix}3 & 2 & 1 \\ 1 & 4 & 1 \\ 2 & 2 & 5 \end{bmatrix} ... ...x}x_1 \\ x_2 \\ x_3 \end{bmatrix} = \begin{bmatrix}10 \\ 12 \\ 21 \end{bmatrix}$$ (39)

式(38)の漸化式に従うプログラム例を以下に示す．

#include <stdio.h>
#include <math.h>
#define N (3)               // 連立方程式の大きさ
#define EPS (1e-15)         // 計算誤差の許容値

int main(void){
  double a[N+1][N+1], x[N+1], b[N+1];
  double dx, absx, sum, new;
  int i,j;

  a[1][1]=3.0;  a[1][2]=2.0;  a[1][3]=1.0;  // 係数行列
  a[2][1]=1.0;  a[2][2]=4.0;  a[2][3]=1.0;
  a[3][1]=2.0;  a[3][2]=2.0;  a[3][3]=5.0;

  b[1]=10.0;                         // 同次項
  b[2]=12.0;
  b[3]=21.0;

  x[1]=0.0;                          // 近似解の初期値
  x[2]=0.0;
  x[3]=0.0;

  do{                                // 反復計算のループ
    dx=0.0;
    absx=0.0;

    for(i=1;i<=N;i++){
      sum=0;
      for(j=1;j<=N;j++){
	if(i != j){
	  sum+=a[i][j]*x[j];
	}
      }

      new=1.0/a[i][i]*(b[i]-sum);   // 反復計算後の近似解
      dx+=fabs(new-x[i]);           // 近似解の変化量を加算
      absx+=fabs(new);              // 近似解の総和計算
      x[i]=new;                     // 新しい近似解を代入
    }
   
  }while(dx/absx > EPS);            // 計算終了条件
 
  for(i=1;i<=N;i++){
    printf("x[%d]=%25.20f\n",i,x[i]);
  }

  return 0;
}

ホームページ: Yamamoto's laboratory
著者: 山本昌志