Python NumPy サンプルコード: SOR 法
今回は、線形方程式の数値計算法の 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]