ガウス•ジョルダン法で連立方程式を解く

プログラム

#include <stdio.h>
#include <stdlib.h>


// 求めたい連立方程式
// 3x1 +2x2 +7x3 +x4 =8
// x1 +5x2 + x3 −x4 =5
// 4x +x +3x −2x =7
// x1 +6x2 +4x3 +3x4 =13

// ↓↓
// x1 x2 x3 x4 = b
// 3  2  7  1  8
// 1  5  1  -1  5
// 4  1  3  -2  7
// 1  6  4  3  13

// 拡大行列 
// 1から4列までは未知数の係数、5列目は定数項
double augmentedMatrix[4][5] = {
    {3, 2, 7, 1, 8},
    {1, 5, 1, -1, 5},
    {4, 1, 3, -2, 7},
    {1, 6, 4, 3, 13}};

//上に同じく
double initialConditions[4][5] = {
    {3, 2, 7, 1, 8},
    {1, 5, 1, -1, 5},
    {4, 1, 3, -2, 7},
    {1, 6, 4, 3, 13}};

double x[4];

//行列表示関数
void printMatrix(double matrix[4][5], int n) {
    int i, j;
    for (i = 0; i < n; i++) {
        printf("|\t");
        for (j = 0; j < n + 1; j++) {
            printf("%f\t", matrix[i][j] + 0);//-0の表示を消すために+0
        }
        printf("|\n");
    }
}

// param: matrix[4][5], n
// matrix[4][5]: 連立方程式の係数行列と定数(?)ｂ
// n: 連立方程式の数
void gauss_jordan(double matrix[4][5], int n) {
    int i, j, k;
    double pivot;
    //階段行列を作る
    for (i = 0; i < n; i++) {//行
        pivot = matrix[i][i];//対角成分をピボットにする
        for (j = 0; j < n + 1; j++) {
            matrix[i][j] = matrix[i][j] / pivot; //対角成分を1にする
            // printf("Pivot .. matrix[%d][%d] = %f\n", j, k, matrix[j][k]);
        }
        for (j = 0; j < n; j++) {
            if (i != j) {
                pivot = matrix[j][i];
                // printf("pivot : %f\n", pivot);
                for (k = 0; k < n + 1; k++) {
                    matrix[j][k] = matrix[j][k] - pivot * matrix[i][k];
                    // printf("matrix[%d][%d] = %f\n", j, k, matrix[j][k]);
                }
            }
        }
    }
    //ガウスジョルダン法を適応した後の拡大行列の表示
    printf("\nガウス・ジョルダン後の拡大行列\n----------------------------------------------------------------------------------------\n");
    printMatrix(matrix, n);
    for (int i = 0; i < 4; i++) {
        x[i] = matrix[i][4]; //解を保存
    }
}

int main() {
    //最初の拡大行列を表示
    printf("\n初期の拡大行列\n----------------------------------------------------------------------------------------\n");
    printMatrix(augmentedMatrix, 4);
    //ガウスジョルダン法で連立方程式を解く
    gauss_jordan(augmentedMatrix, 4);
    //解を表示
    printf("\n未知数x1 x2 x3 x4の解\n----------------------------------------------------------------------------------------\n");
    // printf("x1 = %f, x2 = %f, x3 = %f, x4 = %f", x[0], x[1], x[2], x[3]);
    for (int i = 0; i < 4; i++) {
        printf("- x%d = %f\n", i + 1, x[i]);
    }
		printf("\n検証\n----------------------------------------------------------------------------------------\n");
    int NGCount = 0;
    //未知数が正しいか検算
    for(int i = 0; i < 4; i++){
        double sum = 0;
        int j = 0;
        printf("test%d: ", i+1);
        for(j = 0; j < 4; j++){
            sum += initialConditions[i][j] * x[j];
            printf("%f·x%d", initialConditions[i][j], j + 1);
            if(j != 3)printf(" + ");
            
        }
        printf(" = %f | 理想の定数項 : %f ", sum, initialConditions[i][4]);
        //答えが同じか検証する
        if(sum == initialConditions[i][4]){
            printf("→OK\n");
        }else{
            printf("→NG\n");
            NGCount ++;
        }
    }
    if(NGCount == 0){
        printf("- 全ての検算がOKでした\n");
    }else{
        printf("- 検算にNGが%d件ありました\n", NGCount);
    }
    return 0;
}

出力結果

初期の拡大行列
----------------------------------------------------------------------------------------
|       3.000000        2.000000        7.000000        1.000000        8.000000        |
|       1.000000        5.000000        1.000000        -1.000000       5.000000        |
|       4.000000        1.000000        3.000000        -2.000000       7.000000        |
|       1.000000        6.000000        4.000000        3.000000        13.000000       |

ガウス・ジョルダン後の拡大行列
----------------------------------------------------------------------------------------
|       1.000000        0.000000        0.000000        0.000000        3.500000        |
|       0.000000        1.000000        0.000000        0.000000        1.000000        |
|       0.000000        0.000000        1.000000        0.000000        -1.000000       |
|       0.000000        0.000000        0.000000        1.000000        2.500000        |

未知数x1 x2 x3 x4の解
----------------------------------------------------------------------------------------
- x1 = 3.500000
- x2 = 1.000000
- x3 = -1.000000
- x4 = 2.500000
test1: 3.000000·x1 + 2.000000·x2 + 7.000000·x3 + 1.000000·x4 = 8.000000 | 理想の定数項 : 8.000000 →OK
test2: 1.000000·x1 + 5.000000·x2 + 1.000000·x3 + -1.000000·x4 = 5.000000 | 理想の定数項 : 5.000000 →OK
test3: 4.000000·x1 + 1.000000·x2 + 3.000000·x3 + -2.000000·x4 = 7.000000 | 理想の定数項 : 7.000000 →OK
test4: 1.000000·x1 + 6.000000·x2 + 4.000000·x3 + 3.000000·x4 = 13.000000 | 理想の定数項 : 13.000000 →OK
- 全ての検算がOKでした