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 runge_kutta4_adaptive_method(f, x, y, h, err):
    converged = False
    y_found = None
    h_found = None
    while not converged:
        y_new_full = runge_kutta4_method(f, x, y, h)
        y_new_half = runge_kutta4_method(f, x, runge_kutta4_method(f, x, y, h/2), h/2)
        err_new = abs(y_new_full[0]-y_new_half[0])
        if err_new > err:
            if y_found is None:
                h /= 2
            else:
                converged = True
        else:
            y_found = y_new_half
            h_found = h
            h *= 2
    return y_found, h_found


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_adaptive(err):
    f = lambda x, y: np.array([y[1], -y[0]], dtype=np.float64)
    x = 0
    y = np.array([1, 0], dtype=np.float64)
    x_final = 100 * math.pi
    # initial h, arbitrary
    h = 0.01
    h_list = []
    y_list = []
    while x < x_final:
        y, h = runge_kutta4_adaptive_method(f, x, y, h, err)
        y_list.append(y[0])
        h_list.append(h)
        x += h
    return y_list, h_list

Let's examine the evolution of the step size during the integration:

y_list, h_list = experiment_adaptive(0.0001)

plt.figure()
plt.plot(h_list[:100])
plt.xlabel("steps")
plt.ylabel("step size")

plt.figure()
plt.plot(y_list[:100])
plt.xlabel("steps")
plt.ylabel("y")

plt.show()

The adaptive step size is minimal when the derivative is maximal, showing a good adaptation to the sensitivity of the problem. Let's now examine the mean step size over the whole simulation, while getting increasingly less strict about the desired error:

err_list = np.linspace(0.000001, 0.0001, 128)
n = len(err_list)
mean_step_list = np.zeros((n,), dtype=np.float64)
for i, err in enumerate(err_list):
    y_list, h_list = experiment_adaptive(err)
    mean_step_list[i] = np.mean(h_list)

plt.figure()
plt.plot(err_list, mean_step_list)
plt.xlabel("Desired error")
plt.ylabel("Mean step size")
plt.show()