2017年11月2日木曜日

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

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


§1 ガウス=ザイデル法

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

漸化式の形式で書くと、ヤコビ法は
と反復計算をするのに対し、ガウス=ザイデル法は、収束を早めるために
と反復計算をする。

前回のヤコビ法同様に、
を反復計算の初期値を与えて計算した結果は次の通り。
10回程度の反復計算で(x,y,z)=(1,2,3)に達しており、この計算の場合、ヤコビ法と比較すると、ガウス=ザイデル法は収束速度が約3倍である。




単に収束速度が速いだけではなく、ガウス=ザイデル法は、計算の収束条件が緩和される。
n元連立方程式
に対して、
ヤコビ法が収束する十分な条件が
であるのに対し。ガウス=ザイデル法では
である。

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

STEP1 i=1nの各行に対して次の計算をする
STEP2 収束条件を満たしていなければSTEP1に戻る。満たしていれば計算終了

(1)では、と書いているけれど、プログラムでは
とすればよい。
自動的に値が更新されるので、上のように計算すればよい。


§2 SOR

①より
このように変形すると、(2)式の括弧の中は、の修正量と考えることができる。
そして、次のように適切な値ωを与えることができれば、ガウス=ザイデル法より速く収束することが予想される。
このω>1として計算する方法をSORSuccessive 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]);
    }
}

0 件のコメント:

コメントを投稿