Part 2: Building Better Methods - Runge-Kutta

Learning Objectives¶

By the end of this section, you will be able to:

• Derive explicit midpoint RK2 and explain the stage dependencies of RK4

• Use the mathematical specification to write and verify your own RK implementations

• Connect RK weights to quadrature rules from Module 2

• Trace adaptive step acceptance and rejection using a worked Euler example

• Analyze when higher-order methods help vs harm

• Choose appropriate RK methods based on problem requirements

The Key Insight: Sample Multiple Points¶

Euler extends the initial tangent through the step. Curvature introduces a one-step defect, which motivates sampling additional slopes. Runge-Kutta methods evaluate the derivative at carefully chosen intermediate points, then combine these samples with specific weights to achieve higher-order accuracy.

This connects to quadrature from Module 2: We’re essentially performing numerical integration on:

But now the integrand depends on the unknown solution $x(t)$ itself, requiring us to estimate intermediate values.

RK2 - The Midpoint Method¶

Geometric Intuition¶

Instead of blindly following the tangent from the starting point, we:

1. Take a half-step using the initial derivative

2. Evaluate the derivative at this midpoint

3. Use this midpoint derivative for the full step

This is analogous to the midpoint rule for integration, which we know from Module 2 is second-order accurate!

Complete Mathematical Formulation¶

The RK2 midpoint method is a predictor-corrector method:

Stage 1 - Predictor:

Stage 2 - Corrector:

Error Analysis via Taylor Series¶

Differentiate the ODE along its solution. For a vector state, let $J=\partial f/\partial x$ be the state Jacobian. Then

so the exact solution begins

Try two stages, with $k_1=f(t,x)$ and $k_2=f(t+ch,x+ah f)$, and an update $x+h(b_1k_1+b_2k_2)$. A multivariate expansion gives

and the resulting update is

Matching the independent terms requires

Choose $b_2=1$. Then $b_1=0$ and $a=c=1/2$, giving the explicit midpoint formulas above. The midpoint slope consequently satisfies

Another choice, $b_1=b_2=1/2$ and $a=c=1$, gives Heun’s explicit trapezoidal method. Both are second-order RK methods, but their intermediate states differ. This reading uses explicit midpoint when it says RK2. On a linear autonomous system they have the same one-step map; a nonlinear test can separate them. For $x'=x^2$, $x(0)=1$, their one-step values are $1+h+h^2+h^3/4$ and $1+h+h^2+h^3/2$, respectively. Neither cubic coefficient matches the exact local expansion of $1/(1-h)$; both methods match through second order. Do not mix their stages and weights.

The RK2 update matches the true solution through $O(h^2)$:

• Local truncation error: $O(h^3)$

• Global error: $O(h^2)$

In a resolved truncation regime at fixed duration, halving the timestep reduces the leading global error by approximately four.

From mathematics to your implementation¶

The stage equations are the specification you will translate into Python. Write the RK2 and RK4 implementations yourself. We will use complete Euler examples to demonstrate array handling, stepping, and diagnostics; a working RK routine is not supplied here.

For Python, use independent variable first: f(t, state) returns the rate of change of the NumPy-array state. In the equations above, $x_n$ denotes that state; it can have several components. Addition and scalar multiplication then operate on the entire array. The stage slopes must have the same shape as the state. Their units are state units divided by time units; multiplying a slope by $h$ gives a state increment.

Separate two responsibilities:

First explain the method in mathematical statements, then write your own code. Before writing the loop, answer: which quantities belong to the original state, and which belong to a temporary stage? If a stage accidentally overwrites the starting state, later formulas describe a different method.

Do not judge correctness from a smooth-looking trajectory alone. A constant right-hand side checks basic increments. A time-dependent right-hand side tests whether stage times are used. A nonlinear state-dependent right-hand side tests whether intermediate states are used. Refinement over a fixed interval then checks the claimed global order. These checks probe different mistakes; none is a replacement for deriving the update.

RK4 - The Classical Workhorse¶

RK4 samples the derivative at four carefully chosen points:

Stage 1: $k_1 = f(t_n, x_n)$ (derivative at start)

Stage 2: $k_2 = f(t_n + \frac{h}{2}, x_n + \frac{h}{2}k_1)$ (derivative at midpoint using $k_1$)

Stage 3: $k_3 = f(t_n + \frac{h}{2}, x_n + \frac{h}{2}k_2)$ (derivative at midpoint using $k_2$)

Stage 4: $k_4 = f(t_n + h, x_n + hk_3)$ (derivative at endpoint using $k_3$)

Final update:

The Connection to Simpson’s Rule¶

For a time-only right-hand side $x'=g(t)$, the two midpoint evaluations are identical. The RK4 update then becomes

which is exactly Simpson quadrature over the time interval. Euler similarly reduces to left rectangles and explicit midpoint RK2 to midpoint quadrature.

For $f(t,x)$, however, $k_2$ and $k_3$ generally differ: they evaluate different predicted states at the same time. Simpson’s weights alone do not establish fourth-order ODE accuracy. The stage locations and dependencies must satisfy additional order conditions involving state derivatives of $f$.

Check: for $x'=x$, start at $x=1$, $t=0$ with $h=0.1$. Then $k_1=1$, $k_2=1.05$, $k_3=1.0525$, and $k_4=1.10525$. The unequal midpoint slopes are precisely what the state-dependent problem adds to the quadrature analogy.

Error Analysis¶

Through Taylor series expansion to fourth order:

• Local truncation error: $O(h^5)$

• Global error: $O(h^4)$

In a resolved truncation regime at fixed duration, halving the timestep reduces the leading global error by approximately sixteen.

The General Runge-Kutta Framework¶

All explicit RK methods follow this pattern:

The coefficients $(a_{ij}, b_i, c_i)$ define the method and are displayed in a Butcher tableau:

c₁ |
c₂ | a₂₁
c₃ | a₃₁ a₃₂
⋮  | ⋮
cₛ | aₛ₁ aₛ₂ ... aₛ,ₛ₋₁
---|--------------------
   | b₁  b₂  ...  bₛ

For RK4:

0   |
1/2 | 1/2
1/2 | 0   1/2
1   | 0   0   1
----|----------------
    | 1/6 2/6 2/6 1/6

Follow the dependencies, not just the weights¶

The two midpoint stages occur at the same time but generally at different predicted states. The second midpoint evaluation revises the slope information used to predict the endpoint.

Every stage state in the formulas is built from the original $x_n$ plus a specified increment. The arrows show dependencies; they do not mean that stage states are advanced sequentially by a full ODE step. The two midpoint states are generally not points on the exact trajectory either.

Adaptive Timestep Control¶

Real problems have varying timescales. When a planet is far from the sun, we can take large steps. During close encounters, we need tiny steps. Adaptive methods adjust $h$ dynamically.

The Embedded Runge-Kutta Approach¶

The idea: compute two different order estimates using the same function evaluations, then use their difference to estimate error.

An embedded example: Fehlberg 4(5)¶

The original six-stage Fehlberg pair computes both fourth- and fifth-order solutions. The name RK45 is also used for other pairs: SciPy uses Dormand–Prince, not these six stages. For the Fehlberg example:

1. Compute six $k_i$ values (shared between both estimates)

2. Combine with 4th-order weights $\to$ $x_{(4)}$

3. Combine with 5th-order weights $\to$ $x_{(5)}$

4. Error estimate: $\epsilon = ||x_{(5)} - x_{(4)}||$

5. Suggested next timestep (with a safety factor): $h_{new} = h \cdot 0.9 \cdot (\text{tolerance}/\epsilon)^{1/5}$

For vectors, first scale each component by an allowance such as atol + rtol * abs(reference_state) and use a stated norm. Bound the step-growth factor and handle a zero estimate explicitly. If the scaled estimate is acceptable, accept the step and continue. Otherwise, reject and retry with smaller $h$.

Watch an adaptive controller make a decision¶

We can study acceptance and rejection using Euler steps, while leaving the RK implementations to you. Consider the dimensionless decay problem $u'=-u$, $u(0)=1$, with exact solution $u(t)=e^{-t}$. It gives us a reference for checking the numerical experiment without introducing another physical model.

Take one Euler step of length $h$ and, separately, two Euler steps of length $h/2$, both starting from the same $(t,u)$. If the leading one-step defect is $Ch^2$, the two-half-step result has approximately half the defect of the coarse result. Their difference therefore estimates the two-half-step local error in the smooth, small-step regime:

This is still a first-order method. We accept $u_{\rm fine}$ itself; we do not apply an extrapolation correction. The comparison estimates a local error and does not reveal the exact accumulated error of the trajectory.

For this example, define the scale and normalized estimate explicitly:

Accept when $\eta\le1$. Rejection leaves both the accepted time and accepted state unchanged. The proposed next step uses $h_{\rm new}\propto h\eta^{-1/2}$ because this estimator scales as $h^2$. The exponent would change for another estimator. Growth limits and a safety factor make the controller less abrupt; they do not turn an estimate into a mathematical bound.

Predict before running: Will the initial step of length 1 be accepted? When the solution becomes small, which part of the allowance will eventually matter more? Does every accepted step imply that the final solution is within one tolerance of the exact answer?

import numpy as np
import matplotlib.pyplot as plt


def euler_trial(f, t, u, h):
    """Compare one Euler step with two Euler half-steps; scalar state."""
    coarse = u + h*f(t, u)
    half = u + (h/2)*f(t, u)
    fine = half + (h/2)*f(t + h/2, half)
    return fine, abs(fine-coarse)


def adaptive_euler_decay(rtol, atol, initial_h=1.0, end_time=4.0):
    """A teaching controller for u'=-u; accepts the two-half-step state."""
    if not np.all(np.isfinite([rtol, atol, initial_h, end_time])):
        raise ValueError("Parameters must be finite")
    if rtol < 0 or atol <= 0 or initial_h <= 0 or end_time <= 0:
        raise ValueError("Require rtol >= 0 and positive atol, step, and end time")
    t, u, h = 0.0, 1.0, initial_h
    times, values, trials = [t], [u], []
    for attempt in range(10000):
        if t >= end_time:
            return np.array(times), np.array(values), np.array(trials)
        h = min(h, end_time-t)  # Land on the requested final time.
        if t + h == t:
            raise RuntimeError("The proposed step cannot advance floating-point time")
        fine, estimate = euler_trial(lambda time, value: -value, t, u, h)
        scale = atol + rtol*max(abs(u), abs(fine))
        ratio = estimate/scale
        if not np.isfinite(ratio):
            raise RuntimeError("Non-finite error estimate")
        accepted = ratio <= 1.0
        trials.append((t, h, ratio, accepted))
        if accepted:
            t, u = t+h, fine
            times.append(t)
            values.append(u)
        # Zero estimated error permits bounded growth, not an infinite step.
        factor = 2.0 if ratio == 0 else np.clip(0.9*ratio**(-0.5), 0.2, 2.0)
        h *= factor
    raise RuntimeError("Attempt limit reached")


adaptive_t, adaptive_u, adaptive_trials = adaptive_euler_decay(1e-3, 1e-5)
print("  start t    trial h    scaled estimate   accept?")
for t, h, ratio, accepted in adaptive_trials[:6]:
    print(f"{t:9.4f}  {h:9.4f}  {ratio:15.4f}   {bool(accepted)}")
print("Final actual absolute error:", abs(adaptive_u[-1]-np.exp(-adaptive_t[-1])))

  start t    trial h    scaled estimate   accept?
   0.0000     1.0000         247.5248   False
   0.0000     0.2000           9.9010   False
   0.0000     0.0572           0.8100   True
   0.0572     0.0572           0.8095   True
   0.1144     0.0572           0.8095   True
   0.1716     0.0572           0.8095   True
Final actual absolute error: 0.0010884703398220882

The first three statements in euler_trial are exactly Euler updates. This transparent version evaluates the starting derivative twice; caching it would save one call per trial, without changing the method. The controller records rejected trials as well as accepted ones, so the work cannot be inferred just from the number of stored solution points.

accepted_trials = adaptive_trials[:, 3].astype(bool)
fig, axes = plt.subplots(2, 1, figsize=(8, 6), layout="constrained")
axes[0].plot(adaptive_t, adaptive_u, label="accepted Euler half-steps")
axes[0].plot(adaptive_t, np.exp(-adaptive_t), "--", label="exact decay")
axes[0].set(xlabel="dimensionless time", ylabel="dimensionless u")
axes[0].legend()
for mask, label, marker in ((accepted_trials, "accepted", "o"),
                            (~accepted_trials, "rejected", "x")):
    axes[1].semilogy(adaptive_trials[mask, 0], adaptive_trials[mask, 1],
                     marker, ms=4, label=label)
axes[1].set(xlabel="start time of trial", ylabel="trial step length")
axes[1].legend()
plt.show()

for requested_rtol in (1e-2, 1e-3, 1e-4):
    ts, us, trials = adaptive_euler_decay(requested_rtol, 1e-5)
    actual = abs(us[-1]-np.exp(-ts[-1]))
    final_allowance = 1e-5 + requested_rtol*abs(us[-1])
    print(f"rtol={requested_rtol:.0e}: accepted={len(ts)-1}, attempts={len(trials)}, "
          f"final error={actual:.3e}, error/final allowance={actual/final_allowance:.2f}")

rtol=1e-02: accepted=22, attempts=23, final error=3.226e-03, error/final allowance=20.05
rtol=1e-03: accepted=67, attempts=69, final error=1.088e-03, error/final allowance=39.98
rtol=1e-04: accepted=163, attempts=166, final error=4.780e-04, error/final allowance=40.57

The same controller logic can use an embedded RK error estimate, with the appropriate estimator order and vector scaling. Constructing the RK stages remains your work. A request for evenly spaced output is also different from a request for evenly spaced internal steps: established libraries can interpolate between accepted steps. SciPy’s solver documentation distinguishes these roles and describes component-dependent absolute tolerances.

When to Use Adaptive Methods¶

Good for:

• Problems with varying timescales

• Unknown dynamics

• Achieving specified accuracy efficiently

• Initial exploration of new systems

Bad for:

• Long-term Hamiltonian systems when ordinary state-dependent step selection would compromise a fixed-step symplectic method

• Parallel simulations (variable timesteps complicate synchronization)

Compare accuracy and work with a defined benchmark¶

A comparison needs a defined system, initial state, duration, step size, error metric, and cost metric. We use the dimensionless simple harmonic oscillator (SHO), $q'=p$, $p'=-q$, with $(q(0),p(0))=(1,0)$. Its exact solution is $(\cos t,-\sin t)$ and its period is $2\pi$.

The reference plots below retain the comparisons among Euler, explicit midpoint RK2, classical RK4, and kick–drift–kick leapfrog. The figures were generated separately; no RK implementation is supplied in the page or an accompanying notebook. Use the mathematics to construct your own methods, and use the observables to decide what your implementation needs to establish.

One oscillator, four diagnostics¶

Each method uses 100 steps per period for 20 periods. Before examining the figure, predict which method increases energy, which damps it, and whether bounded energy guarantees accurate phase. The inset enlarges the trajectories near the unit circle; Euler’s growth would otherwise hide their differences.

For this oscillator, the diagnostic definitions are

The phase error is $\phi_n-t_n$. Unwrapping chooses a continuous branch by removing $2\pi$ jumps; it requires sufficiently resolved phase advances. At this resolution, that assumption holds. A point close to the unit circle can still be at the wrong angle, so energy agreement alone is incomplete.

Tools for inspecting a trajectory¶

These diagnostics consume a trajectory; they do not construct RK stages. The next cell demonstrates them using a complete Euler calculation. The state has two dimensionless components; a physical state with mixed units would need specified component scales before taking a combined norm.

def euler_sho(steps_per_period=100, periods=20):
    """A complete Euler example; returns times and two-component states."""
    if any(isinstance(value, (bool, np.bool_)) or not isinstance(value, (int, np.integer))
           or value <= 0 for value in (steps_per_period, periods)):
        raise ValueError("Use positive integer step and period counts")
    n = steps_per_period*periods
    h = 2*np.pi/steps_per_period
    times = h*np.arange(n+1)
    states = np.empty((n+1, 2))
    states[0] = [1.0, 0.0]
    for i in range(n):
        q, p = states[i]  # Both derivatives use the old state.
        states[i+1] = states[i] + h*np.array([p, -q])
    return times, states


def sho_diagnostics(times, states):
    """Measure energy, phase, and Euclidean state error for omega=1."""
    times = np.asarray(times, dtype=float)
    states = np.asarray(states, dtype=float)
    if times.ndim != 1 or states.shape != (len(times), 2):
        raise ValueError("Expected times (n,) and states (n, 2)")
    if not np.all(np.isfinite(times)) or not np.all(np.isfinite(states)):
        raise ValueError("Diagnostics require finite data")
    energy = np.sum(states**2, axis=1)/2
    phase = np.unwrap(np.arctan2(-states[:, 1], states[:, 0]))
    exact = np.column_stack((np.cos(times), -np.sin(times)))
    state_error = np.linalg.norm(states-exact, axis=1)
    return energy, phase-times, state_error


euler_times, euler_states = euler_sho()
euler_energy, euler_phase_error, euler_state_error = sho_diagnostics(euler_times, euler_states)
print("Euler maximum state error:", euler_state_error.max())
print("Euler final relative energy change:", euler_energy[-1]/euler_energy[0]-1)
print("Euler final phase error (rad):", euler_phase_error[-1])

Euler maximum state error: 50.43601922172834
Euler final relative energy change: 2643.2399896798
Euler final phase error (rad): -0.16497620494871512

Apply the definitions to another trajectory only after checking its equations, initial condition, time units, and state ordering. The helper’s exact reference is specific to this SHO setup. A plausible output from a helper used on the wrong problem is not a validation.

Each trajectory has 2000 steps. A force evaluation and a full right-hand-side call need not have equal cost in an application. These counts are not timings. To compare efficiency, select the observable, duration, and error requirement first, then compare work needed to meet that same requirement. Formal order alone does not determine a universal crossover accuracy.

Why high order does not imply energy conservation¶

For the oscillator, put $z=ih$. The explicit RK stability function determines the energy factor per step. Euler gives $|1+ih|^2=1+h^2$; explicit midpoint gives $1+h^4/4$; classical RK4 gives

Thus RK4 damps this oscillator for $0<h<\sqrt8$, while Euler and RK2 increase its energy. Neither the direction nor the size of a Kepler-orbit energy drift follows from this linear calculation. RK4 is globally fourth order and has an $O(h^5)$ one-step defect; those orders say how state error scales at fixed duration, not that a conserved quantity remains exact.

The explicit RK methods discussed here are not symplectic for general Hamiltonian systems. Some implicit RK methods are symplectic, so this is not a defect of every Runge–Kutta method. A measured energy drift, the Jacobian’s volume change, and the full symplectic condition are related diagnostics, but one does not substitute for the others.

Bridge to Part 3: The Need for Geometric Integration¶

The Runge–Kutta family improves finite-time state accuracy by evaluating the derivative at carefully chosen stages. For a short trajectory calculation, adaptive RK can be an excellent choice. For long autonomous Hamiltonian integrations, geometric preservation supplies a different source of control.

In Part 3, we will derive leapfrog and ask exactly what its geometry preserves. Second-order leapfrog does not automatically beat fourth-order RK4: the answer depends on the observable, duration, resolution, and cost. Energy can remain well behaved while orbital phase becomes inaccurate. The comparison experiment keeps these diagnostics separate so you can justify a method choice.

Next: Part 3 - Symplectic Integration