Using RK2 Method to solve the simple harmonic oscillator of a horizontal mass on a spring (1D)

Being new to numerical analysis techniques, in particular RK2, I decided the best way to jump in is by using python to solve the well known mass-spring oscillator using RK2 techniques.

My problem is that the step size seems to influence the period of my solution and I do know why and would appreciate it if someone could help me better understand what I am doing wrong. Commenting on the code is welcomed! Thank you.

The problem is as follows (this is not a textbook or homework problem!): A mass is attached to a spring with a linear spring constant. At time $t=0$ the mass is displaced from equilibrium, $x=0$, to the position $x(t=0) = 1 cm$. The mass is released and left to oscillate about equilibrium. Use the method of RK2 to plot the evolution of the masses position as a function of time. The required first order differential equations for the problem are:

$\dot{x} = -\omega x(t)$

$\frac{dx}{dt} = \dot{x}$

Where the RK2 functions I am using are as follows:

$v(t_{n+1}) = v(t_n) -\omega x(t_n)\delta t$

$x(t_{n+1}) = x(t_n) + v(t_n)\delta t$

$k2_{x} = \delta t (v(t_n) - \omega x(t_n)\frac{\delta t}{2})$

$k2_{v} = \delta t (\omega x(t_n) + v(t_n) \frac{\delta t}{2})$

The last two function are the midpoint values that I think are used for the RK2 approach.

Along with the implemented code :

import numpy as np
import matplotlib.pyplot as plt

print "RK2 Method for Oscillating Spring"

#Spring Constant
k = 1.0
#Mass
m = 0.02
#Angular frequency squared
w = k/m

period = 2*np.pi*(np.sqrt(1/w))

#Initial Conditions
x_start = 1.0
v_start = 0.0

#Set loop variables
x_position = x_start
velocity = v_start

#Step Size
dt = 0.001

#Arrays to store variable values
x_plot = []
v_plot = []
time = []

#Plot the numerically determined position over
# the time interval 1 to 10 seconds in steps of 0.1
for t in np.arange(1,10,0.1):

#The two equations below are used in the Euler Method approach
#Their solution's period, which can be plotted, varies depending on the step size
# and I know it shouldn't
#x_position = x_position + dt*velocity
#velocity = velocity - dt*w*x_position

#The equations below are used for RK2 method. I don't know why I cant get it to work.
velocity = velocity + dt*(-w*x_position)
k2_x = dt*(velocity + (-w*x_position*dt/2.0))
x_position = x_position + k2_x

k2_v = -dt*w*(x_position + velocity*dt/2.0)
velocity = velocity + k2_v

x_plot.append(x_position)
v_plot.append(velocity)
time.append(t)

plt.subplot(3,1,1)
plt.plot(time,x_plot)
plt.subplot(3,1,2)
plt.plot(time,v_plot)
plt.subplot(3,1,3)
plt.plot(x_plot,v_plot) # illustrate divergence over number of steps
plt.show()

My problem is that I don't understand why varying the step size, dt, that the period of the solution, a sinusoidal function, changes.

• Numerical integration schemes have what's called phase and amplitude errors in the produced solution. For your problem, the phase error will affect the period of your solution. As you change $\Delta t$, you will see a change in the phase error and in turn a change in the period, so what you're seeing is likely okay. – spektr May 14 '17 at 4:25
• @Wrzlprmft I cleaned up the post using your suggestions. Hopefully it is easier to read now. Do you have any ideas for why I am experiencing the behavior using the RK2 approach? – kuantumbro May 14 '17 at 14:18
• @C.Howard Interesting I didn't know that was an issue. Is that because of the product of the angular frequency and dt term in the first velocity expression I have written above the segment of code? Could you also further explain or provide a link to a discussion on the phase errors? Thanks for your help! – kuantumbro May 14 '17 at 17:31
• @kuantumbro No, the stuff I describe is a property of any time integration scheme, though some have better error properties than others. I will see if I can provide a link or something later. – spektr May 14 '17 at 17:57
• @kuantumbro: It’s better, but there are still some issues left: 1) Your example code produces a constant solution, presumably because you messed up the variable names: There is an xold and and xtold; also everything is old. 2) The constants $x_{k2}$ and $v_{k2}$ are not used in your abstract description of the algorithm (the second bunch of equations). 3) Please describe the phenomenon in question in detail, so we can easily reproduce it: For what parameters do you observe which period lengths? – Wrzlprmft May 15 '17 at 7:45