自己無撞着法による方程式の解法
プログラムの動作
equation-selfconsistent.py は、自己無撞着法（不動点反復法）を用いて、\(x = f(x)\) の形の方程式の根を数値的に探索するPythonプログラムです。このプログラムは、指定された初期値から反復計算を開始し、連続する反復における \(x\) の値の変化が許容誤差 eps 未満になった時点で収束と判断し、解を出力します。また、最大反復回数 nmaxiter を超えた場合、収束しなかったと判断します。
このプログラムの主な機能は以下の通りです。
特定の非線形方程式 \(x = 0.25 (-x^3 + x^2 - 2)\) の根を探索します。
自己無撞着法の反復プロセスを詳細に標準出力に表示します。
緩和パラメータ
kmixを使用して、反復の収束挙動を調整できます（デフォルトでは通常の不動点反復）。
原理
このプログラムは、自己無撞着法、または不動点反復法として知られる数値計算法に基づいています。この方法は、\(x = f(x)\) の形に変換できる方程式の根を見つけるために使用されます。
不動点反復法
初期値 \(x_0\) から始めて、以下の漸化式を用いて \(x\) の値を繰り返し更新します。
この反復が収束する場合、最終的な \(x\) の値は方程式 \(x = f(x)\) の根、すなわち不動点となります。収束は、連続する反復での値の変化が十分に小さくなったときに判定されます。
緩和法
プログラムでは、kmix という混合パラメータを用いて反復の更新式を一般化しています。
kmix = 1.0の場合、これは標準的な不動点反復 \(x_{i+1} = f(x_i)\) になります。0 < kmix < 1.0の場合、過小緩和（under-relaxation）と呼ばれ、収束を安定化させる効果があります。kmix > 1.0の場合、過大緩和（over-relaxation）と呼ばれ、収束を加速させる効果がありますが、不安定になる可能性もあります。
方程式
func(x) 関数で定義されている方程式は以下の通りです。
したがって、プログラムは \(x = 0.25 (-x^3 + x^2 - 2)\) の根を探索します。
収束判定
反復計算は、現在の値 \(x_i\) と次の反復値 \(x_{i+1}\) との差の絶対値が、あらかじめ設定された許容誤差 eps より小さくなったときに収束したと判断されます。
また、最大反復回数 nmaxiter を超えても収束しない場合は、計算は失敗とみなされます。
必要な非標準ライブラリとインストール方法
このプログラムは、numpy ライブラリをインポートしていますが、現在の実装ではその機能は直接使用されていません。しかし、Python標準ライブラリ以外のモジュールとして明示的にインポートされているため、将来的な拡張や依存関係を考慮してインストールしておくことが推奨されます。
numpy のインストールは、pip コマンドを使用して行います。
pip install numpy
必要な入力ファイル
equation-selfconsistent.py は、外部の入力ファイルを必要としません。初期値 (x0)、緩和パラメータ (kmix)、許容誤差 (eps)、および最大反復回数 (nmaxiter) は、プログラムのソースコード内に直接定義されています。
生成される出力ファイル
equation-selfconsistent.py は、いかなるファイルも生成しません。計算の進行状況と最終的な結果はすべて標準出力に表示されます。
各反復ステップにおける \(x\) の現在値、次の値、およびその差 (
dx) が表示されます。収束した場合は、成功メッセージと収束時の \(dx\) の値が表示されます。
最大反復回数を超えて収束しなかった場合は、失敗メッセージが表示されます。
コマンドラインでの使用例 (Usage)
equation-selfconsistent.py はコマンドライン引数を受け付けません。プログラムの実行はシンプルにPythonインタープリタを介して行います。
python equation-selfconsistent.py
コマンドラインでの具体的な使用例
以下のコマンドは、プログラムをデフォルト設定（x0 = 0.0, kmix = 1.0, eps = 1.0e-10, nmaxiter = 50）で実行します。
python equation-selfconsistent.py
上記のコマンドを実行すると、以下のような出力が標準出力に表示されます（出力の一部を抜粋）。
Solution of an equation by self-consistent method
Iter 0: x: 0.000000000000 => -0.500000000000, dx = -0.5
Iter 1: x: -0.500000000000 => -0.656250000000, dx = -0.1562
Iter 2: x: -0.656250000000 => -0.702930419922, dx = -0.04668
Iter 3: x: -0.702930419922 => -0.714243644917, dx = -0.01131
Iter 4: x: -0.714243644917 => -0.716584218960, dx = -0.002341
Iter 5: x: -0.716584218960 => -0.717079984920, dx = -0.0004958
Iter 6: x: -0.717079984920 => -0.717156947262, dx = -7.696e-05
Iter 7: x: -0.717156947262 => -0.717170562947, dx = -1.362e-05
Iter 8: x: -0.717170562947 => -0.717172778747, dx = -2.216e-06
Iter 9: x: -0.717172778747 => -0.717173169733, dx = -3.909e-07
Iter 10: x: -0.717173169733 => -0.717173238641, dx = -6.891e-08
Iter 11: x: -0.717173238641 => -0.717173250917, dx = -1.228e-08
Iter 12: x: -0.717173250917 => -0.717173253075, dx = -2.158e-09
Iter 13: x: -0.717173253075 => -0.717173253457, dx = -3.824e-10
Success: Convergence reached: dx = -3.824147754676579e-10 < eps = 1e-10
上記の出力では、各反復において \(x\) の値がどのように更新され、その変化量 \(dx\) がどの程度であるかが示されています。この例では、13回の反復で \(dx\) が許容誤差 eps = 1.0e-10 未満になったため、収束に成功し、解 \(x \approx -0.717173253457\) が見つかったことを示しています。