VelocityVerlet双井算法

Question

我正在为双阱势V(x) = x^4-20x^2实施Verlet算法，以创建简单的相位图。生成的相位图具有增强的椭圆形，显然是不正确的。我有一种感觉，我的问题发生在我对x^3的定义中，但我不确定。我还包括了一个经典谐振子的算法，以表明我的代码工作正常

import numpy as np
import matplotlib.pyplot as plt
###Constants
w = 2
m=1
N=500
dt=0.05
t = np.linspace(0, N*dt, N+1)
np.shape(t)
x = np.zeros(N+1)
p = np.zeros(N+1)
p_0 = 0
x_0 = 1
x[0] = x_0
p[0] = p_0

#Velocity Verlet Tuckerman
    #x(dt) = x(0) +p(0)/m*dt + 1/(2m) * F(x(0))
    #p(dt =  p(0) + dt/2[F(x(0)) + F(x(dt))]

#Harmonic Oscillator F(x) = -kx = -mw^2x
for n in range(N):
    x[n+1] = x[n] + (p[n]/m)*dt - (0.5)*w**2*x[n]*dt*dt

    p[n+1] = p[n] - m*(0.5)*w**2*x[n]*dt - m*0.5*w**2*x[n+1]*dt

plt.plot(x,p)

#Symmetric Double Well: F(x) = -4x^3 + 40x
                      #V(x) = x^4 -20x^2
for n in range(N):
    x[n+1] = x[n] + (p[n]/m)*dt +1/(2*m)*( -4*(x[n]*x[n]*x[n])*dt*dt +40*x[n]*dt*dt)

    p[n+1] = p[n]  + (1/2)*(-4*m*(x[n]*x[n]*x[n])*dt +40*m*x[n]*dt - 4*m*(x[n+1]*x[n+1]*x[n+1])*dt +40*m*x[n+1]*dt)


plt.plot(x,p)

谢谢

Answer 1

更准确地说，V在x=+-sqrt(10)处有一个最小值，值为-100，在x=0处的局部最大值给出值0。初始位置x0=1, v0=0将溶液放置在右谷，围绕sqrt(10)振荡

要得到一个8字形，需要一个V(x0)稍大于零的初始点。例如，使用x0=5可以得到V=25*(25-20)=125。或者以x0=4.5 ==> x0^2=20.25 ==> V ~ 5为例