技術系サラリーマン勉強記

数学、物理、プログラミングなど日々勉強した内容を取り扱っていきます。

ガウスの消去法の実装

数値解析の本（数値計算 (理工系の基礎数学 8)）に記載されていたガウスの消去法をpythonで実装してみました。部分ピボットはまだ実装していません。

ガウスの消去法とは

（N元）連立一次方程式の解を求める方法のひとつで、大きく分けて前進消去と後退代入の二つの手続きによって解を導出していきます。

一つ目の操作である前進消去は、方程式を上から下に並べたとして、上の方の一次方程式を用いて他の下の方の一次方程式の変数を一つ消去します。これを繰り返していくと、最終的に変数が一つの方程式が導出されます。

二つ目の後退代入は、下の方にある変数が一つの方程式を解き変数の値を確定し、その値を一つ上の方程式に代入しまた別の変数の値を確定していきます。繰り返すことで連立一次方程式の解を導出することができます。

ガウスの消去法の実装

\begin{align} 6x + 3y + 2z = 5 \\\ 5x + y + 8z = 1 \\\ 9x + 4y + 7z = 4 \end{align}

上の3元連立1次方程式を例にガウスの消去法によって解を求めてみました。コードは下記の通りです。

import numpy as np

# 各係数の設定
A = np.array([[6, 3, 2], 
              [5, 1, 8], 
              [9, 4, 7]], dtype='float')
b = np.array([5, 1, 4], dtype='float')
N = A.shape[0]
x = np.zeros(3)

# 前進消去
for k in range(N-1):
    for i in range(k+1, N):
        t = A[i, k] / A[k, k]
        
        for j in range(k+1, N):
            A[i, j] = A[i, j] - t*A[k, j]
        
        b[i] = b[i] - t*b[k]

# 後退代入
x[N-1] = b[N-1]/A[N-1, N-1]

for i in reversed(range(N-1)):
    u = 0
    
    for k in range(i+1, N):
        u = u + A[i, k]*x[k]
    
    x[i] = (b[i] - u)/A[i,i]

print(x)

これを実行すると、 x = 2.94117647, y=-3.35294118, z=-1.29411765 となります。

この解が正しいかnumpyのlinalg.solve関数を用いて確認してみました。

import numpy as np

A = np.array([[6, 3, 2], 
              [5, 1, 8], 
              [9, 4, 7]])
b = np.array([5, 1, 4])

x2 = np.linalg.solve(A, b)

print(x2)

x = 2.94117647, y=-3.35294118, z=-1.29411765 と、同じ値になっていました。実装に大きな間違いはなさそうでよかったです。

今度は部分ピボット選択をコードに取り入れてみたいと思います。

参考文献