Help us understand the problem. What is going on with this article?

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

More than 1 year has passed since last update.

連立方程式を解く

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

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

ガウス法

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

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

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

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

\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}

左から2つめの行列の形が上三角行列です。この行列に施した操作と同じことを右辺のベクトル(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でのガウス法による連立方程式の求解も見てください。

python_walker
趣味でプログラミングをしています。未熟者ではありますがよろしくお願いします。
Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
Comments
No comments
Sign up for free and join this conversation.
If you already have a Qiita account
Why do not you register as a user and use Qiita more conveniently?
You need to log in to use this function. Qiita can be used more conveniently after logging in.
You seem to be reading articles frequently this month. Qiita can be used more conveniently after logging in.
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
ユーザーは見つかりませんでした