# Simple Harmonic Motion using the leapfrog method

I have to use the leapfrog method to solve the simple harmonic oscillator and I having trouble writing it in code. This is what we were given in class

$$\frac{v_{n+1/2}-v_{n-1/2}}{\Delta t}=-\omega_0^2x_n\\ \frac{x_{n+1}-x_n}{\Delta t}=v_{n+1/2}$$

This is what I have written for the code

import numpy as np
import matplotlib.pyplot as plt

m = 1
k = 1
N = 1000  # steps
x0 = 1  # initial position
w=1

t = np.linspace(0,100,N) #timestep
dt = t[1]-t[0]

v = 0  # initial velocity
x = x0  # initial position

# integration: calculate  v, x for each step
x = np.zeros(N)
v = np.zeros(N)
for i in range(N-1):

v[i+1] = 0.5*v[i]*dt + (-w^2)
x[i+1] = v[i+1] * dt + x[i]

plt.figure()
plt.plot(t, x)

plt.xlabel('time')
plt.ylabel('displacement')
plt.show()

From the above code, this is the graph I get and it is nothing like a simple harmonic motion. The exercise is to get a numerical solution for the simple harmonic oscillator using the leapfrog method from the image at the top. I am new to python and not sure where I went wrong or where to go next so any help would be appreciated

• w^2 in python is a bit-logical operation, you want the square w**2. If you do not initialize $v_{0-\frac12}$ correctly, you will not get the second order Leapfrog Verlet but the first order symplectic Euler method. – Lutz Lehmann Mar 13 at 20:36

Your code has a few issues. First of all, recasting the equations you wrote you have the following relations:

$$v_{n + 1/2} = v_{n - 1/2} - \omega_02 x_n \Delta t\\ x_{n + 1} = x_n + v_{n + 1/2} \Delta t$$

which is not what you implemented. The two lines where you do a step forward should read:

v[i+1] = v[i] + (-w**2) * x[i] * dt
x[i+1] = v[i+1] * dt + x[i]

If you now run the code, you'll obtain a straight line, since you did not set the initial position correctly. Indeed, your overwrite the initial conditions when you set both x and v to zero. The initialisation should read instead:

# integration: calculate  v, x for each step
x = np.zeros(N)
v = np.zeros(N)
x[0] = x0
v[0] = 0