diffeq_euler.py 技術ドキュメント

このドキュメントは、Pythonプログラム diffeq_euler.py の技術的な側面を解説します。

プログラムの動作

diffeq_euler.py は、一階常微分方程式をオイラー法を用いて数値的に解くプログラムです。具体的には、初期値問題 \(dx/dt = -x^2\) を対象とし、その解析解 \(x(t) = 1/(1+t)\) (初期条件 \(x(0)=1.0\) の場合) と比較しながら計算結果を可視化します。

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

  • 微分方程式の数値解法: オイラー法により、\(dx/dt = f(t, x)\) の形式の微分方程式を時間発展させます。

  • 解析解との比較: プログラム内部で定義された解析解と数値計算結果を比較し、その差を確認できます。

  • リアルタイムグラフ表示: matplotlib を使用して、計算の進行とともに数値解と解析解の軌跡をリアルタイムでグラフに表示します。

  • 結果のファイル出力: 計算された時間 \(t\)、数値解 \(x(cal)\)、解析解 \(x(exact)\) のデータをCSV形式でファイルに保存します。

  • コマンドライン引数による設定: 初期値 \(x(0)\)、時間刻み \(dt\)、総ステップ数 \(nt\)、コンソールへの出力間隔 iprint_interval をコマンドライン引数で指定できます。

このプログラムは、数値解析の基本的な手法であるオイラー法の理解と実践、および常微分方程式の挙動の視覚化に役立ちます。

原理

このプログラムは、次の形式の一階常微分方程式を解きます。

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

具体的にこのプログラムが対象とする微分方程式は、\(f(t, x) = -x^2\) であるため、次のようになります。

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

この初期値問題の解析解は、初期条件 \(x(0) = x_0\) の場合、次の式で与えられます。

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

特に、プログラムのデフォルト設定である \(x_0 = 1.0\) の場合、解析解は以下のようになります。

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

数値計算にはオイラー法が用いられます。オイラー法は、微分方程式を次の差分方程式で近似する最も基本的な数値解法の一つです。

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

ここで、\(x_n\) は時刻 \(t_n\) における近似解、\(x_{n+1}\) は時刻 \(t_{n+1} = t_n + h\) における近似解、\(h\) は時間刻み (プログラム中の dt に相当) です。この方法では、\(f(t_n, x_n)\) が点 \((t_n, x_n)\) における曲線の傾きを表し、これを使って次の点の値を予測します。

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

このプログラムは、データ処理とグラフ描画のために以下の非標準ライブラリを使用します。

  • numpy: 数値計算を効率的に行うためのライブラリですが、このプログラムでは直接的な大規模な数値演算には使用されていません。しかし、多くの科学技術計算プログラムで標準的に使用されるため、他の依存関係で必要になる可能性があります。

  • matplotlib: グラフの描画とリアルタイム表示のために使用されます。

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

pip install numpy matplotlib

必要な入力ファイル

diffeq_euler.py は、外部からの入力ファイルを必要としません。全ての初期設定値はプログラム内部で定義されており、必要に応じてコマンドライン引数で上書きすることができます。

生成される出力ファイル

プログラムを実行すると、以下の出力が生成されます。

  1. CSVファイル: 計算された時系列データが diffeq_euler.csv というファイル名で保存されます。

    • ファイル名: diffeq_euler.csv

    • 内容: 各行には、計算時間 t、オイラー法で計算された数値解 x(cal)、および解析解 x(exact) がカンマ区切りで記録されます。ヘッダ行も含まれます。

    ファイルの例：

    t,x(cal),x(exact)
0.5,0.5,0.6666666666666666
1.0,0.375,0.5
1.5,0.3046875,0.4
2.0,0.26048583984375,0.3333333333333333
...

  2. リアルタイムグラフウィンドウ: matplotlib を用いて、計算の進行状況を示すグラフがリアルタイムで表示されます。横軸に時間 \(t\)、縦軸に \(x(t)\) の値をとり、数値解（euler）と解析解（exact）が異なる色でプロットされます。

  3. コンソール出力: 計算の進行状況と、指定された iprint_interval ごとの時間 \(t\)、数値解 \(x(cal)\)、解析解 \(x(exact)\) の値が標準出力に表示されます。

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

diffeq_euler.py は、以下の形式でコマンドラインから実行できます。

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

  • [x0]: 初期値 \(x(0)\) (デフォルト: 1.0)

  • [dt]: 時間刻み \(\Delta t\) (デフォルト: 0.5)

  • [nt]: 総ステップ数 (計算回数) (デフォルト: 501)

  • [iprint_interval]: コンソールに結果を出力する間隔のステップ数 (デフォルト: 20)

引数は全てオプションであり、指定しない場合はデフォルト値が使用されます。引数は左から順に指定する必要があります。例えば、dt を変更したい場合は x0 も指定する必要があります。

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

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

引数を何も指定せずに実行します。

python diffeq_euler.py

  • 実行結果:

    • x0=1.0, dt=0.5, nt=501, iprint_interval=20 が適用されます。

    • コンソールには、初期状態と20ステップごとの計算結果が表示されます。

    • diffeq_euler.csv ファイルが生成され、計算された全データが保存されます。

    • リアルタイムでグラフが表示され、数値解と解析解がプロットされます。

    コンソール出力の例：

    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.076814    0.090909
t=20.00     0.047806    0.047619
t=30.00     0.033785    0.032258
t=40.00     0.025983    0.024390
...
Press ENTER to exit>>

2. 時間刻みを小さくして精度を確認

時間刻み dt を 0.1 に変更し、総ステップ数 nt を 2501 に増やすことで、より長い時間範囲で精度を改善します。

python diffeq_euler.py 1.0 0.1 2501

  • 実行結果:

    • x0=1.0, dt=0.1, nt=2501, iprint_interval=20 が適用されます。

    • dt が小さくなったことで、オイラー法の数値解が解析解に近づき、グラフ上での両者の乖離が減少します。

    • 計算時間が長くなり、グラフの更新も細かくなります。

    • diffeq_euler.csv のデータ行数が増加します。

    コンソール出力の例（一部抜粋）：

    Solve first order diffrential equation by Euler method
Write to [diffeq_euler.csv]
  t       x(cal)      x(exact)
t= 0.10     0.900000    0.909091
t= 2.00     0.320490    0.333333
t= 4.00     0.200927    0.200000
t= 6.00     0.143859    0.142857
t= 8.00     0.111246    0.111111
...
Press ENTER to exit>>

3. 初期値と出力間隔を変更

初期値 x0 を 2.0 に変更し、コンソール出力間隔 iprint_interval を 50 に設定します。

python diffeq_euler.py 2.0 0.5 501 50

  • 実行結果:

    • x0=2.0, dt=0.5, nt=501, iprint_interval=50 が適用されます。

    • 解析解は \(x(t) = 1/(t + 0.5)\) となり、初期値が変更されたグラフが表示されます。

    • コンソールへの出力頻度が減ります。

    コンソール出力の例：

    Solve first order diffrential equation by Euler method
Write to [diffeq_euler.csv]
  t       x(cal)      x(exact)
t= 2.50     0.285854    0.285714
t=25.00     0.038466    0.038462
t=50.00     0.019608    0.019608
t=75.00     0.013245    0.013245
t=100.00     0.009950    0.009950
...
Press ENTER to exit>>