diffeq_heun.py 技術ドキュメント

プログラムの動作

diffeq_heun.py は、Pythonで記述されたプログラムであり、Heun法（修正オイラー法）を用いて1階の常微分方程式を数値的に解きます。このプログラムの主な機能は以下の通りです。

  • 微分方程式の数値解法: dx/dt = f(t, x) の形式で表現される1階常微分方程式に対し、Heun法を適用して数値解を計算します。Heun法は、オイラー法よりも高い精度を持つ予測子・修正子型の数値積分法です。

  • 特定の例題の解決: dx/dt = -x*x という特定の微分方程式を例題として採用しています。この方程式は、初期条件 x(0) = 1.0 の場合、厳密解 x(t) = 1 / (1 + t) を持ちます。

  • 厳密解との比較: 計算された数値解を、定義された厳密解と比較し、その結果を標準出力に出力します。これにより、Heun法の精度と動作を確認できます。

  • パラメータ設定: 初期条件、時間ステップ幅、総計算ステップ数などのパラメータは、プログラムのソースコード内で直接設定されます。

このプログラムは、常微分方程式の数値解法の基本的な実装と、その精度評価を示すことを目的としています。

原理

このプログラムは、1階常微分方程式 $\( \frac{dx}{dt} = f(t, x) \)\( をHeun法（修正オイラー法）で数値的に解きます。Heun法は、現在の時刻 \)t_0\( と状態変数 \)x_0\( から次の時刻 \)t_1 = t_0 + \Delta t\( における状態変数 \)x_1$ を計算する際に、以下のステップを踏みます。

  1. オイラー法による予測 (Predictor): 現在の点 \(f(t_0, x_0)\) での勾配を使用して、次のステップ \(x_{0, \text{pred}}\) をオイラー法で予測します。 $\( k_0 = \Delta t \cdot f(t_0, x_0) \)\( このとき、次の時刻における予測値は \)x_{0, \text{pred}} = x_0 + k_0$ となります。

  2. 予測された点での勾配の計算: 予測された次の点 \((t_0 + \Delta t, x_0 + k_0)\) での勾配を計算します。 $\( k_1 = \Delta t \cdot f(t_0 + \Delta t, x_0 + k_0) \)$

  3. 平均勾配による修正 (Corrector): \(k_0\) と \(k_1\) の平均値を使用して、より正確な次の状態 \(x_1\) を計算します。 $\( x_1 = x_0 + \frac{k_0 + k_1}{2} \)$

この方法により、オイラー法よりも安定性と精度が向上した数値解が得られます。プログラムでは、force(t, x) 関数が \(f(t, x) = -x*x\) を実装し、fsolution(t) 関数が厳密解 \(x(t) = 1 / (1 + t)\) を提供しています。

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

このプログラムは、Pythonの標準ライブラリのみを使用しており、numpy をインポートしていますが、その機能は実際に使用されていません。したがって、特別な非標準ライブラリのインストールは不要です。

必要な入力ファイル

このプログラムは、外部からの入力ファイルを必要としません。 すべての計算パラメータ（初期値 x0、時間ステップ dt、総ステップ数 nt など）は、プログラムのソースコード内に直接記述されています。

生成される出力ファイル

このプログラムは、いかなるファイルも生成しません。 計算結果はすべて標準出力（コンソール）にテキスト形式で表示されます。出力には、各時刻 t におけるHeun法による計算結果 x(cal) と厳密解 x(exact) が含まれます。

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

このプログラムはコマンドライン引数を取らないため、Pythonインタープリタを使用して直接実行します。

python diffeq_heun.py

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

以下のコマンドでプログラムを実行します。

python diffeq_heun.py

実行すると、標準出力に Heun法による計算結果と厳密解が、指定された表示間隔 (iprint_interval = 20) で時刻と共に表示されます。

実行結果の例:

Solve first order diffrential equation by Heun method
  t         x(cal)      x(exact)    
t= 0.00    1.000000    1.000000
t= 0.10    0.909091    0.909091
t= 2.00    0.333333    0.333333
t= 4.00    0.200000    0.200000
t= 6.00    0.142857    0.142857
t= 8.00    0.111111    0.111111
t=10.00    0.090909    0.090909
t=12.00    0.076923    0.076923
t=14.00    0.066667    0.066667
t=16.00    0.058824    0.058824
t=18.00    0.052632    0.052632
t=20.00    0.047619    0.047619
t=22.00    0.043478    0.043478
t=24.00    0.040000    0.040000
t=26.00    0.037037    0.037037
t=28.00    0.034483    0.034483
t=30.00    0.032258    0.032258
t=32.00    0.030303    0.030303
t=34.00    0.028571    0.028571
t=36.00    0.027027    0.027027
t=38.00    0.025641    0.025641
t=40.00    0.024390    0.024390
t=42.00    0.023256    0.023256
t=44.00    0.022222    0.022222
t=46.00    0.021277    0.021277
t=48.00    0.020408    0.020408
t=50.00    0.019608    0.019608

この出力では、x(cal) はHeun法によって計算された各時刻 t における状態変数 x の値、x(exact) はその時刻 t における厳密解を示しています。初期値 x(0)=1.0 から t=50.0 までの計算結果が出力されており、x(cal) と x(exact) が非常に良く一致していることが確認できます。