0
$\begingroup$

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.

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

| cite | improve this question | | | | |
$\endgroup$
  • 1
    $\begingroup$ 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. $\endgroup$ – Lutz Lehmann Mar 13 at 20:36
5
$\begingroup$

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
| cite | improve this answer | | | | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.