`diffeq2nd_euler.py` プログラムの技術ドキュメント

プログラムの動作

diffeq2nd_euler.py は、2階常微分方程式を数値的に解くためのPythonスクリプトです。具体的には、初期値問題として与えられた2階常微分方程式を、最も基本的な数値積分法であるオイラー法（Euler method）を用いて解きます。

このプログラムは以下の機能を提供します。

2階常微分方程式の数値積分: $d^2x/dt^2 = f(t, x)$ の形式の2階常微分方程式を、1階連立微分方程式に変換し、オイラー法でステップごとに時間発展を計算します。
特定の物理現象のシミュレーション: デフォルトでは、$d^2x/dt^2 = -x$ という単振動の方程式を解きます。これは、質量-ばね系や振り子の小振幅振動など、多くの物理現象に見られる基本的なモデルです。
初期条件の設定: $x(0)$ と $v(0) = dx/dt(0)$ の初期条件を設定できます。デフォルトでは $x(0)=0.0, v(0)=1.0$ が使用されます。
厳密解との比較: 計算された数値解と、定義された厳密解（今回の場合は $x(t) = \sin(t)$）を比較し、精度を確認できます。
結果の出力: 計算された時間 $t$、数値解 $x(t)$、厳密解 $x(t)$、および数値解 $v(t)$ をCSV形式のファイルに保存します。また、計算途中の特定ステップでの結果を標準出力にも表示します。
パラメータの柔軟な設定: 時間刻み $dt$、総ステップ数 $nt$、および画面出力間隔 iprint_interval をコマンドライン引数で指定することで、シミュレーション条件を柔軟に変更できます。

このプログラムは、数値計算の入門や、より複雑な物理シミュレーションの基礎を理解するための教育的なツールとして利用できます。

原理

このプログラムは、以下の手順で2階常微分方程式をオイラー法で解きます。

2階微分方程式の定義: プログラムは、以下の形式の2階常微分方程式を扱います。 $$ \frac{d^2x}{dt^2} = f(t, x) $$ 今回のプログラムでは、`force(t, x)` 関数が $f(t, x)$ を実装しており、デフォルトでは $f(t, x) = -x$ と定義されています。
1階連立微分方程式への変換: 2階常微分方程式を直接オイラー法で解くのではなく、まず1階の連立微分方程式に変換します。新しい変数 $v = dx/dt$ を導入すると、元の2階微分方程式は次のように2つの1階微分方程式に分解されます。 $$ \frac{dx}{dt} = v $$ \frac{dv}{dt} = f(t, x) $$
オイラー法による数値積分: 上記で得られた1階の連立微分方程式を、オイラー法を用いて数値的に時間発展させます。オイラー法は、各時間ステップにおいて現在の値と導関数を使って次のステップの値を近似します。時間刻みを $dt$ とすると、現在の時刻 $t_n$ での $x_n$ と $v_n$ から、次の時刻 $t_{n+1} = t_n + dt$ での $x_{n+1}$ と $v_{n+1}$ は次のように計算されます。 $$ x_{n+1} = x_n + dt \cdot v_n $$ v_{n+1} = v_n + dt \cdot f(t_n, x_n) $$ この計算を、所定の総ステップ数 $nt$ まで繰り返します。
厳密解: プログラムがデフォルトで解いている $d^2x/dt^2 = -x$ という方程式は、単振動の方程式であり、その一般解は $x(t) = A\sin(t) + B\cos(t)$ です。このプログラムでは、初期条件として $x(0)=0.0$ および $v(0)=dx/dt(0)=1.0$ が与えられています。これらの初期条件を適用すると、 $x(0) = A\sin(0) + B\cos(0) = B = 0.0$ $v(t) = dx/dt = A\cos(t) - B\sin(t)$ より $v(0) = A\cos(0) - B\sin(0) = A = 1.0$ したがって、この初期条件に対する厳密解は $x(t) = \sin(t)$ となります。プログラム内の xsolution(t) 関数がこの厳密解を実装しています。

プログラム冒頭の「3重引用符」で囲まれたdocstringには、このプログラムがオイラー法を用いて2階常微分方程式を解くことを示しています。

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

このプログラムは、標準ライブラリ以外に numpy を使用しています。

numpy: 数値計算を効率的に行うためのライブラリですが、このプログラムでは直接的な配列操作などには使われず、将来的な拡張を見越してインポートされている可能性が高いです。現状では必須ではありませんが、インストールしておくのが無難です。

numpy のインストールは、以下の pip コマンドで実行できます。

pip install numpy

必要な入力ファイル

このプログラムは、外部の入力ファイルを必要としません。すべての設定（初期値、時間刻み、ステップ数など）はプログラムのソースコード内に直接記述されており、一部はコマンドライン引数として上書きできます。

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

プログラムは、シミュレーション結果を以下のCSVファイルに出力します。

ファイル名: diffeq2nd_euler.csv
内容: 各行がシミュレーションの特定の時間ステップにおけるデータを表します。
- ヘッダー行: t, x(cal), x(exact), v(cal)
- データ行:
  - t: 現在の時間。
  - x(cal): オイラー法によって計算された $x$ の値。
  - x(exact): 厳密解 $x(t) = \sin(t)$ の値。
  - v(cal): オイラー法によって計算された $v = dx/dt$ の値。

このCSVファイルは、計算結果を他のツール（表計算ソフトなど）で分析したり、グラフ化したりするのに便利です。

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

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

python diffeq2nd_euler.py [dt] [nt] [iprint_interval]

[dt] (オプション): 時間刻み幅 (float型)。デフォルト値は 0.01 です。
[nt] (オプション): シミュレーションの総ステップ数 (int型)。デフォルト値は 501 です。
[iprint_interval] (オプション): 画面に結果を出力する間隔 (int型)。例えば 20 を指定すると、20ステップごとに結果が画面に出力されます。デフォルト値は 20 です。

引数を指定しない場合、プログラムはソースコード内に定義されたデフォルト値を使用します。

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

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

引数なしで実行すると、プログラムはデフォルトの設定 (dt=0.01, nt=501, iprint_interval=20) を使用します。

python diffeq2nd_euler.py

実行結果 (標準出力の一部):

Solve second order diffrential equation by Euler method
Write to [diffeq2nd_euler.csv]
  t      x(cal)      x(exact)      v(cal)
t= 0.00    0.000000    0.000000    1.000000
t= 0.20    0.199000    0.198669    0.980000
t= 0.40    0.388080    0.389418    0.920200
t= 0.60    0.557432    0.564642    0.821192
t= 0.80    0.700346    0.717356    0.684524
...
t= 4.40   -0.892095   -0.970597   -0.450379
t= 4.60   -0.963162   -0.999661   -0.291778
t= 4.80   -0.993086   -0.988032   -0.106674
t= 5.00   -0.977421   -0.958924    0.093496

この実行により、diffeq2nd_euler.csv ファイルが生成されます。

diffeq2nd_euler.csv ファイル内容 (一部):

t,x(cal),x(exact),v(cal)
0.0,0.0,0.0,1.0
0.01,0.01,0.009999833333333333,1.0
0.02,0.02,0.0199973336,0.9999
...
5.0, -0.9774209525492212, -0.9589242746631385, 0.09349603099999997

2. 引数を指定した実行例

時間刻みを大きくし ($dt=0.05$)、総ステップ数を少なく ($nt=100$)、画面出力間隔を短く ($iprint_interval=10$) して実行します。

python diffeq2nd_euler.py 0.05 100 10