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

プログラムの動作

diffeq_rungekutta.py は、初期値問題を持つ一次常微分方程式 \(dx/dt = f(t, x)\) の数値解を三段階三次のルンゲクッタ法を用いて計算するPythonプログラムです。このプログラムは、指定された時間ステップと総ステップ数にわたって方程式を積分し、計算された数値解を既知の厳密解と比較して標準出力に表示します。

主な機能は以下の通りです。

常微分方程式の数値積分: \(dx/dt = -x^2\) という形式の一次常微分方程式を数値的に解きます。
三段階三次ルンゲクッタ法: 高精度な数値解を求めるための主要なアルゴリズムとして三段階三次のルンゲクッタ法を採用しています。
初期ステップの予測: ルンゲクッタ法の初期値設定として、ホイン法（修正オイラー法）を使用しています。オイラー法もコメントアウトで用意されています。
厳密解との比較: 対象とする方程式の厳密解がプログラム内に定義されており、計算された数値解との比較をリアルタイムで出力します。
パラメータのカスタマイズ: 初期条件 (\(x_0\))、時間ステップ (\(dt\))、総ステップ数 (\(nt\))、出力間隔 (\(iprint\_interval\)) はプログラムコード内で設定されています。

このプログラムは、解析的に解くのが難しい常微分方程式の近似解を効率的かつ精度よく求める手法を実演し、その結果を厳密解と比較することで数値解析手法の有効性を示します。

原理

本プログラムは、一次常微分方程式 \(dx/dt = f(t, x)\) の初期値問題 \(x(t_0) = x_0\) を数値的に解くために、以下の手法と数式を使用しています。

対象とする常微分方程式と厳密解

本プログラムが対象とする常微分方程式は以下の通りです。

\[ \frac{dx}{dt} = -x^2 \]

初期条件を \(x(0) = 1.0\) とした場合のこの方程式の厳密解は以下のようになります。

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

オイラー法

最も基本的な数値積分法で、次のステップの値を現在の傾きを用いて予測します。

\[ x_{n+1} = x_n + dt \cdot f(t_n, x_n) \]

ホイン法 (修正オイラー法)

オイラー法よりも精度が高い予測子-修正子法の一種です。まずオイラー法で仮の次のステップを予測し、その予測値と現在の値の両方から傾きを計算して平均を取ります。

\[ k_0 = dt \cdot f(t_n, x_n) \]

\[ k_1 = dt \cdot f(t_n + dt, x_n + k_0) \]

\[ x_{n+1} = x_n + \frac{1}{2}(k_0 + k_1) \]

三段階三次ルンゲクッタ法 (Kutta's third-order method)

本プログラムの主要な数値積分アルゴリズムです。ホイン法よりもさらに高い精度で次のステップの値を計算します。プログラムの main 関数内で実装されている形式は、Kuttaの三次法として知られるルンゲクッタ法の一種です。

現在の時刻を \(t_n\)、対応する変数の値を \(x_n\) とすると、次の時刻 \(t_{n+1} = t_n + dt\) での \(x_{n+1}\) は以下のように計算されます。

\[ k_1 = dt \cdot f(t_n, x_n) \]

\[ k_2 = dt \cdot f(t_n + dt, x_n + k_1) \]

\[ k_3 = dt \cdot f(t_n + 2dt, x_n + 4k_2 - 2k_1) \]

\[ x_{n+1} = x_n + \frac{1}{3}(k_1 + 4k_2 + k_3) \]

プログラム内では \(k_1, k_2, k_3\) がそれぞれ k0, k1, k2 という変数名で実装されています。

プログラムの動作としては、まず初期時刻 \(t=0\) での \(x_0\) から \(t=dt\) での \(x_1\) をホイン法で計算します。その後、ループ内で各ステップ \(t_i\) から \(t_{i+1}\) への \(x\) の更新に三段階三次ルンゲクッタ法を適用し、結果を出力します。

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

このプログラムは、数値計算ライブラリ numpy をインポートしていますが、直接その機能を使用している箇所はありません。ただし、標準ライブラリの math モジュールは使用しています。

numpy:
- インストール方法: Pythonのパッケージマネージャ pip を使用してインストールします。
```
pip install numpy
```

必要な入力ファイル

このプログラムは、外部の入力ファイルを必要としません。すべてのシミュレーションパラメータ（初期値 \(x_0\)、時間ステップ \(dt\)、総ステップ数 \(nt\)、出力間隔 \(iprint\_interval\)）は、プログラムコード内に直接ハードコードされています。

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

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

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
...

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

diffeq_rungekutta.py を実行するには、Pythonインタープリタを使用してスクリプトファイルを直接呼び出します。特別なコマンドライン引数は必要ありません。

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

上記は、プログラム内で設定されている以下のパラメータに基づいた実行結果です。

x0 = 1.0 (初期時刻 \(t=0\) における \(x\) の値)
dt = 0.1 (時間ステップ幅)
nt = 501 (総ステップ数、これにより \(t\) は \(0.0\) から \(50.0\) まで変化します)
iprint_interval = 20 (20ステップごとに結果が出力されます)

出力の各行は、時間 t、三段階三次ルンゲクッタ法によって計算された \(x\) の値 x(cal)、および厳密解 x(exact) を示しています。厳密解と計算値がほぼ一致していることから、ルンゲクッタ法の高い精度が確認できます。

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