Chapter 23 — Key Takeaways (Ordinary Differential Equations)

A one-page reference for integrating ODEs in Fortran: the methods, their orders, the code patterns, and the traps.

The initial-value problem

$$\frac{dy}{dt} = f(t, y), \qquad y(t_0) = y_0$$

• $f(t, y)$ is the right-hand-side function (RHS): the rate of change as a function of state and time.
• Solve numerically by marching: produce $y_1, y_2, \dots$ at $t_1, t_2, \dots$ from $t_0$.
• Reduce any higher-order ODE to a first-order system: for $y'' = g(t, y, y')$ set $y_1 = y$, $y_2 = y'$, giving $y_1' = y_2,\ y_2' = g(t, y_1, y_2)$.

The methods at a glance

RK4 — memorize these exactly

$$k_1 = f(t_n, y_n), \quad k_2 = f\!\left(t_n{+}\tfrac h2, y_n{+}\tfrac h2 k_1\right), \quad k_3 = f\!\left(t_n{+}\tfrac h2, y_n{+}\tfrac h2 k_2\right), \quad k_4 = f(t_n{+}h, y_n{+}h k_3)$$ $$y_{n+1} = y_n + \frac{h}{6}\left(k_1 + 2k_2 + 2k_3 + k_4\right)$$

• Stages: left slope, two midpoint slopes, right slope. Weights $\tfrac16,\tfrac13,\tfrac13,\tfrac16$ sum to 1 (they are Simpson's weights).
• The classic bug: $k_3$ accidentally using $k_1$ (or a weight of $\tfrac13$ where $\tfrac16$ belongs). It still runs — and is silently second-order. Test the order to catch it.

Order of accuracy — the numbers to know

Reaching accuracy $\varepsilon$: Euler needs $\sim \varepsilon^{-1}$ steps; RK4 needs $\sim \varepsilon^{-1/4}$. For $\varepsilon = 10^{-6}$, that is ~$10^6$ steps vs ~$30$ — RK4 wins overwhelmingly despite 4 evals/step.

Adaptive stepping (§23.3)

• Estimate error by step doubling: run size $h$ vs two of $h/2$; error of the finer result $\approx (y_{\text{half}} - y_{\text{full}})/(2^{p} - 1)$. For RK4, divide by 15.
• Controller: $h_{\text{new}} = h\,S\,(\text{tol}/\text{err})^{1/(p+1)}$, safety $S \approx 0.9$. err > tol ⇒ reject and shrink; err < tol ⇒ accept and grow.
• Production libraries use embedded pairs (RKF45, Dormand-Prince dopri5 = ode45/SciPy default): two orders from one set of stages.

Method of lines — the bridge to PDEs (§23.4)

Discretize space, keep time continuous → a PDE becomes a large system of ODEs, one per grid point.

• 1D heat: $\dfrac{du_i}{dt} = \dfrac{\alpha}{\Delta x^2}(u_{i-1} - 2u_i + u_{i+1})$, i.e. $\mathbf{u}' = A\mathbf{u}$.
• Euler on this system = the explicit heat scheme of Chapter 24. RK4 drops in for higher-order time.
• Stability limit = ODE absolute-stability limit of the method: for explicit Euler, $\Delta t \le \Delta x^2/(2\alpha)$ (1D) — this is the CFL condition (Chapter 24).

Stiffness and implicit methods (§23.5)

• Stiff = widely separated time scales; explicit methods forced to tiny stability-limited steps.
• Explicit Euler on $y' = \lambda y$ is stable only if $|1 + h\lambda| \le 1$ ⇒ $h \le 2/|\lambda|$ (real $\lambda < 0$). RK4 buys only ~40% (limit $\approx 2.8/|\lambda|$).
• Backward Euler: $y_{n+1} = y_n/(1 - h\lambda)$ — A-stable, decays for any $h > 0$. Cost: a solve each step (linear ⇒ LAPACK dgesv/dgtsv, Chapter 21; nonlinear ⇒ Newton).

The Fortran pattern — one integrator, any equation

abstract interface
  function rhs_sys(t, y) result(dydt)     ! the RHS contract
    import :: dp
    real(dp), intent(in) :: t, y(:)
    real(dp)             :: dydt(size(y))
  end function rhs_sys
end interface
! ...
function rk4_sys(f, t, y, h) result(y_next)
  procedure(rhs_sys)   :: f               ! RHS passed as a procedure argument (Ch 6)
  real(dp), intent(in) :: t, h, y(:)      ! assumed-shape state
  real(dp)             :: y_next(size(y)), 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)   ! whole-array ops (Ch 5)
end function

• Procedure-argument RHS (procedure(rhs_sys) :: f) ⇒ the solver never mentions the equation.
• intent(in) on $t, y$; assumed-shape y(:) ⇒ same routine for 2 or $10^6$ components.
• Stage arithmetic is whole-array — no inner loops, and the compiler vectorizes it.

Intrinsics / features used

Pitfalls

• Mistyped RK4 coefficient → silent order loss. Verify the order numerically.
• Advancing t before using it in a non-autonomous $f(t, y)$ → wrong trajectory.
• Using explicit RK4 on a stiff problem → tiny steps or NaN; switch to implicit.
• Default real over a long integration → round-off swamps the answer; use real(dp) (Chapter 20).
• Reading dt too large for the CFL limit in a MOL/PDE step → blow-up (Chapter 24).

Compile

$ gfortran -std=f2018 -Wall example-02-rk4.f90 -o rk4 && ./rk4

No external libraries — every example is plain, self-contained Fortran.

Heat-solver project increment

The time loop is an ODE integrator: u = u + dt*heat_rhs(u) is Euler's method on the method-of-lines system $\mathbf{u}' = \mathbf{F}(\mathbf{u})$. Saved as heat-solver/heat_mol.f90; Chapter 24 makes it the 2D explicit core (with the CFL condition), and rk4_sys can replace Euler for higher-order time.