For a simple harmonic oscillator $\ddot{y} + y = 0$, with initial conditions $y(0) = 1, \dot{y}(0) = 0$, find $y(t)$ from $t = 0$ to $100 \pi$. 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.

We know that the exact solution of this system is:

$$y(t) = \cos{t} \tag{1}$$

To solve this numerically, we have two variables to integrate: $y$ and $\dot{y}$. We'll use this state vector:

$$z = \begin{bmatrix}y \\ \dot{y}\end{bmatrix}$$

The system can then be formulated in a matrix form:

$$\begin{bmatrix}\dot{y} \\ \ddot{y}\end{bmatrix} = f(t, z(t)) = \begin{bmatrix}0 & 1\\ -1 & 0\end{bmatrix} \begin{bmatrix}y \\ \dot{y}\end{bmatrix}$$

Let's solve this numerically with both an Euler integrator, and a Runge Kutta of order 4. Thanks to Python's callable objects, I can just pass a reference to either function and use it as the integrator:

import numpy as np
import math
import matplotlib.pyplot as plt


def runge_kutta4_method(f, x, y, h):
    k1 = h * f(x, y)
    k2 = h * f(x + h / 2, y + k1 / 2)
    k3 = h * f(x + h / 2, y + k2 / 2)
    k4 = h * f(x + h, y + k3)
    y_new = y + k1 / 6 + k2 / 3 + k3 / 3 + k4 / 6
    return y_new


def euler_method(f, x, y, h):
    return y + h * f(x, y)


def integrate(f, x_axis, y0, func_method):
    n = len(x_axis)
    n_var = len(y0)
    y_all = np.zeros((n, n_var), dtype=np.float64)
    y = y0
    y_all[0, :] = y0
    for i in range(n - 1):
        x = x_axis[i]
        h = x_axis[i + 1] - x
        y = func_method(f, x, y, h)
        y_all[i + 1, :] = y
    return y_all


def experiment(h, func_method):
    f = lambda x, y: np.array([y[1], -y[0]], dtype=np.float64)
    x_axis = np.arange(0, 100 * math.pi, h)
    n = len(x_axis)
    y0 = np.array([1, 0], dtype=np.float64)
    y_all = integrate(f, x_axis, y0, func_method)
    y_true = np.zeros((n, 2))
    y_true[:, 0] = np.cos(x_axis)
    y_true[:, 1] = -np.sin(x_axis)
    mse_err = np.mean(np.square(y_true[:, 0] - y_all[:, 1]))
    rmse_err = math.sqrt(mse_err)
    last_val_err = abs(y_true[-1, 0] - y_all[-1, 0])
    last_slope_err = abs(y_true[-1, 1] - y_all[-1, 1])
    return rmse_err, last_val_err, last_slope_err

h_axis = np.linspace(0.002, 0.15, 64)
n_h = len(h_axis)
rmse_err_euler = np.zeros((n_h,), dtype=np.float64)
val_err_euler = np.zeros((n_h,), dtype=np.float64)
slope_err_euler = np.zeros((n_h,), dtype=np.float64)
rmse_err_rk4 = np.zeros((n_h,), dtype=np.float64)
val_err_rk4 = np.zeros((n_h,), dtype=np.float64)
slope_err_rk4 = np.zeros((n_h,), dtype=np.float64)

for i, h in enumerate(h_axis):
    rmse_err_euler[i], val_err_euler[i], slope_err_euler[i] = experiment(h, euler_method)
    rmse_err_rk4[i], val_err_rk4[i], slope_err_rk4[i] = experiment(h, runge_kutta4_method)

plt.figure()
plt.semilogy(h_axis, rmse_err_euler, "b-")
plt.semilogy(h_axis, rmse_err_rk4, "r-")
plt.grid("on")
plt.xlabel("step size")
plt.ylabel("rmse")
plt.title("Average error")
plt.legend(["Euler", "Runge Kutta 4"])
plt.show()

plt.figure()
plt.semilogy(h_axis, val_err_euler, "b-")
plt.semilogy(h_axis, val_err_rk4, "r-")
plt.grid("on")
plt.xlabel("step size")
plt.ylabel("abs. error")
plt.title("Error in the final value")
plt.legend(["Euler", "Runge Kutta 4"])
plt.show()

plt.figure()
plt.semilogy(h_axis, slope_err_euler, "b-")
plt.semilogy(h_axis, slope_err_rk4, "r-")
plt.grid("on")
plt.xlabel("step size")
plt.ylabel("abs. error")
plt.title("Error in the final slope")
plt.legend(["Euler", "Runge Kutta 4"])
plt.show()