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

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

SOR 法

SOR 法はガウス・ザイデル法の収束速度を改善させた手法です。係数行列 \(\mathbf{A}\) を対角行列 \(\mathbf{D}\), 下三角行列 \(\mathbf{L}\), 上三角行列 \(\mathbf{U}\) の和に分解して計算します。

計算式の導出は

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

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

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

となり、過緩和パラメータ \(\omega\) を用いて方程式は以下のように表されます。

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

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

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

\(\omega\) をあらかじめ最適に定めるられないことが SOR 法の弱点ですが、ヤコビ法の反復行列のスペクトル半径 \(\rho_{jacobi}\) が既知のときは、最適の \(\omega\) は次式となります。

\begin{equation} \omega = \frac{2}{1+\sqrt{1-\rho_{jacobi}^2}} \nonumber \end{equation}

SOR 法は Python では以下のように書くことができます。

import numpy as np


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

    D = np.diag(np.diag(A))
    L = np.tril(A, -1)
    U = A - L - D

    Mj = np.dot(np.linalg.inv(D), -(L+U))
    rho_Mj = max(abs(np.linalg.eigvals(Mj)))
    w = 2/(1+np.sqrt(1-rho_Mj**2))

    T = np.linalg.inv(D+w*L)
    Lw = np.dot(T, -w*U+(1-w)*D)
    c = np.dot(T, w*b)

    while error > tol:
        x = np.dot(Lw, xOld) + c
        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 = SOR(A, b, 1e-9)
print(x)

[-0.80693816  1.11613876 -0.3092006   1.76470588]