モデリングと数値計算(主に化学工学向け)

化学反応、移動現象のモデリングと数値計算の話

Pythonを用いた連立1次方程式の数値解法(ガウス・ザイデル法)

はじめに

連立1次方程式の数値解を得るためのガウス・ザイデル法について、考え方とPythonを用いた計算のやり方を簡単に記す。

ガウス・ザイデル法

次の連立方程式を考える。


\begin{bmatrix}
a_{1,1} & a_{1,2} & a_{1,3} \\
a_{2,1} & a_{2,2} & a_{2,3} \\
a_{3,1} & a_{3,2} & a_{3,3} \\
\end{bmatrix}

\begin{bmatrix}
x_1 \\
x_2 \\
x_3 \\
\end{bmatrix}

=
\begin{bmatrix}
b_1 \\
b_2 \\
b_3 \\
\end{bmatrix}


以前に触れた ヤコビ法 は次の近似式を用い、解の最初の予測値を右辺に代入して x_1, x_2, x_3 を算出し、得られた値を第2の予測値として再び右辺に代入して計算、を繰り返し近似解を得る方法であった。


\begin{aligned}
x_1^{(k+1)} &= \frac{1}{a_{1,1}}(b_1 - a_{1,2}x_2^{(k)} - a_{1, 3}x_3^{(k)}) \\
x_2^{(k+1)} &= \frac{1}{a_{2,2}}(b_2 - a_{2,1}x_1^{(k)} - a_{2, 3}x_3^{(k)}) \\
x_3^{(k+1)} &= \frac{1}{a_{3,3}}(b_3 - a_{3,1}x_1^{(k)} - a_{3, 2}x_2^{(k)}) \\
\end{aligned}

ガウス・ザイデル法はすでに求まった x についてはその場ですぐに使用する。式で書くと次のようになる。1番目の式で求まった x_1 は2番目の式で使用し、1, 2番目の式で求まった x_1, x_2 は両方とも3番目の式で使用する。

f:id:schemer1341:20200809000829p:plain:w250

これら近似式は次のとおり行列の形で比較的簡単に表せる。まずはじめの方程式は


\mathbf{A}=
\begin{bmatrix}
a_{1,1} & a_{1,2} & a_{1,3} \\
a_{2,1} & a_{2,2} & a_{2,3} \\
a_{3,1} & a_{3,2} & a_{3,3} \\
\end{bmatrix}, \qquad

\mathbf{x}=
\begin{bmatrix}
x_1 \\
x_2 \\
x_3 \\
\end{bmatrix}, \qquad

\mathbf{b}=
\begin{bmatrix}
b_1 \\
b_2 \\
b_3 \\
\end{bmatrix}

として、


\mathbf{A} \mathbf{x} = \mathbf{b}

係数行列\mathbf{A}は、その対角成分を取り出した行列 \mathbf{D}、対角と対角より下をゼロとした上三角の行列  \mathbf{U}、対角と対角より上をゼロとした下三角の行列  \mathbf{L}より、 \mathbf{A} = \mathbf{L} + \mathbf{D} + \mathbf{U} と表せるので、


(\mathbf{L} + \mathbf{D} + \mathbf{U}) \mathbf{x} = \mathbf{b}

 \mathbf{U} \mathbf{x} を右辺に移し、両辺に (\mathbf{L} + \mathbf{D})^{-1} をかけると、


\mathbf{x} = (\mathbf{L} + \mathbf{D})^{-1} (\mathbf{b} - \mathbf{U} \mathbf{x})

逐次代入で数値解を得る形を考えると、


\mathbf{x}^{(k+1)} = (\mathbf{L} + \mathbf{D})^{-1} (\mathbf{b} - \mathbf{U} \mathbf{x}^{(k)})

Pythonを用いた解き方は例えば次のようになる。

import numpy as np

def gauss_seidel(A, b):
    x0 = np.zeros_like(b)  # xの初期値
    err = 1e10 # 収束判定用

    D = np.diag(A.diagonal())  # 係数行列Aの対角行列
    U = np.triu(A, k=1)  # 係数行列Aの対角を除く上三角の行列
    L = np.tril(A, k=-1)  # 係数行列Aの対角を除く下三角の行列
    invLD = np.linalg.inv(L + D)  # (L + D) の逆行列

    while err > 1e-5:
        x1 = np.dot(invLD, (b - np.dot(U, x0)))
        err = np.linalg.norm(b - np.dot(A, x1))/np.linalg.norm(b) # 収束判定
        x0 = x1

    return x0

A = np.array([[3, 1, -1],
              [-1, 2, 1],
              [1, -1, 2]])

b = np.array([2,
              6,
              5])

print(gauss_seidel(A, b))

解は\mathbf{x}=[1, 2, 3] であるが、上記コードの計算結果は下記でありよい近似値と言える。(なお上のコードは3x3連立方程式に限らず、nxnについて適用できる。)

[0.9999933  1.99999866 3.00000268]

ガウス・ザイデル法はヤコビ法より早く計算が済む。収束判定条件を同じにしてヤコビ法とガウス・ザイデル法の計算時間を実際に比較してみると(Jupyter Notebook の %timeit で測定)、

ヤコビ法:

631 µs ± 1.57 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

ガウス・ザイデル法:

209 µs ± 725 ns per loop (mean ± std. dev. of 7 runs, 1000 loops each)

となり、この問題については実測としてガウス・ザイデル法のほうが3倍くらい早い結果となった。