連立方程式の解法 反復法１ ヤコビ法

連立方程式の解法３　反復法１　ヤコビ法

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

　　

①をx、②をy、③をzについて解くと次のようになる。

　　

ここで、として上の式の右辺に代入すると、

　　

となる。

このを用いて

　　

さらに、この値を用いて

このように繰り返し計算すると、

　　

と、やがて、３元連立方程式の解である(x,y,z)=(1,2,3)に収束してゆく。

このように反復計算をして連立方程式の解を求める方法をヤコビ法という。

より一般の

　　

つまり、

　　

の場合、次のように計算すればよい。

STEP1　１行からn行まで

　　

を計算する。

STEP2　i=1〜nまでのすべてのをの値に更新する

STEP3　所定の精度に達していなければSTEP1に戻る。達していたら、計算を打ち切る。

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


void Jacovi(double a[][N], double x[], double b[], int n) {
    double xn[n];
    double eps = 1.e-6,res, res_max;
    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);
            xn[i]= x[i]+res/a[i][i];
        }
        if (k%5==0) { // ５回に１度、計算途中の結果を表示
            printf("%2d %10.6f %10.6f %10.6f %10.6f\n", k, xn[0],xn[1], xn[2], res_max);
        }
        for (i=0; i<n; i++)
            x[i]=xn[i]; //xの更新
        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;
   
    Jacovi(a,x,b,n);

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

いま述べたヤコビ法は、直接、連立方程式の解を求める掃き出し法やガウス消去法とは異なり、連立方程式が唯一の解を持っていたとしても、計算が収束しなかったり、発散する場合がある。

優対角条件、つまり、

などの条件が成立しないとき、ヤコビ法などの反復法は振動したり、発散する。

このことは

　　

の②と③を入れ換え

　　

をヤコビ法で解いてみると、このことが確かめられる。