\(\def\bm{\boldsymbol}\)
\(\def\di{\displaystyle}\)
\(\def\ve{\varepsilon_0}\)
\(\def\dd#1#2{\dfrac{\partial #1}{\partial #2}}\)
\(\def\dda#1#2{\dfrac{\partial^2 #1}{\partial #2}}\)
流れ解析の差分法において、ポアソン方程式の解法などで、出てくる連立方程式を
\(\quad\bm{A}\bm{x}=\bm{b}\quad\cdots\,\)(1)
とするとき、直接法では次元が大きいと桁落ちなどで精度で問題が出る。特に対角項の絶対値が他の項より高い場合は、繰り返し法の方が演算回数が低く精度が高い場合がある。
まず、行列\(\,\bm{A}\,\)を対角行列\(\,\bm{D}\,\)と対角成分をゼロとした行列\(\,\bm{A1}\,\)に分けて、反復の式を作成する。
\(\quad\bm{A1}\,\bm{x}^{(n+1)}=\bm{D}\,\bm{x}^{(n)}+\bm{b}\quad\cdots\,\)(2)
式(2)を使って、\(n\,\)番目の近似解\(\,\bm{x}^{(n)}\)から\(\,n+1\,\)番目の近似解\(\,\bm{x}^{(n+1)}\)は、以下のようになる。
\(\quad\bm{x}^{(n+1)}=\bm{D}^{-1}\left(\bm{b}+\bm{A1}\,\bm{x}^{(n)}\right)\quad\cdots\,\)(3)
この方法はヤコビ法と呼ばれるもので、\(\bm{x}^{(n+1)}\)の近似値の計算に全てその前の値\(\,\bm{x}^{(n)}\)を使う。
これとは違い、\(\bm{x}^{(n+1)}\)の各成分の計算が終わると、それを直ちに使うのがガウス・ザイデル法である。
\(\quad x_i^{(n+1)}=a_{ii}^{-1}\left\{b_i-(a_{i1}x_1^{(n+1)}+a_{i2}x_2^{(n+1)}+\cdots+a_{ii+1}x_{i+1}^{(n)}+a_{ii+2}x_{i+2}^{(n)}+\cdots+a_{im}x_m^{(n)})\right\}\quad\cdots\,\)(4)
このガウス・ザイデル法をもう少し改造したのが SOR法 ( Successive Over-Relaxation )である。ガウス・ザイデル法の解の修正は、\(\bm{x}^{(n+1)}-\bm{x}^{(n)}\,\)であったが、これをもっと大きなステップにするものである。
例えば、３行目の計算では加速係数\(\,\omega\,\)を使って、以下のように計算する。
\(\left\{\begin{array}{l}x_3^{(n+1)}=a_{33}^{-1}\left\{b_3-(a_{31}x_1^{(n+1)}+a_{34}x_4^{(n)}+\cdots+a_{3m}x_{m}^{(n)})\right\}\quad\cdots\,(5)\\[2pt]
x_3^{(n+1)}=x_3^{(n)}+\omega(x_3^{(n+1)}-x_3^{(n)})\qquad\qquad\cdots\,(6)\end{array}\right.\)
ここで問題なのは\(\,\omega\,\)の値の取り方である。\(\omega=1.0\,\)だとガウス・ザイデル法となる。流れ解析の差分法のポアソン方程式の解法では、\(\omega=1.2\,\)位を取るケースが多いようである。
import numpy as np
# sor(A,b,w) A は行列 b は定数項 w : 加速係数( 1< w < 2 )
def sor(A,b,w,max_iter=512,eps=1.e-6):
D=np.diag(np.diag(A)) # 対角行列
A1=A-D # 対角行を除いた行列
D1=np.linalg.inv(D) # 逆行列（各値の逆数をとる）
M=D1@A1 # 対角行の逆数で割った行列
c=D1@b # 定数項も対角行の逆数で割る
n=len(b) # 次元
x=np.full(n,0.) # x の初期値 = 0.
for i in range(max_iter):
err=0.
for j in range(n):
s=x[j]
x[j]=(1-w)*s+w*(c[j]-M[j]@x)
err+=abs((x[j]-s)/x[j])
if err<eps: return x,i
・ 実行テスト
a=np.array([[9,1,2],
[1,9,1],
[2,1,9]])
b=np.array([9,18,-5])
x,i=sor(a,b,1.2)
print ('反復回数 = ',i)
print (x)
反復回数 = 10 [ 1.00000001 1.99999999 -0.99999974]