Pythonを用いた連立1次方程式の数値解法(ガウス・ザイデル法)
はじめに
連立1次方程式の数値解を得るためのガウス・ザイデル法について、考え方とPythonを用いた計算のやり方を簡単に記す。
ガウス・ザイデル法
次の連立方程式を考える。
以前に触れた ヤコビ法 は次の近似式を用い、解の最初の予測値を右辺に代入して を算出し、得られた値を第2の予測値として再び右辺に代入して計算、を繰り返し近似解を得る方法であった。
ガウス・ザイデル法はすでに求まった についてはその場ですぐに使用する。式で書くと次のようになる。1番目の式で求まった は2番目の式で使用し、1, 2番目の式で求まった は両方とも3番目の式で使用する。
これら近似式は次のとおり行列の形で比較的簡単に表せる。まずはじめの方程式は
として、
係数行列は、その対角成分を取り出した行列 、対角と対角より下をゼロとした上三角の行列 、対角と対角より上をゼロとした下三角の行列 より、 と表せるので、
を右辺に移し、両辺に をかけると、
逐次代入で数値解を得る形を考えると、
Pythonを用いた解き方は例えば次のようになる。
import numpy as np def gauss_seidel(A, b): x0 = np.zeros_like(b) # xの初期値 err = 1e10 # 収束判定用 D = np.diag(A.diagonal()) # 係数行列Aの対角行列 U = np.triu(A, k=1) # 係数行列Aの対角を除く上三角の行列 L = np.tril(A, k=-1) # 係数行列Aの対角を除く下三角の行列 invLD = np.linalg.inv(L + D) # (L + D) の逆行列 while err > 1e-5: x1 = np.dot(invLD, (b - np.dot(U, x0))) err = np.linalg.norm(b - np.dot(A, x1))/np.linalg.norm(b) # 収束判定 x0 = x1 return x0 A = np.array([[3, 1, -1], [-1, 2, 1], [1, -1, 2]]) b = np.array([2, 6, 5]) print(gauss_seidel(A, b))
解は] であるが、上記コードの計算結果は下記でありよい近似値と言える。(なお上のコードは3x3の連立方程式に限らず、nxnについて適用できる。)
[0.9999933 1.99999866 3.00000268]
ガウス・ザイデル法はヤコビ法より早く計算が済む。収束判定条件を同じにしてヤコビ法とガウス・ザイデル法の計算時間を実際に比較してみると(Jupyter Notebook の %timeit で測定)、
ヤコビ法:
631 µs ± 1.57 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
ガウス・ザイデル法:
209 µs ± 725 ns per loop (mean ± std. dev. of 7 runs, 1000 loops each)
となり、この問題については実測としてガウス・ザイデル法のほうが3倍くらい早い結果となった。