Case Study 23.2: Building a Reusable Adaptive ODE Integrator

"Give me a place to stand and a lever long enough, and I will move the world." — attributed to Archimedes

Executive Summary

In Case Study 23.1 you ported someone's integrator. Here you build one worth reusing: a small Fortran module that integrates any first-order system with automatic step-size control, so a user supplies only an RHS and a tolerance and the module does the rest — growing the step where the solution is smooth, shrinking it where it is sharp. This is the design that turns a one-off script into a tool you reach for on the next problem and the one after.

The engineering content is as important as the numerics. A good integrator module hides its internals behind a clean interface, takes the RHS as a procedure argument so it is model-agnostic, defends itself against a zero error estimate, and returns enough information (the accepted state, the next step, whether the step was accepted) for a driver to build a full adaptive loop. We build it, verify one adaptive step against hand computation, and apply it to the harmonic oscillator, using conserved energy as the diagnostic.

Skills applied: - Designing a model-agnostic API with an abstract interface RHS (§23.1, §23.4, Chapter 6). - Implementing RK4 for systems as whole-array operations (§23.4, Chapter 5). - Step-doubling error estimation and the step-size controller (§23.3). - Defensive numerics: guarding the controller against division by zero (Chapter 13). - Validating with a conserved quantity (§23.4 energy drift, §23.6).

Background

Fixed-step integration wastes work (§23.3): a step small enough for the hardest moment is needlessly small everywhere else. An adaptive integrator estimates each step's local error and steers the step size to keep that error near a tolerance. We use step doubling for the estimate — take the step once at size $h$ and again as two steps of $h/2$, and for a method of order $p$ the difference estimates the finer result's error as $(\mathbf{y}_{\text{half}} - \mathbf{y}_{\text{full}})/(2^{p} - 1)$ — and the standard controller for the next step, $$h_{\text{new}} = h \cdot S \cdot \left(\frac{\text{tol}}{\text{err}}\right)^{1/(p+1)}, \qquad S \approx 0.9.$$ With RK4 as the base method, $p = 4$, so the estimate divides by $2^4 - 1 = 15$ and the controller exponent is $1/5$.

Phase 1 — Design the Interface

Decide the API before writing the body. A reusable integrator must not know what equation it is solving, so the RHS is a procedure argument described by an abstract interface. One adaptive step should report three things: the new state, the recommended next step, and whether this step met tolerance (so the driver can retry a rejected step). That yields a single subroutine:

module ode_adaptive
  implicit none
  integer, parameter :: dp = selected_real_kind(15, 307)

  abstract interface
    function rhs_sys(t, y) result(dydt)
      import :: dp
      real(dp), intent(in) :: t
      real(dp), intent(in) :: y(:)
      real(dp)             :: dydt(size(y))
    end function rhs_sys
  end interface

contains

Phase 2 — The RK4 Core

The base method is the vector RK4 of §23.4, unchanged — a good module reuses its own parts:

  function rk4_sys(f, t, y, h) result(y_next)
    procedure(rhs_sys)   :: f
    real(dp), intent(in) :: t, h, y(:)
    real(dp)             :: y_next(size(y))
    real(dp)             :: k1(size(y)), k2(size(y)), k3(size(y)), k4(size(y))
    k1 = f(t,            y)
    k2 = f(t + 0.5_dp*h, y + 0.5_dp*h*k1)
    k3 = f(t + 0.5_dp*h, y + 0.5_dp*h*k2)
    k4 = f(t + h,        y + h*k3)
    y_next = y + (h/6.0_dp) * (k1 + 2.0_dp*k2 + 2.0_dp*k3 + k4)
  end function rk4_sys

Phase 3 — The Adaptive Step

Now the new part: one controlled step. Take the full step and the two half-steps, form the error estimate as the largest component of their difference over $2^p - 1$, decide acceptance, and propose the next step. Note the max(err, tiny(...)) guard — when the error estimate is zero (an exactly linear stretch) the naive formula would divide by zero, and a robust tool must not crash on its own success.

  subroutine rk4_step_adaptive(f, t, y, h, tol, y_new, h_new, err, accepted)
    procedure(rhs_sys)    :: f
    real(dp), intent(in)  :: t, y(:), h, tol
    real(dp), intent(out) :: y_new(size(y)), h_new, err
    logical,  intent(out) :: accepted
    integer,  parameter   :: p = 4          ! RK4 base order
    real(dp), parameter   :: safety = 0.9_dp
    real(dp) :: y_full(size(y)), y_half(size(y))

    y_full = rk4_sys(f, t, y, h)                       ! one step of size h
    y_half = rk4_sys(f, t,            y,      0.5_dp*h) ! two steps of h/2
    y_half = rk4_sys(f, t + 0.5_dp*h, y_half, 0.5_dp*h)

    err      = maxval(abs(y_half - y_full)) / real(2**p - 1, dp)
    accepted = (err <= tol)
    if (accepted) then
      y_new = y_half                                   ! keep the more accurate result
    else
      y_new = y                                        ! reject: caller retries from y
    end if
    h_new = h * safety * (tol / max(err, tiny(1.0_dp))) ** (1.0_dp/real(p + 1, dp))
  end subroutine rk4_step_adaptive
end module ode_adaptive

Phase 4 — Verify One Adaptive Step by Hand

Test the machinery on the scalar problem $y' = y$, $y(0) = 1$, packaged as a one-component system so it exercises the vector code. With $h = 0.5$ and $\text{tol} = 10^{-4}$: the full RK4 step gives $1.6484375$ (the hand trace from §23.2 with $h = 0.5$), and two half-steps of RK4 give $1.6486995$.

program verify_adaptive
  use ode_adaptive
  implicit none
  real(dp) :: y(1), y_new(1), h_new, err
  logical  :: accepted

  y = [1.0_dp]
  call rk4_step_adaptive(grow, 0.0_dp, y, 0.5_dp, 1.0e-4_dp, y_new, h_new, err, accepted)
  print '(a, f13.7)',  'accepted state y_new    = ', y_new(1)
  print '(a, es14.5)', 'error estimate          = ', err
  print '(a, l3)',     'accepted?               = ', accepted
  print '(a, f13.6)',  'proposed next h         = ', h_new
contains
  function grow(t, y) result(dydt)
    real(dp), intent(in) :: t
    real(dp), intent(in) :: y(:)
    real(dp)             :: dydt(size(y))
    dydt = y
  end function grow
end program verify_adaptive

$ gfortran -std=f2018 -Wall ode_adaptive.f90 verify_adaptive.f90 -o va && ./va
accepted state y_new    =     1.6486995
error estimate          =    1.74646E-05
accepted?               =   T
proposed next h         =      0.637941

By hand: the error estimate is $(1.6486995 - 1.6484375)/15 = 0.0002620/15 \approx 1.746\times10^{-5}$, comfortably under the $10^{-4}$ tolerance, so the step is accepted and the state advances to the more accurate two-half-step value $1.6486995$. Because there is error to spare, the controller grows the step: $h_{\text{new}} = 0.5 \cdot 0.9 \cdot (10^{-4}/1.746\times10^{-5})^{1/5} = 0.637941$. Contrast the chapter's §23.3 example, where a loose Euler step was rejected and shrunk — same controller, opposite verdict, driven entirely by the measured error.

Phase 5 — Apply It, and Diagnose with Energy

Point the finished module at the harmonic oscillator $y_1' = y_2,\ y_2' = -y_1$ from code/example-03-system.f90, wrap rk4_step_adaptive in a do while (t < t_end) loop that retries rejected steps and accepts good ones, and you have a general adaptive solver:

t = 0.0_dp;  y = [1.0_dp, 0.0_dp];  h = 0.5_dp
do while (t < t_end)
  if (t + h > t_end) h = t_end - t                 ! do not overshoot the end
  call rk4_step_adaptive(oscillator, t, y, h, tol, y_new, h_new, err, accepted)
  if (accepted) then
    y = y_new;  t = t + h                          ! commit the step
  end if
  h = min(h_new, 5.0_dp*h)                          ! grow, but clamp runaway growth
end do

Trace the control logic on a rejected step: suppose at some sharp point err comes back at $4\times$ the tolerance. Then accepted is .false., so t and y do not advance — the step is thrown away — and the controller sets $h_{\text{new}} = h \cdot 0.9 \cdot (1/4)^{1/5} \approx 0.68\,h$, shrinking the step for the retry. The loop runs again from the same t with the smaller h, and now err falls under tolerance, the step commits, and h is allowed to grow again. That accept/reject/retry dance, driven entirely by the measured error, is the whole value of adaptivity: the solver spends tiny steps only at the sharp points and long steps through the smooth stretches, with no hand-tuning from you.

The right diagnostic for the run is the conserved energy $E = \tfrac12(y_1^2 + y_2^2)$, which the true oscillator holds at $0.5$ forever. A single RK4 step of $h = 0.5$ already nudged it to $0.499895$ (§23.4); over an adaptive run the tolerance controls how fast that drift accumulates — tighten tol and the energy holds longer, at the cost of more steps. Watching $E$ against wall-clock steps is how you choose a tolerance: the loosest one whose energy drift you can live with over the run you actually need.

The design lesson is what makes this reusable: the module never mentions the oscillator, the test problem, or the heat equation. Any RHS matching the interface — predator-prey from CS-01, a decay chain, the method-of-lines heat system of the Project Checkpoint — drops straight in. That is the power of building the abstraction (an ODE integrator) rather than a script (this integration): one lever, and a place to stand.

A caution for production: step doubling costs three RK4 evaluations per accepted step (one full, two half). The embedded pairs of §23.3 (RKF45, Dormand-Prince) get the same error estimate from a single set of stages and are what libraries ship. Our step-doubling design is the clearest one to build and understand; reach for dopri5 or ODEPACK's lsoda when you need the last drop of efficiency or automatic stiffness handling.

Discussion Questions

1. Why does the adaptive step return accepted as an output rather than deciding internally whether to advance time? What does that buy the calling driver?
2. The controller has a safety factor $S = 0.9$. What goes wrong if you set $S = 1.0$? What if you set it to $0.5$?
3. Step doubling costs three RK4 evaluations per accepted step. Under what circumstances is that overhead worth it compared with a well-chosen fixed step?
4. Why guard the controller with max(err, tiny(1.0_dp))? Construct an input for which the unguarded version divides by zero.

Your Turn: Extensions

• Option A (⭐⭐). Add a relative error control: normalize each component's error by $|\text{tol} \cdot y_i| + \text{atol}$ before taking the max, so large and small components are treated fairly. This is what real solvers do.
• Option B (⭐⭐). Write the full adaptive driver: a do while loop that retries rejected steps, clamps the step growth to a factor of (say) 5 per step, and never overshoots t_end. Integrate the oscillator to $t = 10$ and report the number of accepted versus rejected steps.
• Option C (⭐⭐⭐). Replace step doubling with an embedded RK pair: implement the Bogacki-Shampine 3(2) pair (the method behind SciPy's RK23), which yields a third- and a second-order estimate from the same four stages, and compare its cost per accepted step against step-doubled RK4.

Key Takeaways

• A reusable integrator is defined by its interface: an abstract interface RHS makes it model-agnostic, and returning (y_new, h_new, accepted) lets a driver build any control policy.
• Step doubling gives a local-error estimate for any base method: divide the full-vs-two-halves difference by $2^p - 1$; for RK4, by $15$.
• The step controller $h_{\text{new}} = hS(\text{tol}/\text{err})^{1/(p+1)}$ both grows and shrinks the step from the same formula — the measured error decides.
• Defensive touches (a tiny guard, a growth clamp) separate a toy from a tool.
• Build the abstraction, not the script: the same module integrates predators, oscillators, and the heat equation, changing only the RHS.