見出し画像

＜学習シリーズ＞微分方程式の解法：反復法のガウスザイデル法をPythonで学ぶ

KIYO

１．概要

１－１．緒言

　本記事は”学習シリーズ”として自分の勉強備忘録用になります。

　偏微分方程式の解法の中で反復法は下記手法がありますが、今回はヤコビ（Jacobi）法の改善版であるガウス・ザイデル法を説明します。

  • ヤコビ法：求める変数の計算がすべて完了したらまとめて更新する

  • ガウス・ザイデル法：更新した値を逐次的に使用していく

  • SOR法：ガウス・ザイデル法をさらに修正したモデル

  • 共役勾配法

１－２．サンプル問題

　今回のサンプル問題は3つの変数と式からなる連立方程式を解きます。なお説明のしやすさを考慮して連立方程式と線形代数の２種を記載しました。

$$
a_{11}x_{1} + a_{12}x_{2} + a_{13}x_{3} = b_{1} \\
a_{21}x_{1} + a_{22}x_{2} + a_{23}x_{3} = b_{2} \\
a_{31}x_{1} + a_{32}x_{2} + a_{33}x_{3} = b_{3}
$$

$$
AX=Bより\begin{bmatrix}a_{11}x_{1} + a_{12}x_{2} + a_{13}x_{3} \\a_{21}x_{1} + a_{22}x_{2} + a_{23}x_{3} \\a_{31}x_{1} + a_{32}x_{2} + a_{33}x_{3}\end{bmatrix}=\begin{bmatrix}b_{1} \\ b_{2} \\ b_{3}\end{bmatrix}
$$

$$
A=\begin{bmatrix}a_{11} & a_{12} & a_{13} \\a_{21} & a_{22} & a_{23} \\a_{31} & a_{32} & a_{33}\end{bmatrix},  X=\begin{bmatrix}x_{1} \\ x_{2} \\ x_{3}\end{bmatrix}, B=\begin{bmatrix}b_{1} \\ b_{2} \\ b_{3}\end{bmatrix}
$$

　サンプルの数値は下記を使用します。計算すると解は$${{x1: 1,  x2: 1,  x3: -1}}$$となります。

$$
2x_{1} + x_{2} + x_{3} =2\\
2x_{1} +3x_{2} + x_{3} =4\\
x_{1} + x_{2} +3x_{3} =-1
$$

$$
AX=Bより\begin{bmatrix}2x_{1} + x_{2} + x_{3} \\2x_{1} + 3x_{2} + x_{3} \\x_{1} + x_{2} + 3x_{3}\end{bmatrix}=\begin{bmatrix}2 \\ 4 \\ 1\end{bmatrix}
$$

$$
A=\begin{bmatrix}2 & 1 & 1 \\2 & 3 & 1 \\1 & 1 & 3\end{bmatrix}, X=\begin{bmatrix}x_{1} \\ x_{2} \\ x_{3}\end{bmatrix}, B=\begin{bmatrix}2 \\ 4 \\ 1\end{bmatrix}
$$

[IN]
import sympy as sp
sp.init_printing()

a11, a12, a13, a21, a22, a23, a31, a32, a33 = sp.symbols('a11 a12 a13 a21 a22 a23 a31 a32 a33')
x1, x2, x3 = sp.symbols('x1 x2 x3')
b1, b2, b3 = sp.symbols('b1 b2 b3')

A = sp.Matrix([[a11, a12, a13],
               [a21, a22, a23],
               [a31, a32, a33]])
B = sp.Matrix([[b1],[b2],[b3]])
X = sp.Matrix([[x1],[x2],[x3]])

eq1 = sp.Eq(A@X, B)
display(eq1)

A[0], A[1], A[2], A[3], A[4], A[5], A[6], A[7], A[8] = 2,1,1,2,3,1,1,1,3
B[0], B[1], B[2] = 2,4,-1
display(sp.Eq(A@X, B))

sp.solve(sp.Eq(A@X, B), [x1, x2, x3])

[OUT]

２．ヤコビ法とガウスザイデル法の違い

　ヤコビ法とガウスザイデル法の解法は基本的に同じです。大きな違いとしては下記の通りです。

【解法の違い】
●ヤコビ法：求める解Xをまとめて更新して$${X_{n}->X_{n+1}}$$とする
●ガウスザイデル法：1行ごとに初期値$${x_{n}}$$を使用して$${x_{n+1}}$$を計算する。それ以降の行の解に対してはすでに求めた解${x_{n+1}}$$を使用して計算する。

【ヤコビ法における１ループ（解$${x_{1}->x_{n}}$$までの計算）】

$$
\bm X_{n+1}=\bm D ^{-1}( \bm b-(\bm L+\bm U) \bm X_{n})
$$

　連立方程式で表すと下記の通りです。

$$
x_{1}=\frac{1}{A_{11}}(b_{1}-A_{12}x_{2}-A_{13}x_{3}-\dots-A_{1n}x_{n})
$$

$$
x_{2}=\frac{1}{A_{22}}(b_{2}-A_{21}x_{1}-A_{23}x_{3}-\dots-A_{2n}x_{n})
$$

$$
x_{3}=\frac{1}{A_{33}}(b_{3}-A_{31}x_{1}-A_{32}x_{2}-\dots-A_{3n}x_{n})\\
\vdots
$$

$$
 x_{n}= \frac{1}{A_{nn}}(b_{n}-A_{n1}x_{1}-A_{n2}x_{2}-\dots-A_{nn-1}x_{n-1})
$$

【ガウスザイデル法】
　ガウスザイデル法では更新された解$${x_{n+1}}$$が使用されるため下三角行列Lには更新された$${x}$$、上三角行列Uには更新前の$${x}$$で計算されます。

$$
\bm X_{n+1}=\bm D ^{-1}( \bm b-\bm L \bm X_{n+1}-\bm U \bm X_{n})
$$

　連立方程式で表すと下記の通りです。

$$
x_{1}^{n+1}=\frac{1}{A_{11}}(b_{1}-A_{12}x_{2}-A_{13}x_{3}-\dots-A_{1n}x_{n})
$$

$$
x_{2}=\frac{1}{A_{22}}(b_{2}-A_{21}x_{1}^{n+1}-A_{23}x_{3}-\dots-A_{2n}x_{n})
$$

$$
x_{3}=\frac{1}{A_{33}}(b_{3}-A_{31}x_{1}^{n+1}-A_{32}x_{2}^{n+1}-\dots-A_{3n}x_{n})\\
\vdots
$$

$$
 x_{n}= \frac{1}{A_{nn}}(b_{n}-A_{n1}x_{1}^{n+1}-A_{n2}x_{2}^{n+1}-\dots-A_{nn-1}x_{n-1}^{n+1})
$$

３．ガウスザイデル法の理解１：連立方程式Ver.

　　”直接法（LU分解やガウスの消去法）”では式変形することにより直接的に解Xを算出しました。

　反復法では式変形はせずに求めたい解Xに近似解を代入して逐次的に連立方程式を解くことで近似解を求めていきます。本章では具体的に数式を使用しながら説明します。
　ガウスザイデル法の大きな流れは下記の通りです。

  1. 求めたい解(変数)$${x}$$に適当な初期値を与える

  2. まずは1行目の$${x_{1}}$$を計算する。$${x_{1}}$$以外のxは固定値として$${x_{1_{n+1}}}$$を求める。

  3. 2行目で$${x_{2}}$$を求める場合、$${x_{1}}$$の値は（ヤコビ法のように）初期値$${x_{1_{n}}}$$ではなく先ほど求めた$${x_{1_{n+1}}}$$を使用する。

  4. 3行目の計算では事前に求めた$${x_{1_{n+1}}}$$, $${x_{2_{n+1}}}$$を使用して求める。

  5. i行目は事前に解を求めた行に関しては$${x_{n+1}}$$、i行目移行の解に関しては初期値$${x_{n}}$$を使用して計算する。

  6. 上記の計算を何回も繰り返し計算後の$${x_{i}}$$に変化がなくなった時点で計算を終了する

２－１．連立方程式の概要

　まず初めにわかりやすさを考慮して連立方程式でのガウスザイデル法のイメージを説明します。参考として下記の連立方程式を考えてみます。
　なおaは係数、xは求めたい解、bは定数となります。なお当たり前のことですが連立方程式を解くためには「解の数≦数式の数」ため式の行数は最低でも数式以上の数となります。

$$
A_{11}x_{1}+A_{12}x_{2}+A_{13}x_{3}+\dots+A_{1n}x_{n}=b_{1} \\ A_{21}x_{1}+A_{22}x_{2}+A_{23}x_{3}+\dots+A_{2n}x_{n}=b_{2} \\ A_{31}x_{1}+A_{32}x_{2}+A_{33}x_{3}+\dots+A_{3n}x_{n}=b_{3} \\ \vdots \\ A_{n1}x_{1}+A_{n2}x_{2}+A_{n3}x_{3}+\dots+A_{nn}x_{n}=b_{n}
$$　

　また反復計算をする前に$${x_{i}}$$に適当な(初期)値を与えます。Pythonコードでは正規分布(np.random.randn())を使用しました。（後述で初期値$${c_{i}}$$を与えていますが記号としては$${x_{i}}$$を使用していきます）

$$
\bm x =\begin{bmatrix}   x_{1} & x_{2} & \dots  & x_{n} \end{bmatrix} -> \bm x = \begin{bmatrix}  c_{1} & c_{2} & \dots  & c_{n} \end{bmatrix}
$$

２－２．１行だけで確認

　1行目の式から$${x_{1}}$$を解いてみます。xには既に（適当な初期）値を与えているため$${x_{1}}$$を求めることが出来ます（$${x_{1}}$$も初期値を与えておりますが$${x_{1}}$$を求めるときには使用せず上書きします）。

$$
A_{11}x_{1}+A_{12}x_{2}+A_{13}x_{3}+\dots+A_{1n}x_{n}=b_{1}
$$

$$
x_{1}=\frac{1}{A_{11}}(b_{1}-A_{12}x_{2}-A_{13}x_{3}-\dots-A_{1n}x_{n})
$$

　１行目はヤコビ法と全く同じとなります。上記で求めた解を$${x_{1}^{n+1}}$$として、次に２行目を計算します。

$$
A_{21}x_{1}^{n+1}+A_{22}x_{2}+A_{23}x_{3}+\dots+A_{2n}x_{n}=b_{2}
$$

$$
x_{2}=\frac{1}{A_{22}}(b_{2}-A_{21}x_{1}^{n+1}-A_{23}x_{3}-\dots-A_{2n}x_{n})
$$

　上記の通りヤコビ法では前の行で求めた解を使用しておりませんが、ガウスザイデル法では解を更新するごとに次の行でその値を使用していきます。

２－３．全体で確認

　式全体で考えると下記の通りです。

$$
A_{11}x_{1}+A_{12}x_{2}+A_{13}x_{3}+\dots+A_{1n}x_{n}=b_{1} \\ A_{21}x_{1}+A_{22}x_{2}+A_{23}x_{3}+\dots+A_{2n}x_{n}=b_{2} \\ A_{31}x_{1}+A_{32}x_{2}+A_{33}x_{3}+\dots+A_{3n}x_{n}=b_{3} \\ \vdots \\ A_{n1}x_{1}+A_{n2}x_{2}+A_{n3}x_{3}+\dots+A_{nn}x_{n}=b_{n}
$$

　上式を各行で分解すると下記の通りです。

$$
x_{1}=\frac{1}{A_{11}}(b_{1}-A_{12}x_{2}-A_{13}x_{3}-\dots-A_{1n}x_{n})
$$

$$
x_{2}=\frac{1}{A_{22}}(b_{2}-A_{21}x_{1}^{n+1}-A_{23}x_{3}-\dots-A_{2n}x_{n})
$$

$$
x_{3}=\frac{1}{A_{33}}(b_{3}-A_{31}x_{1}^{n+1}-A_{32}x_{2}^{n+1}-\dots-A_{3n}x_{n})\\
\vdots
$$

$$
 x_{n}= \frac{1}{A_{nn}}(b_{n}-A_{n1}x_{1}^{n+1}-A_{n2}x_{2}^{n+1}-\dots-A_{nn-1}x_{n-1}^{n+1})
$$

　上記より下記が確認できます。

  • i行目にある求めたい解$${x_{i}}$$の分母は$${A}$$成分の対角の値

  • 構造はすべての行で同じであり分子：定数b-対角成分以外の内積の差分、分母：係数Aの対角の値

  • i行目で使用される解xに関してi行目より前の解は更新後x、i行目より後の解は更新前xを使用

４．ガウスザイデル法の理解２：行列Ver.

４－１．行列式の概要

　次に行列で定式化します。$${\bm A}$$は係数行列、$${\bm X}$$は解ベクトル、$${\bm b}$$は定数ベクトルであり求めたいのは$${\bm X}$$です。

$$
\bm A \cdot \bm X = \bm bより  
$$

$$
\begin{bmatrix}A_{11} & A_{12} & A_{13} \dots A_{1n} \\A_{21} & A_{22} & A_{23} \dots A_{2n} \\A_{31} & A_{32} & A_{33} \dots A_{3n} \\ \vdots & \vdots & \vdots \\A_{n1} & A_{n2} & A_{n3} \dots A_{nn} \end{bmatrix}
\begin{bmatrix} x_{1} \\ x_{2}\\x_{3} \\ \vdots \\ x_{n} \end{bmatrix}
= \begin{bmatrix} b_{1} \\ b_{2}\\b_{3}\\ \vdots \\ b_{n} \end{bmatrix}
$$

４－２．LDU分解／余剰行列R

　$${\bm A}$$は係数行列は下記の通り上三角行列U、下三角行列L、対角行列Dに分解できます。

$$
\bm A=\bm U + \bm D + \bm L =
$$

$$
\begin{bmatrix}0 & A_{12} & A_{13} \dots A_{1n} \\ 0 & 0 & A_{23} \dots A_{2n} \\ 0 & 0 & 0 \dots A_{3n} \\ \vdots & \vdots & \vdots \\ 0 & 0 & 0 \dots 0 \end{bmatrix}+
 \begin{bmatrix}A_{11} & 0 & 0 \dots 0 \\ 0 & A_{22} & 0 \dots 0 \\ 0 & 0 & A_{33} \dots 0 \\ \vdots & \vdots & \vdots \\ 0 & 0 & 0 \dots A_{nn} \end{bmatrix}+
\begin{bmatrix}0 & 0 & 0 \dots 0 \\A_{21} & 0 & 0 \dots 0 \\A_{31} & A_{32} & 0 \dots 0 \\ \vdots & \vdots & \vdots \\A_{n1} & A_{n2} & A_{n3} \dots 0 \end{bmatrix}
$$

４－３．数式の解法

　前述の式を考えると行列式は下記の通りに解くことが出来ます。

$$
\bm A \cdot \bm X = \bm bより  (\bm D+(\bm L+\bm U)) \cdot \bm X = \bm b
$$

$$
\bm D \bm X+(\bm L+\bm U) \bm X = \bm b\\
\bm D \bm X= \bm b-(\bm L+\bm U) \bm X \\
\bm X=\bm D ^{-1}( \bm b-(\bm L+\bm U)\bm X) \\
$$

　ガウスザイデル法のルールより「i行目移行の解xは更新後xを使用＝下三角行列Lには更新後xを使用」となるため、初期値$${ \bm X _{n}}$$として次のステップにおける解$${ \bm X _{n+1}}$$は下記の通り計算できます。

$$
\bm X_{n+1}=\bm D ^{-1}( \bm b-\bm L \bm X_{n+1}-\bm U \bm X_{n})
$$

４－４．解法の検算

　検算に関しては下記記事でご確認ください。

５．Pythonによるガウスザイデル法の実装

　サンプル問題は下記の通りであり理論解は$${{x_{1}: 1,  x_{2}: 1,  x_{3}: -1}}$$です。

$$
2x_{1} + x_{2} + x_{3} =2\\
2x_{1} +3x_{2} + x_{3} =4\\
x_{1} + x_{2} +3x_{3} =-1
$$

　実装のポイントは下記の通りです。計算結果よりヤコビ法が147回のループに対してガウスザイデル法では15回で近似解を得ることが出来ました。

  • 行列をLDU分解して余剰行列R=(L+U)を作成できるようにした

  • 解Xを求めるには最初に適当な初期値が必要なため乱数値で初期化

  • ヤコビ法のように行列計算ができないためループ計算で各行で計算したxを更新していく※

  • $${x_{n}->x_{n+1}}$$で更新していき残差が閾値以下になったらループを終了させる（無限ループにならないようにMAXループ数も設定）

※Pythonのループ計算は遅いためCythonやnumba.jit()を使用することで高速化が可能です。

[IN]
import numpy as np

class GaussSeidel:
    def __init__(self, A, b, thresh=1e-7, max_iter=1000):
        self.A = A
        self.b = b
        self.thresh = thresh
        self.max_iter = max_iter
        self.X0 = np.random.randn(len(b))
        
    def LDU_decomp(self, matrix):
        L = np.tril(matrix, -1) #下三角行列
        D = np.diag(np.diag(matrix)) #対角行列
        U = np.triu(matrix, 1) #上三角行列
        return L, D, U
    
    def solve(self):
        L, D, U = self.LDU_decomp(self.A) #L, D, Uを求める
        X = self.X0 #初期値

        for count in range(self.max_iter):
            X0 = X.copy() #更新前の値※X0=Xとすると参照渡しになるので注意
            #ガウスザイデル法による更新
            for i in range(len(X)):
                X[i] = (self.b[i] - (L+U)[i]@X)/D[i, i]

            residual = np.sum(np.sqrt((X - X0)**2))/np.sum(np.sqrt(X**2)) #残差
            
            count += 1
            if residual < self.thresh:
                break
        
        return X, count

np.random.seed(0)

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

b = np.array([2, 4, -1])

gaussseidel = GaussSeidel(A, b)
gaussseidel.solve()

[OUT]
(array([ 1.00000008,  0.99999996, -1.00000001]), 15)

【結果の確認】
　出力を出しながらどういう風に変化するか確認してみました。Xの初期値は$${{x_{1.764}: 1,  x_{2}: 0.4,  x_{3}: 0.979}}$$で計算しています。
　ヤコビ法では更新回数が低い時は理論解$${{x_{1}: 1,  x_{2}: 1,  x_{3}: -1}}$$より遠いところでジグザグ動いていましたがガウスザイデル法では比較的直線的に近似解に向かっております。
　ヤコビ法では１００回以上かかりましたが、ガウスザイデル法では１５回で充分な結果が得られました。

参考資料

あとがき

　これで一度放熱シミュレーションしてみたい・・・・

この記事が気に入ったらサポートをしてみませんか？