Some recent work on physical simulation sent me back to solving ordinary differential equations—and how I had little intuition for how non-trivial numerical solvers are born. After covering some basic properties, this post digs into the quadrature viewpoint for synthesising ‘any’ solvers from the Runge-Kutta family. We illustrate the method by rederiving the celebrated RK4 solver from first principles.

Setting

Motivation

Let $\mathcal{I}$ be some interval of $\mathbb{R}$ and $f:\mathbb{R}^d\times\mathcal{I}\mapsto\mathbb{R}^d$ describe a family of vector fields indexed by some variable $t\in\mathcal{I}$. Let $y:\mathbb{R}\mapsto\mathbb{R}^d$ describe the motion of a ‘particle’ which, at each instant $t\in\mathcal{I}$, follows the vector field $f(\cdot, t)$. It is characterised by a non-autonomous ordinary differential equation (ODE); for any $t\in\mathcal{I}$: $$ \frac{dy}{dt}(t) = f\big(y(t), t\big)\;. $$ When given an initial condition, finding such $y$ is known as a Cauchy problem. Only in some notorious cases (mostly, linear) we know of a closed-form solution for this system. In other cases, we can only be satisfied with numerical approximations—the topic of this post. But first, it will be in good taste to understand the conditions required for a Cauchy problem to be well-posed (i.e. it admits a solution, and only one).

Assumptions

We will work under a simplified setting. Namely, we will assume that the ODE is autonomous: the vector field is constant wrt. the time variable. We will also take $\mathcal{I}=[0, 1]$. (Those assumptions are essentially to limit clutter; all results presented below all easily extend to the non-autonomous case and to arbitrary compact $\mathcal{I}$.) Therefore, throughout this post we are interested in the solution of the Cauchy problem:

The Cauchy-Lipschitz theorem

The Cauchy-Lipschitz theorem guarantees that under some rather mild smoothness assumption, a unique solution exists. Here, we will focus on the ‘strong’ version which requires global Lipschitz continuity—some weaker versions exists, with slightly degraded claims.

“It exists and is unique” is a solid tell that some fixed-point argument is at work here. That’s the path taken in the proof below, which studies the following operator before invoking the Banach fixed-point theorem. $$ \tag{1} \mathcal{T} : y(t) \mapsto y_0 + \int_{0}^t f(y(\tau))d\tau\;. $$

Below, assume this assumption holds, so it makes sense to try and approximate the solution of $(\star)$.

Numerical solvers

Let $y$ be the unique solution of $(\star)$. A numerical solver builds an approximation of $y$ piecemeal, over a regular, discrete mesh $t_1, \ldots, t_N\in[0,1]$ where $t_{i+1} = t_i + h$ and $h>0$. It produces some $y_1, \ldots y_N\in\mathbb{R}^d$ such that $y_i \approx y(t_i)$ for all $i$. Here, we focus on explicit single-step solvers, for which the update writes: $$ \tag{2} y_{i+1} = y_i + h\phi(y_i, h)\;. $$

The goal of this section is to motivate some basic properties we ask of such numerical solvers. Two quantities are especially relevant to our study; the scheme’s global error $E_i := y_i - y(t_i)$ and its truncation error: $$ T_i(h) := \frac{y(t_{i+1}) - y(t_i)}{h} - \phi(y(t_i), h)\;. $$

While the global error measures the total divergence between approximated and true solution after $i$ steps, the truncation error essentially captures the one-step error. Clearly, we want both to be small.

Stability

The scheme detailed in (2) is called stable if it exists $C$ and $H$ such that for all $h<H$ and $x,y\in\mathbb{R}^d$:

$$ \| \phi(x, h) - \phi(y, h) \| \leq C \|x-y\|\;. $$

Stability allows to relate the scheme’s global (accumulated) error to the truncation error. Indeed, if stable: $$ \tag{3} \| E_i \| \leq \blacksquare \sup_{j\leq i} \| T_j(h) \|\; , $$

with $\blacksquare$ some (bad) constant. Given our assumption that $f$ is Lipschitz continuous, stability is hardly a tough condition to meet. Below we assume (2) to be stable so we focus our analysis only on the truncation error.

Consistency and order

A scheme is called consistent if $\lim_{h\to 0} \phi(x, h) = f(x)$ for any $x\in\mathbb{R}^d$. Thanks to (3) we get that a consistent scheme is also a convergent one, as its global error goes to 0 as $h$ does. If consistency captures that the truncation error goes to 0, it doesn’t say how fast—and that would be the role of the scheme’s order. The solver (2) is said to be of order $p\in\mathbb{N}$ if $T_i(h) = O(h^p)$ for all $i\in\{1,\ldots, N\}$. A higher order is desirable, but inevitably comes with extra computations.

Examples

Perhaps the simplest of all solvers is known as the explicit Euler, which writes $ y_{i+1} = y_i + hf(y_i)\;. $ It is trivially stable and consistent, and one easily shows that it is of order 1. Its sibling, the improved Euler: $$ y_{i+1} = y_i + \frac{h}{2}\big(f(y_i)+ f(y_i + hf(y_i))\big)\;, $$ is also stable, consistent, and it enjoys a better $O(h^2)$ rate (the proof being a simple Taylor expansion). Both those solvers are member of the Runge-Kutta family, which the rest of this post studies. Especially, we will be interested in the rationale and mechanisms used to generate them.

Synthesising solvers

The fundamental theorem of calculus essentially translates differential equations into integral ones; $$ \tag{4} y_{i+1} = y_i + \int_{t_i}^{t_{i+1}}f(y(\tau))d\tau\;. $$ Under this perspective, creating a new solver essentially boils down to making some integral approximation choices. For instance, the explicit Euler essentially considers $f$ to be the constant $f(y(t_i))$. This essentially opens up an entire inventory of methods tapping into the many existing quadrature rules. For instance, using the trapezoidal rule yields the (implicit) trapezoidal method: $$ \tag{5} y_{i+1}= y_i + \frac{h}{2}(f(y_{i}) + f(y_{i+1}))\;. $$ Since $y_{i+1}$ stands on both sides of the equal sign, (5) is an implicit solver. We fall back to an explicit one by making a further approximation on the right-hand side, using the explicit Euler quadrature $y_{i+1} = y_i + hf(y_i)$—and recovering the improved Euler method. The following section essentially shows how to generalise and mechanise this approach to higher orders.

Nested quadratures

Let’s rewrite (4) with a simple change of variable: $$ y_{i+1} = y_i + h\int_{z=0}^1 f(y(t_i + hz))dz\;. $$ A quadrature rule for approximating the above integral picks $M\in\mathbb{N}$ anchor points ${z_1, \ldots, z_M}\in[0,1]$ along with some weights $b_1, \ldots, b_M\in\mathbb{R}$ to form a weighted sum: $$ \tag{6} y_{i+1} \approx y_i + h\sum_{j=1}^M b_j f(\tilde{y}_{ij}) \text{ where } \tilde{y}_{ij} := y(t_i + hz_j)\;. $$ Observe that each $\tilde{y}_{ij}$ is itself unknown. They can too be approximated via quadrature rules—an inner integral approximation, which will re-use the anchors points we just defined but aggregate them with different coefficients. Concretely, for $j\in\{1, \ldots, M\}$ we introduce $\{a_{j1}, \ldots, a_{j(j-1)}\}\in\mathbb{R}$ and write: $$ \tag{7} \tilde{y}_{ij} = y_i + h \int_{0}^{z_{j}} f(y(\tau))d\tau \approx y_i + h\sum_{k< j} a_{jk} f(\tilde y_{ik})\;. $$

Combining (6) and (7) yields the entire update rule. At this point we introduced a bunch of free parameters; the anchors $\{z_j\}_j$, their outer weights $\{b_j\}_j$ and inner weights $\{a_{jk}\}_{jk}$. Some constraints can reduce the number of degrees of freedom; for instance, it is typical to require for our quadrature rules to integrate exactly any constant function. This yields the following constraints: $$ \tag{8} \sum_{j=1}^M b_j = 1 \text{ and } \sum_{k<j} a_{jk} = z_j \text{ for all } j\in\{1,\ldots, M\}. $$

Butcher tables

A solver in those likes can be compactly described by a so-called Butcher table: $$ \begin{array}{c|cccc} z_1 & 0 & 0 & \cdots & 0 \\ z_2 & a_{21} & 0 & \cdots & 0 \\ \vdots & \vdots & \ddots & \ddots & \vdots \\ z_M & a_{M1} & a_{M2} & \cdots & 0 \\ \hline & b_1 & b_2 & \cdots & b_M \end{array}\;. $$ The diagonal and upper-diagonal blocks are filled by 0s, as per our requirement for the solver to be explicit. The constraints (8) impose that each row of the inner block sums to its anchor, and that the bottom row sums to one. For instance, the explicit Euler and the improved Euler are respectively described by: $$ \begin{array}{c|c} 0 & 0 \\ \hline & 1 \end{array} \qquad\text{and}\qquad \begin{array}{c|cc} 0 & 0 & 0 \\ 1 & 1 & 0 \\ \hline & 1/2 & 1/2 \end{array}\;. $$

Building RK4

This section uses the quadrature approach to derive the Runge-Kutta 4th-order solver from the ground-up. The starting point is Simpson’s quadrature rule, which writes: $$ \int_0^1 f(x)dx \approx f(0)/6 + 2f(1/2)/3 + f(1)/6\;. $$ To reach a higher order, the RK4 method duplicates the midpoint. Concretely, we will be using:

We can, for completeness, spell out the relationships defined in (7) along with the constraints from (8): $$ \left\{ \begin{aligned} \tilde y_{i1} &= y_i\;,\\ \tilde y_{i2} &= y_i + ha_{21}f(\tilde y_{i1})\;, &\text{ s.t } a_{21}=1/2\\ \tilde y_{i3} &= y_i + ha_{31}f(\tilde y_{i1}) + ha_{32}f(\tilde y_{i2})\;, &\text{ s.t } a_{31}+a_{32}=1/2 \\ \tilde y_{i4} &= y_i + ha_{41}f(\tilde y_{i1}) + ha_{42}f(\tilde y_{i2}) + ha_{43}f(\tilde y_{i3})\;, &\text{ s.t } a_{41}+a_{42}+a_{43}=1\;.\\ \end{aligned} \right. $$

It stands out that we still have $3$ unresolved degrees of freedom. To settle them, we resort to Taylor expansion matching: we choose the remaining degrees of freedom so the scheme can be of the highest possible order. Concretely, we will try to set the $\{a_{jk}\}_{jk}$ so that they respect the aforementioned constraints and allow to write $T_i(h) = O(h^p)$ with the highest possible value of $p$. The computation is rather tedious as it involves somewhat long Taylor expansions; we omit it here to directly reveal the result:

$$ \begin{aligned} \frac{y(t_{i+1}) - y_{i+1}}{h} = & h^2 \left[ \frac{1}{12} \left( 2 - 2a_{32} - a_{42} - a_{43} \right) \right] f’^2 f \\ & + h^3 \left[ \frac{1}{48} \left( 8 - 6a_{32} - 5(a_{42} + a_{43}) \right) \right] f’f’’ f^2 \\ & + h^3 \left[ \frac{1}{24} \left( 1 - 2 a_{43} a_{32} \right) \right] f’^3 f \\ & + O(h^4)\;, \end{aligned} $$ where $f$ and its derivatives are all evaluated at $y_i$. Now we can retrieve have our remaining constraints, that ensure that all $O(h^2)$ and $O(h^3)$ term are zeroes. That’s three of them, one for each unresolved degree of freedom. Combining them with previous ones, we get the system: $$ \left\{ \begin{aligned} &a_{21} = 1/2\;, \\ &a_{31} + a_{32} = 1/2\;, \\ &a_{41} + a_{42} + a_{43} = 1\;, \\ &2a_{32} + a_{42} + a_{43} = 2\;, \\ &6a_{32} + 5(a_{42} + a_{43}) = 8\;, \\ &2a_{43}\,a_{32} = 1\;. \end{aligned} \right. $$ This system admits a unique solution: $a_{21} = 1/2$, $a_{31} = 0$, $a_{32} = 1/2$, $a_{41} = a_{42} = 0$ and $a_{43} = 1$. From this we recover the Butcher table of the celebrated RK4 solver:

or in a form that stands close to the actual implementation:

As a by-product, we proved that $T_i(h)=O(h^4)$ for any $i\in\{1,\ldots, M\}$. In other words, RK4 is of order 4.

Implicit schemes

This post so far focused essentially on explicit single-step solvers. While powerful, many (most) advanced applications require more sophistication. There is a very large zoology out there; symplectic solvers, multistep solvers, adaptive-step schemes, etc. Covering them is out of scope here; instead, we hereinafter briefly introduce implicit solvers. To understand their relevance, we first need to discuss stiffness.

Stiffness

Consider the one-dimensional Cauchy problem: $$ \frac{dy}{dt}(t) = - \beta y(t) \text{ and } y(0)=1\;, \tag{9} $$ with $\beta\gg 1$. Its unique solution is the quickly decaying $y(t) = \exp(-\beta t)$. It is fairly easy to show that the sequence of estimates produced by the explicit Euler solver is, for any $i\in\mathbb{N}$: $$ y_i = (1 - h\beta)^{i}\;. $$ For any $h>2/\beta$ this sequence quickly diverges, even though the explicit Euler scheme is stable and consistent—and hence convergent (observe that this does not contradict earlier statements). Essentially, the requirement for the scheme to behave is that the discretisation step is small—at least compared to $\beta$.

This becomes an issue for multidimensional systems, where some quickly decaying components impose an absurdly small global step-size for approximating other components. For example, the system described by: $$ \frac{dy}{dt}(t) = \begin{pmatrix} -\beta & 0 \\ 0 & -\beta^{-1}\end{pmatrix}y(t)\;, $$ is stiff. Solving it (naively) via an explicit numerical solver requires $h\lessapprox 1/\beta$ because of its first component; but this is risible (as in, unnecessarily small) step-size to approximate the second component which evolves slowly. Explicit solvers ultimately all suffer from this defect (formally, one says that they are not A-stable).

An implicit scheme

Single-step implicit solvers write as: $$ y_{i+1} = y_i + h\phi(y_i, y_{i+1}, h)\;. $$ In general, they require their own solver (e.g., a root finder) to materialise $y_{i+1}$; hence, they will typically ask for more computations. In return, they generally pair much better with stiff equations. There is not (to my knowledge) any general demonstration of this claim (perhaps because there is no universally agreed upon definition of stiffness) and this is mostly folk knowledge. It is educative to look at the behaviour of the implicit Euler scheme $y_{i+1} = y_i + hf(y_{i+1})$ on the example from (9); one can easily show that then: $$ y_{i} = (1+h\beta)^{-i}\;, $$ which provides reasonable a estimate regardless of the value of $h/\beta$ (it is A-stable). Implicit solvers also have a Butcher table representation, with a non-zero entries on and/or above the diagonal. For instance, the implicit Euler being a diagonally implicit scheme has the following representation: $$ \begin{array}{c|c} 1 & 1 \\ \hline & 1 \end{array}\;. $$