"""
概要: ヤコビ法による連立一次方程式の数値解法を提供するモジュール。
詳細説明:
このモジュールは、連立一次方程式 Ax = b をヤコビ法を用いて数値的に解くための関数を定義しています。
係数行列 A は `Hij`、右辺ベクトル b は `yi` に対応します。
解は反復計算により求められ、指定された収束条件または最大反復回数に達するまで計算が続行されます。
関連リンク: :doc:`Jacobi_usage`
"""
[ドキュメント]
def Jacobi(Hij, yi, nmaxiter, eps):
"""
概要: ヤコビ法を用いて連立一次方程式の解を計算します。
詳細説明:
連立一次方程式 Hx = y をヤコビ法で反復的に解きます。
Hは`Hij`、yは`yi`に対応します。
各反復ステップで、現在の近似解 `x` を用いて次の近似解 `x2` を計算します。
この計算では、`i`番目の未知数の更新に、現在の `i` 以外の未知数の値として、
インデックスが `(i+1)%n` および `(i+2)%n` のものを使用します。
前後ステップの解の差が許容誤差`eps`以下になった場合に収束と判断します。
収束しないまま最大反復回数`nmaxiter`に達した場合、未収束として処理を終了します。
計算の途中経過と最終的な収束状態は標準出力に出力されます。
対角要素 `Hij[i][i]` が0の場合、ゼロ除算エラーが発生する可能性があります。
:param Hij: list[list[float]]: 連立一次方程式の係数行列 H。H[i][j] は i行j列の要素を示します。
:param yi: list[float]: 連立一次方程式の右辺ベクトル y。yi[i] は i番目の要素を示します。
:param nmaxiter: int: 最大反復回数。この回数を超えても収束しない場合は、計算を打ち切ります。
:param eps: float: 収束判定に使用する許容誤差。各要素の絶対差がこの値より小さければ収束と判断します。
:returns: None: 計算結果（近似解および収束状況）は標準出力に出力され、関数自体は値を返しません。
"""
n = len(yi)
x = [0] * n
x2 = [0] * n
for it in range(nmaxiter):
for i in range(n):
i1 = (i+1) % n
i2 = (i+2) % n
x2[i] = (yi[i] - Hij[i][i1] * x[i1] - Hij[i][i2] * x[i2]) / Hij[i][i]
print(str(it+1) + "-th iteration")
print(" x=", x2)
IsPassed = True
for i in range(n):
if abs(x2[i] - x[i]) > eps:
IsPassed = False
break
if IsPassed:
print("Converged")
break
for i in range(n):
x[i] = x2[i]
else:
print("Not converged")
if __name__ == "__main__":
Hij = [[3, 2, 1], [1, 4, 1], [2, 2, 5]]
yi = [10, 12, 21]
eps = 1.0e-16
nmaxiter = 300
Jacobi(Hij, yi, nmaxiter, eps)