Python NumPy サンプルコード: ヤコビ法

今回は、線形方程式の数値計算法の ヤコビ法 を紹介します。

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

大きな連立方程式を計算する場合、LU 分解などによる直接法よりも計算時間が早い 反復法 がよく用いられます。

ヤコビ法

反復法の一つであるヤコビ法は、係数行列 \(\mathbf{A}\) を対角行列 \(\mathbf{D}\) とその残余 \(\mathbf{R}\) の和に分解して計算します。 \(\boldsymbol{x}\) に適当な初期値を与え、解が収束するまで反復計算していきます。 解が収束する十分条件は、係数行列の対角要素の絶対値が非対角要素の絶対値よりも相対的に大きいことです。

計算式の導出は

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

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

\begin{equation} (\mathbf{D}+\mathbf{R})\boldsymbol{x} = \boldsymbol{b} \nonumber \end{equation}

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

\begin{equation} \mathbf{D}\boldsymbol{x} = \boldsymbol{b}-\mathbf{R}\boldsymbol{x} \nonumber \end{equation}

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

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

ここで、\(\mathbf{D}\) は対角行列なので、その逆行列は対角成分の逆数を並べた対角行列であり、計算が簡単なのがポイントです。

Python にヤコビ法を解くメソッドはないので、自前で定義すると以下のように書くことができます。

import numpy as np


def jacobi(A, b, tol):
    # 線形連立方程式をヤコビ法で解く
    xOld = np.empty_like(b)
    error = 1e12

    D = np.diag(A)
    R = A - np.diagflat(D)

    while error > tol:
        x = (b-np.dot(R, xOld))/D
        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 = jacobi(A, b, 1e-9)
print(x)
[-0.80693816  1.11613876 -0.3092006   1.76470588]

コメント

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