[C言語]ガウスの消去法を用いた掃き出し法のプログラムを組んでみた[コード付き]

どうもシュモクザメです。
突然ですが一次方程式ってあるじゃないですか、

大学に入ってからもたくさん解く機会はありますが
クソめんどくさいですよね

特にこのくらい変数が多いのだと

ただ時間がかかるだけだしなぁ..
大学生にもなってこんな計算にやりたくないわ…

てことで今回はプログラミングの力を使ってこの問題を解決しましょう。

流れの確認
準備
コードの実装
実際に動かしてみた
全体のコード

流れの確認

前提として変数の数と同じ数方程式がある場合のみで考えてます

今回の記事ではガウスの消去法を用います。
まず先ほどの連立方程式をもう一度どうぞ

ここから

行列形式表示への変換
係数行列の上三角化
後退代入法による解の導出

という流れでやっていきます。

準備

今回連立方程式はテキストファイルに格納した状態から始めるので、
ファイルを扱う関数を作っときましょう。

問題を生成してファイルに書き込む関数

まずは問題を作ることからいきましょう。
問題を作るって具体的に何をするかっていうと、上の連立方程式の各変数の係数と定数項を決めるってことですね。
係数が決まればそれに対する答えがありますから、それを求めることにします。

係数の決め方ですが、正直なんでも良いのでrand関数を用いて無作為に決めていきます。
ただ流石に範囲がないと煩雑になるので、
係数も定数項もすべて-1~1の間の値をランダムに持たせることにします。
一旦コードです

void Q(int N, char s[NMAX]){
    int i,j;
    FILE *fp;
    
    if((fp = fopen(s,"w"))==NULL)
        printf("ファイルを開けません\n");
    else{
        for(i=0; i<N; i++){
            for(j=0; j<N+1; j++){
                double r = 2*(double)(rand()/(double)RAND_MAX)-1;
                fprintf(fp,"%f ",r);
            }
            fprintf(fp,"\n");
        }
        fclose(fp);
    }
}

仮引数のNは、係数行列がN×Nであるという意味です。
char sは問題を書き込むファイルの名前です

ファイルのオープンはよくある形 if(fp=fopen) の形ですね
fopenはrモードで開きます

１つ目のforループの回数はNですが、
２つ目のforループの回数はN+1です。
これは定数項の部分が増えるからですね。

こんな感じでファイルに書き込んでます⬇︎

またループの最も内側ではfprintfを使って係数や定数項を開いたファイルに書き込んでいます。
このとき-1~1の範囲に乱数を収めるため、

 double r = 2*(double)(rand()/(double)RAND_MAX)-1;

このようにしています。

rand()/(double)RAND_MAXで0~1の範囲にして、
それを二倍して0~2、
そこから-1することで-1~1ということです。

問題が入ったファイルを読み込む関数

上の関数で作った問題を、
main関数で扱うために係数と定数項に分けてそれぞれ配列に格納する関数を作ります。

void SCAN(int N, double A[N][N], double B[N], char s[NMAX]){
    double b = 0;
    FILE *fp;
    
    
    if((fp = fopen(s,"r"))==NULL)
        printf("ファイルを開けません\n");
    else{
        for(int i=0; i<N; i++){
            double a[N];
            for(int j=0; j<N; j++){
                fscanf(fp,"%lf",&a[j]);
                A[i][j] = a[j];
            }
            fscanf(fp,"%lf",&b);
                  B[i] = b;
        }
        fclose(fp);
    }
}

仮引数のA[N][N]が係数行列が入る配列、
B[N]が定数項が入る配列で残りは上の関数と同じです。

fopenをwモードで開きます。

外側のforループは行を表してますね。N回まわります。
次にその内側ですが、
まずN回forループを回して係数をすべてfscanfで配列Aに書き込みます。

内側のループを抜けたあとは外側のループを抜けずに、一度fscanfを行います。これでN+1列目に書かれている定数項を回収できるということです。
なのでここで回収する値は配列Bに書き込みます。

ちなみにファイルの読み書き操作についてはこの記事で詳しく扱っているのでよければ是非参考にしてください！

配列の中身を表示する関数

mainで中身を表示することがあるので作ります。
2次元配列用と1次元配列用で分けます。

void SHOW_2(int N, double A[N][N]){
    int i, j;
    for(i=0; i<N; i++){
        for(j=0; j<N; j++){
            printf("%f ",A[i][j]);
        }
        printf("\n");
    }
}

void SHOW_1(int N, double B[N]){
    int i;
    for(i=0; i<N; i++){
        printf("%f\n",B[i]);
    }
}

これはそのままですね、

てことで準備は完了です。

コードの実装

行列形式表示

まず行列形式表示について。
連立方程式を以下のように行列の計算に直すことを指します

まあこれは具体的な例で、一般化すると

こんな感じ。
Aが係数が入っている係数行列、
ｂが定数項が入ってる定数項ベクトル、
ｘは求める解が入っている解ベクトルです。

上で問題を生成してそれをファイルに格納して、そのファイルを配列に格納するところまでやりました。
つまり上のファイル読み込みの時点で係数行列表示ができています。
その係数行列表示の状態からA、b別々の配列に格納もできているのでここはおkですね。

三角化（ガウスの消去法）

三角化とは？

まず三角化について。
以下のようにNxN係数行列の左下半分をすべて0に変形することを三角化と言います。

ピボット選択の関数

ここからはややこしいので一旦整理してから取り組みます。

係数行列において、要素のことをピボットと呼びます。
ガウスの消去法を調べていただければわかると思うのですが、
三角化の際に対角上に残すピボットはその列で絶対値が最も大きい値を持つ要素にしなければなりません。

そうでないと途中で計算がバグる可能性があります。
なので

列で最も絶対値が大きい要素を探す
その要素が対角上に来るように行の入れ替え
それより下の列の要素の消去（０にする）

に分けて考えます

int PIBOT(int N, double A[N][N], int k){
    double D[N];
    int p = 0;
    
  
    for(int l=k; l<N; l++){
        D[l] = A[l][k]*A[l][k];//絶対値で最大値を探すので二乗する
    }
    
    //最大値を探して保存
    double m = D[0];
    for(int a=1; a<N; a++)
        if(m<D[a]){
            m=D[a];
        }
    
    if(m==0){
          printf("解が一意に定まりません");
          exit(0);//exit関数を用いてきプログラムを終了する
      }
      
    
    for(int b=0; b<N; b++)
        if(m==D[b])
            p = b;//保存した最大値と同じ値をもつ添字を保存する
    
    
    return p;
}

仮引数のAが係数行列、kは操作する列を表してます。
絶対値が一番大きいものを選択するので一旦二乗して一番大きいものを考えます。

二乗したピボット候補を入れる配列Dを用意してfor文ですべてに0を代入。
そのあとに

for(l=k; l<N; l++){
        D[l] = A[l][k]*A[l][k];//絶対値で最大値を探すので二乗する
    }

k番目からN-１番目の要素だけに操作をします。

k番目より前は関係ないですからね。

その後はDの中で最も大きい値をmに代入。

それが0ならランク(階数)が足りないので解が定まりません。
なのでプログラムを終了させます（exit(0)）
（ランクについてわからない人はググって！）

for(int b=0; b<N; b++)
        if(m==D[b])
            p = b;//保存した最大値と同じ値をもつ添字を保存する

この部分は最大値と同じ値を持つ配列Dの要素の添字を走査してます。
それを関数の返り値にします。

選択したピボットの列をk行目と交換する関数

上で返ってきた添字の行を、k行目の列と総入れ替えします。

void CHANGE(int N, double A[N][N], double B[N], int k, int p){
    int i;
    double S;
    double T;
    
    T = B[k];
    B[k] = B[p];
    B[p] = T;
    
    for(i=0; i<N; i++){
        S = A[k][i];
        A[k][i] = A[p][i];
        A[p][i] = S;
    }
    
}

仮引数はAが係数行列、Bが定数項
kは処理を行っている列の番号で、
pが上の関数で返ってきた選択したピボットがある行の番号です。

中身自体はただ入れ替えてるだけです。

列の、k行目のピボットより下を0にする関数

あるピボットを0にするためには、
そのピボットと同じ値で引いてあげればよいですね。

今、上の関数の時点でk列目の選択ピボットはk行目に来ています。
この選択ピボットの行をx倍したものを、0にしたいピボットがある行から引くことを考えます。

0にしたいピボットがi行目にあるとすると

double r = A[i][k]/A[k][k];

このようにまずrを計算します。

i行から、k行にrを掛けたものを引けば求めている操作になりますね。

void DASH(int N, double A[N][N], double B[N], int k, int i){
    int l;
    double r = A[i][k]/A[k][k];
    
    for(l=k; l<N; l++){
        A[i][l] = A[i][l] - r*A[k][l];
    }
    
    B[i] = B[i] - r*B[k];
}

仮引数のkが処理している列の番号、
iは今0にしようとしているピボットの行の番号です。

for文でAの i行 k~N-1列番目の要素に処理を行っています。

Bもしっかり処理してます。

これらの関数を制御する関数

void DELETE(int N, double A[N][N], double B[N]){
    int k;
    int p=0;
    for(k=0; k<N-1; k++){
        p = PIBOT(N,A,k);
        CHANGE(N, A, B, k, p);
        for(int i=k+1; i<N; i++){
            DASH(N, A, B, k, i);
        }
    }
}

３つの関数を制御して三角化を行う関数です。

外側のforで処理する列の番号kを決めて、

ピボット選択
列交換

としたあと、
係数行列Aの　k行　k+1〜N-1列目、
定数項Bの　k行目に
DASH関数の処理を行わせてます。

つまりBには無駄に何回も処理を行わせてしまっていることですね。（まあ結果には影響しないのでとりあえずおkってことで…）

とりあえずこれで、この関数を通せば
三角化が完了したということになります。

解の導出

ここまでくれば簡単です。
最初に述べた通り、後退代入法によって導出します。

具体的には以下の式です。

三角化をしたAとbの要素を代入するだけです

このタイプの計算ならプログラミングで一瞬ですね

void ANS(int N, double A[N][N], double B[N], double X[N]){
    int i,j;
    X[N-1] = B[N-1]/A[N-1][N-1];
    
    for(i=N-2; i>=0; i--){
        
        double S=0;
        for(j=i+1; j<N; j++)
            S = S + A[i][j] * X[j];
        
        X[i] = (B[i]-S)/A[i][i];
    }
    
}

仮引数のXは解を格納する配列です。

公式の通りに計算して解をXに入れていきます。

実際に動かしてみた

適当にmain関数も実装して動かしてみました。

おｋ

全体のコード

//問題生成
void Q(int N, char s[NMAX]){
    int i,j;
    FILE *fp;
    
    if((fp = fopen(s,"w"))==NULL)
        printf("ファイルを開けません\n");
    else{
        for(i=0; i<N; i++){
            for(j=0; j<N+1; j++){
                double r = 2*(double)(rand()/(double)RAND_MAX)-1;//-1~1までにするために、randを最大値で割った後２を掛けてそこから１を引いた
                fprintf(fp,"%f ",r);
            }
            fprintf(fp,"\n");
        }
        fclose(fp);
    }
}

//２次元配列の中身を表示する
void SHOW_2(int N, double A[N][N]){
    int i, j;
    for(i=0; i<N; i++){
        for(j=0; j<N; j++){
            printf("%f ",A[i][j]);
        }
        printf("\n");
    }
}

//1次元配列の中身を表示する
void SHOW_1(int N, double B[N]){
    int i;
    for(i=0; i<N; i++){
        printf("%f\n",B[i]);
    }
}

//ピボットの列の中で最も大きい値を持つ場所の添字を返す
int PIBOT(int N, double A[N][N], int k){
    double D[N];
    int p=0;
    int l;
    
    for(l=0; l<N; l++)
        D[l] = 0;
    
    
    for(l=k; l<N; l++){
        D[l] = A[l][k]*A[l][k];//絶対値で最大値を探すので二乗する
    }
    
    //最大値を探して保存
    double m = D[0];
    for(int a=1; a<N; a++)
        if(m<D[a]){
            m=D[a];
        }
    
    if(m==0){
          printf("解が一意に定まりません");
          exit(0);//exit関数を用いてきプログラムを終了する
      }
      
    
    for(int b=0; b<N; b++)
        if(m==D[b])
            p = b;//保存した最大値と同じ値をもつ添字を保存する
    
    
    return p;
}

//k列において、k行とPIBOTで得た行の交換
void CHANGE(int N, double A[N][N], double B[N], int k, int p){
    int i;
    double S;
    double T;
    
    T = B[k];
    B[k] = B[p];
    B[p] = T;
    
    for(i=0; i<N; i++){
        S = A[k][i];
        A[k][i] = A[p][i];
        A[p][i] = S;
    }
    
}

void DASH(int N, double A[N][N], double B[N], int k, int i){
    int l;
    double r = A[i][k]/A[k][k];
    
    for(l=k; l<N; l++){
        A[i][l] = A[i][l] - r*A[k][l];
    }
    
    B[i] = B[i] - r*B[k];
}


//ガウスの消去法による三角化
void DELETE(int N, double A[N][N], double B[N]){
    int k;
    int p=0;
    for(k=0; k<N-1; k++){
        p = PIBOT(N,A,k);
        CHANGE(N, A, B, k, p);
        for(int i=k+1; i<N; i++){
            DASH(N, A, B, k, i);
        }
    }
}

//ファイルから行列AとBを取り出す
void SCAN(int N, double A[N][N], double B[N], char s[NMAX]){
    double b = 0;
    FILE *fp;
    
    
    if((fp = fopen(s,"r"))==NULL)
        printf("ファイルを開けません\n");
    else{
        for(int i=0; i<N; i++){
            double a[N];
            for(int j=0; j<N; j++){
                fscanf(fp,"%lf",&a[j]);
                A[i][j] = a[j];
            }
            fscanf(fp,"%lf",&b);
                  B[i] = b;
        }
        fclose(fp);
    }
}

//解を公式から求めて配列に格納
void ANS(int N, double A[N][N], double B[N], double X[N]){
    int i,j;
    X[N-1] = B[N-1]/A[N-1][N-1];
    
    for(i=N-2; i>=0; i--){
        
        double S=0;
        for(j=i+1; j<N; j++)
            S = S + A[i][j] * X[j];
        
        X[i] = (B[i]-S)/A[i][i];
    }
    
}

int main(void) {
    int N;
    printf("行列のサイズを教えて下さい:");
    scanf("%d",&N);
     
    double A[N][N];
    double B[N];
    double X[N];
     
    
    char s[NMAX];
    printf("使用するファイルを入力してください:");
    scanf("%s",s);
    
    Q(N,s);
    SCAN(N,A,B,s);
     
    SHOW_2(N,A);
    printf("\n");
    SHOW_1(N,B);
    printf("\n");
    printf("三角化\n");
    DELETE(N,A,B);
     
    SHOW_2(N,A);
    printf("\n");
    SHOW_1(N,B);
    printf("\n");
    
    ANS(N,A,B,X);
    printf("解\n");
    SHOW_1(N,X);
  
    return 0;
}

また今回のガウスの消去法プログラムを応用して最小二乗法をとくプログラムの記事も書いたのでそちらもよければご覧ください