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

プログラムの動作

diffeq_simpson.py は、1階常微分方程式 $dx/dt = -x^2$ を数値的に解くためのPythonプログラムです。特に、このプログラムはシンプソン則の重み付けを応用したRunge-Kutta法の一種（ここでは「シンプソン法」と呼称）を用いて、与えられた初期条件 $x(0) = 1.0$ のもとで微分方程式の解を時間的に追跡します。

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

微分方程式の定義: $dx/dt = -x^2$ の右辺を関数 dxdt で定義します。
解析解の提供: この微分方程式の解析解 $x(t) = 1/(1+t)$ を関数 fsolution で計算します。
数値解法の実装:
- diffeq_euler: オイラー法による1ステップの計算。
- diffeq_heun: ホイン法（修正オイラー法）による1ステップの計算。
- diffeq_simpson: シンプソン則の重み付けを応用したRunge-Kutta法による1ステップの計算。
シミュレーション実行: main 関数は、定義された初期値と時間刻み幅を用いて、「シンプソン法」によって微分方程式を数値的に積分します。
結果の比較と出力: 計算された数値解と解析解を、指定された間隔で標準出力に整形して表示し、両者を比較できるようにします。

本プログラムは、数値積分法の原理を理解し、その精度を解析解と比較する教育的な目的にも使用できます。

原理

このプログラムは、以下の1階常微分方程式を対象としています。

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

初期条件は $x(0) = 1.0$ と設定されています。

この微分方程式の解析解は、変数分離法を用いて求めることができ、次のようになります。

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

\[ -\frac{1}{x} = -t + C \]

\[ \frac{1}{x} = t - C \]

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

初期条件 $x(0) = 1.0$ を代入すると、$1 = 1/(0 - C)$ となり、$C = -1$ が導かれます。したがって、解析解は以下の通りです。

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

プログラムでは、この解析解と数値解を比較します。数値解法としては、以下の手法が実装されています。

オイラー法 (Euler's Method): 最も基本的な数値積分法です。現在の時刻 $t_i$ における傾き $f(t_i, x_i)$ を用いて次の時刻 $t_{i+1}$ の値を予測します。 $$ x_{i+1} = x_i + \Delta t \cdot f(t_i, x_i) $$
ホイン法 (Heun's Method): 修正オイラー法とも呼ばれ、オイラー法よりも高い精度を持ちます。現在の傾きと、オイラー法で予測した次の点の傾きの平均を取ることで、より正確な予測を行います。 $$ k_1 = \Delta t \cdot f(t_i, x_i) $$ k_2 = \Delta t \cdot f(t_i + \Delta t, x_i + k_1) $$ x_{i+1} = x_i + \frac{k_1 + k_2}{2} $$
シンプソン法（類似のRunge-Kutta法）: プログラムの diffeq_simpson 関数は、Runge-Kutta法の一種として実装されていますが、その重み付けはシンプソン則に似ています。通常、シンプソン則は積分範囲を3点（両端と中央）で評価する求積法ですが、ここでは常微分方程式の数値解法にそのアイデアを適用しています。具体的には、3つの傾き $k_0, k_1, k_2$ を計算し、それらを $1/3, 4/3, 1/3$ の重みで加重平均して $x$ を更新します。

\[ k_0 = \Delta t \cdot f(t_i, x_i) \]

\[ k_1 = \Delta t \cdot f(t_i + \Delta t, x_i + k_0) \]

\[ k_2 = \Delta t \cdot f(t_i + 2\Delta t, x_i + k_0 + k_1) \]

\[ x_{i+1} = x_i + \frac{1}{3} (k_0 + 4k_1 + k_2) \]

これは標準的な4次Runge-Kutta法とは異なる形ですが、シンプソン則の重み $1, 4, 1$ を持つことから、この名称が付けられていると考えられます。

main 関数では、初期値 $x_0 = 1.0$、時間刻み幅 $dt = 0.1$、総ステップ数 $nt = 501$ が設定され、diffeq_simpson 関数を使って数値積分が実行されます。結果は iprint_interval で指定された20ステップごとに標準出力に表示されます。

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

本プログラムは以下のPythonライブラリを使用しています。

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

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

pip install numpy

標準ライブラリである math は、追加のインストールなしで使用できます。

必要な入力ファイル

このプログラムは、外部の入力ファイルを必要としません。シミュレーションの初期値 ($x0$), 時間刻み幅 ($dt$), 総ステップ数 ($nt$), および出力間隔 ($iprint_interval$) は、プログラムのソースコード内に直接定義されています。

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

このプログラムは、いかなるファイルも生成しません。すべてのシミュレーション結果は、標準出力（コンソール）に直接表示されます。出力には、時刻 t、シンプソン法による計算値 x(cal)、および解析解 x(exact) の3つの列が含まれます。

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

diffeq_simpson.py は、コマンドライン引数を取らずに直接実行することを想定しています。

基本的な実行コマンドは以下の通りです。

python diffeq_simpson.py

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

以下のコマンドを実行します。

python diffeq_simpson.py

実行結果の説明:

プログラムを実行すると、標準出力に微分方程式 $dx/dt = -x^2$ の数値解と解析解が整形された形式で表示されます。出力は、t=0.0 から t=50.0 までの範囲で、20ステップ（dt=0.1 なので $20 \times 0.1 = 2.0$ 刻み）ごとに表示されます。最初の行と、その後の iprint_interval で指定された間隔（この場合は20ステップごと）で、現在の時刻 t、シンプソン法による数値解 x(cal)、および対応する解析解 x(exact) が小数点以下6桁まで表示されます。これにより、シンプソン法の精度を解析解と比較することができます。

具体的な出力の最初の数行は以下のようになります（実際の出力はもっと長く続きます）。

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