I am implemententing the Verlet algorithm for a double well potential V(x) = x^4-20x^2, as to create a simple phase portrait.The generated phase portrait has an augmented oval shape and is clearly incorrect. I have a feeling that my problem is occurring in my definition of the of x^3 but I am not sure. I have also included the algorithm for a classical harmonic oscillator to show that my code works correctly.
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)
Thanks!
To be more precise, V has a minimum at x=+-sqrt(10) with value -100, the local maximum at x=0 gives value 0. The initial position x0=1, v0=0 places the solution in the right valley, oscillating around sqrt(10).
To get a figure-8 shape you need an initial point with V(x0) slightly larger than zero. For instance with x0=5 one gets V=25*(25-20)=125. Or take x0=4.5 ==> x0^2=20.25 ==> V ~ 5.