(7.2)

For a simple harmonic oscillator y''+y = 0, with initial conditions y(0) = 1; y'(0) = 0, find y(t) from t = 0 to 100pi(). Use an Euler method and a fixed-step fourth-order Runge-Kutta method. For each method check how the average local error, and the error in the final value and slope, depend on the step size.

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
    

h = 0.01
t = np.arange(0, 100*np.pi, h)

# Euler method
yt = [1]  # y(t) = 1, when t=0: inital value
ydt = [0]  # y''(t) = 0, when t=0: inital value

E_error = 0
for i in t:
    # compute new y(x) and new y'(x)
    yt_n = yt[-1] + (h * ydt[-1])
    ydt_n = ydt[-1] - (h * yt[-1])
    yt.append(yt_n)
    ydt.append(ydt_n)
    E_error += np.abs(np.cos(i) - yt[-1])

plt.plot(t, np.cos(t))
plt.plot(t, yt[:-1])
plt.show()

     This is error accumulated plot

# Runge-Kutta method method
y0 = [1]
y1 = [0]
k1 = [0,0]
k2 = [0,0]
k3 = [0,0]
k4 = [0,0]

RK_error = 0
for i 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) )
    RK_error += np.abs(np.cos(i) - y0[-1])

plt.plot(t, np.cos(t))
plt.plot(t, y0[:-1])
plt.show()
print("     This is error accumulated plot")

     This is error accumulated plot