First order Runge-Kutta methods - Executable Engineering: Interactive Numerical Methods in Python

import numpy as np
import plotly.graph_objects as go

try:
    import executable_engineering as exe
except ImportError:
    %pip install -q executable_engineering
    import executable_engineering as exe

Forward Euler Method¶

The Forward Euler method (also known as the Euler-Cauchy / point-slope method, or Explicit Euler) is the simplest Runge-Kutta method. One simply takes the slope at the current point, given by the differential equation, and assumes it remains constant over the step:

Given

\frac{\mathrm{d}y}{\mathrm{d}t} = f(t, y)

(1)

the next time step is:

y_{i+1} = y_i + f(t_i, y_i) h

(2)

which is equivalent to the left Riemann sum, carrying $O(h^2)$ local error over the step and $O(h)$ global error over the full solution.

differential_equations_module.py

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
def forward_euler(f, t0, y0, t_end, h):
    t_points = np.arange(t0, t_end + h / 2, h)
    y_points = np.zeros_like(t_points, dtype=float)
    y_points[0] = y0

    for i in range(len(t_points) - 1):
        y_points[i + 1] = y_points[i] + h * f(t_points[i], y_points[i])

    return {
        "t_points": t_points,
        "y_points": y_points,
        "h": h,
        "method": "Forward Euler",
    }

exe.display_source(exe.forward_euler)

def f(t, y):
    return 4.0 * np.exp(0.8 * t) - 0.5 * y

def y_true(t):
    return (4.0 / 1.3) * (np.exp(0.8 * t) - np.exp(-0.5 * t)) + 2.0 * np.exp(-0.5 * t)

res_euler = exe.forward_euler(f, t0=0.0, y0=2.0, t_end=4.0, h=1.0)
fig = exe.plot_ode_comparison((0, 4), y_true, {"Forward Euler (h=1)": res_euler})
fig.show()

Heun’s Method¶

Now that we have a prediction for $y_{i+1}$ , can we use this as an initial guess and correct our slope estimate? Heun’s method does exactly this and is called a predictor-corrector algorithm.

Consider $\frac{\mathrm{d}y}{\mathrm{d}t} = f(t, y)$ and make a Forward Euler prediction:

y^0_{i+1} = y_i + f(t_i, y_i) h

(3)

Now take the average of the slopes at the start and end of the interval:

\begin{aligned} \bar{y}' &= \frac{y'_i + y'_{i+1}}{2} \\ &= \frac{f(t_i, y_i) + f(t_{i+1}, y^0_{i+1})}{2} \end{aligned}

(4)

which is used in the corrector:

y_{i+1} = y_i + \frac{f(t_i, y_i) + f(t_{i+1}, y^0_{i+1})}{2} h

(5)

differential_equations_module.py

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
def heun_method(f, t0, y0, t_end, h, iterations=1):
    t_points = np.arange(t0, t_end + h / 2, h)
    y_points = np.zeros_like(t_points, dtype=float)
    y_points[0] = y0
    cycle_history = []

    for i in range(len(t_points) - 1):
        t_i, y_i = t_points[i], y_points[i]
        k1 = f(t_i, y_i)
        y_pred = y_i + h * k1
        step_cycles = [y_pred]

        for _ in range(iterations):
            k2 = f(t_i + h, y_pred)
            y_pred = y_i + h * (k1 + k2) / 2.0
            step_cycles.append(y_pred)

        y_points[i + 1] = y_pred
        cycle_history.append(step_cycles)

    return {
        "t_points": t_points,
        "y_points": y_points,
        "cycle_history": cycle_history,
        "h": h,
        "iterations": iterations,
        "method": f"Heun's Method ({iterations} iter)",
    }

exe.display_source(exe.heun_method)

Example: Heun’s Method with Multiple Predictor-Corrector Cycles¶

One could repeat the predictor-corrector cycle iteratively with the intention of converging toward the exact answer, but we will see why this does not eliminate truncation error.

Integrate

y' = 4 e^{0.8 t} - 0.5 y, \quad y(0) = 2

(6)

using Heun’s method with 1 and 15 predictor-corrector cycles ( $h = 1$ ). Unfortunately, iterating Heun’s corrector converges to the fixed point of the trapezoidal discretization, not the exact analytical solution!

res_heun_1 = exe.heun_method(f, t0=0.0, y0=2.0, t_end=4.0, h=1.0, iterations=1)
res_heun_15 = exe.heun_method(f, t0=0.0, y0=2.0, t_end=4.0, h=1.0, iterations=15)

for t_val, y1, y15 in zip(res_heun_1["t_points"][1:], res_heun_1["y_points"][1:], res_heun_15["y_points"][1:]):
    print(f"t = {t_val:.0f} | Heun (1 iter) = {y1:.4f} | Heun (15 iter) = {y15:.4f} | Exact = {y_true(t_val):.4f}")

fig = exe.plot_heun_iterations(f, y_true, t0=0.0, y0=2.0, t_end=4.0, h=1.0, max_iter=15)
fig.show()

Special case: connection to the Trapezoid Rule¶

In the special case where the slope does not depend on the dependent variable, $y'(t) = f(t)$ , the endpoint slope is known directly without a predictor cycle:

y_{i+1} = y_i + \frac{f(t_i) + f(t_{i+1})}{2} h

(7)

which is simply the composite trapezoid rule, carrying $O(h^3)$ local accuracy and $O(h^2)$ global accuracy.

The Midpoint Method¶

Recall from differentiation and integration that sampling at the midpoint is often superior to either endpoint in isolation because leading over- and under-estimates cancel.

As an alternative to Heun’s predictor-corrector method, we subdivide the interval and estimate the state at the midpoint via a half-step:

y_{i+1/2} = y_i + f(t_i, y_i) \frac{h}{2}

(8)

We then evaluate the slope at the midpoint, $y'_{i+1/2} = f(t_{i+1/2}, y_{i+1/2})$ , to advance the full step:

y_{i+1} = y_i + f(t_{i+1/2}, y_{i+1/2}) h

(9)

which achieves $O(h^3)$ local and $O(h^2)$ global error. The midpoint method demonstrates that we can calculate and combine substeps to achieve higher-order predictions.

Optimal construction of substeps will be the foundation for the general Runge-Kutta methods.

differential_equations_module.py

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
def midpoint_ode(f, t0, y0, t_end, h):
    t_points = np.arange(t0, t_end + h / 2, h)
    y_points = np.zeros_like(t_points, dtype=float)
    y_points[0] = y0

    for i in range(len(t_points) - 1):
        t_i, y_i = t_points[i], y_points[i]
        k1 = f(t_i, y_i)
        k2 = f(t_i + 0.5 * h, y_i + 0.5 * h * k1)
        y_points[i + 1] = y_points[i] + h * k2

    return {
        "t_points": t_points,
        "y_points": y_points,
        "h": h,
        "method": "Midpoint Method (RK2)",
    }

exe.display_source(exe.midpoint_ode)

res_mid = exe.midpoint_ode(f, t0=0.0, y0=2.0, t_end=4.0, h=1.0)
fig = exe.plot_ode_comparison(
    (0, 4),
    y_true,
    {"Forward Euler": res_euler, "Heun's Method": res_heun_1, "Midpoint Method": res_mid},
)
fig.show()