diffeq_euler_heun.py 技術ドキュメント

プログラムの動作

diffeq_euler_heun.py は、1階常微分方程式を数値的に解くためのPythonプログラムです。特に、与えられた常微分方程式 \(\frac{dx}{dt} = f(t, x)\) に対して、オイラー法とホイン法（Heun法）という2つの異なる数値積分手法を用いて数値解を計算します。また、厳密解が既知の場合には、それと比較して各数値解法の精度を評価します。

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

  • 常微分方程式の右辺 \(f(t, x)\) および厳密解を関数として定義。

  • オイラー法とホイン法による数値積分の実装。

  • シミュレーション結果（時刻、厳密解、オイラー解、ホイン解、各手法の誤差）をCSVファイルに出力。

  • 計算の進行状況と結果をコンソールに表示。

  • Matplotlibライブラリを使用し、時刻に対する厳密解、数値解、および各数値解法の誤差をリアルタイムでグラフ表示。

これにより、常微分方程式の解析解が複雑または存在しない場合に、効率的に近似解を求め、その精度を視覚的に評価することができます。

原理

このプログラムは、以下の1階常微分方程式の初期値問題を数値的に解きます。

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

具体的な例として、プログラムでは以下の微分方程式を扱います。

\[f(t, x) = -x^2\]

この微分方程式に対する初期条件 \(x(0) = 1.0\) の厳密解は以下の通りです。

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

プログラムは、この厳密解と、以下の2つの数値積分手法で得られる近似解を比較します。

オイラー法 (Euler method)

オイラー法は、現在の点 \((t_n, x_n)\) における微分係数 \(f(t_n, x_n)\) を一定と仮定して、次の点の値 \(x_{n+1}\) を予測する最も単純な数値積分法です。

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

ここで \(dt\) は時間ステップ幅です。

ホイン法 (Heun method)

ホイン法は、予測子-修正子法の一種であり、オイラー法よりも高い精度を持ちます。まずオイラー法で次の点の値を予測し（予測子）、その予測値と現在の点の微分係数の平均を用いて次の点を計算します（修正子）。

  1. 予測子: オイラー法で一時的な次の値を計算します。 $\(x'_{n+1} = x_n + dt \cdot f(t_n, x_n)\)$

  2. 修正子: 現在の微分係数 \(f(t_n, x_n)\) と予測された点での微分係数 \(f(t_n + dt, x'_{n+1})\) の平均を用いて、最終的な \(x_{n+1}\) を計算します。 $\(x_{n+1} = x_n + \frac{dt}{2} \cdot [f(t_n, x_n) + f(t_n + dt, x'_{n+1})]\)\( プログラムの実装では、予測子の部分を \)k_0 = dt \cdot f(t_n, x_n)\( とし、\)x_n + k_0\( を使って \)k_1 = dt \cdot f(t_n + dt, x_n + k_0)\( を計算しています。最終的な更新は \)x_{n+1} = x_n + \frac{k_0 + k_1}{2.0}$ となります。

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

本プログラムは、以下の非標準ライブラリを使用します。

  • numpy: 数値計算を効率的に行うためのライブラリです。

  • matplotlib: グラフ描画のためのライブラリです。

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

pip install numpy matplotlib

必要な入力ファイル

このプログラムは、外部からの入力ファイルを必要としません。 すべてのパラメータはプログラムのソースコード内にデフォルト値として設定されており、コマンドライン引数によってそれらの値を変更することが可能です。

生成される出力ファイル

プログラムは以下のファイルとウィンドウを生成します。

  • diffeq_euler_heun.csv:

    • ファイル形式: CSV (Comma Separated Values)

    • 内容: 各タイムステップにおける数値積分結果と厳密解を記録します。

    • ヘッダ行: t, x(cal), x(Euler), x(Heun)

      • t: 時刻

      • x(cal): 微分方程式の厳密解 \(x(t)\)

      • x(Euler): オイラー法による数値解 \(x(t)\)

      • x(Heun): ホイン法による数値解 \(x(t)\)

  • Matplotlibによるリアルタイムプロットウィンドウ:

    • 内容: シミュレーションの進行と同時に、以下の3つのグラフを動的に更新して表示します。

      1. 上段グラフ: 時刻 \(t\) に対する \(x(t)\) の変化。厳密解、オイラー法による解、ホイン法による解がそれぞれ異なる色と凡例で示されます。

      2. 中段グラフ: 時刻 \(t\) に対するオイラー法の誤差。誤差は \((x_{\text{Euler}} - x_{\text{exact}})\) として計算されます。

      3. 下段グラフ: 時刻 \(t\) に対するホイン法の誤差。誤差は \((x_{\text{Heun}} - x_{\text{exact}})\) として計算されます。

    • 終了方法: シミュレーション終了後、このウィンドウがアクティブな状態でコンソールに表示される「Press ENTER to exit>>」のプロンプトに従い、Enterキーを押すとプログラムが終了し、ウィンドウが閉じられます。

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

プログラムはコマンドライン引数を受け入れ、初期条件、ステップ幅、総ステップ数、コンソール出力間隔をカスタマイズできます。

python diffeq_euler_heun.py [x0] [dt] [nt] [iprint_interval]

引数の説明:

  • x0 (オプション): 初期条件 \(x(0)\) の値。浮動小数点数。デフォルト値は 1.0。

  • dt (オプション): 時間ステップ幅 \(\Delta t\) の値。浮動小数点数。デフォルト値は 0.1。

  • nt (オプション): 計算の総ステップ数。整数。デフォルト値は 501。

  • iprint_interval (オプション): コンソールに結果を出力するステップの間隔。整数。デフォルト値は 20。

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

例1: デフォルトパラメータでの実行

python diffeq_euler_heun.py

実行結果の説明: プログラムはデフォルトのパラメータ（\(x_0 = 1.0\), \(dt = 0.1\), \(nt = 501\), \(iprint\_interval = 20\)）で実行されます。 コンソールには、時刻 t、厳密解 x(cal)、オイラー法による解 x(euler)、ホイン法による解 x(heun) が20ステップごとに出力されます。 diffeq_euler_heun.csv というCSVファイルが生成され、すべてのステップの結果が記録されます。 同時にMatplotlibのウィンドウが開き、3つのグラフ（厳密解と数値解の比較、オイラー法の誤差、ホイン法の誤差）がリアルタイムで更新されながら表示されます。

コンソール出力例 (一部):

Solve first order diffrential equation by Heun method
  t          x(cal)      x(euler)      x(heun)
t= 0.00    1.000000      1.000000      1.000000
t= 2.00    0.333333      0.281858      0.323568
t= 4.00    0.200000      0.155702      0.189582
t= 6.00    0.142857      0.108420      0.134583
t= 8.00    0.111111      0.082710      0.104278
t=10.00    0.090909      0.066224      0.084814
...
Press ENTER to exit>>

例2: 特定のパラメータを指定して実行

初期条件 \(x_0 = 0.5\)、時間ステップ幅 \(dt = 0.05\)、総ステップ数 \(nt = 1000\) で実行する場合。

python diffeq_euler_heun.py 0.5 0.05 1000

実行結果の説明: プログラムは指定されたパラメータ（\(x_0 = 0.5\), \(dt = 0.05\), \(nt = 1000\)）で実行されます。iprint_interval はデフォルトの 20 のままです。 コンソールには同様に20ステップごとに結果が出力され、diffeq_euler_heun.csv ファイルには指定されたステップ幅と総ステップ数に応じたデータが記録されます。 Matplotlibウィンドウも同様に、指定された条件下でのシミュレーション結果をリアルタイムで表示します。

コンソール出力例 (一部):

Solve first order diffrential equation by Heun method
  t          x(cal)      x(euler)      x(heun)
t= 0.00    0.500000      0.500000      0.500000
t= 1.00    0.333333      0.299538      0.326280
t= 2.00    0.250000      0.202359      0.241517
t= 3.00    0.200000      0.157835      0.191563
t= 4.00    0.166667      0.130138      0.158730
t= 5.00    0.142857      0.110757      0.136054
...
Press ENTER to exit>>