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

プログラムの動作

diffeq_rungekutta.py は、1階の常微分方程式 $dx/dt = f(t, x)$ を数値的に解くPythonスクリプトです。具体的には、このプログラムは $dx/dt = -x^2$ という特定の方程式を対象とし、初期条件 $x(0) = 1.0$ のもとで3段3次のルンゲクッタ法を使用して数値解を求めます。計算された数値解は、同じ初期条件での厳密解 $x(t) = 1 / (1 + t)$ と比較され、その結果が標準出力に表示されます。

プログラムの主な機能は以下の通りです。

特定の常微分方程式 $dx/dt = -x^2$ を定義します。
その厳密解 $x(t) = 1 / (1 + t)$ を定義します。
初期ステップではHeun法（2次ルンゲクッタ法の一種）を使用して次の点の近似値を計算します。
以降のステップでは、3段3次のルンゲクッタ法を使用して、時間発展をシミュレーションします。
指定された時間間隔で、計算された数値解と厳密解を並べて標準出力に表示し、比較を容易にします。

原理

本プログラムは、以下の1階常微分方程式を解きます。

\[ \frac{dx}{dt} = f(t, x) \]

具体的には、$f(t, x) = -x^2$ とし、初期条件 $x(0) = 1.0$ で積分を行います。この特定の方程式の厳密解は以下の通りです。

\[ x(t) = \frac{1}{1 + t} \]

数値計算には、3段3次のルンゲクッタ法が使用されます。時間 $t_n$ における値 $x_n$ から、次の時間 $t_{n+1} = t_n + \Delta t$ における値 $x_{n+1}$ を計算するアルゴリズムは以下の通りです。

最初のステップ: プログラムの実行開始時、最初の1ステップはHeun法（修正オイラー法とも呼ばれる2次ルンゲクッタ法）で計算されます。 $$ k_A = \Delta t \cdot f(t_n, x_n) $$ k_B = \Delta t \cdot f(t_n + \Delta t, x_n + k_A) $$ x_{n+1}^* = x_n + \frac{1}{2} (k_A + k_B) $$ ここで $x_{n+1}^*$ はHeun法による1ステップ後の近似値です。
以降のステップ: ループ内の主要な計算では、以下の3段3次ルンゲクッタ法が使用されます。 $$ k_0 = \Delta t \cdot f(t_n, x_n) $$ k_1 = \Delta t \cdot f(t_n + \Delta t, x_n + k_0) $$ k_2 = \Delta t \cdot f(t_n + 2\Delta t, x_n + 4k_1 - 2k_0) $$ x_{n+1} = x_n + \frac{1}{3} (k_0 + 4k_1 + k_2) $$ ここで、$x_n$ は現在の時刻での計算値、$x_{n+1}$ は次の時刻での計算値、$k_0, k_1, k_2$ は傾きの推定値です。この式は、特定の組み合わせで3次精度を達成するよう設計されたルンゲクッタ法の一種です。

プログラムは、時間刻み幅 $\Delta t = 0.1$ で $nt = 501$ ステップの計算を実行します。計算された数値解は、厳密解と比較され、指定された出力間隔 (iprint_interval = 20 ステップごと) で標準出力に表示されます。

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

本プログラムは、数値計算のために以下のPythonライブラリを使用します。

numpy: 配列操作や高度な数値計算機能を提供します。本プログラムではインポートされていますが、直接的な関数呼び出しは行われていません。将来の拡張のために含まれている可能性があります。

numpy ライブラリは、以下の pip コマンドでインストールできます。

pip install numpy

math モジュールはPythonの標準ライブラリであるため、別途インストールする必要はありません。

必要な入力ファイル

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

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

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

出力される情報は以下の形式で、各行が時間 $t$、数値解 $x(cal)$、厳密解 $x(exact)$ を示します。

    t       x(cal)          x(exact)
t= 0.00  1.0000000000    1.0000000000
t= 2.00  0.3333333333    0.3333333333
t= 4.00  0.2000000000    0.2000000000
...

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

diffeq_rungekutta.py は、引数なしで実行されます。

python diffeq_rungekutta.py

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

以下のコマンドを実行することで、プログラムを実行し、常微分方程式の数値解と厳密解の比較結果を標準出力に表示できます。

python diffeq_rungekutta.py

上記のコマンドを実行すると、以下のような出力がコンソールに表示されます（出力の一部を抜粋）。

Solve first order diffrential equation by three-stage third-order Runge-kutta method
  t       x(cal)          x(exact)    
t= 0.00  1.0000000000    1.0000000000
t= 2.00  0.3333333333    0.3333333333
t= 4.00  0.2000000000    0.2000000000
t= 6.00  0.1428571429    0.1428571429
t= 8.00  0.1111111111    0.1111111111
t=10.00  0.0909090909    0.0909090909
t=12.00  0.0769230769    0.0769230769
t=14.00  0.0666666667    0.0666666667
t=16.00  0.0588235294    0.0588235294
t=18.00  0.0526315789    0.0526315789
t=20.00  0.0476190476    0.0476190476
t=22.00  0.0434782609    0.0434782609
t=24.00  0.0400000000    0.0400000000
t=26.00  0.0370370370    0.0370370370
t=28.00  0.0344827586    0.0344827586
t=30.00  0.0322580645    0.0322580645
t=32.00  0.0303030303    0.0303030303
t=34.00  0.0285714286    0.0285714286
t=36.00  0.0270270270    0.0270270270
t=38.00  0.0256410256    0.0256410256
t=40.00  0.0243902439    0.0243902439
t=42.00  0.0232558140    0.0232558140
t=44.00  0.0222222222    0.0222222222
t=46.00  0.0212765957    0.0212765957
t=48.00  0.0204081633    0.0204081633
t=50.00  0.0196078431    0.0196078431

出力結果を見ると、計算された数値解 x(cal) が厳密解 x(exact) と非常に高い精度で一致していることがわかります。これは、プログラムが実装している3段3次のルンゲクッタ法が適切に機能していることを示しています。