Jacobi.py 技術ドキュメント

プログラムの動作

Jacobi.py は、連立一次方程式 \(H \mathbf{x} = \mathbf{y}\) を、ヤコビ反復法（Jacobi iteration method）を用いて数値的に解くためのPythonプログラムです。特に、このプログラムの更新式は、連立方程式の次元 \(n\) が3の場合に一般的なヤコビ反復法と等価となります。

主な機能は以下の通りです。

係数行列 \(H\) と右辺ベクトル \(\mathbf{y}\) を入力として受け取ります。
指定された最大反復回数 nmaxiter または収束条件が満たされるまで、反復計算を実行します。
各反復計算のステップで、現在の解ベクトル \(\mathbf{x}\) の値を標準出力に表示します。
隣接する反復ステップ間での解ベクトルの変化が、許容誤差 eps 以下になった場合に収束と判断し、計算を終了します。
最大反復回数に達しても収束しない場合は、その旨を通知します。

このプログラムは、線形代数方程式の近似解を効率的に求めることを目的としています。

原理

本プログラムは、連立一次方程式 \(H \mathbf{x} = \mathbf{y}\) をヤコビ反復法で解きます。ここで、\(H\) は \(n \times n\) の係数行列、\( \mathbf{x} \) は \(n\) 次元の未知数ベクトル、\( \mathbf{y} \) は \(n\) 次元の右辺ベクトルです。

ヤコビ反復法では、各未知数 \(x_i\) を他の未知数の前回の反復値を用いて独立に更新します。プログラムで実装されている更新式は以下の通りです。

\[x_i^{(k+1)} = \frac{1}{H_{ii}} \left( y_i - H_{i,(i+1)\%n} x_{(i+1)\%n}^{(k)} - H_{i,(i+2)\%n} x_{(i+2)\%n}^{(k)} \right)\]

ここで、\(x_i^{(k+1)}\) は \(k+1\) 回目の反復における \(i\) 番目の未知数、\(x_j^{(k)}\) は \(k\) 回目の反復における \(j\) 番目の未知数、\(H_{ij}\) は行列 \(H\) の \(i\) 行 \(j\) 列の要素、\(y_i\) はベクトル \(\mathbf{y}\) の \(i\) 番目の要素です。添字の (i+1)%n と (i+2)%n は、リストの範囲内で循環するインデックスを意味します。

この更新式は、特に連立方程式の次元 \(n=3\) の場合に、一般的なヤコビ反復法の更新式

\[x_i^{(k+1)} = \frac{1}{H_{ii}} \left( y_i - \sum_{j \neq i} H_{ij} x_j^{(k)} \right)\]

と一致します。プログラムは初期解ベクトル \(\mathbf{x}^{(0)}\) をゼロベクトルとし、指定された最大反復回数までこの更新を繰り返します。

収束判定は、各未知数の現在の値と前回の値との絶対差が、指定された許容誤差 eps (\(|x_i^{(k+1)} - x_i^{(k)}| < \text{eps}\)) を下回る場合に収束したと判断します。全ての未知数についてこの条件が満たされた場合、反復計算は終了します。

必要な非標準ライブラリとインストール方法

本プログラム Jacobi.py はPythonの標準ライブラリのみを使用しており、特別な非標準ライブラリを必要としません。そのため、追加のインストール作業は不要です。

必要な入力ファイル

本プログラムは、外部の入力ファイルを必要としません。連立方程式の係数行列 \(H\) (Hij)、右辺ベクトル \(\mathbf{y}\) (yi)、収束判定の許容誤差 eps、および最大反復回数 nmaxiter は、プログラムのソースコード内に直接記述されています。

Hij: 係数行列 \(H\) を表す二次元リスト。Hij[i][j] が \(H_{ij}\) に対応します。例： Hij = [[3, 2, 1], [1, 4, 1], [2, 2, 5]]
yi: 右辺ベクトル \(\mathbf{y}\) を表す一次元リスト。yi[i] が \(y_i\) に対応します。例： yi = [10, 12, 21]
eps: 収束判定の許容誤差を表す浮動小数点数。例： eps = 1.0e-16
nmaxiter: 最大反復回数を表す整数。例： nmaxiter = 300

これらの値は、if __name__ == "__main__": ブロック内で定義されており、プログラム実行時に直接使用されます。

生成される出力ファイル

本プログラム Jacobi.py は、いかなる出力ファイルも生成しません。計算結果はすべて標準出力（コンソール）に直接表示されます。

出力される主な内容は以下の通りです。

各反復ステップ（イテレーション）の回数。
その反復ステップで計算された現在の解ベクトル \(\mathbf{x}\) の値。
収束条件が満たされた場合に表示される「Converged」メッセージ。
最大反復回数に達しても収束しなかった場合に表示される「Not converged」メッセージ。

コマンドラインでの使用例 (Usage)

Jacobi.py プログラムは、コマンドライン引数を受け付けません。直接Pythonインタープリタで実行するだけです。

基本的な実行コマンドは以下の通りです。

python Jacobi.py

コマンドラインでの具体的な使用例

以下のコマンドで Jacobi.py を実行します。

python Jacobi.py

プログラムに設定されているデフォルト値は以下の通りです。

係数行列 Hij: [[3, 2, 1], [1, 4, 1], [2, 2, 5]]
右辺ベクトル yi: [10, 12, 21]
許容誤差 eps: 1.0e-16
最大反復回数 nmaxiter: 300

この実行により、ヤコビ反復法による解の計算過程と、最終的な収束結果が標準出力に表示されます。

実行結果例:

1-th iteration
  x= [3.3333333333333335, 3.0, 4.2]
2-th iteration
  x= [1.6, 1.875, 2.533333333333333]
3-th iteration
  x= [2.8583333333333335, 1.9666666666666666, 2.37]
4-th iteration
  x= [2.011111111111111, 2.059722222222222, 2.158333333333333]
5-th iteration
  x= [2.6075925925925927, 2.000347222222222, 2.08375]
6-th iteration
  x= [2.2185185185185187, 2.013449074074074, 2.028046296296296]
7-th iteration
  x= [2.4939012345679013, 2.001962962962963, 2.0084537037037037]
8-th iteration
  x= [2.304020987654321, 2.003185185185185, 2.0029302469135804]
9-th iteration
  x= [2.428549382716049, 2.0009953703703703, 2.0011506172839507]
10-th iteration
  x= [2.3456249999999997, 2.001099537037037, 2.000454320987654]
...
27-th iteration
  x= [2.0000000000000004, 2.0000000000000004, 2.0000000000000004]
28-th iteration
  x= [2.0, 2.0, 2.0]
Converged

上記の出力例では、28回目の反復で収束条件を満たし、最終的な解ベクトル \(\mathbf{x} \approx [2.0, 2.0, 2.0]^T\) が得られたことを示しています。