# Euler's Method for Simple Harmonic Motion¶

In [22]:
import numpy as np
import matplotlib.pyplot as plt

def euler_shm():
h = 0.01
t = np.arange(0, 100*np.pi, h)
# init values for y(x) and y'(x)
yx = [1]
ypx = [0]
cum_error = 0

for ts in t:
# compute new y(x) and new y'(x)
yx_new = yx[-1] + (h * ypx[-1])
ypx_new = ypx[-1] - (h * yx[-1])
yx.append(yx_new)
ypx.append(ypx_new)
cum_error += np.abs(np.cos(ts) - yx[-1])

plt.plot(t, np.cos(t))
plt.plot(t, yx[:-1])
print "Cumulative Error", cum_error
plt.show()

if __name__ == "__main__":
euler_shm()
Cumulative Error 28518.200244

# Runge-Kutta Method for Simple Harmonic Motion¶

In [20]:
# yx -> y0
# ypx -> y1; Neil's notation easier to follow

def rk_shm():
h = 0.01
t = np.arange(0, 100*np.pi, h)
# init values
y0 = [1]
y1 = [0]
k1 = [0,0]
k2 = [0,0]
k3 = [0,0]
k4 = [0,0]
cum_error = 0

for ts in t:
k1[0] = h * y1[-1]
k1[1] = -h * y0[-1]
k2[0] = h * (y1[-1] + k1[1]/2.0)
k2[1] = -h * (y0[-1] + k1[0]/2.0)
k3[0] = h * (y1[-1] + k2[1]/2.0)
k3[1] =-h * (y0[-1] + k2[0]/2.0)
k4[0] = h * (y1[-1] + k3[1]/2.0)
k4[1] =-h * (y0[-1] + k3[0]/2.0)
y0.append(y0[-1] + (k1[0] / 6.0) + (k2[0] / 3.0) + (k2[0] / 3.0) + (k4[0] / 6.0) )
y1.append(y1[-1] + (k1[1] / 6.0) + (k2[1] / 3.0) + (k2[1] / 3.0) + (k4[1] / 6.0) )
cum_error += np.abs(np.cos(ts) - y0[-1])

plt.plot(t, np.cos(t))
plt.plot(t, y0[:-1])
print "Cumulative Error", cum_error
plt.show()

if __name__ == "__main__":
rk_shm()
Cumulative Error 2881.64536239

# Runge-Kutta 4th Order with Adaptive Stepper¶

In [24]:
# Incorrect - the y value diverging from cos solution?

def make_step(si, sf, s):
k1 = [0,0]
k2 = [0,0]
k3 = [0,0]
k4 = [0,0]

k1[0] = s * si[1]
k1[1] = -s * si[0]
k2[0] = s * (si[1] + k1[1]/2.0)
k2[1] = -s * (si[0] + k1[0]/2.0)
k3[0] = s * (si[1] + k2[1]/2.0)
k3[1] =-s * (si[0] + k2[0]/2.0)
k4[0] = s * (si[1] + k3[1])
k4[1] =-s * (si[0] + k3[0])
sf[0] = (si[0] + (k1[0] / 6.0) + (k2[0] / 3.0) + (k2[0] / 3.0) + (k4[0] / 6.0) )
sf[1] = (si[1] + (k1[0] / 6.0) + (k2[0] / 3.0) + (k2[0] / 3.0) + (k4[0] / 6.0) )
return (si, sf)

allowed_error = 0.1
scale_down = 1.3
scale_up = 1.2
h = 0.1
t = 2 * np.pi
# init values
y0 = [1,0]
y1 = [0,0]
y2 = [0,0]
y3 = [0,0]
y4 = [0,0]

error = 0
out = []
time = []
i = 0
while i < t:
#         print i
y0, y1 = make_step(y0, y1, h/2.0)
y1, y2 = make_step(y1, y2, h/2.0)
y0, y3 = make_step(y0, y3, h)
error = np.abs(y3[0] - y2[0])
#         print error
if error > allowed_error:
h = h / scale_down
else:
y0[0] = y2[0]
y0[1] = y2[1]
print y0[0]
out.append(y0[0])
time.append(i)
i+=h
if error < allowed_error * 0.20:
h = h * scale_up

plt.plot(time, np.cos(time))
print len(time)
print len(out)
plt.plot(time, out)
print "Error", error
plt.show()

if __name__ == "__main__":
0.997439608721
0.993429383975
0.987124534322
0.977140725213
0.961162529291
0.935223965544
0.892344427356
0.819854219283
0.721861210113
0.554928437003
0.260204286158
-0.281125798592
-0.934302533198
-1.66577715277
-2.42968163753
-3.18053233579
-4.09306334051
-4.94621571296
-5.94107146687
-7.10116666562
-8.14213903223
-9.31632798052
-10.6407816649
-12.1347299469
-13.4314704687
-14.8559870599
26
26
Error 0.0669775850878

# Numerically solve pendulum problem¶

In [ ]:

def rk_shm_pendulum():
h = 0.01
t = np.arange(0, 100*np.pi, h)
l = 2
g = 9.81

# init values
#     y0 = [1]
#     y1 = [0]

# initial values should be a function of l, g, and z?

k1 = [0,0]
k2 = [0,0]
k3 = [0,0]
k4 = [0,0]
cum_error = 0

for ts in t:
k1[0] = h * y1[-1]
k1[1] = -h * y0[-1]
k2[0] = h * (y1[-1] + k1[1]/2.0)
k2[1] = -h * (y0[-1] + k1[0]/2.0)
k3[0] = h * (y1[-1] + k2[1]/2.0)
k3[1] =-h * (y0[-1] + k2[0]/2.0)
k4[0] = h * (y1[-1] + k3[1]/2.0)
k4[1] =-h * (y0[-1] + k3[0]/2.0)
y0.append(y0[-1] + (k1[0] / 6.0) + (k2[0] / 3.0) + (k2[0] / 3.0) + (k4[0] / 6.0) )
y1.append(y1[-1] + (k1[1] / 6.0) + (k2[1] / 3.0) + (k2[1] / 3.0) + (k4[1] / 6.0) )
cum_error += np.abs(np.cos(ts) - y0[-1])

plt.plot(t, np.cos(t))
plt.plot(t, y0[:-1])
print "Cumulative Error", cum_error
plt.show()

if __name__ == "__main__":
rk_shm()