速度ベルレー(Verlet)法により一次元調和振動子の運動方程式，$ M d^2x/dt^2 = F(x,t)=- k x $を初期条件 $x = 0$, $dx/dt =1.0$ の下で解く。質点の質量M = 1 kg，バネ定数 $k = 1.0$ N/mとする。

速度ベルレ法では座標xと速度v($=dx/dt$)を以下のように更新する[3]。

$x_{n+1} = x_n +v_n \ \delta t +0.5 \ F(x_n,t_n)\ (\delta t)^2/M $
$v_{n+1} = v_n+ 0.5 \ \delta t\ (\ F(x_n, t_n) + F(x_{n+1}, t_{n+1}))/M$

ルンゲ-クッタ法よりも全エネルギー保存性が高くなる点が特徴である(シンプレクティック性を持つ)[1,2]。また，速度ベルレ法は通常のベルレ法と数学的に等価であるが丸め誤差に強く数値的により堅牢な表式となっている。これらの理由から，分子動力学シミュレーションでしばしば速度ベルレ法が用いられる。

import numpy as np
import matplotlib.pyplot as plt

"""
速度ベルレー法による解法
"""



def f(x, t): # 
    return -k*x

M = 1.0 # 質量 mass [Kg]
k=1.0  # ばね定数 spring constant [N/m]

t0 = 0.0 # シミュレーション時間の初期値
t1 = 20.0 # シミュレーション時間の最大値
N = 200 # t1からt0までの時間刻み数: 
del_t = (t1-t0)/N # 時間刻み time step

tpoints = np.arange(t0, t1, del_t) # t0からt1までの時間点をdel_t刻みで生成
xpoints = [] # 時間にx座標を格納するためのリスト
vpoints = []# 時間毎に速度v座標を格納するためのリスト
Etot=[] # 時間毎に力学的エネルギーを格納するためのリスト

# initial condition
x0 = 0.0 # xの初期条件
v0 = 1.0 #速度vの初期条件

x, v = x0, v0


for t in tpoints:
    xpoints.append(x)
    vpoints.append(v)
    Etot.append((M*v**2)/2+(k*x**2)/2) # 力学的エネルギー (mv^2/2 + k x^2/2)

    xx=x  # xの前ステップの格納
    x +=  v*del_t + 0.5 * f(x,t)*(del_t**2) #速度ベルレー法による座標更新
    v += 0.5*del_t*(f(xx , t) + f(x, t+del_t)) # 速度ベルレー法による速度更新
     
# for plot  #厳密解 x = sin(t)との比較
plt.plot (tpoints, xpoints, 'o',label='Velocity Verlet')
plt.xlabel("t",  fontsize=24)
plt.ylabel("x(t)",  fontsize=24)

tt=np.arange(t0, t1, 0.01*del_t)
plt.plot (tt, np.sin(tt), '-',label='Exact: Sin (t)',color = 'Green')
plt.legend(loc='upper left')

plt.show()

##厳密解 v = cos(t)との比較
plt.plot (tpoints,vpoints, 'o',label='Velocity Verlet')
plt.xlabel("t", fontsize=24)
plt.ylabel("v(t)", fontsize=24)

tt=np.arange(t0, t1, 0.01*del_t)
plt.plot (tt, np.cos(tt), '-',label='Exact: Cos (t)', color = 'red')
plt.legend(loc='upper left')

plt.show()


# 全エネルギー Etotの描画。厳密解はEtot=0.5
plt.plot (tpoints, Etot, '-',label='Velocity Verlet')
plt.xlabel("t",  fontsize=24)
plt.ylabel("Etot(t)",  fontsize=24)


plt.legend(loc='upper left')

plt.show()

##補遺: シンプレクティック法の長所と短所 [4]
長所 : 保存系の力学の高精度計算が可能， (ニュートン方程式が持つ)時間反転対称性を持つ，長期的な力学的エネルギー保存性を有する，など。

短所: 時間刻みの変更には不便 (自動精度制御不可能)，非保存系では特色が生かされない（定温・定圧手法、速度に依存する力がでてくるとき，など。

##参考文献

[1] 遠藤理平, 「常微分方程式の解法」

[2]陰山聡, 「シンプレクティック積分法について」

[3] webではココがわかりやすい。本では，たとえば，ハーベイ著計算物理学入門やクーニン著計算機物理学などに初等的で分かりやすい説明がある。

[4] 能勢修一, [「分子動力学シミュレーション・数値積分法」] (http://www.phys.keio.ac.jp/guidance/labs/riron/pdf/nosenote1.pdf)