C言語でのガウス法による線形連立方程式の求解

#連立方程式を解く

学校などで数値計算の授業があると、「方程式を以下にして解くか」ということが一大テーマになっています。今回はその中で最も基本的な、「線形連立方程式を解く」ということに着目します。

線形の連立方程式は一般に、下のように係数行列Aと定数ベクトルbを使って下のように書くことができることはよく知られた事実です。
$$
Ax = b
$$
この形式で書くと、解xを求めるには、Aの逆行列を求めれば良いことがわかります。
$$
x=A^{-1} b
$$

#ガウス法
今回は、ガウス法を用いて線形連立方程式を解く方法を紹介します。これは明示的に逆行列を計算して解くような方法ではありませんが、本質的には逆行列を用いて解く方法と同じです。

ガウス法では、行変形を行って係数行列を上三角行列に変形します。行変形とは、ある行を定数倍したり、ある行に別の行を足したり引いたりする操作のことです。

上三角行列に変形したあとには、下の行から、対角項を1に、それ以外を0にしていきます。詳しくはこちらのHPを見るとわかりやすいと思います。

例として、２元連立方程式

\begin{align}
\begin{cases}
2x+y=4 \\
4x-3y=-2
\end{cases}
\end{align}

を解くことを考えます。係数行列をガウス法で計算していくと、

\begin{pmatrix}
2 & 1 \\ 4& 3
\end{pmatrix}\rightarrow
\begin{pmatrix}
2 & 1 \\ 0& -5
\end{pmatrix}\rightarrow
\begin{pmatrix}
2 & 1 \\ 0& 1
\end{pmatrix}\rightarrow
\begin{pmatrix}
1 & 0 \\ 0& 1
\end{pmatrix}

左から２つめの行列の形が上三角行列です。この行列に施した操作と同じことを右辺のベクトル(4, -2)についてもやっていくと、

\begin{pmatrix}
4 \\ -2
\end{pmatrix}\rightarrow
\begin{pmatrix}
4 \\ -10
\end{pmatrix}\rightarrow
\begin{pmatrix}
4 \\ 2
\end{pmatrix}\rightarrow
\begin{pmatrix}
1 \\ 2
\end{pmatrix}

となります。この結果から、方程式の解は

x=1, y=2

と求めることができます。複雑そうなことをしているように見えますが、やっていることは中学で習った加減法と一緒です。上三角行列を作ることは、加減法で言うxの消去に対応します。

C言語でのガウス法

上で書いたことを素直にプログラムに書き起こすと下のようになります。

#include <stdio.h>

#define N  3   //the dimension of equation

void vec_diff(float a[N], float b[N]){
    /* Calcurate the difference of two vectors. Be caution that b[N] changes.*/
    for (int i = 0; i < N; i++){
        b[i] -= a[i];
    }
}

int main(){
    float m[N][N] = {{5,-1,-1}, {2,1,-3},{1,1,1}};    // The matrix
    float b[N] = {0,-5,6};

    printf("The coefficient matrix is : \n");
    for (int i = 0; i < N; i++){
        for (int j = 0; j < N; j++){
            printf("%1.f ", m[i][j]);
            if (j == N-1){
                printf("\n");
            }
        }
    }

    printf("\nUse Gauss method to solve equations : \n");
    for (int i = 0; i < N; i++){
        for (int j = i+1; j < N; j++){
            float coef = m[j][i] / m[i][i];
            float del[N];

            for (int k = 0; k < N; k++){
                del[k] = m[i][k] * coef;
            }
            vec_diff(del, m[j]);
            b[j] -= b[i] * coef;
        }
    }

    for (int i = N -1; i >= 0; i--){
        float x = 1. / m[i][i];
        m[i][i] *= x;
        b[i] *= x;

        for (int j = 0; j < i; j++){
            b[j] -= b[i]*m[j][i];
            m[j][i] = 0;
        }
    }

    for (int i = 0; i < N; i++){
        for (int j = 0; j < N; j++){
            printf("%1.f ", m[i][j]);
            if (j == N - 1){
                printf("\n");
            }
        }
    }

    for (int i = 0; i < N; i++){
        printf("%f ", b[i]);
    }

    return 0;
}

このプログラムでは方程式、

\begin{align}
\begin{cases}
5x-y-z=0 \\
2x+y-3z = -5 \\
x+y+z=6
\end{cases}
\end{align}

を解いています。これを実行すると解

x=1, \ \ y=2, \ \ z=3

が得られます。上のプログラムで

#define N 3

のところを変えれば多元線形連立方程式を解くことができます。

このプログラムは簡単なものであり、桁落ち等、数値解析で考えなくてはならないところはちゃんと考えていません。もっとしっかり考えるためには桁落ちを防ぐために計算の順序を考慮に入れたプログラムを書く必要があります。

ちなみに、このコードはPythonで書くと非常にシンプルになります。興味のある方はPythonでのガウス法による連立方程式の求解も見てください。