Diffeq2nd Verlet Method ドキュメント
プログラムの動作
diffeq2nd_verlet.py は、2階常微分方程式をVerlet法を用いて数値的に解くPythonプログラムです。特に、単振動の方程式 \(d^2x/dt^2 = -x\) を初期条件 \(x(0)=0.0\), \(v(0)=1.0\) の下で解き、その結果を解析解 \(x(t)=\sin(t)\) と比較します。
主な機能は以下の通りです。
指定された時間ステップ
dtと総ステップ数ntに基づいて、時刻tにおける位置xと速度vを計算します。計算された数値解を解析解と比較し、その誤差を評価します。
計算結果を標準出力に表示するとともに、CSVファイルに保存します。
このプログラムは、物理シミュレーションにおける基本的な数値積分法であるVerlet法の動作と精度を、解析解と比較することで理解することを目的としています。
原理
このプログラムは、以下の2階常微分方程式を対象とします。
本プログラムでは、具体的な関数 \(f(t, x)\) として \(f(t, x) = -x\) を採用しています。これは単振動を表す方程式です。
初期条件は \(x(0)=0.0\) および \(v(0)=1.0\) です。この初期条件における解析解は以下の通りです。
この方程式を解くために、Verlet法（特にVelocity Verletではない、古典的なVerlet法）が用いられます。Verlet法は、位置の更新を以下のように行います。
ここで、\(\Delta t\) は時間ステップです。速度 \(v(t)\) は、中心差分近似を用いて以下のように導出されます。
Verlet法を適用するには、\(x(t)\) と \(x(t-\Delta t)\) の2つの時刻における位置が必要になります。しかし、初期状態では \(x(0)\) と \(v(0)\) しか与えられていません。そのため、プログラムでは最初の1ステップ（\(x(0)\) から \(x(\Delta t)\) へ）を計算するために、以下のような2次近似の手法を採用しています。
時刻 \(t=0\) での速度 \(v(0)\) と位置 \(x(0)\) を用いて、時刻 \(\Delta t\) での速度 \(v(\Delta t)\) を計算します。これはHeun法（改良Euler法）に類似しています。
\(\Delta v_0 = \Delta t \cdot f(0, x(0))\)
\(\Delta v_1 = \Delta t \cdot f(\Delta t, x(0) + \Delta v_0)\) (速度予測子の位置を使用)
\(v(\Delta t) = v(0) + (\Delta v_0 + \Delta v_1) / 2.0\)
次に、計算された \(v(0)\) と \(v(\Delta t)\) を用いて、時刻 \(\Delta t\) での位置 \(x(\Delta t)\) を計算します。これは台形公式に類似しています。
\(\Delta x_0 = \Delta t \cdot v(0)\)
\(\Delta x_1 = \Delta t \cdot v(\Delta t)\) (上記で計算された \(v(\Delta t)\) を使用)
\(x(\Delta t) = x(0) + (\Delta x_0 + \Delta x_1) / 2.0\)
これらの計算により、\(x(0)\) と \(x(\Delta t)\) が得られ、その後のステップで通常のVerlet法の位置更新が適用可能となります。
必要な非標準ライブラリとインストール方法
このプログラムは、以下のPythonライブラリをインポートしています。
csv: Python標準ライブラリであり、追加インストールは不要です。numpy: 高度な数値計算を行うためのライブラリです。このプログラムでは直接的にNumpyの関数は使用されていませんが、将来的な拡張のためインポートされています。math: Python標準ライブラリであり、追加インストールは不要です。
numpy は非標準ライブラリのため、インストールが必要です。以下の pip コマンドでインストールできます。
pip install numpy
必要な入力ファイル
このプログラムは、実行に必要な外部入力ファイルを指定しません。すべてのパラメータ（初期位置、初期速度、時間ステップ、総ステップ数、出力ファイル名など）はプログラムのソースコード内にハードコードされています。
生成される出力ファイル
プログラムを実行すると、以下のCSVファイルが生成されます。
diffeq2nd_verlet.csv:内容: 各時刻における計算された位置、解析解、計算された速度がCSV形式で保存されます。
ファイル形式: カンマ区切り（CSV）。各行の終端は改行コード（
\n）です。ヘッダ行: 以下の4つの列名を含みます。
t: 時刻x(cal): Verlet法で計算された位置x(exact): 解析解 \(x(t)=\sin(t)\) による位置v(cal): Verlet法で計算された速度
ファイルの内容の例（最初の数行）：
t,x(cal),x(exact),v(cal)
0.0,0.000000,0.000000,1.000000
0.01,0.010000,0.010000,0.999950
0.02,0.020000,0.020000,0.999800
...
コマンドラインでの使用例 (Usage)
このプログラムは、コマンドライン引数を受け付けません。すべての設定はソースコード内に記述されています。
基本的な実行コマンドは以下の通りです。
python diffeq2nd_verlet.py
コマンドラインでの具体的な使用例
以下のコマンドを実行することで、プログラムを起動できます。
python diffeq2nd_verlet.py
実行結果の説明:
プログラムが実行されると、標準出力には以下のような計算の進捗と結果の一部が表示されます。
Solve second order diffrential equation by Verlet method
Write to [diffeq2nd_verlet.csv]
t x(cal) x(exact) v(cal)
t= 0.00 0.000000 0.000000 1.000000
t= 0.01 0.010000 0.010000 0.999950
t= 0.20 0.198669 0.198669 0.980067
t= 0.40 0.389418 0.389418 0.921057
...
t= 5.00 -0.958925 -0.958924 0.283662
同時に、現在のディレクトリに diffeq2nd_verlet.csv というファイルが生成されます。このファイルには、すべての時刻ステップにおける計算結果と解析解がCSV形式で詳細に記録されています。例えば、このファイルの最初の数行と最後の数行は以下のようになるでしょう。
t,x(cal),x(exact),v(cal)
0.0,0.000000,0.000000,1.000000
0.01,0.010000,0.010000,0.999950
0.02,0.020000,0.020000,0.999800
0.03,0.030000,0.030000,0.999550
...
5.00, -0.958925, -0.958924, 0.283662
5.01, -0.961759, -0.961757, 0.277853