[C言語] 連立1次方程式を解く方法(吐き出し法/クラメルの公式/反復法)

#include <stdio.h>
#define N 3  // 方程式の数
void gaussElimination(float matrix[N][N+1]) {
    for (int i = 0; i < N; i++) {
        // ピボット選択
        for (int j = i + 1; j < N; j++) {
            if (matrix[j][i] > matrix[i][i]) {
                for (int k = 0; k <= N; k++) {
                    float temp = matrix[i][k];
                    matrix[i][k] = matrix[j][k];
                    matrix[j][k] = temp;
                }
            }
        }
        // 行の操作
        for (int j = i + 1; j < N; j++) {
            float ratio = matrix[j][i] / matrix[i][i];
            for (int k = i; k <= N; k++) {
                matrix[j][k] -= ratio * matrix[i][k];
            }
        }
    }
}
void backSubstitution(float matrix[N][N+1], float result[N]) {
    for (int i = N - 1; i >= 0; i--) {
        result[i] = matrix[i][N];
        for (int j = i + 1; j < N; j++) {
            result[i] -= matrix[i][j] * result[j];
        }
        result[i] /= matrix[i][i];
    }
}
int main() {
    float matrix[N][N+1] = {
        {2, 1, -1, 8},
        {-3, -1, 2, -11},
        {-2, 1, 2, -3}
    };
    float result[N];
    gaussElimination(matrix);
    backSubstitution(matrix, result);
    printf("解: \n");
    for (int i = 0; i < N; i++) {
        printf("x%d = %f\n", i + 1, result[i]);
    }
    return 0;
}

解: 
x1 = 2.000000
x2 = 3.000000
x3 = 1.000000

#include <stdio.h>
#define N 3 // 方程式の数

float determinant(float matrix[N][N]) {
    float det = 0;
    det = matrix[0][0] *
        (matrix[1][1] * matrix[2][2] - matrix[1][2] * matrix[2][1]) -
        matrix[0][1] *
        (matrix[1][0] * matrix[2][2] - matrix[1][2] * matrix[2][0]) +
        matrix[0][2] *
        (matrix[1][0] * matrix[2][1] - matrix[1][1] * matrix[2][0]);
    return det;
}

void cramer(float matrix[N][N + 1], float result[N]) {
    float detA = determinant((float[N][N]){{matrix[0][0], matrix[0][1], matrix[0][2]},
                                           {matrix[1][0], matrix[1][1], matrix[1][2]},
                                           {matrix[2][0], matrix[2][1], matrix[2][2]}});
    for (int i = 0; i < N; i++) {
        float tempMatrix[N][N];
        for (int j = 0; j < N; j++) {
            for (int k = 0; k < N; k++) {
                if (k == i) {
                    tempMatrix[j][k] = matrix[j][N]; // 右辺の列を置き換える
                }
                else {
                    tempMatrix[j][k] = matrix[j][k];
                }
            }
        }
        result[i] = determinant(tempMatrix) / detA; // 解を求める
    }
}

int main() {
    float matrix[N][N + 1] = {
        {2,  1,  -1, 8  },
        {-3, -1, 2,  -11},
        {-2, 1,  2,  -3 }
    };
    float result[N];
    cramer(matrix, result);
    printf("解: \n");
    for (int i = 0; i < N; i++) {
        printf("x%d = %f\n", i + 1, result[i]);
    }
    return 0;
}

解: 
x1 = 2.000000
x2 = 3.000000
x3 = -1.000000

#include <stdio.h>
#include <math.h>
#define N 3  // 方程式の数
#define EPSILON 0.0001  // 収束条件

void jacobi(float matrix[N][N + 1], float result[N]) {
    float temp[N];
    float old_result[N];
    do {
        for (int i = 0; i < N; i++) {
            old_result[i] = result[i];  // 現在の結果を保存
        }
        for (int i = 0; i < N; i++) {
            temp[i] = matrix[i][N];  // 右辺の値を初期化
            for (int j = 0; j < N; j++) {
                if (j != i) {
                    temp[i] -= matrix[i][j] * old_result[j];
                }
            }
            temp[i] /= matrix[i][i];  // 新しい値を計算
        }
        for (int i = 0; i < N; i++) {
            result[i] = temp[i];  // 結果を更新
        }
    } while (fabs(result[0] - old_result[0]) > EPSILON || fabs(result[1] - old_result[1]) > EPSILON || fabs(result[2] - old_result[2]) > EPSILON);
}

void gaussSeidel(float matrix[N][N + 1], float result[N]) {
    float temp[N];
    do {
        for (int i = 0; i < N; i++) {
            temp[i] = result[i];  // 現在の結果を保存
        }
        for (int i = 0; i < N; i++) {
            result[i] = matrix[i][N];  // 右辺の値を初期化
            for (int j = 0; j < N; j++) {
                if (j != i) {
                    result[i] -= matrix[i][j] * result[j];
                }
            }
            result[i] /= matrix[i][i];  // 新しい値を計算
        }
    } while (fabs(result[0] - temp[0]) > EPSILON || fabs(result[1] - temp[1]) > EPSILON || fabs(result[2] - temp[2]) > EPSILON);
}

int main() {
    float matrix[N][N + 1] = {
        {4, -1, 0, 3},
        {-1, 4, -1, 3},
        {0, -1, 4, 3}
    };
    float result[N] = { 0, 0, 0 };  // 初期値
    printf("ヤコビ法の解:\n");
    jacobi(matrix, result);
    for (int i = 0; i < N; i++) {
        printf("x%d = %f\n", i + 1, result[i]);
    }
    // 初期値をリセット
    for (int i = 0; i < N; i++) {
        result[i] = 0;
    }
    printf("ガウス＝ザイデル法の解:\n");
    gaussSeidel(matrix, result);
    for (int i = 0; i < N; i++) {
        printf("x%d = %f\n", i + 1, result[i]);
    }
    return 0;
}

ヤコビ法の解:
x1 = 1.071396
x2 = 1.285675
x3 = 1.071396
ガウス＝ザイデル法の解:
x1 = 1.071426
x2 = 1.285713
x3 = 1.071428