diffeq2nd_2d_verlet.py 技術ドキュメント

プログラムの動作

diffeq2nd_2d_verlet.py は、Verlet法を用いて2次元における連立2階常微分方程式を数値的に解くPythonプログラムです。特に、以下のような単振動を表す運動方程式を扱います。

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

このプログラムは、指定された初期条件 (\(x_0, y_0, v_{x0}, v_{y0}\)) と時間刻み幅 (\(\Delta t\)) に基づいて、これらの微分方程式の数値解を計算します。計算結果はコンソールに定期的に表示されるほか、CSVファイルに記録されます。CSVファイルには、各時間ステップにおける計算された座標 (\(x_{cal}, y_{cal}\)) と、対応する厳密解 (\(x_{exact}, y_{exact}\)) が含まれ、数値解の精度を評価できるようになっています。

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

  • 2次元の2階常微分方程式をVerlet法で数値積分します。

  • 最初の1ステップのみ、Verlet法の初期化のためにHeun法を使用します。

  • 厳密解が既知であるため、計算結果と厳密解を比較して出力します。

  • 計算結果はCSVファイルに保存され、後で分析できます。

原理

このプログラムは、主にVerlet法を用いて2階常微分方程式を数値的に解きます。Verlet法は、力学系の時間発展を計算する際に、エネルギー保存則などの物理的な制約を比較的よく満たすという特徴があります。

対象とする方程式

プログラムが解く対象となる方程式は、質量の影響を無視した場合の単純な調和振動子に対応する以下の形です。

\[ \frac{d^2x}{dt^2} = F_x(t, x, y) \]
\[ \frac{d^2y}{dt^2} = F_y(t, x, y) \]

ここで、force(t, x, y) 関数は \(F_x = -x\) および \(F_y = -y\) を返します。

厳密解

初期条件 \(x(0)=0.0, y(0)=2.0, v_x(0)=1.0, v_y(0)=0.0\) に対して、厳密解は以下のようになります。

\[ x(t) = \sin(t) \]
\[ y(t) = 2\cos(t) \]

これは xysolution(t) 関数で実装されています。

Verlet法

Verlet法は、位置 \(\mathbf{r}\) と時間 \(t\) の関数である力 \(\mathbf{F}(\mathbf{r}, t)\) を用いて、次の時間ステップ \(\mathbf{r}(t+\Delta t)\) での位置を計算します。

\[ \mathbf{r}(t+\Delta t) = 2\mathbf{r}(t) - \mathbf{r}(t-\Delta t) + \mathbf{F}(\mathbf{r}(t), t) (\Delta t)^2 \]

プログラムでは、これを2次元座標 \(x, y\) に適用しています。

\[ x(t+\Delta t) = 2x(t) - x(t-\Delta t) + F_x(t, x(t), y(t)) (\Delta t)^2 \]
\[ y(t+\Delta t) = 2y(t) - y(t-\Delta t) + F_y(t, x(t), y(t)) (\Delta t)^2 \]

Verlet法は直接速度を計算しませんが、後処理として現在の速度を以下のように概算できます。

\[ v_x(t) = \frac{x(t+\Delta t) - x(t-\Delta t)}{2\Delta t} \]
\[ v_y(t) = \frac{y(t+\Delta t) - y(t-\Delta t)}{2\Delta t} \]

Heun法による初期ステップの計算

Verlet法は、現在の位置 \(\mathbf{r}(t)\) と一つ前の位置 \(\mathbf{r}(t-\Delta t)\) の2点が必要なため、最初のステップ (\(t=0\) から \(t=\Delta t\)) の計算には別の方法を用いる必要があります。このプログラムでは、Heun法（修正オイラー法）に似た手法を用いて最初の \(\Delta t\) での位置 \(x_1, y_1\) および速度 \(v_{x1}, v_{y1}\) を計算しています。

プログラムの動作は以下の通りです。

  1. 現在の時刻 \(t_0=0\) での力 \(F_x(t_0, x_0, y_0)\) と \(F_y(t_0, x_0, y_0)\) を計算します (\(kvx0, kvy0\) に格納)。

  2. 次の時刻 \(t_1=\Delta t\) での力を予測します。この予測計算は force(0.0+dt, x0+kvy0, y0+kvy0) として実装されており、x0 と y0 には、それぞれ現在の \(F_y(t_0, x_0, y_0)\) の値が加算されたものを次ステップの位置の予測値として利用しています。この計算方法は一般的なHeun法の位置予測とは異なります。 一般的なHeun法での速度 \(v\) の更新は次のようになります。 $\( v_{p}(t+\Delta t) = v(t) + \Delta t \cdot F(t, x(t), y(t)) \)\( \)\( v(t+\Delta t) = v(t) + \frac{\Delta t}{2} [F(t, x(t), y(t)) + F(t+\Delta t, x(t)+\Delta t v(t), y(t)+\Delta t v(t))] \)\( 位置 \)x\( の更新は、同様に平均速度を用いて行われます。 \)\( x(t+\Delta t) = x(t) + \frac{\Delta t}{2} [v(t) + v(t+\Delta t)] \)$

プログラムでは、上記のkvx0, kvy0、kvx1, kvy1を利用して、速度vx1, vy1を計算し、さらにこれらを使って位置x1, y1を計算しています。

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

このプログラムは、以下のライブラリに依存しています。

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

csv および math はPythonの標準ライブラリであり、追加のインストールは不要です。

numpy がインストールされていない場合は、以下のpipコマンドでインストールできます。

pip install numpy

必要な入力ファイル

このプログラムは、コマンドライン引数や外部ファイルからの入力を必要としません。すべての初期条件とパラメーターは、プログラムのソースコード内にハードコードされています。

具体的には、以下のパラメーターが設定されています。

  • outfile: 'diffeq2nd_2D_Verlet.csv'

  • x0, vx0: 初期位置 \(x=0.0\), 初期速度 \(v_x=1.0\)

  • y0, vy0: 初期位置 \(y=2.0\), 初期速度 \(v_y=0.0\)

  • dt: 時間刻み幅 \(0.01\)

  • nt: 総ステップ数 \(501\)

  • iprint_interval: コンソール出力の頻度 \(20\) ステップごと

生成される出力ファイル

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

  • ファイル名: diffeq2nd_2D_Verlet.csv

ファイルの内容: CSVファイルはヘッダー行を持ち、その後に各時間ステップにおけるデータが続きます。

列名

内容

t

現在の時刻

x(cal)

Verlet法によって計算された \(x\) 座標

x(exact)

厳密解による \(x\) 座標

y(cal)

Verlet法によって計算された \(y\) 座標

y(exact)

厳密解による \(y\) 座標

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

このプログラムは引数を取らずに直接実行されます。

python diffeq2nd_2d_verlet.py

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

以下のコマンドでプログラムを実行します。

python diffeq2nd_2d_verlet.py

実行結果 (コンソール出力の例):

Solve simulataneous second order diffrential equations by Verlet method
Write to [diffeq2nd_2D_Verlet.csv]
  t        x(cal)      x(exact)      y(cal)      y(exact)  
t= 0.00    0.000000    0.000000    2.000000    2.000000
t= 0.01    0.010000    0.010000    1.999900    1.999900
t= 0.20    0.198669    0.198669    1.980066    1.980067
t= 0.40    0.389417    0.389418    1.921061    1.921061
t= 0.60    0.564642    0.564642    1.821878    1.821878
t= 0.80    0.717356    0.717356    1.684940    1.684940
t= 1.00    0.841470    0.841471    1.513476    1.513476
... (中略) ...
t= 4.80   -0.996165   -0.996165   -0.199465   -0.199466
t= 5.00   -0.958924   -0.958924   -0.573932   -0.573933

生成される diffeq2nd_2D_Verlet.csv ファイルの内容 (例):

t,x(cal),x(exact),y(cal),y(exact)
0.0,0.0,0.0,2.0,2.0
0.01,0.010000000000000002,0.009999833334166665,1.9999,1.9999000000000001
0.02,0.01999999333333333,0.019999200000000004,1.9996000000000002,1.9996000000000002
... (中略) ...
5.0, -0.9589242746631385, -0.9589242746631385, -0.5739324021759048, -0.5739324021759048