Euler's Method for Simple Harmonic Motion

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

# 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

# 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)

def rk_adaptive():
    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__":
    rk_adaptive()

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

    
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()