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

この記事は後編です

前回の記事はこちら

今回はガウスの消去法と掃き出し法の実際の計算部分を扱います
ではやっていきましょう



どういうアルゴリズムか

簡単に説明します。
とりあえず連立方程式を行列式に直します。

上の図は具体的な例です。
一般化して

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

これを三角化を行って

このように変形します。

この状態まで持っていくと、

このように計算すれば解が求まります。(不思議ですね)
これを計算するプログラムを組んでいきます。



コードの実装

前回の記事で問題を生成してそれをファイルに格納して、そのファイルを配列に格納するところまでやりました。
なので行列式は配列に入っている前提で進めていきます。

とりあえずアルゴリズムから、今回の方針は

  • 行列に三角化を行う関数
  • 解を公式から求める関数


が主な部分ですね。
三角化はややこしいので、その中でも関数を分けて実装します。



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

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

係数行列において、要素のことをピボットと呼びます。

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

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

  • 列で最も絶対値が大きい要素を探す
  • その要素が対角上に来るように行の入れ替え
  • それより下の列の要素の消去(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には無駄に何回も処理を行わせてしまっている)

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


解を求める関数

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


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

よかったら御覧ください


今回はここまで


ではまた、

syumokuzame

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です

コメントする