Ordinary Differential Equations

General framework

Let $\boldsymbol{ f} \in C(\mathbb{ R} \times \mathbb{ R}^{d}, \mathbb{ R}^{d})$ be continuous in the first variable and Lipschitz continuous in the second variable. For all $\boldsymbol{ p} \in \mathbb{ R}^{d}$, the Picard–Lindelöf theorem ensures the existence of an open interval $I \subset \mathbb{ R}$ containing $0$ and a unique $\boldsymbol{ u}: I \to \mathbb{ R}^{d}$ such that

Besides, $\boldsymbol{ u}$ depends continuously on $\boldsymbol{ u}$; in other words, this initial value problem is well-posed.

We are willing to approximate the flow of $\boldsymbol{ f}$ in the interval $[0, T]$ for some $T \in I \cap \mathbb{ R}_{+}$. For simplicity, we construct an approximate flow using constant time steps. Let $n \geq 1$, define $h = T / n$ and $t_{k} = k h$ for $k \in \{0, \ldots, n\}$. We consider a general marching procedure described by the following framework:

• A starting map $\mathsf{ S}_{h}: \mathbb{ R}^{d} \to \mathbb{ S}$ embeds the initial condition into a state space $\mathbb{ S}$,

• A marching map $\mathsf{ M}_{h}: \mathbb{ S} \to \mathbb{ S}$ evolves the state for one time step,

• A finishing map $\mathsf{ F}_{h}: \mathbb{ S} \to \mathbb{ R}^{d}$ converts the state space into the solution space.

We define the recurrences $(\boldsymbol{ s}_{k})_{0 \leq k \leq n}$ and $(\boldsymbol{ u}_{k})_{0 \leq k \leq n}$ by

In other words, $\boldsymbol{ s}_{k} = \mathsf{ M}_{h}^{k} \circ \mathsf{ S}_{h}(\boldsymbol{ u}_{0})$ and $\boldsymbol{ u}_{k} = \mathsf{ F}_{h} \circ \mathsf{ M}_{h}^{k} \circ \mathsf{ S}_{h}(\boldsymbol{ u}_{0})$. We want to design the scheme $(\mathsf{ S}_{h}, \mathsf{ M}_{h}, \mathsf{ F}_{h})$ so that $\boldsymbol{ u}_{k}$ approximates $\boldsymbol{ u}(t_{k})$.

Let $\mathsf{ E}_{h}: \mathbb{ R}^{d} \to \mathbb{ R}^{d}$ denote the exact flow of $\boldsymbol{ f}$ for a time step $h$, that is $\mathsf{ E}(\boldsymbol{ u}(t)) = \boldsymbol{ u}(t + h)$. A first idea is to compare the approximate and continuous solutions at time $T$. The corresponding error is called the global truncation error and is defined as

We are only in the dependence of $\mathsf{ GTE}$ on $h$ and the derivatives of $\boldsymbol{ f}$; not on $\boldsymbol{ p}$. It is easier to look at the error after a single step. This is called the local truncation error,

We can now state the consistency and convergence of a scheme.

We now move to the linear stability of a scheme. It is defined in terms of the behaviour of the scheme when applied to the test problem $\boldsymbol{ f}(t, \boldsymbol{ u}) = \lambda \boldsymbol{ u}$ for $\lambda \in \mathbb{ C}$. The exact solution is $\boldsymbol{ u}(t) = \exp(\lambda t) \boldsymbol{ p}$; in particular $\boldsymbol{ u}$ vanishes at infinity when $\Re \lambda < 0$. We are interested in schemes that mimic this property numerically.

We now move to some nonlinear stability properties. We are interested in dissipative problems where

In that case, if $\boldsymbol{ u}_{0}, \boldsymbol{ v}_{0} \in \mathbb{ R}^{d}$ are two initial conditions, then the corresponding solutions $\boldsymbol{ u}, \boldsymbol{ v}$ satisfy

so $\boldsymbol{ u}$ and $\boldsymbol{ v}$ merge at infinity. It is desirable to reproduce this behaviour numerically.

Similar stability and order notions can be defined in the state space or in some internal process of the marching procedure.

In the following, we briefly introduce the $\theta$-method, Runge-Kutta schemes and linear multi-step. We then provide a detailed analysis of general linear methods, which are a generalisation of theses schemes.

$\theta$-method

The $\theta$-method is a simple family of schemes that includes the explicit and implicit Euler schemes, as well as the trapezoidal rule as special cases. We simply consider $\mathbb{ S} = \mathbb{ R}^{d}$, $\mathsf{ S} = \mathsf{ F} = \operatorname{Id}_{\mathbb{ R}^{d}}$ and $\mathsf{ M}(\boldsymbol{ s}_{k}) = \boldsymbol{ s}_{k+1}$ is defined implicitly through

where $\theta \in \mathbb{ R}$ characterises the scheme. The choice $\theta = 0$ is the only explicit member of the family and is known as explicit Euler. The choice $\theta = 1$ corresponds to implicit Euler, and $\theta = 1/2$ is the trapezoidal rule.

The third-order Taylor expansion of the exact and approximate solutions are

Observe that regardless of $\theta$, the expansions $\boldsymbol{ u}(t_{k+1})$ and $\boldsymbol{ u}_{k+1}$ coincide up to $O(h)$, so the $\theta$-method is at least of order $1$. They coincide up to $O(h^{2})$ if and only if $\theta = 1/2$.

In the case of a linear first-order ODE, the definition of $\boldsymbol{ u}_{k+1}$ is

so we infer

Observe that

We distinguish three classes:

• If $\theta < 1/2$, the stability region is the circle of radius $1 / (1 - 2 \theta)$ centred at $(-1/(1 - 2 \theta), 0)$. This scheme is not $A$-stable.

• If $\theta = 1/2$, the stability region is the whole left complex plane. The scheme is $A$-stable.

• If $\theta > 1/2$, the stability region includes the whole complex plane except the circle of radius $1 / (2 \theta - 1)$ and centred at $(1 /(2 \theta - 1), 0)$. In particular, this scheme is $A$-stable.

The scheme is only $L$-stable when $\theta = 1$.

Runge-Kutta

The Runge-Kutta method is a multi-stage method consisting in enforcing a zero residual at intermediate times between $t_{k}$ and $t_{k+1}$, and forming $u_{k+1}$ from linear combinations of these intermediate solutions. Here again, $\mathbb{ S} = \mathbb{ R}^{d}$ and $\mathsf{ S} = \mathsf{ F} = \operatorname{Id}_{\mathbb{ R}^{d}}$. A RK scheme depends on $s \geq 1$, a matrix $\boldsymbol{ A} \in \mathbb{ R}^{s \times s}$ and vectors $\boldsymbol{ b}, \boldsymbol{ c} \in \mathbb{ R}^{s}$. They are conveniently arranged in a Butcher tableau

The marching procedure is defined as

where $\boldsymbol{ w}_{1}, \ldots, \boldsymbol{ w}_{s} \in \mathbb{ R}^{d}$ are such that

Linear stability

Applied to a linear ODE with constant step size $h$, letting $\boldsymbol{ w} = (\boldsymbol{ w}_{1}, \ldots, \boldsymbol{ w}_{s}) \in \mathbb{ R}^{s}$, this recursion becomes

We infer $\boldsymbol{ w} = \lambda (\boldsymbol{ I} - z \boldsymbol{ A})^{-1} \boldsymbol{ 1}_{s} \boldsymbol{ u}_{k}$ and $\boldsymbol{ u}_{k+1} = (1 + z \boldsymbol{ b}^{\ast} (\boldsymbol{ I} - z \boldsymbol{ A})^{-1} \boldsymbol{ 1}) \boldsymbol{ u}_{k}$, that is

Zero-stability, pre-consistency and consistency

All RK schemes are zero-stable and pre-consistent. The consistency condition is $\boldsymbol{ b}^{\ast} \boldsymbol{ 1} = 1$. One usually assumes stage consistency, that is $\boldsymbol{ A} \boldsymbol{ 1} = \boldsymbol{ c}$.

Order conditions

Higher order conditions can be studied using B-series: introduce $\boldsymbol{ \xi} \in \mathsf{ B}^{s}$ for each stage ($\boldsymbol{ \xi} \sim \boldsymbol{ u}(t + \boldsymbol{ c} h)$). The order conditions read

where $q$ is the stage-order and $p$ is the order.

Linear multi-step

Linear multi-step methods store previous iterates and predict the next. As notable examples, this class of methods includes Adams-Moulton, Adams-Bashford and Backward Differentiation Formula (BDF) methods.

In this setting, $\mathbb{ S} = (\mathbb{ R}^{d})^{r}$ for some $r \geq 1$, $\mathsf{ S}$ can be a Runge-Kutta method and represents $\boldsymbol{ u}(k h)$ for $k \in \{0, \ldots, r - 1\}$. In that case $\mathsf{ F}$ simply extracts the first component of the state, that is $\mathsf{ F}(\boldsymbol{ s}_{1}, \ldots, \boldsymbol{ s}_{r}) = \boldsymbol{ s}_{1}$. The marching procedure is defined as

where $\boldsymbol{ s}_{0} \in \mathbb{ R}^{d}$ solves

Here $\boldsymbol{ \alpha} \in \mathbb{ R}^{r}$ and $\boldsymbol{ \beta} \in \mathbb{ R}^{r+1}$ are the parameters of the method.

Linear stability

The stability matrix is the companion matrix of the polynomial $\pi(w, z) = \rho(w) - z \sigma(w)$, so the stability condition at a given $z \in \mathbb{ C}$ is that the roots of $\pi(w, z)$ are in the unit circle.

Zero-stability, pre-consistency and consistency

It is common to introduce the polynomials

The zero-stability requires the roots of $\rho$ to lie in the unit circle. The pre-consistency condition is then $\sum_{i = 1}^{r} \boldsymbol{ \alpha}_{i} = 1$, that is $\rho(1) = 0$. The consistency condition is $\sum_{i = 0}^{r} \boldsymbol{ \beta}_{i} = r - \sum_{i = 1}^{r} \boldsymbol{ \alpha}_{i} (r - i)$, that is $\rho'(1) - \sigma(1) = 0$.

Order conditions

The order conditions are $\pi(\exp(z), z) = \rho(\exp(z)) - z \sigma(\exp(z)) = O(z^{p+1})$. Alternatively, using the formalism of B-series, the order conditions are written

General linear methods

Runge-Kutta and linear multi-step methods can be examined under the broad class of general linear methods (GLM). In this class of methods, we allow for $r \geq 1$ reused steps and $s \geq 1$ internal stages. Runge-Kutta methods are recovered for $r = 1$, while linear multi-step methods correspond to $s = 1$. GLM were introduced to offer more flexibility than RK and LM methods, provide a better understanding of order conditions, and break order barriers.

The state space is thus $\mathbb{ S} = (\mathbb{ R}^{d})^{r}$. The starting and finishing procedures depend on the particular scheme; we describe them below. The marching procedure has the form

where $\boldsymbol{ w}_{1}, \ldots, \boldsymbol{ w}_{s} \in \mathbb{ R}^{d}$ are defined implicitly via

These schemes are characterised by the matrices $\boldsymbol{ A} \in \mathbb{ R}^{s \times s}$, $\boldsymbol{ B} \in \mathbb{ R}^{r \times s}$, $\boldsymbol{ U} \in \mathbb{ R}^{s \times r}$ and $\boldsymbol{ V} \in \mathbb{ R}^{r \times r}$, as well as assumptions on the interpretation of the state variables, the shape of the stage variables and the finishing procedure. One usually assume that the starting procedure is such that the state at time $t$ represents approximations of $\boldsymbol{ u}(t)$ involving higher-order derivatives, that is

where $\boldsymbol{ \alpha}_{0}, \ldots, \boldsymbol{ \alpha}_{p} \in \mathbb{ R}^{r}$. We do not focus on the design of the starting procedure for now. The internal stages are assumed to represent

for some $\boldsymbol{ c} \in \mathbb{ R}^{s}$. Finally, we assume that the finishing procedure is just a linear combination of the state vectors, that is

for some $\boldsymbol{ \beta} \in \mathbb{ R}^{r}$.

Linear stability

Applied to the linear ODE $\boldsymbol{ f}(t, \boldsymbol{ u}) = \lambda \boldsymbol{ u}$, the recursions defining $\boldsymbol{ w}$ and $\boldsymbol{ s}_{k}$ read

We readily infer $\boldsymbol{ s}_{k+1} = (\boldsymbol{ V} + z \boldsymbol{ B} (\boldsymbol{ I} - z \boldsymbol{ A})^{-1} \boldsymbol{ U}) \boldsymbol{ s}_{k}$, and thus

Strictly speaking, the linear stability analysis requires the knowledge of $\mathsf{ S}_{h}$ and $\mathsf{ F}_{h}$, but one usually directly relates it to the stability in the state space, and introduces

Zero-stability, pre-consistency and consistency

Applying a GLM to $\boldsymbol{ f} = \boldsymbol{ 0}$, we obtain

In particular for $k = 0$, one must have $\mathsf{ F}_{h} \circ \mathsf{ S}_{h} = \operatorname{Id}_{\mathbb{ R}^{d}}$, so $\mathsf{ S}_{h}$ must be left-invertible and $\mathsf{ F}_{h}$ has to be a left-inverse of $\mathsf{ S}_{h}$. Besides, $\boldsymbol{ V}$ must be power-bounded for $\boldsymbol{ u}_{k}$ to remain bounded.

The pre-consistency and consistency conditions are $\boldsymbol{ \beta}^{\ast} \boldsymbol{ V} \boldsymbol{ \alpha}_{0} = 1$ and $\boldsymbol{ \beta}^{\ast} (\boldsymbol{ B} \boldsymbol{ 1} + \boldsymbol{ V} \boldsymbol{ \alpha}_{1}) = 1$. In fact, these are projections on the output space of the pre-consistency and consistency conditions in the state space, $\boldsymbol{ V} \boldsymbol{ \alpha}_{0} = \boldsymbol{ \alpha}_{0}$ and $\boldsymbol{ \alpha}_{0} + \boldsymbol{ \alpha}_{1} = \boldsymbol{ B} \boldsymbol{ 1} + \boldsymbol{ V} \boldsymbol{ \alpha}_{1}$, respectively. It is often of interest to look for the pre-consistency and consistency at the stage level. These conditions are $\boldsymbol{ U} \boldsymbol{ \alpha}_{0} = \boldsymbol{ 1}$ and $\boldsymbol{ c} = \boldsymbol{ A} \boldsymbol{ 1} + \boldsymbol{ U} \boldsymbol{ \alpha}_{1}$, respectively.

Order conditions

The analysis of GLM is simplified through the formalism of B-series. We introduce $\boldsymbol{ \xi} \in \mathsf{ B}^{s}$ to represent the stage vectors (corresponding to $\boldsymbol{ u}(t + \boldsymbol{ c} h)$) and $\boldsymbol{ \eta} \in \mathsf{ B}^{r}$ for the state vectors (corresponding to $\boldsymbol{ \alpha}_{k} \boldsymbol{ u}^{(k)}(t)$). They are related by

Here $q$ is the stage order, $p$ is the state order and $o$ is the order of the method.