1

我正在实现双阱电位的 Verlet 算法V(x) = x^4-20x^2，以创建简单的相图。生成的相图具有增强的椭圆形，显然是不正确的。我有一种感觉，我的问题出现在我对 of 的定义中，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)

谢谢！

4

1 回答 1

0

更准确地说，在 valueV处有最小值，在x=+-sqrt(10)value处有-100局部最大值。初始位置将解决方案置于右谷，围绕 振荡。x=00x0=1, v0=0sqrt(10)

要获得 8 字形，您需要一个V(x0)略大于零的初始点。例如用x0=5一个得到V=25*(25-20)=125。或采取x0=4.5 ==> x0^2=20.25 ==> V ~ 5。

于 2021-01-25T10:18:32.907 回答