diffeq_rungekutta4.py テクニカルドキュメント

プログラムの動作

diffeq_rungekutta4.py は、1階常微分方程式を数値的に解くPythonプログラムです。具体的には、初期値問題として与えられた常微分方程式 \(dx/dt = f(t, x)\) を、4段4次ルンゲ＝クッタ法（のバリアントである予測子-修正子型多段階法）を用いて近似的に解きます。

このプログラムは以下の機能を提供します。

  • 特定の微分方程式 \(dx/dt = -x^2\) を対象とします。この方程式には解析解 \(x(t) = 1/(1+t)\) が存在し、初期条件 \(x(0) = 1.0\) の下でこの解析解と比較しながら数値計算の精度を確認できます。

  • 時間ステップごとに近似解を計算し、解析解と比較した結果を標準出力（コンソール）に表示します。

  • 計算の初期ステップではヘウン法を用いて次ステップの値を予測し、その後のステップでは多段階法に似たルンゲ＝クッタ型の計算式で次ステップの値を計算します。

  • 計算は指定された時間ステップ数だけ実行され、定期的に計算結果がコンソールに出力されます。

原理

このプログラムは、以下の原理に基づいて1階常微分方程式を数値的に解きます。

1階常微分方程式

一般に、1階常微分方程式は以下の形式で表されます。

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

このプログラムでは、\(f(t, x) = -x^2\) という具体的な関数を対象としています。

解析解

初期条件 \(x(0) = 1.0\) の場合、\(dx/dt = -x^2\) の解析解は以下のようになります。

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

プログラムはこの解析解と数値計算結果を比較することで、近似の精度を評価します。

数値積分アルゴリズム

プログラムは、以下に示す数値積分アルゴリズムの組み合わせを使用しています。

1. 初期ステップの予測（ヘウン法）

メインループに入る前に、最初のステップ \(x_1\) はヘウン法（Heun's method）によって予測されます。ヘウン法は、オイラー法を発展させた予測子-修正子型アルゴリズムであり、次のステップ \(x_{n+1}\) を以下の式で計算します。

\[k_0 = \Delta t \cdot f(t_n, x_n)\]
\[k_1 = \Delta t \cdot f(t_n + \Delta t, x_n + k_0)\]
\[x_{n+1} = x_n + \frac{1}{2} (k_0 + k_1)\]

コードでは x1 = diffeq_heun(dxdt, t0, x0, dt) として使用されています。

2. 多段階予測子-修正子型ルンゲ＝クッタ法

プログラムの主要な数値計算はループ内で行われます。コードのコメントでは「four-stage fourth-order Runge-kutta method」とされていますが、その実装は標準的な古典的4次ルンゲ＝クッタ法（RK4）とは異なる形式であり、多段階予測子-修正子法のような構造を持っています。具体的には、\(x_{n+2}\) (コード内の x2) を計算するために、以下の一連の式が用いられます。

\[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 + \Delta t, x_n + k_1)\]
\[k_3 = \Delta t \cdot f(t_n + 2\Delta t, x_n + 2k_2)\]
\[x_{n+2} = x_n + \frac{1}{3} (k_0 + 2k_1 + 2k_2 + k_3)\]

ここで \(t_n\) は現在の時刻 t0、\(x_n\) は現在の近似解 x0 に対応します。計算された \(x_{n+2}\) は次のステップで \(x_{n+1}\) (コード内の x1) の値として更新され、さらにその次のステップでは \(x_n\) (コード内の x0) の値として利用されます。このように、過去のステップの情報 (x0, x1 のシフト) を用いて次のステップを計算する形式は、多段階法の特徴を示しています。

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

このプログラムは以下の非標準ライブラリに依存しています。

  • numpy: 数値計算を効率的に行うためのライブラリです。本プログラムではインポートされていますが、実際には使用されていません。

インストール方法

numpy がお使いの環境にインストールされていない場合は、以下の pip コマンドでインストールできます。

pip install numpy

必要な入力ファイル

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

生成される出力ファイル

diffeq_rungekutta4.py は、実行時にいかなるファイルも生成しません。 すべての計算結果は標準出力（コンソール）に出力されます。出力形式は以下の通りです。

       t            x(cal)        x(exact)
t=    0.00  1.0000000000e+00  1.0000000000e+00
t=100000.00  XXXXXXXXXXXXe+XX  9.9999000000e-06
...

  • t: 現在の時刻。

  • x(cal): プログラムによって計算された近似解。

  • x(exact): 微分方程式の解析解。

これらの値は、iprint_interval で指定された間隔ごとに表示されます。

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

diffeq_rungekutta4.py は引数を取らないため、単純にPythonインタープリタで実行します。

python diffeq_rungekutta4.py

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

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

python diffeq_rungekutta4.py

実行すると、プログラムは常微分方程式の数値積分を開始し、計算の進行状況と結果を標準出力に表示します。 初期設定では dt = 0.1、nt = 10000001、iprint_interval = 1000000 となっています。

実行結果の例:

Solve first order diffrential equation by four-stage fourth-order Runge-kutta method
       t            x(cal)        x(exact)
t=    0.00  1.0000000000e+00  1.0000000000e+00
t=100000.00  9.9999000000e-06  9.9999000000e-06
t=200000.00  4.9999750001e-06  4.9999750001e-06
t=300000.00  3.3333222222e-06  3.3333222222e-06
t=400000.00  2.5000000000e-06  2.5000000000e-06
t=500000.00  2.0000000000e-06  2.0000000000e-06
t=600000.00  1.6666666667e-06  1.6666666667e-06
t=700000.00  1.4285714286e-06  1.4285714286e-06
t=800000.00  1.2500000000e-06  1.2500000000e-06
t=900000.00  1.1111111111e-06  1.1111111111e-06
t=1000000.00  1.0000000000e-06  1.0000000000e-06

上記はプログラムの初期出力行と、iprint_interval (1000000) ごとに計算される結果の一部を示しています。最終的な時刻 t は (nt - 1) * dt = (10000001 - 1) * 0.1 = 1000000.0 となり、その時点での計算値と解析解が表示されます。この例では、計算された近似解 x(cal) が解析解 x(exact) と非常に良い一致を示していることが分かります。