連立方程式の解法 ガウス＝ザイデル法

連立方程式の解法　ガウス＝ザイデル法

§１　ガウス＝ザイデル法

次の連立方程式があるとする。

まず、この連立方程式を、ヤコビ法と同様に次のように変形する。

漸化式の形式で書くと、ヤコビ法は

と反復計算をするのに対し、ガウス＝ザイデル法は、収束を早めるために

と反復計算をする。

前回のヤコビ法同様に、

を反復計算の初期値を与えて計算した結果は次の通り。

１０回程度の反復計算で(x,y,z)=(1,2,3)に達しており、この計算の場合、ヤコビ法と比較すると、ガウス＝ザイデル法は収束速度が約３倍である。

単に収束速度が速いだけではなく、ガウス＝ザイデル法は、計算の収束条件が緩和される。

n元連立方程式

に対して、

ヤコビ法が収束する十分な条件が

であるのに対し。ガウス＝ザイデル法では

である。

ガウス＝ザイデル法の計算法、アルゴリズムは次の通り。

STEP１　i=1〜nの各行に対して次の計算をする

STEP2　収束条件を満たしていなければSTEP１に戻る。満たしていれば計算終了

（１）では、と書いているけれど、プログラムでは

とすればよい。

自動的に値が更新されるので、上のように計算すればよい。

§２　SOR法

①より

このように変形すると、（２）式の括弧の中は、との修正量と考えることができる。

そして、次のように適切な値ωを与えることができれば、ガウス＝ザイデル法より速く収束することが予想される。

このω>1として計算する方法をSOR（Successive Over Relaxion 逐次加速緩和法）と呼ぶ。

経験的にω=1.5くらいにとると速く収束することが知られている。

#include <stdio.h>
#include <math.h>
#define N 100


void Gauss_Seidel(double a[][N], double x[], double b[], int n) {
    double eps = 1.e-6,res=0,res_max;
    double omega = 1.0;// omega>1ならばSOR法
    int i, j, k;
    int limit = 100;

    printf("反復回数 x(1)      x(2)      x(3)      残差\n");
    for (k = 1 ; k <= limit; k++) {
        res_max = 0;
        for (i=0; i < n; i++) {
            res = b[i];
            for (j=0; j < n; j++) {
                res = res - a[i][j]*x[j];
            }
            if(res_max<fabs(res)) res_max = fabs(res);
            x[i]= x[i]+omega*res/a[i][i];
        }
        printf("%2d %10.6f %10.6f %10.6f %10.6f\n", k, x[0],x[1], x[2], res_max);
        if (eps>res_max) break;
    }
}

main() {
    int i, n;
    double a[N][N] = {{3,1,1}, {1,3,1}, {1,1,3}};
    double x[N]={0.};
    double b[N]= {8., 10., 12.};
   
    n=3;
   
    Gauss_Seidel(a,x,b,n);

    for (i = 0; i < n; i++) {
        printf("x[%d] = %f\n", i+1,x[i]);
    }
}