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

プログラムの動作

diffeq_simpson.py は、1階常微分方程式 $dx/dt = f(t, x)$ を数値的に解くためのPythonプログラムです。特に、このプログラムは $dx/dt = -x^2$ という具体的な方程式を対象とし、その初期条件 $x(0) = 1.0$ のもとでの解を求めます。計算された数値解は、厳密解 $x(t) = 1/(1+t)$ と比較され、その結果が標準出力に表形式で表示されます。

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

常微分方程式の定義: $f(t, x) = -x^2$ を dxdt 関数として定義します。
厳密解の定義: $x(t) = 1/(1+t)$ を fsolution 関数として定義します。
数値積分法の選択:
- オイラー法 (diffeq_euler)
- ホイン法 (diffeq_heun - コメントアウトされており、デフォルトでは使用されません)
- シンプソン法（と名付けられたルンゲ・クッタ型の手法） (diffeq_simpson) が実装されています。
シミュレーション実行: main 関数内で、最初のステップはオイラー法で計算され、その後は diffeq_simpson 関数によって各時間ステップの数値解が計算されます。
結果の出力: 計算された時間 $t$ 、数値解 $x(cal)$ 、および厳密解 $x(exact)$ が、指定された時間間隔で標準出力に整形されて表示されます。

このプログラムは、数値積分法の基本的な概念と、その精度を厳密解と比較して評価するデモンストレーションとして機能します。

原理

このプログラムは、1階常微分方程式の初期値問題 $$ \frac{dx}{dt} = f(t, x) $$ を数値的に解きます。ここでは、具体的な方程式として $f(t, x) = -x^2$ を使用し、初期条件 $x(0) = 1.0$ を与えています。この方程式の厳密解は $$ x(t) = \frac{1}{1+t} $$ です。

プログラムでは、いくつかの異なる数値積分法が実装されています。

オイラー法 (Euler's Method)

最も基本的な数値積分法です。現在の点 $(t_n, x_n)$ における傾き $f(t_n, x_n)$ を用いて次のステップ $x_{n+1}$ を近似します。 $$ x_{n+1} = x_n + \Delta t \cdot f(t_n, x_n) $$ ここで $\Delta t$ は時間刻み幅です。

ホイン法 (Heun's Method)

オイラー法を改良した予測子-修正子法の一種で、より高い精度を持ちます。 $$ 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) $$ このプログラムでは実装されていますが、デフォルトでは使用されません。

シンプソン法 (Simpson Method - 実装されている形式)

このプログラムの diffeq_simpson 関数で実装されている手法は、ルンゲ・クッタ法の一種として構築されており、次のように定義されます。 $$ 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 + k_0 + k_1) $$ x_{n+1} = x_n + \frac{1}{3}(k_0 + 4k_1 + k_2) $$ この形式は、一般的なシンプソン積分公式を直接適用する予測子-修正子法とは異なり、ルンゲ・クッタ法における $k$ の定義に基づいて次のステップの値を計算します。プログラムの実行では、最初のステップのみオイラー法が使われ、その後の全てのステップでこの「シンプソン法」が適用されます。

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

このプログラムはPythonの標準ライブラリのみを使用しており、特別な非標準ライブラリのインストールは不要です。

numpy: プログラム内で import numpy as np がありますが、np の機能は実際には使用されていません。したがって、numpy はこのプログラムの実行には必須ではありません。
math: Pythonの標準ライブラリであり、exp, sqrt, sin, cos, pi がインポートされています。fsolution 関数で exp が使用されています。標準ライブラリであるため、追加のインストールは不要です。

必要な入力ファイル

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

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

diffeq_simpson.py は、いかなるファイルも生成しません。計算結果は全て標準出力（コンソール）に直接表示されます。出力は、時間 $t$、計算された数値解 $x(cal)$、および厳密解 $x(exact)$ の3列からなる表形式です。

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

diffeq_simpson.py は、引数なしで直接Pythonインタプリタによって実行されます。

python diffeq_simpson.py

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

以下のコマンドを実行すると、プログラムが定義された常微分方程式の数値解を計算し、その結果を標準出力に出力します。

python diffeq_simpson.py

実行結果の例:

Solve first order diffrential equation by Simpson method
  t       x(cal)      x(exact)  
t= 0.00   1.000000    1.000000
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

この出力は、プログラム内で設定された iprint_interval (20ステップごと) に基づいて、時間 $t$ 、計算された $x$ の値、および厳密解の $x$ の値が並んで表示されていることを示しています。この例では、シンプソン法（と名付けられた手法）が厳密解と非常に高い精度で一致していることがわかります。