Pythonでのガウス法を用いた連立方程式の求解

#環境

• Linux Mint MATE 18.2
• Python3.6.3
• Atom エディター

コーディング

線形の連立方程式をガウス法と呼ばれる最も簡単な方法を用いて解きたいと思います。

連立方程式をプログラムで解くという話はC言語でのガウス法による線形連立方程式の求解で一度しています。このアルゴリズムの詳細についてはこの記事とそのリンクを見てください。この記事ではC言語を用いてコードを書きましたが、このような行列計算はPythonが得意としています。ですので、今回はPython3を用いて記事のコードを書きなおしてみたいと思います。

Cで書いたコードは、Pythonでnumpyを使うことでかなり簡潔に書くことができます。

import numpy as np

# solve "mat * x = b" about x
mat = np.array([[1., 2., 4.], [3., 4., 8.], [-1., 4., 2.]])
b = np.array([7., 15., 5.])

print(mat)  #original matrix
for i in range(len(mat)-1):
    for j in range(i+1, len(mat)):
        coef = mat[j][i] / mat[i][i]
        mat[j] -= mat[i] * coef
        b[j] -= b[i] * coef

for i in range(len(mat)-1, 0, -1):
    b[i] /= mat[i][i]
    mat[i] /= mat[i][i]
    for j in range(i):
        b[j] -= b[i] * mat[j][i]
        mat[j][i] = 0

print(b) #Solution

このコードでは、例として

\begin{align}
\begin{cases}
x+2y+4z=7 \\
3x+4y+8z=15 \\
-x+4y+2z = 5
\end{cases}
\end{align}

を解いています。これを実行すると[1.,1.,1.]が得られます。