C言語で解く連立1次方程式：ガウスの消去法実装例と注意点を解説

本稿では、C言語を使って連立1次方程式を解くために、ガウスの消去法の実装方法をわかりやすく紹介します。

前進消去から後退代入までの流れを具体例を交えて解説し、ゼロ除算など注意すべき点についても触れています。

初心者でも理解しやすい内容となっております。

ガウスの消去法の基本

ガウスの消去法は、$Ax=b$ の形の連立一次方程式を解くための手法です。

行列を段階的に変形し、上三角行列へ変換することで、後退代入により各未知数を求める方法です。

アルゴリズムの流れ

ガウスの消去法は大きく分けて「前進消去」と「後退代入」の二段階で進みます。

まず、前進消去で係数行列を上三角行列に変換し、次に後退代入で各未知数を順に解いていきます。

以下の数式で示すように、前進消去では各行の係数を利用して他の行の係数をゼロにしていきます。

$$\text{factor}=\frac{A_{k,i}}{A_{i,i}},A_{k,j}=A_{k,j}–\text{factor}\times A_{i,j}$$

この処理により、すべての行で主要対角成分以外がゼロとなった上三角行列が得られます。

前進消去の手順

前進消去では、上から順に各行の係数を利用し、下の行に対して以下の操作を行います。

1. 各行について、現在のピボット要素 $A_{i,i}$ がゼロでないかを確認します。
2. ゼロの場合は、ピボット選択の手法を用いて、下位の行と交換し、ゼロでない要素をピボットにします。
3. 各下位行 $k$ に対して、係数 $\text{factor}=\frac{A_{k,i}}{A_{i,i}}$ を計算し、該当行の全ての要素を更新します。

これにより行列全体が上三角行列となり、後続の計算が容易になります。

後退代入の手順

前進消去により得られた上三角行列を用いて、後退代入で解を求めます。

1. 最終行から始め、$x_i$ を次の式で計算します。

$$x_i=\frac{1}{A_{i,i}}\left(b_i–\sum _{j=i+1}^n{A}_{i,j}{x}_j\right)$$

2. この処理を上方向に繰り返し、すべての未知数 $x_n,x_{n-1},\dots ,x_1$ を求めます。

各ステップで既に求めた $x_j$ を用いるため、計算順序が非常に重要となります。

C言語での実装方法

C言語では、配列やループを活用してガウスの消去法を効率的に実装することができます。

ここでは、連立方程式を解くための基本的な構造とサンプルコードをご紹介します。

基本構造と配列の利用

C言語での実装では、連立方程式の係数と定数項をまとめた拡大係数行列を2次元配列として取り扱います。

配列の各行は一つの方程式に対応し、最後の列には定数項が格納されます。

たとえば、3元連立方程式の場合、行列のサイズは $3\times 4$ となります。

2次元配列の活用

2次元配列を利用することで、各行・各列へのアクセスが容易になり、行の入れ替えや要素の更新といった操作も直感的に行うことができます。

また、配列のサイズは定数やマクロで定義することで、コードの可読性とメンテナンス性が向上します。

インデックス管理の方法

C言語においては、配列のインデックスは0から始まるため、ループを用いる際は0からサイズ未満の範囲で管理する必要があります。

たとえば、行数を定義したマクロを使用して

for (i = 0; i < N; i++)

という形でループを構築し、正確な範囲で処理を実施します。

サンプルコードの紹介

以下に、ガウスの消去法を実装するサンプルコードを紹介します。

コード中のコメントや文字列リテラルは日本語で記述しており、変数名や関数名は英語表記となっています。

前進消去の実装例

コード中の forwardElimination関数は、前進消去の手順を実装しています。

各行ごとにピボット要素のチェックと、下位行への消去処理を行っています。

ゼロ除算回避の処理

前進消去の際、ピボット要素がゼロの場合の対処として、下位の行との入れ替え処理を実装しています。

これにより、ゼロ除算の発生を未然に防ぐとともに、計算の安定性が保たれます。

ピボット選択の方法

コード内では、現在のピボット $A[i][i]$ がゼロの場合に、下の行から非ゼロのピボット候補を探して入れ替える方法を用いています。

この方法により、計算の途中で不適切な除算が避けられ、正確な解が得られます。

後退代入の実装例

backwardSubstitution関数では、前進消去により得られた上三角行列を元に、後退代入を実施しています。

最終行から順に計算を進め、未知数 $x$ を求めています。

解の計算手順

最終的に main関数内で、前進消去と後退代入の両方を実行し、得られた解を出力しています。

以下のサンプルコードは、3元連立方程式の例を用いて、実際に連立方程式を解く流れを示しています。

#include <stdio.h>
#include <stdlib.h>
#define N 3  // 連立方程式の数
// 前進消去を行う関数
void forwardElimination(double A[N][N+1]) {
    int i, j, k;
    for (i = 0; i < N - 1; i++) {
        // ピボットがゼロの場合、下の行と入れ替えてゼロ除算を回避
        if (A[i][i] == 0) {
            for (k = i + 1; k < N; k++) {
                if (A[k][i] != 0) {
                    // 行の入れ替え
                    for (j = i; j < N + 1; j++) {
                        double temp = A[i][j];
                        A[i][j] = A[k][j];
                        A[k][j] = temp;
                    }
                    break;
                }
            }
        }
        // 下位行に対して消去処理を実施
        for (k = i + 1; k < N; k++) {
            // 今回もゼロ除算が発生する場合はスキップ
            if (A[i][i] == 0)
                continue;
            double factor = A[k][i] / A[i][i];
            for (j = i; j < N + 1; j++) {
                A[k][j] -= factor * A[i][j];
            }
        }
    }
}
// 後退代入により解を求める関数
void backwardSubstitution(double A[N][N+1], double x[N]) {
    int i, j;
    for (i = N - 1; i >= 0; i--) {
        x[i] = A[i][N];
        for (j = i + 1; j < N; j++) {
            x[i] -= A[i][j] * x[j];
        }
        // ゼロ除算のチェック(通常は前進消去で回避済み)
        if (A[i][i] == 0) {
            printf("解が一意に定まりません\n");
            exit(EXIT_FAILURE);
        }
        x[i] /= A[i][i];
    }
}
int main(void) {
    // 拡大係数行列の定義
    double A[N][N+1] = {
        {2, 1, -1, 8},    // 方程式1: 2x + y - z = 8
        {-3, -1, 2, -11}, // 方程式2: -3x - y + 2z = -11
        {-2, 1, 2, -3}    // 方程式3: -2x + y + 2z = -3
    };
    double x[N]; // 解を格納する配列
    // 前進消去による上三角行列への変換
    forwardElimination(A);
    // 後退代入による解の計算
    backwardSubstitution(A, x);
    // 解の出力
    printf("解は：\n");
    for (int i = 0; i < N; i++) {
        printf("x[%d] = %f\n", i, x[i]);
    }
    return 0;
}

解は：
x[0] = 2.000000
x[1] = 3.000000
x[2] = -1.000000

実装時の注意点

C言語でガウスの消去法を実装する際には、数値演算の精度やエラー処理など、いくつか注意すべきポイントがあります。

数値演算の精度管理

浮動小数点数(double型など)を使用することで、計算の精度をある程度確保することができます。

ただし、丸め誤差や計算結果の誤差が積み重なることがあるため、精度が要求される場合は特別な対策を検討する必要があります。

ゼロ除算と例外処理

前述の通り、ピボット要素がゼロとなるとゼロ除算が発生するため、必ずゼロ除算回避の処理を実装します。

また、例外的なケース(例えば、解が一意に定まらない場合)には、適切なメッセージ出力やプログラムの終了処理を加えることで、後続の処理に悪影響が出ないように工夫します。

デバッグ時のエラーチェック

実装したコードが期待通りの動作をするかどうか、デバッグ時には各段階で中間結果の出力や配列の状態を確認すると良いです。

特に、前進消去での行の入れ替え処理や後退代入の計算過程は、細かいエラーが混入しやすいため、デバッグ用のログを一時的に出力して確認するとトラブルシュートが容易になります。

まとめ

この記事では、行列を上三角行列に変換し未知数を解くガウスの消去法の基本と手順を解説しました。

前進消去で係数操作を行い、後退代入で解を求める流れを詳述し、C言語での具体的な実装例をサンプルコードを交えて説明しました。

さらに、ゼロ除算の回避やデバッグ時の注意点も触れ、実装時のポイントについてまとめています。