ガウスの消去法

こんばんは

スキちゃんスキちゃんガウスがスキちゃん

先月ガウスの消去法を用いて連立方程式を解きましたが、今回はそのプログラムを改良してきました。

具体的には係数の入力をファイルではなく手入力にし、解の表示をmain関数内で行うようにしました。

また、前回のプログラムの先頭でmath.hをインクルードしていたのですが、あれは必要ありません。

以下ソースコードになります。

#include <stdio.h>
#define N 4                     //未知数の数
void gauss(double a[N][N+1]);

int main(void)
{
    double x[N][N+1],d[N*(N+1)];
    int i,j;
    
    /*方程式の入力*/
    printf("未知数が%d個の連立方程式の係数を入力してください\n",N);
    for( i = 0; i < N; i++ ){
        for( j = 0; j < N+1; j++ ){
            printf("B[%d,%d]=", i + 1,j + 1 );
            scanf("%lf", &x[i][j]);
        }
    }
    
    /*行列xの表示*/
    printf("行列\n");
    for( i = 0; i < N; i++ ){
        for( j = 0; j < (N+1); j++ ){
            printf("%10.5f", x[i][j] );
        }
        printf("\n");
    }
    
    gauss(x);/*行列x[N][N+1]を関数gaussに渡す*/
    
    /*解の表示*/
    printf("解は\n");
    for(i=0;i<N;i++) {
        printf("x%d = %f\n",i+1,x[i][N]);
    }
    return 0;
}

void gauss(double a[N][N+1]){
    int i,j,k;
    
    /*前進消去ステップ*/
    for(k=0;k<N;k++) {
        double p,q;
        
        /*step1*/
        p = a[k][k];         /*対角要素を保存*/
        for(j=k+1;j<N+1;j++) {
            a[k][j] /= p;
        }
        a[k][k] = 1;         /*対角要素は1になることがわかっているので直接代入*/
        
        /*step2*/
        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];
            }
            a[i][k]=0;       /*0となることがわかっているので直接代入*/
        }
    }
    
    /*後退代入ステップ*/
    for(i=N-2;i>=0;i--) {
        for(j=N-1;j>i;j--) {
            a[i][N] -= a[i][j] * a[j][N];
        }
    }
}
 

また改良点があれば治したいと思います。

それではおやすみなさい