Python NumPy サンプルコード: SOR 法

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

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

SOR 法

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

計算式の導出は

15060401.png

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

15060402.png

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

15060403.png

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

15060304.png

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

15060305.png

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

# coding: utf-8
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]

コメント

Comments powered by Disqus
書籍更新情報
2017-01-16
Pythonによる科学技術計算 基礎編
1.3版への更新が可能になりました。
サポートページはこちら
電子書籍
Pythonによる科学技術計算 基礎編
Kindle ストアで販売中です
Pythonによる科学技術計算 基礎編
同人誌
技術書典(2016/6/25)
Emacs/org-modeのPDF作成術
電子版をBOOTHで販売中です
Emacs/org-modeのPDF作成術
Share