オイラー法による常微分方程式の数値解法プログラム `diffeq_euler.py` の技術ドキュメント

プログラムの動作

diffeq_euler.py は、1階の常微分方程式 $dx/dt = force(t, x)$ をオイラー法を用いて数値的に解くためのPythonプログラムです。このプログラムの主な目的は、数値解法の結果を厳密解と比較し、その過程を可視化することにあります。

主な機能:

常微分方程式の数値解法: $dx/dt = -x^2$ という形式の1階常微分方程式をオイラー法で数値的に解きます。
厳密解との比較: 解かれた微分方程式の厳密解 $x(t) = 1/(1+t)$ (初期条件 $x(0)=1.0$ の場合) を計算し、数値解と比較します。
結果の出力:
- 計算された時間、数値解、厳密解をCSVファイル (diffeq_euler.csv) に保存します。
- 計算の進行状況をコンソールに定期的に出力します。
リアルタイム可視化: matplotlib ライブラリを使用して、数値解と厳密解の比較グラフをリアルタイムで表示し、計算が進むにつれてグラフを更新します。
柔軟な設定: 初期条件 x0、時間刻み幅 dt、総ステップ数 nt、コンソール出力間隔 iprint_interval をコマンドライン引数で指定できます。

解決する課題:

本プログラムは、数値計算の基本的なアルゴリズムであるオイラー法の理解を助け、実際の数値解が厳密解とどのように異なるか、また時間刻み幅が結果にどのように影響するかを視覚的に学ぶためのツールとして機能します。教育目的や、より複雑な微分方程式の数値解法への導入として有用です。

原理

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

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

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

この微分方程式は、初期条件 $x(0) = 1.0$ が与えられた場合、以下の厳密解を持ちます。

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

プログラム内の force(t, x) 関数が $dx/dt$ の右辺を定義し、fsolution(t) 関数がこの厳密解を計算します。

オイラー法

オイラー法は、常微分方程式の数値解法で最も単純な手法の一つです。与えられた微分方程式 $dx/dt = f(t, x)$ を、以下の差分近似式で解きます。

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

ここで：

$x_n$ は時刻 $t_n$ における $x$ の値。
$x_{n+1}$ は時刻 $t_{n+1} = t_n + \Delta t$ における $x$ の推定値。
$\Delta t$ は時間刻み幅（プログラムでは dt）。
$f(t_n, x_n)$ は時刻 $t_n$ と $x_n$ における $dx/dt$ の値（プログラムでは force(t_n, x_n)）。

この式は、現在の点 $(t_n, x_n)$ における接線の傾き $f(t_n, x_n)$ を用いて、次の点の $x$ の値を線形に近似するという考えに基づいています。プログラム内の diffeq_euler 関数がこの1ステップのオイラー法を実装しています。

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

本プログラムは、数値計算のための numpy と、グラフ描画のための matplotlib を利用します。numpy は直接的には使われていませんが、matplotlib が依存している場合があるため、両方をインストールすることをお勧めします。

これらのライブラリは、Pythonのパッケージマネージャー pip を使ってインストールできます。コマンドライン（ターミナルやコマンドプロンプト）で以下のコマンドを実行してください。

pip install numpy matplotlib

必要な入力ファイル

diffeq_euler.py は、実行時に特定の入力ファイルを必要としません。すべての設定（初期値、時間刻み幅など）は、プログラム内のデフォルト値を使用するか、コマンドライン引数として直接渡すことで行われます。

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

プログラムの実行により、以下のファイルと出力が生成されます。

1. CSVファイル

ファイル名: diffeq_euler.csv
内容: 計算された時間、数値解、厳密解のデータがカンマ区切りで保存されます。1行目はヘッダ行です。

ファイル形式の例:
```
t,x(cal),x(exact)
0.0,1.0,1.0
0.5,0.5,0.666666
1.0,0.333333,0.5
1.5,0.25,0.4
...
```
- t: 時刻。
- x(cal): オイラー法によって計算された $x$ の値。
- x(exact): 定義された微分方程式の厳密解 $x$ の値。

2. リアルタイムグラフ表示

形式: matplotlib によるインタラクティブなグラフウィンドウが表示されます。
内容: 時間 t を横軸、x(t) を縦軸として、オイラー法による数値解（euler）と厳密解（exact）が同時にプロットされます。計算が進むにつれてグラフはリアルタイムで更新され、それぞれの曲線の変化を視覚的に追うことができます。
終了方法: プログラムの計算がすべて完了した後もグラフウィンドウは開いたままで、ユーザーがコンソールで Enter キーを押すまで表示を維持します。

3. コンソール出力

計算の途中経過が、iprint_interval で指定されたステップ間隔でコンソールに出力されます。

出力形式の例:

Solve first order diffrential equation by Euler method
Write to [diffeq_euler.csv]
  t        x(cal)      x(exact)
t= 0.50    0.500000    0.666667
t= 10.00    0.090909    0.090909
t= 20.00    0.047619    0.047619
...

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

diffeq_euler.py は、以下の形式でコマンドラインから実行できます。引数はオプションであり、指定しない場合はプログラム内のデフォルト値が使用されます。

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

引数の説明:

x0 (float, オプション):
- 初期条件 $x(0)$ の値。
- デフォルト値: 1.0
dt (float, オプション):
- 時間刻み幅 $\Delta t$ 。オイラー法の精度に影響します。
- デフォルト値: 0.5
nt (int, オプション):
- 計算する総ステップ数。
- デフォルト値: 501 (0からnt-1までループするため、ntステップ)
iprint_interval (int, オプション):
- コンソールに計算経過を出力する間隔のステップ数。
- デフォルト値: 20

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

1. デフォルト設定での実行

引数を何も指定せずに実行します。この場合、x0=1.0, dt=0.5, nt=501, iprint_interval=20 のデフォルト値が使用されます。

python diffeq_euler.py

実行結果の概要:

コンソールには、約20ステップごとに計算された $t$ , $x(cal)$ , $x(exact)$ の値が出力されます。
diffeq_euler.csv ファイルが生成され、すべての計算ステップにおける $t$ , $x(cal)$ , $x(exact)$ の値が記録されます。
matplotlib のグラフウィンドウが開き、オイラー法による計算値と厳密解がリアルタイムでプロットされ、更新されていきます。初期条件 $x(0)=1.0$ の場合、厳密解と数値解が比較されます。

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

Solve first order diffrential equation by Euler method
Write to [diffeq_euler.csv]
  t        x(cal)      x(exact)
t= 0.50    0.500000    0.666667
t= 10.00    0.090909    0.090909
t= 20.00    0.047619    0.047619
t= 250.00    0.003984    0.003984
Press ENTER to exit>>

diffeq_euler.csv の内容例 (一部):

t,x(cal),x(exact)
0,1.0,1.0
5,0.5,0.6666666666666666
0,0.3333333333333333,0.5
5,0.25,0.4
0,0.2,0.3333333333333333
...

2. 特定の引数を指定して実行

初期値を $x0=2.0$ , 時間刻み幅を $dt=0.1$ , 総ステップ数を $nt=100$ に設定して実行します。

python diffeq_euler.py 2.0 0.1 100