Python NumPy サンプルコード: ガウス・ザイデル法

今回は、線形方程式の数値計算法の ガウス・ザイデル法 を紹介します。

Python で線形連立方程式を反復法で解く

ガウス・ザイデル法

ガウス・ザイデル法は 先日のヤコビ法 と同じで反復法の一種です。 係数行列 \(\mathbf{A}\) を下三角行列 \(\mathbf{L}\) と上三角行列 \(\mathbf{U}\) の和に分解して計算します。 ヤコビ法と同じで \(\boldsymbol{x}\) に適当な初期値を与え、解が収束するまで反復計算していきます。 解が収束する十分条件も同じで、係数行列の対角要素の絶対値が非対角要素の絶対値よりも相対的に大きいことです。 直列計算の計算速度ではガウス・ザイデル法が有利ですが、ヤコビ法は容易に計算を並列化できる点がメリットだそうです。

計算式の導出は

\begin{equation} \mathbf{A}\boldsymbol{x} = \boldsymbol{b} \nonumber \end{equation}

に \(\mathbf{A} = \mathbf{L}+\mathbf{U}\) を代入すると

\begin{equation} (\mathbf{L}+\mathbf{U})\boldsymbol{x} = \boldsymbol{b} \nonumber \end{equation}

となり、\(\mathbf{U}\boldsymbol{x}\) を右辺に移動させて以下のようになります。

\begin{equation} \mathbf{L}\boldsymbol{x} = \boldsymbol{b}-\mathbf{U}\boldsymbol{x} \nonumber \end{equation}

n 回目の反復で得られた \(\boldsymbol{x}\) の値を \(\boldsymbol{x}_n\) すると、 以下のような反復法の漸化式が得られます。

\begin{equation} \boldsymbol{x}_{n+1} = \mathbf{L}^{-1}(\boldsymbol{b}-\mathbf{U}\boldsymbol{x}_{n}) \nonumber \end{equation}

Python にガウス・ザイデル法のメソッドはないので、自前で定義すると以下のように書くことができます。

import numpy as np


def gaussSeidel(A, b, tol):
    # 線形連立方程式をガウス・ザイデル法で解く
    xOld = np.empty_like(b)
    error = 1e12

    L = np.tril(A)
    U = A - L
    LInv = np.linalg.inv(L)

    while error  > tol:
        x = np.dot(LInv, b-np.dot(U, xOld))
        error = np.linalg.norm(x-xOld)/np.linalg.norm(x)
        xOld = x
    return x

A = np.array([[7, -1, 0, 1],
              [-1, 9, -2, 2],
              [0, -2, 8, -3],
              [1, 2, -3, 10]])
b = np.array([-5, 15, -10, 20])

x = gaussSeidel(A, b, 1e-9)
print(x)
[-0.80693816  1.11613876 -0.3092006   1.76470588]

コメント

Comments powered by Disqus
書籍更新情報
2017-02-18
Pythonによる科学技術計算 基礎編
1.4版への更新が可能になりました。
サポートページはこちら
電子書籍
Pythonによる科学技術計算 基礎編
電子書籍
線形代数(1): Pythonによる科学技術計算 実践編