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

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

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}

これら方程式を次の近似解を得る形に書き直し、


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

解の最初の予測値を右辺に代入して計算し、得られた値を第2の予測値として再び右辺に代入して計算、を繰り返し近似解を得る数値解法がヤコビ法である。 新しい解と1つ前の解の差が予め決めた許容値を下回るまで計算を繰り返す。

なおこれら近似式は行列の形で表すと以下のとおりとなる。まずはじめの方程式は


\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{R} より、 \mathbf{A} = \mathbf{D} + \mathbf{R} と表せるので、


(\mathbf{D}+\mathbf{R}) \mathbf{x} = \mathbf{b}

更に式変形していくと、次の形になる。


\mathbf{x} = \mathbf{D}^{-1} (\mathbf{b} - \mathbf{R} \mathbf{x})

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


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

これは、上述の3つ並べた近似式と同内容である。

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

import numpy as np

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

    D = np.diag(A.diagonal())  # 係数行列Aの対角成分を取り出した行列
    R = A - D  # 係数行列Aの対角成分を0にした行列
    invD = np.linalg.inv(D)  # Dの逆行列

    while err > 1e-5:
        x1 = np.dot(invD, (b - np.dot(R, 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(jacobi(A, b))

解は\mathbf{x}=[1, 2, 3] であるが、上記コードの計算結果は下記でありよい近似値と言える。

[0.99999306 2.00002081 2.99998959]