diffeq_heun.py の技術ドキュメント

プログラムの動作

diffeq_heun.py は、Heun法（改良Euler法または予測子修正子法とも呼ばれる）を用いて1階常微分方程式を数値的に解くPythonプログラムです。このプログラムは、特定の常微分方程式 $dx/dt = -x^2$ を対象とし、初期値 $x(0)=1.0$ からの時間発展を計算します。主な機能は以下の通りです。

Heun法による数値積分: 指定された微分方程式をHeun法でステップバイステップで数値積分します。
厳密解との比較: 対象とする微分方程式の厳密解 $x(t) = 1/(1+t)$ を計算し、数値解と比較してその精度を評価します。
結果の出力: 各時間ステップ（または指定された間隔）における時間 t、Heun法による計算値 x(cal)、および厳密解 x(exact) を標準出力に表示します。

このプログラムは、Heun法による常微分方程式の数値解法の実装例を提供し、その動作と精度を確認することを目的としています。

原理

1. 対象となる微分方程式と厳密解

プログラムが解いている微分方程式は以下の通りです。 $$ \frac{dx}{dt} = -x^2 $$ 初期条件は $x(0) = 1.0$ です。この微分方程式の厳密解は、変数分離法を用いて求めることができます。 $$ \frac{dx}{-x^2} = dt $$ \int \frac{dx}{-x^2} = \int dt $$ \frac{1}{x} = t + C $$ 初期条件 $x(0) = 1.0$ を代入すると、$1/1 = 0 + C \implies C = 1$ となります。したがって、厳密解は以下の通りです。 $$ x(t) = \frac{1}{1+t} $$ プログラム内の fsolution(t) 関数がこの厳密解を計算します。

2. Heun法 (改良Euler法)

Heun法は、常微分方程式 $dx/dt = f(t, x)$ を数値的に解くための予測子修正子法の一種です。各ステップで以下の計算を行います。

予測子 (Predictor): Euler法を用いて、次のステップの仮の値を予測します。 $$ x_p = x_n + \Delta t \cdot f(t_n, x_n) $$ ここで $t_n$ と $x_n$ は現在の時間と値、$\Delta t$ は時間刻み幅です。プログラム内の `diffeq_heun` 関数では、`k0 = dt * force(t0, x0)` が $\Delta t \cdot f(t_n, x_n)$ に相当し、`x0 + k0` が予測子 $x_p$ に相当します。
修正子 (Corrector): 現在の点と予測された次の点の勾配の平均を用いて、次のステップの値を修正します。これは台形公式の考え方に基づいています。 $$ x_{n+1} = x_n + \frac{\Delta t}{2} [f(t_n, x_n) + f(t_{n+1}, x_p)] $$ プログラム内の `diffeq_heun` 関数では、`k1 = dt * force(t0+dt, x0+k0)` が $\Delta t \cdot f(t_{n+1}, x_p)$ に相当します。最終的な次のステップの値 `x1` は以下の式で計算されます。 $$ x_1 = x_0 + \frac{k_0 + k_1}{2.0} $$ この計算により、より安定で正確な数値解が得られます。

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

このプログラムは標準ライブラリのみを使用しており、numpy をインポートしていますが、現在のコードでは numpy の機能は直接利用されていません。したがって、実行に必要な非標準ライブラリはありません。

必要な入力ファイル

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

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

このプログラムは、外部にファイルを出力しません。計算結果は全て標準出力（コンソール）に表示されます。出力形式は以下の通りです。

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: 現在の時間。
x(cal): Heun法によって計算された $x$ の値。
x(exact): 微分方程式の厳密解 $x(t) = 1/(1+t)$ の値。

出力は、初期値 ($t=0.00$)、最初のステップ ($t=0.10$)、および iprint_interval（デフォルトでは20）ステップごとに表示されます。

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

本プログラムは引数を取らずに直接実行されます。

python diffeq_heun.py

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

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

python diffeq_heun.py

実行結果の説明

プログラムを実行すると、Heun法による常微分方程式の数値解が時間ステップごとに標準出力に表示されます。以下は実行結果の冒頭部分の抜粋です。

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
 ... (計算がnt=501ステップ、つまりt=50.00まで続く)
 t=50.00    0.019608    0.019608

t= 0.00 の行は初期値を示しており、Heun法による計算値と厳密解が一致していることがわかります。
t= 0.10 の行は最初の時間ステップでの結果です。
その後の行は、iprint_interval (デフォルトでは20) ステップごとに結果が出力されています。例えば、dt=0.1 なので、t=2.00 は i=20 のステップに相当し、t=4.00 は i=40 のステップに相当します。
各行において、x(cal) の値と x(exact) の値が非常に近いことが確認でき、Heun法がこの微分方程式に対して高い精度で解を求めていることを示しています。