ガウスの消去法

ガウスの消去法とは、連立1次方程式を解く際によく使われる方法の1つです。ガウスの消去法では、連立1次方程式にある道の係数と、右辺の値を行列で並べ、ある操作を加えて三角行列の形にして、解を求めます。

未知数 ｎ 個を含む、ｎ 個の方程式からなる連立１次方程式
の係数と右辺の値を並べた行列を
と書くとわかりやすいです。しかし、ここでは一部の表記方法を変えて
と書きます。このように書いたほうがプログラムを書きやすいです。

ここらはこの行列を三角行列に変えていく方法を説明します。まず第1行をで割り、第2行から第1行を倍したものを引きます。以後第3行以降も同様の計算を行い、未知数 の係数を第1行以外全て0にします。一般化すると第 k 行を で割り、第 i 行( i > k ) から第 k 行を 倍したものを引くという事です。そうすると、行列は次のような形になります。ただし、それぞれの行列の成分は計算されて先ほどの行列の成分とは違う値に変わっている事に注意してください。
よって、
となり、さらに
と続けて、まで解を求めることができます。

しかし、この計算方法で、もしが0だと、発散してしまいます。が非常に小さい値の場合でもコンピューターは非常に大きな値を扱う事になり、困るケースとなる事もあります。そこで、
とおいて、もし p = 0 だったら、第 k 行から第 n 行までの間のどれかの行と交換して p = 0 とならないようにします。これを部分ピボット選択法と呼びます。また、さらに場合によって列交換まで行う方法を完全ピボット選択法と言います。

このプログラムでは、ガウスの消去法の部分ピボットを用いています。一応任意の ｎ 個の未知数と方程式から成る連立方程式なら解ける事になっていますが、最小限の事しか行わないプログラムなので、全ての問題が解けるわけではありません。問題によっては、誤差が大きくなったり、三角行列にできなかったりして、解くことができない場合があります。

このプログラムで解く場合、解きたい連立1次方程式の各係数を配列 a[N][N+1] に入れます。それと未知数の個数 N の設定もする必要があります。

#include <stdio.h>
#include <math.h>

#define N 4                     //未知数の数
/*計算結果を表示する際に、三角行列になったか
確認したいときは1を、解のみ表示した場合は0にしておく*/
#define CHECK 1 

void gauss(double a[N][N+1])
{
        int i,j,k,l,pivot;
        double x[N];
        double p,q,m,b[1][N+1];

        for(i=0;i<N;i++) {
                m=0;
                pivot=i;

                for(l=i;l<N;l++) {
/*i列の中で一番値が大きい行を選ぶ*/
                        if(fabs(a[l][i])>m) {
                                m=fabs(a[l][i]);
                                pivot=l;
                        }
                }

/*pivotがiと違えば、行の入れ替え*/
                if(pivot!=i) {  
                        for(j=0;j<N+1;j++) {
                                b[0][j]=a[i][j];        
                                a[i][j]=a[pivot][j];
                                a[pivot][j]=b[0][j];
                        }
                }
        }

        for(k=0;k<N;k++) {
                p=a[k][k];              //対角要素を保存
/*対角要素は1になることがわかっているので直接代入*/
                a[k][k]=1;      

                for(j=k+1;j<N+1;j++) {
                        a[k][j]/=p;
                }

                for(i=k+1;i<N;i++) {
                        q=a[i][k];

                        for(j=k+1;j<N+1;j++) {
                                a[i][j]-=q*a[k][j];
                        }
/*0となることがわかっているので直接代入*/
                a[i][k]=0;
                }
        }

/*解の計算*/
        for(i=N-1;i>=0;i--) {
                x[i]=a[i][N];
                for(j=N-1;j>i;j--) {
                        x[i]-=a[i][j]*x[j];
                }
        }

/*行列が最後どうなったか見たいときに実行*/
#if CHECK==1
        for(i=0;i<N;i++) {
                for(j=0;j<N+1;j++) {
                        printf("%10.3f",a[i][j]);
                }
                printf("\n");
                
        }
#endif

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

}

int main(void)
{
/*テストデータ用(プログラムが正常に動いているかの確認用)
3つの解1,2,3がでる。未知数の個数Nは3にする必要あり*/
        /*double a[N][N+1]={{4.0,1.0,1.0,9.0},
                                        {1.0,3.0,1.0,10.0},
                                        {2.0,1.0,5.0,19.0}
                                        };*/
        double a[N][N+1]={{1.0,1.0,-3.0,-4.0,-1.0},
        {2.0,1.0,5.0,1.0,5.0},
        {3.0,6.0,-2.0,1.0,8.0},
        {2.0,2.0,2.0,-3.0,2.0}};

/*解きたい連立方程式から行列をつくり配列a[N][N+1]に入れ、
関数gaussに渡す*/
        gauss(a);

        return 0;
}

上記のプログラムでCHECKを1にした場合は、三角行列となった行列の値の一覧と解が表示されます。CHECKを0にした場合は解のみ表示されます。

     1.000     2.000    -0.667     0.333     2.667
     0.000     1.000    -1.667     1.833     1.667
     0.000     0.000     1.000     4.375     3.500
     0.000     0.000     0.000     1.000     0.800
解は
2.000000
0.200000
0.000000
0.800000
Press any key to continue