Pythonを用いた連立1次方程式の数値解法(ヤコビ法)
はじめに
連立1次方程式の数値解を得るためのヤコビ法について、考え方とPythonを用いた計算のやり方を簡単に記す。
ヤコビ法
例えば次の連立方程式を考える。
これら方程式を次の近似解を得る形に書き直し、
解の最初の予測値を右辺に代入して計算し、得られた値を第2の予測値として再び右辺に代入して計算、を繰り返し近似解を得る数値解法がヤコビ法である。 新しい解と1つ前の解の差が予め決めた許容値を下回るまで計算を繰り返す。
なおこれら近似式は行列の形で表すと以下のとおりとなる。まずはじめの方程式は
として、
係数行列は、その対角成分を取り出した行列と対角成分をゼロとした行列 より、 と表せるので、
更に式変形していくと、次の形になる。
逐次代入で数値解を得る形を考えると、
これは、上述の3つ並べた近似式と同内容である。
Pythonを用いた解き方は例えば次のようになる。
import numpy as np def jacobi(A, b): x0 = np.zeros_like(b) # xの初期値 err = 1e10 # 収束判定用 D = np.diag(A.diagonal()) # 係数行列Aの対角成分を取り出した行列 R = A - D # 係数行列Aの対角成分を0にした行列 invD = np.linalg.inv(D) # Dの逆行列 while err > 1e-5: x1 = np.dot(invD, (b - np.dot(R, 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(jacobi(A, b))
解は] であるが、上記コードの計算結果は下記でありよい近似値と言える。
[0.99999306 2.00002081 2.99998959]