比较分子动力学模拟中Verlet算法、Speed Verlet算法及辛算法的精度

比较分子动力学模拟中Verlet算法、Speed Verlet算法及辛算法的精度

原理

分子动力学模拟是一种研究凝聚态系统的有力方法，对于经典体系，分子动力学假定所有粒子的运动服从牛顿运动方程，只要给出粒子的初始位置、速度（初始条件），边界条件，以及粒子之间相互作用的势函数，就能得到系统在未来任意时刻的状态信息。
由于多原子体系的相互作用非常复杂，难以用解析法求解体系的运动方程，因此用数值方法求解运动方程的近似解。这个过程就不可避免地引入了误差，包括截断误差和舍入误差。前者与不同算法有关，如有限差分法采用TaylorTaylorTaylor展开后截断某些高阶项会引入截断误差，后者来源于计算机本身的精度。两种误差都能通过减小时间步长δt\delta tδt来减小，当δt\delta tδt较大时，截断误差起主要作用，但随着δt\delta tδt的减小会迅速降低，舍入误差随δt\delta tδt变化不大，并在δt\delta tδt较小时起主要作用。
进行数值求解的算法有很多种，对于差分法，其共同点是将粒子的位置，速度，加速度展开为TaylorTaylorTaylor级数。下面我们来分析常用的VerletVerletVerlet算法和更精确一些的SpeedSpeedSpeed VerletVerletVerlet算法的具体步骤以及它们的截断误差。

VerletVerletVerlet AlgorithmAlgorithmAlgorithm

VerletVerletVerlet算法在分子动力学模拟中应用极为广泛，而且原理也很简明。这种方法利用粒子在ttt时刻的位置和加速度，以及t−δtt-\delta tt−δt时刻的位置，计算出t+δtt+\delta tt+δt时刻的位置。
将粒子的位置用TaylorTaylorTaylor公式展开，忽略四阶以上小量：
r(t+δt)=r(t)+ddtr(t)δt+12!d2dt2r(t)(δt)2+13!d3dt3r(t)(δt)3+o[(δt)4]r(t+\delta t)=r(t)+\frac{d}{dt} r(t) \delta t+\frac{1}{2!} \frac{d^2}{dt^2}r(t)(\delta t)^2+\frac{1}{3!} \frac{d^3}{dt^3}r(t)(\delta t)^3+o[(\delta t)^4]r(t+δt)=r(t)+dtdr(t)δt+2!1dt2d2r(t)(δt)2+3!1dt3d3r(t)(δt)3+o[(δt)4]
将式中的δt\delta tδt换为−δt-\delta t−δt得：
r(t−δt)=r(t)−ddtr(t)δt+12!d2dt2r(t)(δt)2−13!d3dt3r(t)(δt)3+o[(δt)4]r(t-\delta t)=r(t)-\frac{d}{dt} r(t) \delta t+\frac{1}{2!} \frac{d^2}{dt^2}r(t)(\delta t)^2-\frac{1}{3!} \frac{d^3}{dt^3}r(t)(\delta t)^3+o[(\delta t)^4]r(t−δt)=r(t)−dtdr(t)δt+2!1dt2d2r(t)(δt)2−3!1dt3d3r(t)(δt)3+o[(δt)4]
以上两式相加，并忽略o[(δt)4]o[(\delta t)^4]o[(δt)4]小量得：
r(t+δt)=2r(t)−r(t−δt)+d2dt2r(t)(δt)2r(t+\delta t)=2r(t)-r(t-\delta t)+\frac{d^2}{dt^2} r(t) (\delta t)^2r(t+δt)=2r(t)−r(t−δt)+dt2d2r(t)(δt)2
两式相减，并忽略o[(δt)2]o[(\delta t)^2]o[(δt)2]小量得：
v(t)=drdt=12δt[r(t+δt)−r(t−δt)]v(t)=\frac{dr}{dt}=\frac{1}{2\delta t}[r(t+\delta t)-r(t-\delta t)]v(t)=dtdr=2δt1[r(t+δt)−r(t−δt)]
至此，我们得到了t+δtt+\delta tt+δt时刻粒子位置和ttt时刻速度的表达式，并且知道了位置的误差正比于o[(δt)4]o[(\delta t)^4]o[(δt)4]，速度的误差正比于o[(δt)2]o[(\delta t)^2]o[(δt)2]。VerletVerletVerlet算法的数值稳定性比简单的欧拉方法高很多，并保持了物理系统中的时间可逆性与相空间体积元体积守恒的性质。然而，这种算法得出的轨迹与速度无关（无法与热浴耦联），并且不能给出同一时刻的位置和速度。

SpeedSpeedSpeed VerletVerletVerlet AlgorithmAlgorithmAlgorithm

在VerletVerletVerlet AlgorithmAlgorithmAlgorithm中，我们发现其无法给出同一时刻的速度和位置，SpeedSpeedSpeed VerletVerletVerlet（有些地方也叫VelocityVelocityVelocity VerletVerletVerlet）正是用来克服这个缺陷的。
将位置展开到二阶小量：
r(t+δt)=r(t)+ddtr(t)δt+12!d2dt2r(t)(δt)2r(t+\delta t)=r(t)+\frac{d}{dt} r(t) \delta t+\frac{1}{2!} \frac{d^2}{dt^2}r(t)(\delta t)^2r(t+δt)=r(t)+dtdr(t)δt+2!1dt2d2r(t)(δt)2
引入t+δt2t+\frac{\delta t}{2}t+2δt时刻的速度：
v(t+δt2)=v(t)+12ddtv(t)δtv(t+\frac{\delta t}{2})=v(t)+\frac{1}{2}\frac{d}{dt} v(t) \delta tv(t+2δt)=v(t)+21dtdv(t)δt
则t+δtt+\delta tt+δt时刻的速度可以表示为
v(t+δt)=v(t+δt2)+12ddtv(t+δt2)δtv(t+\delta t)=v(t+\frac{\delta t}{2})+\frac{1}{2}\frac{d}{dt} v(t+\frac{\delta t}{2}) \delta tv(t+δt)=v(t+2δt)+21dtdv(t+2δt)δt
因此，t+δtt+\delta tt+δt时刻的速度又可表示为：
v(t+δt)=v(t)+12ddtv(t)δt+12ddtv(t+δt2)δtv(t+\delta t)=v(t)+\frac{1}{2}\frac{d}{dt} v(t) \delta t+\frac{1}{2}\frac{d}{dt} v(t+\frac{\delta t}{2}) \delta tv(t+δt)=v(t)+21dtdv(t)δt+21dtdv(t+2δt)δt
其中加速度项通过a(t)=ddtv(t)=−1m∇V(r(t))a(t)=\frac{d}{dt} v(t)=-\frac{1}{m}\nabla V(r(t))a(t)=dtdv(t)=−m1∇V(r(t))得出。因此，SpeedSpeedSpeed VerletVerletVerlet能够给出同一时刻的速度和位置。

SymplecticSymplecticSymplectic AlgorithmAlgorithmAlgorithm

经典动力学方程的求解算法主要包括各级各阶的RK算法，中心差分法，WilsonWilsonWilson θ\thetaθ方法和NewmarkNewmarkNewmark方法等。然而，这些算法本身耗散系统的能量，并使得动态响应的相位滞后，因此其长期跟踪能力有限。
冯康等[1][1][1]建立了哈密顿系统的辛几何算法，从理论上阐明了传统数值方法的耗散来自于其截断项，而辛算法对应的差分格式严格保持哈密顿系统的辛结构，有限阶辛算法的截断部分不会导致系统能量发生线性变化。
辛格式有许多种。Kane等证明了隐式中点格式和NewmarkNewmarkNewmark格式都是变分积分格式，对某些参数的NewmarkNewmarkNewmark格式为辛形式。邢誉峰[2][2][2]等证明了δ=0.5\delta=0.5δ=0.5和α=0.25\alpha=0.25α=0.25的NewmarkNewmarkNewmark算法就是EulerEulerEuler中点隐式差分格式。
由于哈密顿方程可写为q.=Tp\stackrel{.}{q}=Tpq.=Tpp.=−Vq\stackrel{.}{p}=-Vqp.=−Vq对其进行差分可得到1h(qm+3/2−qm+1/2)=TPm+1\frac{1}{h}\Big(q_{m+3/2}-q_{m+1/2}\Big)=TP_{m+1}h1(qm+3/2−qm+1/2)=TPm+11h(pm+1−qm)=−Vqm+1/2\frac{1}{h}\Big(p_{m+1}-q_{m}\Big)=-Vq_{m+1/2}h1(pm+1−qm)=−Vqm+1/2可以证明，上述表示是一种辛格式，即可分线性哈密顿体系的交叉显式辛格式。选用上式进行模拟。

设计思路

一维双谐振子体系

首先通过对一维双谐振子体系的模拟，比较VerletVerletVerlet AlgorithmAlgorithmAlgorithm及SpeedSpeedSpeed VerletVerletVerlet AlgorithmAlgorithmAlgorithm及辛算法的速度及位置精度。
对单个一维线性谐振子，有V(x)=12kx2V(x)=\frac{1}{2}kx^2V(x)=21kx2
md2xdt2=−kxm\frac{d^2x}{dt^2}=-kxmdt2d2x=−kx

对于两个谐振子的情形，如下图所示，给予它们共线的初速度，考虑它们运动的情形

假设两球质量mA=mB=1kgm_A=m_B=1kgmA=mB=1kg，弹簧劲度系数k=0.5N/mk=0.5N/mk=0.5N/m,初始位置xA=−5m,xB=5mx_A=-5m,x_B=5mxA=−5m,xB=5m,初始速度xA=1m/s,xB=−1m/sx_A=1m/s,x_B=-1m/sxA=1m/s,xB=−1m/s，初始加速度aA=aB=0m/s2a_A=a_B=0 m/s^2aA=aB=0m/s2，时间步长δt=0.01\delta t=0.01δt=0.01，步数N=1000−10000N=1000-10000N=1000−10000（谐振子周期在此条件下为T=2πT=2\piT=2π,时间步长应取特征时间的10−310^{-3}10−3~10−210^{-2}10−2,因此取δt=0.01\delta t=0.01δt=0.01是在范围内的）。

运算结果

对于双振子中的一个，比较VerletVerletVerlet AlgorithmAlgorithmAlgorithm，SpeedSpeedSpeed VerletVerletVerlet AlgorithmAlgorithmAlgorithm，辛算法与精确解的轨迹图有（pythonpythonpython实现）：

将它们放在一起比较：

放大局部区域可见

可见VerletVerletVerlet AlgorithmAlgorithmAlgorithm与SpeedSpeedSpeed VerletVerletVerlet AlgorithmAlgorithmAlgorithm以及SymplecticSymplecticSymplectic AlgorithmAlgorithmAlgorithm与精确解的轨迹都重合得很好，这是可以预期的，因为前两种算法位置的误差都是正比于o[(δt)4]o[(\delta t)^4]o[(δt)4]的，而辛算法的耗散在理论上应该非常小。
为了比较三种算法速度的精度，作x−px-px−p相空间轨迹图如下：

放大局部区域可见

可见SpeedSpeedSpeed VerletVerletVerlet AlgorithmAlgorithmAlgorithm以及辛算法与精确解的相空间轨迹图重合得很好。
为了更细致地比较，列出N取不同值时各个算法得到的最后时刻的位置表（δt=0.01\delta t =0.01δt=0.01，保留8位小数）：

N	Verlet Algorithm	Speed Verlet Algorithm	Sympletic Algorithm	Exact Solution
1000	4.17077401	4.17077401	4.17077401	4.17312046
10000	5.33561231	5.33561231	5.33561231	5.29607788

最后时刻的速度表：

N	Verlet Algorithm	Speed Verlet Algorithm	Sympletic Algorithm	Exact Solution
1000	-0.55892883	-0.55892883	-0.55478270	-0.56237908
10000	0.94200169	0.94200169	0.940323631	0.95516380
对比可得，对于位置，三种算法得到的结果相当一致；对于速度，VerletVerletVerlet AlgorithmAlgorithmAlgorithm与SpeedSpeedSpeed VerletVerletVerlet AlgorithmAlgorithmAlgorithm得到的结果相当一致，而辛算法与它们略有差别，差别在0.0010.0010.001量级。
三种算法得到的速度和位置都和精确解略有差别，在δt=0.01,N=1000\delta t =0.01,N=1000δt=0.01,N=1000的条件下，三种算法的速度与精确解相差在0.010.010.01量级，位置相差在0.0010.0010.001量级。

二维三谐振子体系

除一维双谐振子体系之外，还进行了二维三谐振子体系的模拟。但由于这种体系的精确解难以求解，就没有进行几种算法与精确解的对照，只采用VerletVerletVerlet AlgorithmAlgorithmAlgorithm以及SpeedSpeedSpeed VerletVerletVerlet AlgorithmAlgorithmAlgorithm进行了模拟。
假设三球质量mA=mB=mC=1kgm_A=m_B=m_C=1kgmA=mB=mC=1kg，弹簧劲度系数k=1N/mk=1N/mk=1N/m,
初始位置(xA,yA)=(5,0.5)m,(xB,yB)=(−5,0)m,(xC,yC)=(0,1)m(x_A,y_A)=(5,0.5)m,(x_B,y_B)=(-5,0)m,(x_C,y_C)=(0,1)m(xA,yA)=(5,0.5)m,(xB,yB)=(−5,0)m,(xC,yC)=(0,1)m,
初始速度（xA，yA)=(−1,0)m/s,(xB,yB)=(0,0)m/s,(xC,yC)=(0.5,−0.5)m/s（x_A，y_A)=(-1,0)m/s,(x_B,y_B)=(0,0)m/s,(x_C,y_C)=(0.5,-0.5)m/s（xA，yA)=(−1,0)m/s,(xB,yB)=(0,0)m/s,(xC,yC)=(0.5,−0.5)m/s，
初始加速度全为零，时间步长δt=0.01\delta t=0.01δt=0.01，步数N=100N=100N=100。这里取较小的步数是为了得到清晰易于观察的图像。
下面三幅图从上到下依次是VerletVerletVerlet AlgorithmAlgorithmAlgorithm与SpeedSpeedSpeed VerletVerletVerlet AlgorithmAlgorithmAlgorithm模拟三谐振子体系的x−t,y−t,x−yx-t,y-t,x-yx−t,y−t,x−y图。

下图是两种算法模拟三谐振子体系的相空间轨迹图，选取了其中一个谐振子的x−px-px−p轨迹。

可见在上述模拟参数下，两种算法得到的位置和速度没有明显区别。
为了进一步验证VerletVerletVerlet AlgorithmAlgorithmAlgorithm与SpeedSpeedSpeed VerletVerletVerlet AlgorithmAlgorithmAlgorithm得到的速度和位置的一致性，取N=1000N=1000N=1000时最后一帧的速度和位置如下表（δt=0.01,\delta t =0.01,δt=0.01,保留8位小数）：

—	Speed Verlet Algorithm	Verlet Algorithm
DisplacementDisplacementDisplacement	172.25260099	172.25260099
VelocityVelocityVelocity	0.74219359	0.74219359
可见，两种算法得到的位置和速度相当一致。

小结

通过对一维双谐振子体系的模拟，比较了VerletVerletVerlet AlgorithmAlgorithmAlgorithm，SpeedSpeedSpeed VerletVerletVerlet AlgorithmAlgorithmAlgorithm以及辛算法的精度，三者和精确解重合度都很高，而且在速度和位置上的精度在测试的模拟参数下没有明显的区别。辛算法由于具有低耗散的优势，其优势在长时段的模拟中应该能够体现出来。
通过对二维三谐振子体系的模拟，验证了VerletVerletVerlet AlgorithmAlgorithmAlgorithm与SpeedSpeedSpeed VerletVerletVerlet AlgorithmAlgorithmAlgorithm得到的速度及位置是相当一致的。

参考文献

[1]冯康,秦孟兆.哈密顿系统的辛几何算法.[M].杭州：浙江科学技术出版社,2003
[2]邢誉峰,杨蓉.动力学平衡方程的Euler中点辛差分求解格式[J].力学学报,2007(01):100-105.

python 代码

一维双谐振子体系

from pylab import *
import osdt = 0.01
N=1000
x0=5
v0=1
a0=0
m1=1
m2=1
k=0.5
# 用verlet算法模拟谐振子的运动轨迹
def trj_v():t = [0,dt]a1 = [a0]a2 = [a0]vx1 = [-v0]vx2 = [v0]x1=[x0,x0+vx1[0]*dt+ k/m1*a1[0]*1/2*dt**2]x2=[-x0,-x0+vx2[0]*dt+ k/m2*a2[0]*1/2*dt** 2]while(t[-1]<N):a1.append(-k / m1 * (x1[-1] - x2[-1] - (x1[0] - x2[0])))a2.append(-k / m2 * (x2[-1] - x1[-1] - (x2[0] - x1[0])))x1.append(2*x1[-1] -x1[-2] + ( a1[-1]*dt**2))x2.append(2*x2[-1] -x2[-2] + (a2[-1]*dt**2))vx1.append((x1[-1]-x1[-3])/(2*dt))vx2.append((x2[-1]-x2[-3])/(2*dt))t.append(t[-1]+dt)return x1[:-1],x2[:-1],t[:-1],vx1,vx2#删去最后一个元素，便于和speed verlet的结果相比较#输出坐标文件
x1_v=trj_v()[0]
vx1=trj_v()[3]
t_v=trj_v()[2]
coordinate = open("verlet", 'w')
for i in range(len(t_v)):coordinate.write(str(t_v[i]) + "\t" + str(x1_v[i]) + "\t" + str(vx1[i]) + "\n")
coordinate.close()#画verlet算法相应的轨迹图
def plot_v_trj():plot(t_v,x1_v,'yellowgreen',label='verlet algorithm')grid(True)#画verlet算法的相空间轨迹图
def plot_v_xp():x1 = trj_v()[0]p1 = m1*trj_v()[3]plot(x1, p1, 'yellowgreen',label='verlet algorithm',marker="<")grid(True)#####################
#精确解的轨迹
def plot_exact():t=trj_v()[2]x=-sin(array(t))+5 #适用于以上参数情形v = -cos(array(t))plot(t,x,'orange',label='exact solution')return x,t,vdef plot_exact_xp():t = trj_v()[2]x = -sin(array(t))+ 5p = -cos(array(t))plot(x, p, 'orange', label='exact solution',marker="x")t_ex=plot_exact()[1]
x_ex=plot_exact()[0]
v_ex=plot_exact()[2]
#输出坐标文件
coordinate = open('Exact Solution', 'w')
for i in range(len(t_ex)):coordinate.write(str(t_ex[i]) + "\t" + str(x_ex[i]) + "\t" + str(v_ex[i]) + "\n")
coordinate.close()###################### 用 speed verlet算法模拟谐振子的运动轨迹
def trj_sv():t = [0]x1=[x0]x2=[-x0]vx1 = [-v0]vx2=[v0]a1=[a0]a2=[a0]k=0.5while(t[-1]<N):x1.append(x1[-1] + vx1[-1]*dt + a1[-1]*1/2*dt**2)x2.append(x2[-1] + vx2[-1]*dt + a2[-1]*1/2*dt**2)a1.append(-k / m1 * (x1[-1] - x2[-1]-(x1[0]-x2[0])))a2.append(-k / m2 * (x2[-1] - x1[-1]-(x2[0]-x1[0])))# 定义加速度时要小心vx1.append(vx1[-1] + 1 / 2 * dt * (a1[-1] + a1[-2]))vx2.append(vx2[-1] + 1 / 2 * dt * (a2[-1] + a2[-2]))t.append(t[-1]+dt)return x1,x2,t,vx1,vx2#画speed verlet算法相应的轨迹图
def plot_sv_trj():plot(t,x1,'steelblue',label='speed verlet algorithm')grid(True)#画speed_verlet算法的相空间轨迹图
def plot_sv_xp():x1 = trj_sv()[0]p1 = trj_sv()[3]plot(x1, p1, 'r', label='speed verlet algorithm',marker=">")grid(True)#输出坐标文件
x1=trj_sv()[0]
x2=trj_sv()[1]
t=trj_sv()[2]
vx1=trj_sv()[3]
coordinate = open("speed_verlet", 'w')
for i in range(len(t)):coordinate.write(str(t[i]) + "\t" + str(x1[i]) + "\t" + str(vx1[i]) + "\n")
coordinate.close()###################辛算法,可分线性哈密顿体系的交叉显式辛格式
def trj_e2():t = [0]vx1 = [-v0]vx2 = [v0]x1 = [x0]x2 = [-x0]#w = 2 ** (1 / 2)k=1while (t[-1] < N):x1.append(x1[-1] +dt*vx1[-1])x2.append(x2[-1] +dt*vx2[-1])# a1.append(-k / m1 *(x1[-1] - x2[-1] - (x1[0] - x2[0])))# a2.append(-k / m2 * (x2[-1] - x1[-1] - (x2[0] - x1[0])))vx1.append(vx1[-1]-dt*k*(x1[-1]-x1[0])/m1)vx2.append(vx1[-1]-dt*k*(x2[-1]-x2[0])/m2)t.append(t[-1] + dt)#print(x1)return x1, x2, t, vx1, vx2def plot_e2_trj():plot(t_e2,x1_e2,'r',label='Sympletic Algorithm')grid(True)def plot_e2_xp():x1 = trj_e2()[0]p1 = m1*trj_e2()[3]plot(x1, p1, 'g', label='Sympletic Algorithm',marker="*")grid(True)#输出坐标文件
x1_e2=trj_e2()[0]
x2_e2=trj_e2()[1]
t_e2=trj_e2()[2]
vx1_e2=trj_e2()[3]
coordinate = open('Sympletic Algorithm', 'w')
for i in range(len(t_e2)):coordinate.write(str(t_e2[i]) + "\t" + str(x1_e2[i]) + "\t" + str(vx1_e2[i]) + "\n")
coordinate.close()################################ 使得图像可以随鼠标滚轮而缩放
#思路：
#使用fig.canvas.mpl_connect()函数来绑定相关fig的滚轮事件
#利用事件event的inaxes属性获取当前鼠标所在坐标系ax
#使用get_xlim()函数获取坐标系ax的x/y轴坐标刻度范围
#使用set()函数对坐标系ax进行放大/缩小import numpy as np
fig = figure()
def call_back(event):axtemp = event.inaxesx_min, x_max = axtemp.get_xlim()increment = (x_max - x_min) / 10if event.button == 'up':axtemp.set(xlim=(x_min + increment, x_max - increment))print('up')elif event.button == 'down':axtemp.set(xlim=(x_min - increment, x_max + increment))print('down')fig.canvas.draw_idle()  # 绘图动作实时反映在图像上fig.canvas.mpl_connect('scroll_event', call_back)
fig.canvas.mpl_connect('button_press_event', call_back)# plot_v_trj()
# plot_sv_trj()
# plot_e2_trj()
# plot_exact()
# legend(loc="right")
# xlabel('t',fontsize="16")
# ylabel('x',fontsize="16")
# title('Twins Harmonic Oscillators',fontsize="16")
# show()subplot(4,1,1)
plot_v_trj()
title('Twins Harmonic Oscillators',fontsize="16")
legend(loc="right")
ylabel('x',fontsize="16")subplot(4,1,2)
plot_sv_trj()
ylabel('x',fontsize="16")
legend(loc="right")subplot(4,1,3)
plot_e2_trj()
legend(loc="right")
ylabel('x',fontsize="16")subplot(4,1,4)
plot_exact()
legend(loc="right")
xlabel('t',fontsize="16")
ylabel('x',fontsize="16")show()fig = figure()
def call_back(event):axtemp = event.inaxesx_min, x_max = axtemp.get_xlim()increment = (x_max - x_min) / 10if event.button == 'up':axtemp.set(xlim=(x_min + increment, x_max - increment))print('up')elif event.button == 'down':axtemp.set(xlim=(x_min - increment, x_max + increment))print('down')fig.canvas.draw_idle()  # 绘图动作实时反映在图像上fig.canvas.mpl_connect('scroll_event', call_back)
fig.canvas.mpl_connect('button_press_event', call_back)
plot_v_xp()
plot_sv_xp()
plot_e2_xp()
plot_exact_xp()
legend(loc="right")
xlabel('x',fontsize="16")
ylabel('p',fontsize="16")
title('Twins Harmonic Oscillators',fontsize="16")
show()

二维三谐振子体系

from matplotlib.pyplot import *
from vpython import *class Harmonics:def __init__(self,_m1=1,_x1=5,_y1=0.5,_vx1=-1,_vy1=0,_ax1=0,_ay1=0,_m2=1,_x2=-5,_y2=0,_vx2=1,_vy2=0,_ax2=0,_ay2=0,_m3=1,_x3=0,_y3=1,_vx3=0.5,_vy3=-0.5,_ax3=0,_ay3=0,_t=0,_dt=0.01,_k=1,_n=1000):self.x1=[]self.x1.append(_x1)self.vx1=[]self.vx1.append(_vx1)self.ax1=[]self.ax1.append(_ax1)self.y1=[]self.y1.append(_y1)self.vy1=[]self.vy1.append(_vy1)self.ay1 = []self.ay1.append(_ay1)self.x2=[]self.x2.append(_x2)self.vx2=[]self.vx2.append(_vx2)self.ax2 = []self.ax2.append(_ax2)self.y2=[]self.y2.append(_y2)self.vy2=[]self.vy2.append(_vy2)self.ay2 = []self.ay2.append(_ay2)self.x3=[]self.x3.append(_x3)self.vx3=[]self.vx3.append(_vx3)self.ax3 = []self.ax3.append(_ax3)self.y3=[]self.y3.append(_y3)self.vy3=[]self.vy3.append(_vy3)self.ay3 = []self.ay3.append(_ay3)self.m1=_m1self.m2=_m2self.m3=_m3self.dt=_dtself.n=_nself.t=[]self.t.append(_t)self.k=_kreturndef calculate_speedverlet(self):while self.t[-1]<self.n:self.x1.append(self.x1[-1] + self.vx1[-1] * self.dt + self.ax1[-1] * 1 / 2 * self.dt ** 2)self.x2.append(self.x2[-1] + self.vx2[-1] * self.dt + self.ax2[-1] * 1 / 2 * self.dt ** 2)self.x3.append(self.x3[-1] + self.vx3[-1] * self.dt + self.ax3[-1] * 1 / 2 * self.dt ** 2)self.y1.append(self.y1[-1] + self.vy1[-1] * self.dt + self.ay1[-1] * 1 / 2 * self.dt ** 2)self.y2.append(self.y2[-1] + self.vy2[-1] * self.dt + self.ay2[-1] * 1 / 2 * self.dt ** 2)self.y3.append(self.y3[-1] + self.vy3[-1] * self.dt + self.ay3[-1] * 1 / 2 * self.dt ** 2)self.ax1.append(-self.k / self.m1 * (self.x1[-1] - self.x2[-1] - (self.x1[0] - self.x2[0]))-self.k / self.m1 * (self.x1[-1] - self.x3[-1] - (self.x1[0] - self.x3[0])))self.ax2.append(-self.k / self.m2 * (self.x2[-1] - self.x1[-1] - (self.x2[0] - self.x1[0]))-self.k / self.m2 * (self.x2[-1] - self.x3[-1] - (self.x2[0] - self.x3[0])))self.ax3.append(-self.k / self.m3 * (self.x3[-1] - self.x1[-1] - (self.x3[0] - self.x1[0]))-self.k / self.m3 * (self.x3[-1] - self.x2[-1] - (self.x3[0] - self.x2[0])))self.ay1.append(-self.k / self.m1 * (self.y1[-1] - self.y2[-1] - (self.y1[0] - self.y2[0]))- self.k / self.m1 * (self.y1[-1] - self.y3[-1] - (self.y1[0] - self.y3[0])))self.ay2.append(-self.k / self.m2 * (self.y2[-1] - self.y1[-1] - (self.y2[0] - self.y1[0]))- self.k / self.m2 * (self.y2[-1] - self.y3[-1] - (self.y2[0] - self.y3[0])))self.ay3.append(-self.k / self.m3 * (self.y3[-1] - self.y1[-1] - (self.y3[0] - self.y1[0]))- self.k / self.m3 * (self.y3[-1] - self.y2[-1] - (self.y3[0] - self.y2[0])))self.vx1.append(self.vx1[-1] + 1 / 2 * self.dt * (self.ax1[-1] + self.ax1[-2]))self.vx2.append(self.vx2[-1] + 1 / 2 * self.dt * (self.ax2[-1] + self.ax2[-2]))self.vx3.append(self.vx3[-1] + 1 / 2 * self.dt * (self.ax3[-1] + self.ax3[-2]))self.vy1.append(self.vy1[-1] + 1 / 2 * self.dt * (self.ay1[-1] + self.ay1[-2]))self.vy2.append(self.vy2[-1] + 1 / 2 * self.dt * (self.ay2[-1] + self.ay2[-2]))self.vy3.append(self.vy3[-1] + 1 / 2 * self.dt * (self.ay3[-1] + self.ay3[-2]))self.t.append(self.t[-1] + self.dt)return self.t , self.x1,self.vx1def plot_tx(self):plot(self.t,self.x1,'r',lw=1,label='Speed Verlet Algorithm')plot(self.t,self.x2,'yellowgreen',lw=1.5,label='Speed Verlet Algorithm')plot(self.t,self.x3,'orange',lw=1,label='Speed Verlet Algorithm')xlabel('t')ylabel('x')grid(True)legend(loc='right')returndef plot_ty(self):plot(self.t,self.y1,'r',lw=1,label='Speed Verlet Algorithm')plot(self.t,self.y2,'yellowgreen',lw=1.5,label='Speed Verlet Algorithm')plot(self.t,self.y3,'orange',lw=1,label='Speed Verlet Algorithm')xlabel('t')ylabel('y')grid(True)legend(loc='right')returndef plot_2d(self):plot(self.x1,self.y1,'r',lw=1,label='Speed Verlet Algorithm')plot(self.x2,self.y2,'yellowgreen',lw=1.5,label='Speed Verlet Algorithm')plot(self.x3,self.y3,'orange',lw=1,label='Speed Verlet Algorithm')xlabel('x')ylabel('y')grid(True)legend(loc='right')returndef plot_xp(self):plot(self.x1, self.m1*self.vx1, 'steelblue', label='speed verlet algorithm',lw="2")grid(True)class Harmonics_2(Harmonics):def calculate_verlet(self):self.t.append(self.dt)self.x1.append(5 + self.vx1[0]*self.dt+ self.k/self.m1*self.ax1[0]*1/2*self.dt**2)self.y1.append(0.5 + self.vy1[0] * self.dt + self.k / self.m1 * self.ay1[0] * 1 / 2 * self.dt ** 2)self.x2.append(-5 + self.vx2[0] * self.dt + self.k / self.m2 * self.ax2[0] * 1 / 2 * self.dt ** 2)self.y2.append(0 + self.vy2[0] * self.dt + self.k / self.m2 * self.ay2[0] * 1 / 2 * self.dt ** 2)self.x3.append(0 + self.vx3[0] * self.dt + self.k / self.m3 * self.ax3[0] * 1 / 2 * self.dt ** 2)self.y3.append(1 + self.vy3[0] * self.dt + self.k / self.m3 * self.ay3[0] * 1 / 2 * self.dt ** 2)while self.t[-1]<self.n:self.ax1.append(-self.k / self.m1 * (self.x1[-1] - self.x2[-1] - (self.x1[0] - self.x2[0]))- self.k / self.m1 * (self.x1[-1] - self.x3[-1] - (self.x1[0] - self.x3[0])))self.ax2.append(-self.k / self.m2 * (self.x2[-1] - self.x1[-1] - (self.x2[0] - self.x1[0]))- self.k / self.m2 * (self.x2[-1] - self.x3[-1] - (self.x2[0] - self.x3[0])))self.ax3.append(-self.k / self.m3 * (self.x3[-1] - self.x1[-1] - (self.x3[0] - self.x1[0]))- self.k / self.m3 * (self.x3[-1] - self.x2[-1] - (self.x3[0] - self.x2[0])))self.ay1.append(-self.k / self.m1 * (self.y1[-1] - self.y2[-1] - (self.y1[0] - self.y2[0]))- self.k / self.m1 * (self.y1[-1] - self.y3[-1] - (self.y1[0] - self.y3[0])))self.ay2.append(-self.k / self.m2 * (self.y2[-1] - self.y1[-1] - (self.y2[0] - self.y1[0]))- self.k / self.m2 * (self.y2[-1] - self.y3[-1] - (self.y2[0] - self.y3[0])))self.ay3.append(-self.k / self.m3 * (self.y3[-1] - self.y1[-1] - (self.y3[0] - self.y1[0]))- self.k / self.m3 * (self.y3[-1] - self.y2[-1] - (self.y3[0] - self.y2[0])))self.x1.append(2 * self.x1[-1] - self.x1[-2] + (self.ax1[-1] * self.dt ** 2))self.x2.append(2 * self.x2[-1] - self.x2[-2] + (self.ax2[-1] * self.dt ** 2))self.x3.append(2 * self.x3[-1] - self.x3[-2] + (self.ax3[-1] * self.dt ** 2))self.y1.append(2 * self.y1[-1] - self.y1[-2] + (self.ay1[-1] * self.dt ** 2))self.y2.append(2 * self.y2[-1] - self.y2[-2] + (self.ay2[-1] * self.dt ** 2))self.y3.append(2 * self.y3[-1] - self.y3[-2] + (self.ay3[-1] * self.dt ** 2))self.vx1.append((self.x1[-1] - self.x1[-3]) / (2 * self.dt))self.vx2.append((self.x2[-1] - self.x2[-3]) / (2 * self.dt))self.vx3.append((self.x3[-1] - self.x3[-3]) / (2 * self.dt))self.vy1.append((self.y1[-1] - self.y1[-3]) / (2 * self.dt))self.vy2.append((self.y2[-1] - self.y2[-3]) / (2 * self.dt))self.vy3.append((self.y3[-1] - self.y3[-3]) / (2 * self.dt))self.t.append(self.t[-1] + self.dt)return self.t ,self.x1,self.vx1def plot_tx(self):plot(self.t,self.x1,'r',label='Verlet Algorithm',linestyle=':',lw=2.5)plot(self.t,self.x2,'yellowgreen',label='Verlet Algorithm',linestyle=':',lw=2.5)plot(self.t,self.x3,'orange',label='Verlet Algorithm',linestyle=':',lw=2.5)xlabel('t')ylabel('x')grid(True)legend(loc='right')returndef plot_ty(self):plot(self.t,self.y1,'r',label='Verlet Algorithm',ls=':',lw=2.5)plot(self.t,self.y2,'yellowgreen',label='Verlet Algorithm',linestyle=':',lw=2.5)plot(self.t,self.y3,'orange',label='Verlet Algorithm',linestyle=':',lw=2.5)xlabel('t')ylabel('y')grid(True)legend(loc='right')returndef plot_2d(self):plot(self.x1,self.y1,'r',label='Verlet Algorithm',ls=':',lw=2.5)plot(self.x2,self.y2,'yellowgreen',label='Verlet Algorithm',ls=':',lw=2.5)plot(self.x3,self.y3,'orange',label='Verlet Algorithm',ls=':',lw=2.5)xlabel('x')ylabel('y')grid(True)legend(loc='right')returndef plot_xp(self):x1_cut = self.x1[:-1]plot(x1_cut, self.m1*self.vx1, 'r', label='verlet algorithm',ls=':',lw=2.5)grid(True)a=Harmonics()
a.calculate_speedverlet()fig = figure()
def call_back(event):axtemp = event.inaxesx_min, x_max = axtemp.get_xlim()increment = (x_max - x_min) / 10if event.button == 'up':axtemp.set(xlim=(x_min + increment, x_max - increment))print('up')elif event.button == 'down':axtemp.set(xlim=(x_min - increment, x_max + increment))print('down')fig.canvas.draw_idle()  # 绘图动作实时反映在图像上fig.canvas.mpl_connect('scroll_event', call_back)
fig.canvas.mpl_connect('button_press_event', call_back)subplot(3,1,1)
a.plot_tx()
title("Three Harmonic Oscillators")
#title("speed verlet algorithm")
subplot(3,1,2)
a.plot_ty()
subplot(3,1,3)
a.plot_2d()##############################
b=Harmonics_2()
b.calculate_verlet()subplot(3,1,1)
b.plot_tx()
#title("verlet algorithm")
subplot(3,1,2)
b.plot_ty()
subplot(3,1,3)
b.plot_2d()
show()###############################fig = figure()
def call_back(event):axtemp = event.inaxesx_min, x_max = axtemp.get_xlim()increment = (x_max - x_min) / 10if event.button == 'up':axtemp.set(xlim=(x_min + increment, x_max - increment))print('up')elif event.button == 'down':axtemp.set(xlim=(x_min - increment, x_max + increment))print('down')fig.canvas.draw_idle()  # 绘图动作实时反映在图像上fig.canvas.mpl_connect('scroll_event', call_back)
fig.canvas.mpl_connect('button_press_event', call_back)a.plot_xp()
b.plot_xp()
xlabel('x',fontsize="15")
ylabel('p',fontsize="15")
title("Three Harmonic Oscillators")
legend()
show()#输出坐标文件
t=a.calculate_speedverlet()[0]
x1=a.calculate_speedverlet()[1]
vx1=a.calculate_speedverlet()[2]
coordinate = open("speed_verlet_three", 'w')
for i in range(len(t)):coordinate.write(str(t[i]) + "\t" + str(x1[i]) + "\t" + str(vx1[i]) + "\n")
coordinate.close()#输出坐标文件
t_v=b.calculate_verlet()[0]
x1_v=b.calculate_verlet()[1]
vx1_v=b.calculate_verlet()[2]
coordinate = open("verlet_three", 'w')
for i in range(len(t)):coordinate.write(str(t_v[i]) + "\t" + str(x1_v[i]) +  "\n")
coordinate.close()coordinate = open("verlet_three_v", 'w')
for i in range(len(vx1)):coordinate.write( str(vx1_v[i])+ "\n")
coordinate.close()