第2-2回　ガウスの消去法[python]

$$\bf Ax=b$$

今回は上記のような行列を含む方程式を解く代表的な方法であるガウスの消去法について行います．A,bがわかっており、xを求めます。

まずは具体的な3×3の行列で考えます．

$$ \left[\begin{array}{ccc}a_{00}&a_{01}&a_{02} \\ a_{10}&a_{11}&a_{12} \\ a_{20}&a_{21}&a_{22} \end{array}\right] \bf x = \left[\begin{array}{c}b_{0} \\ b_{1} \\ b_{2} \end{array}\right] $$

一行目の式を①，二行目を②，三行目を③として計算を行います．

“前進消去”

対角成分の下側を全て0にすることを目標に進めます。
まず、\(a_{10}\)と\(a_{20}\)を\(0\)にするために次の計算を行います。

$$②式-①式\frac{a_{10}}{a_{00}}$$

$$③式-①式\frac{a_{20}}{a_{00}}$$

すると下記の様になります．

$$ \left[\begin{array}{ccc}a_{00}&a_{01}&a_{02} \\ 0&a_{11}-a_{01}\frac{a_{10}}{a_{00}}&a_{12} -a_{02}\frac{a_{10}}{a_{00}} \\ 0&a_{21} -a_{01}\frac{a_{20}}{a_{00}} &a_{22} -a_{02}\frac{a_{20}}{a_{00}} \end{array}\right] \bf x = \left[\begin{array}{c}b_{0} \\ b_{1}-
b_{0}\frac{a_{10}}{a_{00}} \\ b_{2} -b_{0}\frac{a_{20}}{a_{00}} \end{array}\right] $$

改めて上記の式を添え字を用いて表します。この添え字は計算し値が変化した回数を示します。

$$ \left[\begin{array}{ccc}a_{00}&a_{01}&a_{02} \\ 0&a^{[1]}_{11}&a^{[1]} _{12} \\ 0&a ^{[1]} _{21}&a ^{[1]} _{22} \end{array}\right] \bf x = \left[\begin{array}{c}b_{0} \\ b ^{[1]} _{1} \\ b ^{[1]} _{2} \end{array}\right] $$

さらに次の計算を行います。

$$③式-②式\frac{a ^{[1]} _{21}}{a ^{[1]} _{11}}$$

$$ \left[\begin{array}{ccc}a_{00}&a_{01}&a_{02} \\ 0&a ^{[1]} _{11}&a ^{[1]} _{12} \\ 0&0&a ^{[2]} _{22} \end{array}\right] \bf x = \left[\begin{array}{c}b_{0} \\ b ^{[1]} _{1} \\ b ^{[2]} _{2} \end{array}\right] $$

“後退消去”

対角成分の上側を全て0にすることを目標に計算を行います。
今度は\(a ^{[1]} _{12}\)と \(a_{02}\)を\(0\)にします。

$$②式-③式\frac{a ^{[1]} _{12}}{a ^{[2]} _{22}}$$

$$①式-③式\frac{a_{02}}{a ^{[2]} _{22}}$$

$$ \left[\begin{array}{ccc}a_{00}&a_{01}&0 \\ 0&a ^{[1]} _{11}&0 \\ 0&0&a ^{[2]} _{22} \end{array}\right] \bf x = \left[\begin{array}{c}b ^{[1]} _{0} \\ b ^{[2]} _{1} \\ b ^{[2]} _{2} \end{array}\right] $$

さらに次の計算を行います。

$$①式-②式\frac{a_{01}}{a ^{[1]} _{11}}$$

$$ \left[\begin{array}{ccc}a_{00}&0&0 \\ 0&a ^{[1]} _{11}&0 \\ 0&0&a ^{[2]} _{22} \end{array}\right] \bf x = \left[\begin{array}{c}b ^{[2]} _{0} \\ b ^{[2]} _{1} \\ b ^{[2]} _{2} \end{array}\right] $$

以上のことからxが導かれます。

$$\bf x= \left[\begin{array}{c} b^{[2]}_{0}/a_{00} \\ b ^{[2]} _{1}/a ^{[1]} _{11} \\ b ^{[2]} _{2}/a ^{[2]} _{22} \end{array}\right] $$

このように解を求めることができます。

ガウスの消去法をpythonの関数として作り、テストケースとして下記の例題を解いてみました。
また、テストケースのコードでは今後も使用できるように、list型のA行列,b行列を与えると、解を返すガウスの消去法を関数として定義しました。

$$\bf A = \left[\begin{array}{ccc}3&2&1 \\ 1&3&1 \\ 4&1&3 \end{array}\right], \bf b = \left[\begin{array}{c}10\\ 10 \\ 15 \end{array}\right] $$

解析解は\(\bf x = \left[\begin{array}{c}1\\ 2 \\ 3 \end{array}\right] \)となっています。

結果はこのようになりました．浮動小数を使っているのでわずかな誤差はありますが，一致しました．

[1.0000000000000002, 1.9999999999999996, 3.0000000000000004]

次回はこの方法で定常拡散方程式を解いてみたいと思います。

by hide

import numpy


def gaussian_elimination(A_list, b_list):
    A_array = numpy.array(A_list)
    b_array = numpy.array(b_list)

    """ 前進消去 """
    for i in range(1, A_array.shape[0]):
        for j in range(i, A_array.shape[0]):
            hi = A_array[j, i-1] / A_array[i-1, i-1]
            A_array[j] = A_array[j] - A_array[i-1] * hi
            b_array[j] = b_array[j] - b_array[i-1] * hi

    """ 後退代入 """
    for i in range(A_array.shape[0]-2, -1, -1):
        for j in range(i, -1, -1):
            hi = A_array[j, i+1] / A_array[i+1, i+1]
            A_array[j] = A_array[j] - A_array[i+1] * hi
            b_array[j] = b_array[j] - b_array[i+1] * hi

    """ 解の算出 """
    true_x = []
    for i in range(A_array.shape[0]):
        true_x.append(b_array[i] / A_array[i, i])

    return true_x


# ------------------------
# -- main               --
# ------------------------

A = [[3., 2., 1.], [1., 3., 1.], [4., 1., 3.]]
b = [10., 10., 15.]

print(gaussian_elimination(A, b))

# 真の解x=[1.0, 2.0, 3.0]