Adaptive Time Stepping Methods - Executable Engineering: Interactive Numerical Methods in Python

It is common for engineering systems to exhibit a combination of fast transients alongside otherwise slow, smooth behaviour. In general, we do not know in advance when a sharp transient will occur or how rapidly it will evolve!

For example, the function $y(t) = \tanh\!\left(\frac{t - 5}{0.1}\right) - t$ combines a steady linear drift with a sudden jump near $t = 5$ . If a uniform step size $h$ is chosen very small, the jump is resolved accurately, but thousands of evaluations are wasted along the quasi-linear regions. Conversely, if $h$ is too large, the transient can be missed or distorted.

import numpy as np
from scipy.integrate import solve_ivp

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

Adaptive Step-Size Control¶

Adaptive time steppers automatically adjust their step size $h$ based on a local estimate of the truncation error.

If the current step of size $h_{\text{current}}$ produces an estimated error $E_{\text{current}}$ and we aim for a target error tolerance $E_{\text{goal}}$ , we can update the step size via:

h_{\text{new}} = h_{\text{current}} \left|\frac{E_{\text{goal}}}{E_{\text{current}}}\right|^{\alpha}

(1)

where $\alpha$ is a scaling exponent tied to the method’s order ( $\alpha \approx 0.2$ is a standard choice for fourth/fifth-order schemes).

Two practical design questions arise in any adaptive solver:

Should we reject the current time step whenever $E_{\text{current}} > E_{\text{tolerance}}$ and recompute from $t_i$ with a smaller $h$ ?
Should we include a meta-analysis of the solver to detect non-convergent or stiff behaviour?

To estimate $E_{\text{current}}$ without excessive computational burden, we typically use one of three strategies:

Step-halving: compare a coarse step to two half-steps.
Embedded Runge-Kutta pairs: compare two integration orders over the same step.
Direct comparison: monitor solution increments or solver convergence metrics.

exe.plot_transient_function(t_range=(0.0, 10.0), n_coarse=12)

Step-Halving Methods¶

Error analysis of particular schemes tell you how the error scales with step size $h$ , but we need an estimate of the error itself.

Step-halving (Richardson extrapolation for ODEs) compares the result of a single full step of size $h$ ( $y_{\text{single step}}$ ) against two successive half-steps of size $h/2$ ( $y_{\text{double step}}$ ).

For the fourth-order Classical RK4 scheme, the difference between the two approximations should scale as $\mathcal{O}(h^5)$ , the local truncation error of the more accurate two-step result is approximately $E \approx \frac{y_{\text{double step}} - y_{\text{single step}}}{15}$

Example: Integrate $\frac{\mathrm{d}y}{\mathrm{d}t} = 4e^{0.8t} - 0.5y$ with $y(0) = 2$ from $t = 0$ to $t = 2$ (exact value $y(2) \approx 14.84392$ ). Note that evaluating one full step and two half-steps requires $4 + 2 \times 4 - 1 = 11$ (or 12 unshared) evaluations of $f(t, y)$ .

differential_equations_module.py

def rk4_step(f, t, y, h):
    k1 = h * f(t, y)
    k2 = h * f(t + 0.5 * h, y + 0.5 * k1)
    k3 = h * f(t + 0.5 * h, y + 0.5 * k2)
    k4 = h * f(t + h, y + k3)
    return y + (k1 + 2.0 * k2 + 2.0 * k3 + k4) / 6.0

exe.display_source(exe.rk4_step)

f = lambda t, y: 4.0 * np.exp(0.8 * t) - 0.5 * y
y_exact_2 = 14.84392190704581

y_onestep = exe.rk4_step(f, 0.0, 2.0, h=2.0)
y_mid = exe.rk4_step(f, 0.0, 2.0, h=1.0)
y_twostep = exe.rk4_step(f, 1.0, y_mid, h=1.0)

est_error = (y_twostep - y_onestep) / 15.0
true_error = y_exact_2 - y_twostep

print(f"One step (h = 2):   {y_onestep:.6f}")
print(f"Two steps (h = 1):  {y_twostep:.6f}")
print(f"Estimated error:    {est_error:.6f}  vs.  True error: {true_error:.6f}")

Runge-Kutta-Fehlberg (RKF45)¶

A more efficient strategy is to compare two different integration orders over the same step by reusing intermediate stage evaluations $k_1, \dots, k_s$ . The classic Runge-Kutta-Fehlberg (RKF45) method uses 6 stages to construct both a 4th-order and a 5th-order approximation:

\begin{array}{c|cccccc} 0 & & & & & & \\ \frac{1}{5} & \frac{1}{5} & & & & & \\ \frac{3}{10} & \frac{3}{40} & \frac{9}{40} & & & & \\ \frac{3}{5} & \frac{3}{10} & -\frac{9}{10} & \frac{6}{5} & & & \\ 1 & -\frac{11}{54} & \frac{5}{2} & -\frac{70}{27} & \frac{35}{27} & & \\ \frac{7}{8} & \frac{1631}{55296} & \frac{175}{512} & \frac{575}{13824} & \frac{44275}{110592} & \frac{253}{4096} & \\ \hline y^{(4)} & \frac{37}{378} & 0 & \frac{250}{621} & \frac{125}{594} & 0 & \frac{512}{1771} \\ y^{(5)} & \frac{2825}{27648} & 0 & \frac{18575}{48384} & \frac{13525}{55296} & \frac{277}{14336} & \frac{1}{4} \end{array}

(2)

The two bottom rows yield the 4th-order and 5th-order updates:

\begin{aligned} y_{i+1}^{(4)} &= y_i + h \left( \frac{37}{378} k_1 + \frac{250}{621} k_3 + \frac{125}{594} k_4 + \frac{512}{1771} k_6 \right) \\ y_{i+1}^{(5)} &= y_i + h \left( \frac{2825}{27648} k_1 + \frac{18575}{48384} k_3 + \frac{13525}{55296} k_4 + \frac{277}{14336} k_5 + \frac{1}{4} k_6 \right) \end{aligned}

(3)

The local error estimate is simply $E = \left|y_{i+1}^{(5)} - y_{i+1}^{(4)}\right|$ . Embedded RK45 pairs are the default workhorse in scientific libraries such as scipy.integrate.solve_ivp.

f = lambda t, y: 4.0 * np.exp(0.8 * t) - 0.5 * y
y_true = lambda t: (4.0 / 1.3) * (np.exp(0.8 * t) - np.exp(-0.5 * t)) + 2.0 * np.exp(-0.5 * t)

sol = solve_ivp(f, (0.0, 2.0), [2.0], method="RK45")
print(f"Adaptive time nodes: {np.round(sol.t, 4)} ({sol.nfev} function evaluations)")
print(f"RK45 y(2) = {sol.y[0, -1]:.6f}  |  Absolute Error = { abs(sol.y[0, -1] - y_true(2.0)):.6e}")

exe.plot_rk45_adaptive(sol, y_true=y_true)

Direct Comparison¶

A less formal, yet often pragmatic, approach uses empirical solution metrics to regulate the step size.

After advancing a time step, the finite-difference rate of change can be evaluated:

\frac{\mathrm{d}y}{\mathrm{d}t} \approx \frac{y_{i+1} - y_i}{h}

(4)

from which a bound on $\left\|\frac{\mathrm{d}y}{\mathrm{d}t}\right\|$ (or on the relative increment $\|y_{i+1} - y_i\|$ ) can be enforced using the Euclidean or $\infty$ -norm. Note that this heuristic bounds the change in the state variable rather than the mathematical truncation error.

In implicit solvers, an even more direct indicator is the convergence rate of the nonlinear root finder at each time step. Because Newton-Raphson converges quadratically when the initial guess $y_i$ is close to $y_{i+1}$ , taking too many Newton iterations signals that the new state has changed too drastically and the step size $h$ should be reduced.