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

C言語コード

この記事は後編です

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

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

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

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

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

↓前回の記事はこちら

流れの確認

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

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

ここから

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

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

コードの実装

行列形式表示

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

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

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

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

三角化(ガウスの消去法)

三角化とは?

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

ピボット選択の関数

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

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

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

  • 列で最も絶対値が大きい要素を探す
  • その要素が対角上に来るように行の入れ替え
  • それより下の列の要素の消去(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-1番目の要素だけに操作をします。

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);
        }
    }
}

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

外側の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関数も実装して動かしてみました。

おk

全体のコード

//問題生成
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を最大値で割った後2を掛けてそこから1を引いた
                fprintf(fp,"%f ",r);
            }
            fprintf(fp,"\n");
        }
        fclose(fp);
    }
}

//2次元配列の中身を表示する
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;
}


この全体のコード内で使われている関数で、
今回の記事で解説してないものはこちらで解説してます

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

今回はここまで
ではまた〜

コメント

タイトルとURLをコピーしました