How to implement a velocity Verlet integrator which works for the harmonic oscillator in python?

Question

I am new to python and i am trying to implement a velocity Verlet integrator which works for the harmonic oscillator. As you can see from my notebook below (taken from: http://hplgit.github.io/prog4comp/doc/pub/._p4c-solarized-Python022.html), the Euler's method works, but not the verlet integrator. What am i missing?

## Euler's method

from numpy import zeros, linspace, pi, cos, array
import matplotlib.pyplot as plt
import numpy as np

omega = 1
m=1
N=500
dt=0.08
t = linspace(0, N*dt, N+1)

u = zeros(N+1)
vi = zeros(N+1)

X_0 = 2
u[0] = X_0
vi[0] = 0

for n in range(N):
    u[n+1] = u[n] + dt*vi[n]
    vi[n+1] = vi[n] - dt*omega**2*u[n]

fig = plt.figure()
l1, l2 = plt.plot(t, u, 'b-', t, X_0*cos(omega*t), 'r--')
fig.legend((l1, l2), ('numerical', 'exact'), 'upper left')
plt.xlabel('t')
plt.show()

## Velocity Verlet

for n in range(N):
    v_next = vi[n] - 0.5*omega**2*u[n]*dt
    u[n+1]= u[n]+ dt*v_next
    vi[n+1] = v_next - 0.5*omega**2*u[n]*dt
    np.append(u, x_next)

plt.plot(t, X_0*cos(omega*t), label = 'Analytical')
plt.plot(t, u, 'r--', label = 'Verlet')
plt.title('Verlet', fontweight = 'bold', fontsize = 1)
plt.xlabel('t', fontweight = 'bold', fontsize = 14)
plt.ylabel('X', fontweight = 'bold', fontsize = 14)
plt.legend()
plt.show()

Answer 1

You did not implement the velocity Verlet step correctly, the second velocity update uses the newly computed position, not the old one,

vi[n+1] = v_next - 0.5*omega**2*u[n+1]*dt

This small change should restore second order and energy/amplitude preservation.

---

Also remove the remains of a previous change, there is no x_next to append.