Velocity-Verlet 双井算法 Python

Question

我正在为双阱势 V(x) = x^4-20x^2 实施 Verlet 算法，以创建一个简单的相 portrait.The 生成的相图具有增强的椭圆形状，这显然是不正确的。我感觉我的问题出现在我对 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.

Velocity-Verlet 双井算法 Python

Velocity-Verlet Double-Well Algorithm Python

python

numpy

physics

verlet-integration