連立方程式を計算機に解かせる上で情報系や物理系の学生は必ず目にするのが、ヤコビ法[Jacobi]とガウス・ザイデル法[Gauss-Seidel]かなと思います。今回は大学院の課題でやったのを備忘録的に残しておきます。
忙しい人はgitへどうぞ
ヤコビ法[Jacobi]とガウス・ザイデル法[Gauss-Seidel]とは
ヤコビ法[Jacobi]とガウス・ザイデル法[Gauss-Seidel]両方ともに連立方程式を計算機に解かせるためのアルゴリズムです。ガウス・ザイデル法は計算結果をすぐ使うような工夫がされているため、ヤコビ方に比べて計算がとんでもなく早いです。
直説法と反復法の違い
直説法は普通に行列で連立方程式を解く際に用いられる方法で、下記の図のように連立方程式を行列で手計算で解く時に使われる方法です。
直説法だとコンピューターはあまり嬉しくない
直説法は人間にとっては早いかもしれないですが、ゼロでない部分が多く発生してしまい、その分計算を高速化することができません。
また、単調な繰り返しを得意とするコンピューターの性格上からしても複雑であまり得意ではないです。
そこで、ヤコビ法をはじめとする反復法では、対象の行列Aを対角行列D、 下三角行列Lと上三角行列Uに分解して考えます。
コード
以下にヤコビ法[Jacobi]とガウス・ザイデル法[Gauss-Seidel]の比較実験した際のコードを記載しますので、参考にしてみてください。
解いてみた方程式
下記のコードで解いてみている方程式は以下のようなものです。初期設定値を変えることでいくらでも遊べると思うので、遊んでみてください。
\left\{
\begin{array}{ll}
10x_1 - x_2 - 2x_3 &= 72 \\
-1x_1 + 10x_2 - 2x_3 &= 83 \\
-1x_1 - x_2 + 5x_3 &= 42 \\
\end{array}
\right.
$x_0$を初期値として。
x_0 = (0, 0, 0)
初期値の設定
A = np.array([[10, -1, -2],
[-1, 10, -2],
[-1, -1, 5]])
b = np.array([72, 83, 42])
ヤコビ法[Jacobi]
import numpy as np
from matplotlib import pyplot as plt
def jacobi(A, b, tolerance):
''' 連立方程式をヤコビ法（反復法）で解く '''
# 初期化（適当な解と残差）
x0 = [0, 0, 0]
residual = 1e20 #残差
# A = D(対角行列) + R(R := L + U)
D = np.diag(A)
R = A - np.diag(D)
# 反復計算=>残差がトレランスより小さくなったら終了
i = 0
print('----------Start iteration----------')
res = []
while residual > tolerance:
# 解と残差を計算
x = (b - (R @ x0)) / D
# xがbでx0がAxのつもりで、どれだけ解が欲しいターゲットの値と乖離しているかを計算している。
denominator = np.sqrt((x - x0) ** 2)
numerator = np.sqrt(x ** 2)
residual = np.sum(denominator) / np.sum(numerator)
res.append(residual)
# 100回ごとの誤差を出力
if i % 100 == 0:
print('Iteration=', i)
print('Residual=', residual)
print('x=', x)
i += 1
x0 = x
print('----------End iteration----------')
print('Iteration=', i)
print('Residual=', residual)
print('x=', x, '\n')
return x, res
# ヤコビ法で解を反復計算
x_jacobi, res1 = jacobi(A, b, 1e-6)
print('Jacobi method x=', x_jacobi, 'Validated b=', A @ x_jacobi)
ガウス・ザイデル法[Gauss-Seidel]
def gauss_seidel(A, b, tolerance):
''' 連立方程式をガウス・ザイデル法（反復法）で解く '''
# 初期化（適当な解と残差）
x = np.random.random(len(b))
residual = 1e20
# 係数行列を分解:D=対角成分ベクトル、R=剰余行列
D = np.diag(A)
R = A - np.diag(D)
# 反復計算（逐次求めた解を次のループで使うため、forを使っている）
iteration = 0
print('----------Start iteration----------')
res = []
while residual > tolerance:
# 収束判定用の1ステップ前の解を保持
x0 = x.copy()
# 逐次計算
for i in range(len(x)):
x[i] = (b[i] - R[i] @ x) / D[i]
# 収束判定
residual = np.sum(np.sqrt((x - x0) ** 2)) / np.sum(np.sqrt(x ** 2))
res.append(residual)
if iteration % 100 == 0:
print('Iteration=', iteration)
print('Residual=', residual)
print('x=', x)
iteration += 1
print('----------End iteration----------')
print('Iteration=', iteration)
print('Residual=', residual)
print('x=', x, '\n')
return x, res
# ガウス・ザイデル法で解を反復計算
x_gs, res2 = gauss_seidel(A, b, 1e-6)
print('Gauss-Seidel method x=', x_gs, 'Validated b=', A @ x_gs)
グラフの描画
def plot(r_jacobi, r_gauss_seidel):
''' Residualの変化をプロットする '''
# フォントの種類とサイズを設定する。
plt.rcParams['font.size'] = 14
plt.rcParams['font.family'] = 'DejaVu Serif'
# 目盛を内側にする。
plt.rcParams['xtick.direction'] = 'in'
plt.rcParams['ytick.direction'] = 'in'
# グラフの上下左右に目盛線を付ける。
fig = plt.figure()
ax1 = fig.add_subplot(111)
ax1.yaxis.set_ticks_position('both')
ax1.xaxis.set_ticks_position('both')
# 軸のラベルを設定する。
ax1.set_xlabel('Step')
ax1.set_ylabel('Residual')
# スケールの設定をする。
ax1.set_xlim(0, 120)
#ax1.set_ylim(-0.2, 1.2)
ax1.set_yscale('log')
# プロットを行う。
x1 = np.arange(0, len(r_jacobi), 1)
x2 = np.arange(0, len(r_gauss_seidel), 1)
ax1.plot(x1, r_jacobi, label='Jacobi', lw=1, color='red')
ax1.plot(x2, r_gauss_seidel, label='Gauss-Seidel', lw=1, color='blue')
ax1.legend()
# レイアウト設定
fig.tight_layout()
# グラフを表示する。
plt.show()
plt.close()
return
plot(res1, res2)
参考
コードは一般的な数式の表現に合うように改変して、コメントを追加したのみで下記とほとんど変わりませんので下記も一緒に参照してください。
数式などより詳しい解説は下記の方の方がいいと思います。