2.1 Preliminaries

We are concerned with autonomous systems of differential equations in $\mathbb{R}^{D}$ ,

(2.1) $$\begin{eqnarray}\displaystyle \frac{\text{d}}{\text{d}t}x=f(x);\end{eqnarray}$$

the function (vector field) $f$ is assumed throughout to be defined in the whole of $\mathbb{R}^{D}$ and to be sufficiently smooth. An important role is played by the particular case

(2.2) $$\begin{eqnarray}\displaystyle \frac{\text{d}}{\text{d}t}q=M^{-1}p,\quad \displaystyle \frac{\text{d}}{\text{d}t}p=F(q),\end{eqnarray}$$

where $x=(q,p)\in \mathbb{R}^{D}$ , $D=2d$ , $q\in \mathbb{R}^{d}$ , $p\in \mathbb{R}^{d}$ and $M$ is a constant, invertible matrix, so that $f=(M^{-1}p,F(q))$ . By eliminating $p$ , (2.2) is seen to be equivalent to

$$\begin{eqnarray}M\displaystyle \frac{\text{d}^{2}}{\text{d}t^{2}}q=F(q);\end{eqnarray}$$

this is not the most general autonomous system of second-order differential equations in $\mathbb{R}^{d}$ because the derivative $(\text{d}/\text{d}t)q$ does not appear on the right-hand side. When the forces depend only on the positions, Newton’s second law for a mechanical system gives rise to differential equations of the form (2.2); then $q$ , $(\text{d}/\text{d}t)q$ , $p$ and $F$ are respectively the vectors of coordinates, velocities, momenta and forces, and $M$ is the matrix of masses.

We let $\unicode[STIX]{x1D711}_{t}$ denote the $t$ -flow of the system (2.1) under consideration. By definition, for fixed real $t$ , $\unicode[STIX]{x1D711}_{t}:\mathbb{R}^{D}\rightarrow \mathbb{R}^{D}$ is the map that associates with each $\unicode[STIX]{x1D6FC}\in \mathbb{R}^{D}$ the value at time $t$ of the solution of (2.1) that at the initial time 0 takes the initial value $\unicode[STIX]{x1D6FC}$ .

Example 2.1. As a very simple but important example, we consider the standard harmonic oscillator, the system in $\mathbb{R}^{2}$ of the special form (2.2) given by

(2.3) $$\begin{eqnarray}\displaystyle \frac{\text{d}q}{\text{d}t}=p,\quad \displaystyle \frac{\text{d}p}{\text{d}t}=-q.\end{eqnarray}$$

For future reference, we note that, with matrix notation, the solutions satisfy

(2.4) $$\begin{eqnarray}\left[\begin{array}{@{}c@{}}q(t)\\ p(t)\end{array}\right]=M_{t}\left[\begin{array}{@{}c@{}}q(0)\\ p(0)\end{array}\right]\!,\quad M_{t}=\left[\begin{array}{@{}cc@{}}\phantom{-}\cos t & \sin t\\ -\sin t & \cos t\end{array}\right]\!.\end{eqnarray}$$

Thus the flow has the expression

(2.5) $$\begin{eqnarray}\unicode[STIX]{x1D711}_{t}(\unicode[STIX]{x1D709},\hspace{2.22198pt}\unicode[STIX]{x1D702})=(\unicode[STIX]{x1D709}\cos t+\unicode[STIX]{x1D702}\sin t,\hspace{2.22198pt}-\unicode[STIX]{x1D709}\sin t+\unicode[STIX]{x1D702}\cos t).\end{eqnarray}$$

When $\unicode[STIX]{x1D6FC}=(\unicode[STIX]{x1D709},\unicode[STIX]{x1D702})\in \mathbb{R}^{2}$ is fixed and $t$ varies, the right-hand side of (2.5) yields the solution that at $t=0$ takes the initial value $(\unicode[STIX]{x1D709},\unicode[STIX]{x1D702})$ . The notation $\unicode[STIX]{x1D711}_{t}(\unicode[STIX]{x1D709},\hspace{2.22198pt}\unicode[STIX]{x1D702})$ emphasizes that, in the flow, it is the parameter $t$ that is seen as fixed, while $(\unicode[STIX]{x1D709},\unicode[STIX]{x1D702})\in \mathbb{R}^{2}$ is regarded as a variable. Geometrically, $\unicode[STIX]{x1D711}_{t}$ is the clockwise rotation of angle $t$ around the origin of the $(q,p)$ -plane.

For a given system (2.1), it is well possible that for some choices of $\unicode[STIX]{x1D6FC}$ and $t$ , the vector $\unicode[STIX]{x1D711}_{t}(\unicode[STIX]{x1D6FC})\in \mathbb{R}^{D}$ is not defined; this will happen if $t$ is outside the interval in which the solution of (2.1) with initial value $\unicode[STIX]{x1D6FC}$ exists. For simplicity in the statements, we shall assume hereafter that $\unicode[STIX]{x1D711}_{t}(\unicode[STIX]{x1D6FC})$ is always defined.

Flows possess the group property: $\unicode[STIX]{x1D711}_{0}$ is the identity map in $\mathbb{R}^{D}$ and, for arbitrary real $s$ and $t$ ,

(2.6) $$\begin{eqnarray}\unicode[STIX]{x1D711}_{t}\circ \unicode[STIX]{x1D711}_{s}=\unicode[STIX]{x1D711}_{s+t}.\end{eqnarray}$$

In particular, for each $t$ ,

(2.7) $$\begin{eqnarray}(\unicode[STIX]{x1D711}_{t})^{-1}=\unicode[STIX]{x1D711}_{-t},\end{eqnarray}$$

that is, $\unicode[STIX]{x1D711}_{-t}$ is the inverse of the map $\unicode[STIX]{x1D711}_{t}$ . For the harmonic oscillator example, the group property simply states that a rotation of angle $s$ followed by a rotation of angle $t$ has the same effect as a rotation of angle $s+t$ .

2.2 Hamiltonian systems

The Hamiltonian formalism (Marsden and Ratiu Reference Marsden and Ratiu1999, Arnol’d Reference Arnol’d1989) is essential for understanding HMC algorithms.

2.2.1 Hamiltonian vector fields

Assume that the dimension $D$ of (2.1) is even, $D=2d$ , and write $x=(q,p)$ with $q,p\in \mathbb{R}^{d}$ . Then the system (2.1) is said to be Hamiltonian if there is a function $H:\mathbb{R}^{2d}\rightarrow \mathbb{R}$ such that, for $i=1,\ldots ,d$ , the scalar components $f^{i}$ of $f$ are given by

$$\begin{eqnarray}f^{i}(q,p)=+\displaystyle \frac{\unicode[STIX]{x2202}H}{\unicode[STIX]{x2202}q^{i}}(q,p),\quad f^{d+i}(q,p)=-\displaystyle \frac{\unicode[STIX]{x2202}H}{\unicode[STIX]{x2202}p^{i}}(q,p).\end{eqnarray}$$

Thus, the system is

$$\begin{eqnarray}\displaystyle \frac{\text{d}q^{i}}{\text{d}t}=+\displaystyle \frac{\unicode[STIX]{x2202}H}{\unicode[STIX]{x2202}q^{i}}(q,p),\quad \displaystyle \frac{\text{d}p^{i}}{\text{d}t}=-\displaystyle \frac{\unicode[STIX]{x2202}H}{\unicode[STIX]{x2202}p^{i}}(q,p),\end{eqnarray}$$

or, in vector notation (Sanz-Serna and Calvo Reference Sanz-Serna and Calvo1994),

(2.8) $$\begin{eqnarray}\displaystyle \frac{\text{d}}{\text{d}t}\left[\begin{array}{@{}c@{}}q\\ p\end{array}\right]=J^{-1}\unicode[STIX]{x1D6FB}H(q,p),\end{eqnarray}$$

where

$$\begin{eqnarray}\unicode[STIX]{x1D6FB}H=\biggl[\displaystyle \frac{\unicode[STIX]{x2202}H}{\unicode[STIX]{x2202}q^{1}},\ldots ,\displaystyle \frac{\unicode[STIX]{x2202}H}{\unicode[STIX]{x2202}q^{d}},\displaystyle \frac{\unicode[STIX]{x2202}H}{\unicode[STIX]{x2202}p^{1}},\ldots ,\displaystyle \frac{\unicode[STIX]{x2202}H}{\unicode[STIX]{x2202}p^{d}}\biggr]^{T}\end{eqnarray}$$

and

$$\begin{eqnarray}J=\left[\begin{array}{@{}cc@{}}0_{d\times d} & -I_{d\times d}\\ I_{d\times d} & \phantom{-}0_{d\times d}\end{array}\right]\!.\end{eqnarray}$$

The function $H$ is called the Hamiltonian, $\mathbb{R}^{2d}$ is the phase space, and $d$ is the number of degrees of freedom.

A system of the special form (2.2) is Hamiltonian if and only if $F=-\unicode[STIX]{x1D6FB}{\mathcal{U}}(q)$ for a suitable real-valued function ${\mathcal{U}}$ , that is,

(2.9) $$\begin{eqnarray}\displaystyle \frac{\text{d}}{\text{d}t}q=M^{-1}p,\quad \displaystyle \frac{\text{d}}{\text{d}t}p=-\unicode[STIX]{x1D6FB}{\mathcal{U}}(q);\end{eqnarray}$$

when that is the case,

(2.10) $$\begin{eqnarray}H(q,p)={\mathcal{T}}(p)+{\mathcal{U}}(q),\quad \text{with }{\mathcal{T}}(p)=\displaystyle \frac{1}{2}p^{T}M^{-1}p.\end{eqnarray}$$

In applications to mechanics, ${\mathcal{T}}$ and ${\mathcal{U}}$ are respectively the potential and kinetic energy and $H$ represents the total energy in the system. The harmonic oscillator (2.3) provides the simplest example; there ${\mathcal{T}}=(1/2)p^{2}$ , ${\mathcal{U}}=(1/2)q^{2}$ .

2.2.2 Symplecticness and preservation of oriented volume

A mapping $\unicode[STIX]{x1D6F7}:\mathbb{R}^{2d}\rightarrow \mathbb{R}^{2d}$ is said to be symplectic or canonical if, at each point $(q,p)\in \mathbb{R}^{2d}$ ,

(2.11) $$\begin{eqnarray}\unicode[STIX]{x1D6F7}^{\prime }(q,p)^{T}J\unicode[STIX]{x1D6F7}^{\prime }(q,p)=J\end{eqnarray}$$

( $\unicode[STIX]{x1D6F7}^{\prime }(q,p)$ denotes the $2d\times 2d$ Jacobian matrix of $\unicode[STIX]{x1D6F7}$ ). The (analytic) condition (2.11) has a geometric interpretation in terms of preservation of two-dimensional areas (Arnol’d Reference Arnol’d1989); this interpretation is not required to understand the rest of the paper.

When $d=1$ , if we set $\unicode[STIX]{x1D6F7}(q,p)=(q^{\ast },p^{\ast })$ , the condition (2.11), after multiplying the matrices on the left-hand side, is seen to be equivalent to

$$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}q^{\ast }}{\unicode[STIX]{x2202}q}\displaystyle \frac{\unicode[STIX]{x2202}p^{\ast }}{\unicode[STIX]{x2202}p}-\displaystyle \frac{\unicode[STIX]{x2202}q^{\ast }}{\unicode[STIX]{x2202}p}\displaystyle \frac{\unicode[STIX]{x2202}p^{\ast }}{\unicode[STIX]{x2202}q}=1.\end{eqnarray}$$

The left-hand side is the Jacobian determinant of $\unicode[STIX]{x1D6F7}$ and therefore the transformation $\unicode[STIX]{x1D6F7}$ is symplectic if and only if the mapping $(q,p)\mapsto (q^{\ast },p^{\ast })$ preserves oriented area in the $(q,p)$ -plane, that is, for any domain $D$ the oriented area of the image $\unicode[STIX]{x1D6F7}(D)\subset \mathbb{R}^{2}$ coincides with the oriented area of $D$ .Footnote ¹ For instance, for each $t$ , the rotation in (2.5) is a symplectic transformation in  $\mathbb{R}^{2}$ .

For general $d$ the following result holds (Arnol’d Reference Arnol’d1989, Section 38).

Proposition 2.2. For a symplectic transformation the determinant of $\unicode[STIX]{x1D6F7}^{\prime }$ equals 1. Therefore symplectic transformations preserve the oriented volume in $\mathbb{R}^{2d}$ , that is, $\unicode[STIX]{x1D6F7}(D)$ and $D$ have the same oriented volume for each domain $D\subset \mathbb{R}^{2d}$ .

For $d>1$ , preservation of oriented volume is a strictly weaker property than symplecticness.

The proof of the following two results is easy using (2.11).

Proposition 2.3. The composition $\unicode[STIX]{x1D6F7}_{1}\circ \unicode[STIX]{x1D6F7}_{2}$ of two symplectic mappings is itself symplectic.

Proposition 2.4. The change of variables $(q,p)=\unicode[STIX]{x1D6F7}(\bar{q},\bar{p})$ with $\unicode[STIX]{x1D6F7}$ symplectic transforms the Hamiltonian system of differential equations (2.8) into a system for $(\bar{q},\bar{p})$ that is also Hamiltonian. Moreover, the Hamiltonian function $\bar{H}$ of the transformed system is the result of changing variables in $H$ , that is, $\bar{H}=H\circ \unicode[STIX]{x1D6F7}$ .

In view of the following important general result (Marsden and Ratiu Reference Marsden and Ratiu1999, Proposition 2.6.2), the symplecticness of the rotation (2.5) noted above is a manifestation of the Hamiltonian character of the harmonic oscillator.

Theorem 2.5. Let $D=2d$ . The system (2.1) with flow $\unicode[STIX]{x1D711}_{t}$ is Hamiltonian if and only if, for each real $t$ , $\unicode[STIX]{x1D711}_{t}$ is a symplectic mapping.

Thus symplecticness is a characteristic property that allows us to decide whether a differential system is Hamiltonian or otherwise in terms of its flow, without knowing the vector field (right-hand side) $f$ of the equation.

We recall that a flow preserves oriented volume if and only if the corresponding vector field $f$ is divergence-free ( $\unicode[STIX]{x1D6FB}\cdot f=\sum _{i}(\unicode[STIX]{x2202}/\unicode[STIX]{x2202}x^{i})f^{i}=0$ ). If $d>1$ , there are divergence-free differential systems in $\mathbb{R}^{2d}$ that are not Hamiltonian; their flows preserve oriented volume but are not symplectic.

The behaviour of the solutions of Hamiltonian problems is very different from that encountered in ‘general’ systems; some features that are ‘the rule’ in Hamiltonian systems are exceptional in non-Hamiltonian systems. Such special behaviour of Hamiltonian solutions may always be traced back to the symplecticness of the flow. As a very simple example, we consider once more the harmonic oscillator (2.3). The origin is a centre: a neutrally stable equilibrium surrounded by periodic trajectories. Small perturbations of the right-hand side of (2.3) generically destroy the centre; after perturbation the trajectories become spirals and the origin becomes either an asymptotically stable node (trajectories spiral in) or an unstable node (trajectories spiral out). However, if the perturbation is such that the perturbed system is also Hamiltonian, then the centre will not disappear under small perturbations.

2.2.3 Preservation of energy

If $(q(t),p(t))$ is a solution of (2.8), then

$$\begin{eqnarray}\displaystyle \frac{\text{d}}{\text{d}t}H(q(t),p(t))=\unicode[STIX]{x1D6FB}H(q(t),p(t))^{T}J^{-1}\unicode[STIX]{x1D6FB}H(q(t),p(t))=0,\end{eqnarray}$$

because $J^{-1}$ is skew-symmetric. Therefore we may state the following.

Theorem 2.6. The value of the Hamiltonian function is preserved by the flow of the corresponding Hamiltonian system, that is, $H\circ \unicode[STIX]{x1D711}_{t}=H$ for each real  $t$ .

In applications to the physical sciences, this result is usually the mathematical expression of the principle of conservation of energy. Unlike symplecticness, which is a characteristic property, conservation of energy on its own does not ensure that the underlying system is Hamiltonian. There are many examples of non-Hamiltonian systems whose solutions preserve the value of a suitable energy function.

2.2.4 Preservation of the canonical probability measure

Let $\unicode[STIX]{x1D6FD}$ denote a positive constant and assume that $H$ is such that

$$\begin{eqnarray}Z=\int _{\mathbb{R}^{2d}}\exp (-\unicode[STIX]{x1D6FD}H(q,p))\,\text{d}q\,\text{d}p<\infty .\end{eqnarray}$$

Then we have the following result, which is a direct consequence of the fact that $\unicode[STIX]{x1D711}_{t}$ preserves both the volume element $\text{d}q\,\text{d}p$ (Proposition 2.2) and the value of $\exp (-\unicode[STIX]{x1D6FD}H)$ (because, according to Theorem 2.6, it preserves the value of $H$ ).

Theorem 2.7. The probability measure $\unicode[STIX]{x1D707}$ in $\mathbb{R}^{2d}$ with density (with respect to the ordinary Lebesgue measure) $Z^{-1}\exp (-\unicode[STIX]{x1D6FD}H(q,p))$ is preserved by the flow of the Hamiltonian system (2.8), that is, $\unicode[STIX]{x1D707}(\unicode[STIX]{x1D711}_{t}(D))=\unicode[STIX]{x1D707}(D)$ for each domain $D\subset \mathbb{R}^{2d}$ and each real $t$ .

In statistical mechanics (Allen and Tildesley Reference Allen and Tildesley1987, Frenkel and Smit Reference Frenkel and Smit2002, Tuckerman Reference Tuckerman2010), if (2.8) describes the dynamics of a physical system and $\unicode[STIX]{x1D6FD}$ is the inverse of the absolute temperature, then $\unicode[STIX]{x1D707}$ is the canonical measure that governs the distribution of $(q,p)$ over an ensemble of many copies of the given system when the system is in contact with a heat bath at constant temperature, that is, $Z^{-1}\exp (-\unicode[STIX]{x1D6FD}H(q,p))\,\text{d}q\,\text{d}p$ represents the fraction of copies with momenta between $p$ and $p+\text{d}p$ , and configuration between $q$ and $q+\,\text{d}q$ . Note that under the canonical distribution, (local) minima of the energy $H$ correspond to (local) maxima of the probability density function, that is, to modes of the distribution. Also, if the temperature decreases ( $\unicode[STIX]{x1D6FD}$ increases), it is less likely to find the system at a location $(q,p)$ with high energy.

For Hamiltonian functions of the particular form in (2.10), we may factorize

$$\begin{eqnarray}\exp (-\unicode[STIX]{x1D6FD}H(q,p))=\exp \biggl(-\displaystyle \frac{1}{2}\unicode[STIX]{x1D6FD}\,p^{T}M^{-1}p\biggr)\times \exp (-\unicode[STIX]{x1D6FD}\,{\mathcal{U}}(q))\end{eqnarray}$$

and therefore, under the canonical distribution, $q$ and $p$ are stochastically independent. The (marginal) distribution of the configuration variables $q$ has probability density function proportional to $\exp (-\unicode[STIX]{x1D6FD}\,{\mathcal{U}}(q))$ . The momenta $p$ possess a Gaussian distribution with zero mean and covariance matrix equal to $M$ . These distributions are associated with the names of Boltzmann, Gibbs and Maxwell. Hereafter we refer to the canonical measure as the Boltzmann–Gibbs distribution.

2.3 Reversible systems

Assume now that $S$ is a linear involution in $\mathbb{R}^{D}$ , that is, a linear map such that $S(S(x))=x$ for each $x$ . A mapping $\unicode[STIX]{x1D6F7}:\mathbb{R}^{D}\rightarrow \mathbb{R}^{D}$ is said to be reversible (with respect to $S$ ) if, for each $x$ , $S(\unicode[STIX]{x1D6F7}(x))=\unicode[STIX]{x1D6F7}^{-1}(S(x))$ or, more compactly,

(2.12) $$\begin{eqnarray}S\circ \unicode[STIX]{x1D6F7}=\unicode[STIX]{x1D6F7}^{-1}\circ S.\end{eqnarray}$$

The following results have easy proofs.

Proposition 2.8. If $\unicode[STIX]{x1D6F7}$ is $S$ reversible, then

$$\begin{eqnarray}|\text{det}\,\unicode[STIX]{x1D6F7}^{\prime }(S(\unicode[STIX]{x1D6F7}(x)))|=|\text{det}\,\unicode[STIX]{x1D6F7}^{\prime }(x)|^{-1},\end{eqnarray}$$

for each $x$ .

Systems of differential equations that satisfy the conditions in the theorem are said to be reversible (with respect to $S$ ). Systems of the particular form (2.2) are reversible with respect to the momentum flip involution

(2.13) $$\begin{eqnarray}S(q,p)=(q,-p).\end{eqnarray}$$

If (2.2) describes a mechanical system, then (2.12) expresses the well-known time-reversibility of mechanics: if $(q_{in},p_{in})$ is the initial state of a system and $(q_{f},p_{f})$ the final state after $t$ units of time have elapsed, then the state $(q_{f},-p_{f})$ evolves in $t$ units of time to the state $(q_{in},-p_{in})$ .

Figure 6.1 (see page 160) illustrates the reversibility of the Hamiltonian flow corresponding to a one-degree-of-freedom double-well potential.

As is the case for Hamiltonian systems, reversible systems have flows with special geometric properties, not shared by ‘general’ systems (Lamb and Roberts Reference Lamb and Roberts1998).

2.4 The Lie bracket

If $\unicode[STIX]{x1D711}_{t}^{(f)}$ and $\unicode[STIX]{x1D711}_{t}^{(g)}$ denote respectively the flows of the $D$ -dimensional systems

$$\begin{eqnarray}\displaystyle \frac{\text{d}}{\text{d}t}x=f(x),\quad \displaystyle \frac{\text{d}}{\text{d}t}x=g(x),\end{eqnarray}$$

in general $\unicode[STIX]{x1D711}_{s}^{(g)}\circ \unicode[STIX]{x1D711}_{t}^{(f)}\neq \unicode[STIX]{x1D711}_{t}^{(f)}\circ \unicode[STIX]{x1D711}_{s}^{(g)}$ ; in fact, a Taylor expansion shows that, as $t$ , $s$ approach $0$ ,

$$\begin{eqnarray}\unicode[STIX]{x1D711}_{s}^{(g)}(\unicode[STIX]{x1D711}_{t}^{(f)}(x))-\unicode[STIX]{x1D711}_{t}^{(f)}(\unicode[STIX]{x1D711}_{s}^{(g)}(x))=st\,[f,g](x)+O(t^{3}+s^{3}),\end{eqnarray}$$

where the Lie bracket (or Lie–Jacobi bracket or commutator) $[f,g]$ of $f$ and $g$ (Arnol’d Reference Arnol’d1989, Section 39) is the mapping $\mathbb{R}^{D}\rightarrow \mathbb{R}^{D}$ that at $x\in \mathbb{R}^{D}$ takes the value

(2.14) $$\begin{eqnarray}[f,g](x)=g^{\prime }(x)f(x)-f^{\prime }(x)g(x).\end{eqnarray}$$

Thus the magnitude of $[f,g]$ measures the lack of commutativity of the corresponding flows. The following result holds.

Theorem 2.12. With the preceding notation, $\unicode[STIX]{x1D711}_{s}^{(g)}\circ \unicode[STIX]{x1D711}_{t}^{(f)}=\unicode[STIX]{x1D711}_{t}^{(f)}\circ \unicode[STIX]{x1D711}_{s}^{(g)}$ for arbitrary $t$ and $s$ if and only if the Lie bracket $[f,g](x)$ vanishes at each $x\in \mathbb{R}^{D}$ . When these conditions hold, we say that $f$ and $g$ commute.

For commuting $f$ and $g$ , $\unicode[STIX]{x1D711}_{t}^{(g)}\circ \unicode[STIX]{x1D711}_{t}^{(f)}=\unicode[STIX]{x1D711}_{t}^{(f)}\circ \unicode[STIX]{x1D711}_{s}^{(g)}$ provides the $t$ flow of the system

$$\begin{eqnarray}\displaystyle \frac{\text{d}}{\text{d}t}x=f(x)+g(x).\end{eqnarray}$$

The Lie bracket is skew-symmetric: $[f,g]=-[g,f]$ for arbitrary $f$ and $g$ . In addition, it satisfies the Jacobi identity,

$$\begin{eqnarray}[f_{1},[f_{2},f_{3}]]+[f_{2},[f_{3},f_{1}]]+[f_{3},[f_{1},f_{2}]]=0,\end{eqnarray}$$

for any $f_{1}$ , $f_{2}$ , $f_{3}$ . In this way the vector space of all vector fields together with the operation $[\cdot ,\cdot ]$ is a Lie algebra.

For Hamiltonian systems it is possible to work in terms of the so-called Poisson bracket of the Hamiltonian functions (Arnol’d Reference Arnol’d1989, Section 40), rather than in terms of the Lie bracket of the fields.

Theorem 2.13. If the fields $f$ and $g$ are Hamiltonian with Hamiltonian functions $H$ and $K$ respectively, that is, $f=J^{-1}\unicode[STIX]{x1D6FB}H$ , $g=J^{-1}\unicode[STIX]{x1D6FB}K$ , then $[f,g]$ is also a Hamiltonian vector field. Moreover, the Hamiltonian function of $[f,g]$ is given by $-\{H,K\}$ ,Footnote ² where $\{H,K\}$ is the Poisson bracket of the functions $H$ and $K$ , defined as

$$\begin{eqnarray}\{H,K\}=(\unicode[STIX]{x1D6FB}H)^{T}J^{-1}\unicode[STIX]{x1D6FB}K.\end{eqnarray}$$

For real-valued functions in $\mathbb{R}^{2d}$ , the Poisson bracket operation is skew-symmetric, that is, $\{H,K\}=-\{K,H\}$ , and satisfies the Jacobi identity

$$\begin{eqnarray}\{H_{1},\{H_{2},H_{3}\}\}+\{H_{2},\{H_{3},H_{1}\}\}+\{H_{3},\{H_{1},H_{2}\}\}=0.\end{eqnarray}$$

Let $f$ and $g$ be reversible vector fields. Differentiation in $g(S(x))=-S(g(x))$ implies for the Jacobian that $g^{\prime }(S(x))S=-S(g^{\prime }(x))$ , and it follows that $g^{\prime }(S(x))f(S(x)))=S(g^{\prime }(x)f(x))$ (no minus sign!). Then the Lie bracket of two reversible fields is not reversible, but rather satisfies the following property.

Proposition 2.14. If two vector fields $f$ and $g$ are $S$ -reversible, then

$$\begin{eqnarray}[f,g](S(x))=S([f,g](x)).\end{eqnarray}$$

For three $S$ -reversible vector fields, the iterated commutator $[f_{1},[f_{2},f_{3}]]$ is $S$ -reversible.

3.1 Preliminaries

Each one-step numerical method or one-step integrator for (2.1) is described by a map $\unicode[STIX]{x1D713}_{h}:\mathbb{R}^{D}\rightarrow \mathbb{R}^{D}$ that depends on a real parameter $h$ , the step size. Given the initial value $\unicode[STIX]{x1D6FC}$ , and a value of $h$ ( $h\neq 0$ ), the integrator generates a numerical trajectory, $x_{0},x_{1},x_{2},\ldots ,$ defined by $x_{0}=\unicode[STIX]{x1D6FC}$ and, iteratively,

(3.1) $$\begin{eqnarray}x_{n+1}=\unicode[STIX]{x1D713}_{h}(x_{n}),\quad n=0,1,2,\ldots .\end{eqnarray}$$

To compute $x_{n+1}$ when $x_{n}$ has already been found is to perform a (time) step. For each $n$ , the vector $x_{n}$ is an approximation to the value at time $t_{n}=nh$ of the solution $x(t)$ of (2.1) with initial condition $x(0)=\unicode[STIX]{x1D6FC}$ , i.e. to $\unicode[STIX]{x1D711}_{t_{n}}(\unicode[STIX]{x1D6FC})$ . Typically $h$ is positive and then the integration is forward in time $0=t_{0}<t_{1}<t_{2}<\cdots \,$ , but in some applications it may be of interest to use $h<0$ so as to get $0=t_{0}>t_{1}>t_{2}>\cdots \,.$

The simplest and best known integrator, Euler’s rule, with

(3.2) $$\begin{eqnarray}x_{n+1}=x_{n}+hf(x_{n}),\end{eqnarray}$$

corresponds to the mapping $\unicode[STIX]{x1D713}_{h}(x)=x+hf(x)$ . It uses one evaluation of $f$ per step. Explicit $s$ -stage Runge–Kutta formulas use $s$ evaluations of $f$ per step, $s=1,2,\ldots ,$ and are therefore $s$ times more expensive per step than Euler’s rule; examples include Runge’s method

$$\begin{eqnarray}\unicode[STIX]{x1D713}_{h}(x)=x+hf\biggl(x+\displaystyle \frac{h}{2}f(x)\biggr)\end{eqnarray}$$

(with two stages), several well-known formulas of Kutta with four stagesFootnote ³ and the formulas within the popular MATLAB function ode45. A method with $s$ stages will be competitive with Euler’s rule only if it gives more accurate approximations than Euler’s rule when this is operated with a step size $s$ times shorter, so as to equalize computational costs.

Implicit Runge–Kutta integrators are also used in practice, and then $\unicode[STIX]{x1D713}_{h}$ is defined by means of algebraic equations. For instance, the midpoint rule has

$$\begin{eqnarray}\unicode[STIX]{x1D713}_{h}(x)=x+hf\biggl(\displaystyle \frac{1}{2}(x+\unicode[STIX]{x1D713}_{h}(x))\biggr),\end{eqnarray}$$

that is,

(3.3) $$\begin{eqnarray}x_{n+1}=x_{n}+hf\biggl(\displaystyle \frac{1}{2}(x_{n}+x_{n+1})\biggr).\end{eqnarray}$$

Here, computing $x_{n+1}$ when $x_{n}$ is already known requires us to solve a system of algebraic equations in $\mathbb{R}^{D}$ . There are many more useful examples of implicit Runge–Kutta methods, including the so-called Gauss methods.

A one-step integrator is called symmetric or self-adjoint Footnote ⁴ if

(3.4) $$\begin{eqnarray}(\unicode[STIX]{x1D713}_{h})^{-1}=\unicode[STIX]{x1D713}_{-h},\end{eqnarray}$$

so as to mimic the property (2.7) of the exact solution flow. The midpoint rule (3.3) provides an example. Explicit Runge–Kutta methods are never symmetric.

Multistep integrators, including the well-known Adams formulas, where the computation of $x_{n+1}$ requires knowledge of $k\geq 2$ past values $x_{n},x_{n-1},$ $\ldots ,x_{n-k+1}$ , may be very efficient, but will not be considered in this paper; they have seldom been applied in the context of sampling algorithms.

3.2 Order

Since, for the exact solution, the sequence $x(0)=\unicode[STIX]{x1D6FC},x(t_{1}),x(t_{2}),\ldots$ satisfies $x(t_{n+1})=\unicode[STIX]{x1D711}_{h}(x(t_{n}))$ , $n=0,1,\ldots ,$ rather than (3.1), for the numerical integrator to make sense it is necessary that $\unicode[STIX]{x1D713}_{h}$ be an approximation to $\unicode[STIX]{x1D711}_{h}$ . The integrator is said to be consistent if, at each fixed $x\in \mathbb{R}^{D}$ , $\unicode[STIX]{x1D713}_{h}(x)-\unicode[STIX]{x1D711}_{h}(x)=O(h^{2})$ as $h\rightarrow 0$ . If $\unicode[STIX]{x1D713}_{h}(x)-\unicode[STIX]{x1D711}_{h}(x)=O(h^{\unicode[STIX]{x1D708}+1})$ , $\unicode[STIX]{x1D708}$ a positive integer, then the integrator is (consistent) of order $\geq \unicode[STIX]{x1D708}$ . A method of order $\geq \unicode[STIX]{x1D708}$ that is not of order $\geq \unicode[STIX]{x1D708}+1$ is said to have order $\unicode[STIX]{x1D708}$ . Euler’s rule (3.2) has order 1, the four-stage formulas of Kutta have order 4 and the midpoint rule (3.3) has order 2. We shall see later (Theorem 4.7) that the order of a symmetric integrator is an even integer. All integrators to be considered hereafter are assumed to be consistent.

The vector $\unicode[STIX]{x1D713}_{h}(x)-\unicode[STIX]{x1D711}_{h}(x)$ is called the local error at $x$ : it is the difference between the result of a single time-step of the numerical method starting from $x$ and the result of the application of the $h$ -flow to the same point $x$ . For future reference, we note that the expansion of the local error for Euler’s method is given by

(3.5) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D713}_{h}(x)-\unicode[STIX]{x1D711}_{h}(x) & = & \displaystyle (x+hf(x))-\biggl(x+hf(x)+\displaystyle \frac{h^{2}}{2}f^{\prime }(x)f(x)+O(h^{3})\biggr)\nonumber\\ \displaystyle & = & \displaystyle -\displaystyle \frac{h^{2}}{2}f^{\prime }(x)f(x)+O(h^{3}).\end{eqnarray}$$

In the expansion of the true solution flow we have used that $(\text{d}/\text{d}t)x=f(x)$ and $(\text{d}^{2}/\text{d}t^{2})x=f^{\prime }(x)(\text{d}/\text{d}t)x=f^{\prime }(x)f(x)$ .

The local error does not per se give information on the global error at $t_{n}$ , i.e. on the difference $x_{n}-x(t_{n})$ . This is because $x_{n}=\unicode[STIX]{x1D713}_{h}(x_{n-1})$ while $x(t_{n})$ is the result of the application of the $h$ -flow to $x(t_{n-1})$ : $\unicode[STIX]{x1D713}_{h}$ and $\unicode[STIX]{x1D711}_{h}$ are not applied to the same point.Footnote ⁵ However, the following result holds.

Theorem 3.2. If the (one-step) integrator is consistent of order $\unicode[STIX]{x1D708}$ , then, for each fixed initial value $x_{0}=x(0)$ and $T>0$ ,

(3.6) $$\begin{eqnarray}\max _{0\leq t_{n}\leq T}|x_{n}-x(t_{n})|=O(h^{\unicode[STIX]{x1D708}}),\quad h\rightarrow 0+.\end{eqnarray}$$

This expresses the fact that the integrator is convergent of order $\unicode[STIX]{x1D708}$ . Of course, the corresponding result holds if the integration is carried out backward in time over $-T\leq t\leq 0$ .

Remark 3.3. The global errors $x_{n}-x(t_{n})$ , in addition to possessing an $O(h^{\unicode[STIX]{x1D708}})$ bound as above, have an asymptotic expansion in powers of $h$ . For instance, if we restrict the attention to the leading $h^{\unicode[STIX]{x1D708}}$ -term in the expansion, we have

$$\begin{eqnarray}x_{n}-x(t_{n})=h^{\unicode[STIX]{x1D708}}a(x(0),t_{n})+h^{\unicode[STIX]{x1D708}}r(x(0),t_{n};h),\end{eqnarray}$$

where the function $a$ is independent of $h$ and, for each fixed initial value $x(0)$ , $\max _{0\leq t\leq T}|r(x(0),t,h)|$ tends to $0$ as $h\rightarrow 0\!+$ . Thus, for $h$ sufficiently small, the global error is approximately equal to $h^{\unicode[STIX]{x1D708}}a(x(0),t_{n})$ : halving $h$ divides the error by a factor $2^{\unicode[STIX]{x1D708}}$ .

3.3 Splitting methods

Euler’s method and the other integrators mentioned above may be applied to any given system (2.1). Other techniques only make sense for particular classes of systems and cannot therefore be incorporated into general-purpose software; however, they have gained much popularity in recent decades due, among other things, to their role in geometric integration (see Section 4). Of special importance to us is the class of splitting integrators that we consider next. Blanes and Casas (Reference Blanes and Casas2016) provide a very good source of information.

3.3.1 The Lie–Trotter formula

Splitting methods are applicable to cases where the right-hand side of (2.1) may be split into two parts,

(3.7) $$\begin{eqnarray}\displaystyle \frac{\text{d}}{\text{d}t}x=f(x)=f^{(\text{A})}(x)+f^{(\text{B})}(x),\end{eqnarray}$$

in such a way that the flows $\unicode[STIX]{x1D711}_{t}^{(\text{A})}$ and $\unicode[STIX]{x1D711}_{t}^{(\text{B})}$ of the split systems

(3.8) $$\begin{eqnarray}\displaystyle \frac{\text{d}}{\text{d}t}x=f^{(\text{A})}(x),\quad \displaystyle \frac{\text{d}}{\text{d}t}x=f^{(\text{B})}(x),\end{eqnarray}$$

are available analytically. To avoid trivial cases, we shall hereafter assume that the Lie bracket $[f^{(\text{A})},f^{(\text{B})}]$ does not vanish identically, because otherwise Theorem 2.12 shows that (3.7) may also be solved analytically without resorting to numerical approximations.

Systems of the particular form (2.2) provide an important example. When they are split by taking

(3.9) $$\begin{eqnarray}\displaystyle (\text{A})\qquad \displaystyle \frac{\text{d}}{\text{d}t}q & = & \displaystyle M^{-1}p,\quad \displaystyle \frac{\text{d}}{\text{d}t}p=0\end{eqnarray}$$

and

(3.10) $$\begin{eqnarray}\displaystyle (\text{B})\qquad \displaystyle \frac{\text{d}}{\text{d}t}q & =0,\quad \displaystyle \frac{\text{d}}{\text{d}t}p=F(q), & \displaystyle\end{eqnarray}$$

the flows are explicitly given by

$$\begin{eqnarray}\unicode[STIX]{x1D711}_{t}^{(\text{A})}(q,p)=(q+tM^{-1}p,p)\end{eqnarray}$$

and

$$\begin{eqnarray}\unicode[STIX]{x1D711}_{t}^{(\text{B})}(q,p)=(q,p+tF(q)).\end{eqnarray}$$

In molecular dynamics (Schlick Reference Schlick2002) these mappings are respectively called a drift ( $q$ advances with constant speed and the momentum $p$ remains constant) and a kick (the system stays in its current configuration, and the momentum is incremented by the action of the force). A direct computation of the Lie bracket shows that, except in the trivial case where $F$ vanishes identically, the vector fields in (3.9)–(3.10) do not commute.

A simple Taylor expansion proves that the Lie–Trotter formula

(3.11) $$\begin{eqnarray}\unicode[STIX]{x1D713}_{h}=\unicode[STIX]{x1D711}_{h}^{(\text{B})}\circ \unicode[STIX]{x1D711}_{h}^{(\text{A})}\end{eqnarray}$$

defines a first-order integrator for (3.7):

(3.12) $$\begin{eqnarray}\unicode[STIX]{x1D713}_{h}(x)-\unicode[STIX]{x1D711}_{h}(x)=\displaystyle \frac{h^{2}}{2}[f^{(\text{A})},f^{(\text{B})}](x)+O(h^{3}).\end{eqnarray}$$

(Observe that, not unexpectedly, the coefficient of the leading power of $h$ is proportional to the Lie bracket.) While $f^{(\text{A})}$ and $f^{(\text{B})}$ contribute simultaneously to the change of $x$ in (3.7), they do so successively in the Lie–Trotter integrator. By swapping the roles of $\text{A}$ and $\text{B}$ , we have the alternative integrator

$$\begin{eqnarray}\unicode[STIX]{x1D713}_{h}=\unicode[STIX]{x1D711}_{h}^{(\text{A})}\circ \unicode[STIX]{x1D711}_{h}^{(\text{B})}.\end{eqnarray}$$

3.3.2 Strang’s formula

The most popular splitting integrator for (3.7) corresponds to Strang’s formula (Strang Reference Strang1963)

(3.13) $$\begin{eqnarray}\unicode[STIX]{x1D713}_{h}=\unicode[STIX]{x1D711}_{(1/2)h}^{(\text{B})}\circ \unicode[STIX]{x1D711}_{h}^{(\text{A})}\circ \unicode[STIX]{x1D711}_{(1/2)h}^{(\text{B})}\end{eqnarray}$$

and, as a Taylor expansion shows, has second-order accuracy, $\unicode[STIX]{x1D713}_{h}(x)-\unicode[STIX]{x1D711}_{h}(x)=O(h^{2})$ as $h\rightarrow 0$ . More precisely,

(3.14) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D713}_{h}(x)-\unicode[STIX]{x1D711}_{h}(x) & = & \displaystyle \displaystyle \frac{h^{3}}{12}[f^{(\text{A})},[f^{(\text{A})},f^{(\text{B})}]](x)\nonumber\\ \displaystyle & & \displaystyle \quad +\,\displaystyle \frac{h^{3}}{24}[f^{(\text{B})},[f^{(\text{A})},f^{(\text{B})}]](x)+O(h^{4}),\end{eqnarray}$$

so that the leading term of the local error is a linear combination of two so-called iterated Lie brackets.

When applied to the particular case (2.2), the formula (3.13) yields the well-known velocity Verlet integrator, the method of choice in molecular dynamics; the step $n\rightarrow n+1$ comprises two kicks of duration $h/2$ separated by a drift of duration $h$ :

$$\begin{eqnarray}\displaystyle p_{n+1/2} & = & \displaystyle p_{n}+\displaystyle \frac{h}{2}F(q_{n}),\nonumber\\ \displaystyle q_{n+1} & = & \displaystyle q_{n}+hM^{-1}p_{n+1/2},\nonumber\\ \displaystyle p_{n+1} & = & \displaystyle p_{n+1/2}+\displaystyle \frac{h}{2}F(q_{n+1}).\nonumber\end{eqnarray}$$

The evaluations of $F$ represent the bulk of the computational cost of the algorithm. The value $F(q_{n+1})$ to be used in the first kick of the next step, $n+1\rightarrow n+2$ , is the same used in the second kick of the current step. In this way, while the very first step requires two evaluations of $F$ , subsequent steps only need one.

When $N$ steps of the method (3.13) are taken, the map that advances the numerical solution from $x_{0}$ to $x_{N}$ , that is,

$$\begin{eqnarray}\unicode[STIX]{x1D713}_{h}^{N}=\overbrace{(\unicode[STIX]{x1D711}_{(1/2)h}^{(\text{B})}\circ \unicode[STIX]{x1D711}_{h}^{(\text{A})}\circ \unicode[STIX]{x1D711}_{(1/2)h}^{(\text{B})})\circ \cdots \circ (\unicode[STIX]{x1D711}_{(1/2)h}^{(\text{B})}\circ \unicode[STIX]{x1D711}_{h}^{(\text{A})}\circ \unicode[STIX]{x1D711}_{(1/2)h}^{(\text{B})})}^{N\text{ times}},\end{eqnarray}$$

may be rewritten with the help of the group property (2.6) in the leapfrog form

$$\begin{eqnarray}\unicode[STIX]{x1D713}_{h}^{N}=\unicode[STIX]{x1D711}_{(1/2)h}^{(\text{B})}\circ \overbrace{(\unicode[STIX]{x1D711}_{h}^{(\text{A})}\circ \unicode[STIX]{x1D711}_{h}^{(\text{B})})\circ \cdots \circ (\unicode[STIX]{x1D711}_{h}^{(\text{A})}\circ \unicode[STIX]{x1D711}_{h}^{(\text{B})})}^{N-1\text{ times}}\circ \,\unicode[STIX]{x1D711}_{h}^{(\text{A})}\circ \unicode[STIX]{x1D711}_{(1/2)h}^{(\text{B})};\end{eqnarray}$$

now the right-hand side only uses $N+1$ times the flow $\unicode[STIX]{x1D711}_{t}^{(\text{B})}$ . In the particular case (3.9)–(3.10), the combination $\unicode[STIX]{x1D711}_{h}^{(\text{A})}\circ \unicode[STIX]{x1D711}_{h}^{(\text{B})}$ corresponds to the following formulas to advance the numerical solution, $n=1,\ldots ,N-1$ :

$$\begin{eqnarray}\displaystyle p_{n+1/2} & =p_{n-1/2}+hF(q_{n}), & \displaystyle \nonumber\\ \displaystyle q_{n+1} & =q_{n}+hM^{-1}p_{n+1/2}; & \displaystyle \nonumber\end{eqnarray}$$

$p$ jumps over $q$ and then $q$ jumps over $p$ like children playing leapfrog. The leapfrog implementation makes apparent the truth of an earlier observation: $N$ steps of the velocity Verlet integrator may be implemented with $N+1$ evaluations of $F$ .

Strang’s method is symmetric:

$$\begin{eqnarray}\displaystyle (\unicode[STIX]{x1D713}_{h})^{-1} & = & \displaystyle (\unicode[STIX]{x1D711}_{(1/2)h}^{(\text{B})})^{-1}\!\circ (\unicode[STIX]{x1D711}_{h}^{(\text{A})})^{-1}\!\circ (\unicode[STIX]{x1D711}_{(1/2)h}^{(\text{B})})^{-1}\nonumber\\ \displaystyle & = & \displaystyle \unicode[STIX]{x1D711}_{-(1/2)h}^{(\text{B})}\circ \unicode[STIX]{x1D711}_{-h}^{(\text{A})}\circ \unicode[STIX]{x1D711}_{-(1/2)h}^{(\text{B})}\nonumber\\ \displaystyle & = & \displaystyle \unicode[STIX]{x1D713}_{-h}.\nonumber\end{eqnarray}$$

It is clear that the symmetry is a consequence of the palindromic structure of (3.13), that is, the formula reads the same from left to right as from right to left.

The roles of $\text{A}$ and $\text{B}$ in (3.13) may be interchanged:

(3.15) $$\begin{eqnarray}\unicode[STIX]{x1D713}_{h}=\unicode[STIX]{x1D711}_{(1/2)h}^{(\text{A})}\circ \unicode[STIX]{x1D711}_{h}^{(\text{B})}\circ \unicode[STIX]{x1D711}_{(1/2)h}^{(\text{A})}.\end{eqnarray}$$

For (3.9)–(3.10) one then obtains the position Verlet integrator: one step comprises two drifts of duration $h/2$ and one kick of duration $h$ .

3.3.3 Splitting formulas with more stages

It is of course possible to use splitting formulas more sophisticated than Strang’s. For instance, for any choice of the real parameter $b$ , we may consider the integrator

(3.16) $$\begin{eqnarray}\unicode[STIX]{x1D713}_{h}=\unicode[STIX]{x1D711}_{bh}^{(\text{B})}\circ \unicode[STIX]{x1D711}_{(1/2)h}^{(\text{A})}\circ \unicode[STIX]{x1D711}_{(1-2b)h}^{(\text{B})}\circ \unicode[STIX]{x1D711}_{(1/2)h}^{(\text{A})}\circ \unicode[STIX]{x1D711}_{bh}^{(\text{B})},\end{eqnarray}$$

where we observe that, in one step, the $\text{A}$ and $\text{B}$ flows of the systems (3.8) act for a total duration of $h$ units of time each to ensure consistency. Due to its palindromic structure, this integrator is symmetric, and from Theorem 4.7 its order is at least two. It turns out that the order is exactly two for all choices of $b$ (Blanes, Casas and Sanz-Serna Reference Blanes, Casas and Sanz-Serna2014). The order of splitting integrators is discussed later in relation to the concept of modified equations.

Even though three $\text{B}$ flows and two $\text{A}$ flows feature in (3.16), $N$ steps of the integrator only require the computation of $2N+1$ $\text{B}$ flows and $2N$ $\text{A}$ flows; this is seen by combining flows as we did above for the Strang case. We say that (3.16) is a palindromic two-stage integrator.Footnote ⁶ The method (3.16) may be denoted by

(3.17) $$\begin{eqnarray}(b,1/2,(1-2b),1/2,b).\end{eqnarray}$$

Similarly, one may consider the two-parameter family of palindromic three-stage splittings

(3.18) $$\begin{eqnarray}(b,a,1/2-b,1-2a,1/2-b,a,b).\end{eqnarray}$$

A full description of this family is given by Campos and Sanz-Serna (Reference Campos and Sanz-Serna2017), who suggest parameter choices for various applications. There is a unique choice of $a$ and $b$ resulting in a fourth-order method often associated with Yoshida’s name (Yoshida Reference Yoshida1990); for all other choices, the order is $\unicode[STIX]{x1D708}=2$ .

The family of palindromic $s$ -stage splitting formulas is given by

(3.19) $$\begin{eqnarray}(b_{1},a_{1},b_{2},a_{2},\ldots ,a_{s^{\prime }},b_{s^{\prime }+1},a_{s^{\prime }},\ldots ,a_{2},b_{2},a_{1},b_{1})\end{eqnarray}$$

if $s=2s^{\prime }$ is even, and by

(3.20) $$\begin{eqnarray}(b_{1},a_{1},b_{2},a_{2},\ldots ,b_{s^{\prime }},a_{s^{\prime }},b_{s^{\prime }},\ldots ,a_{2},b_{2},a_{1},b_{1})\end{eqnarray}$$

if $s=2s^{\prime }-1$ . After imposing the consistency requirement that at each step the $\text{A}$ and $\text{B}$ flows act during $h$ units of time each, the family has $s-1$ parameters left. By taking $s$ sufficiently high it is possible to achieve any desired order (Sanz-Serna and Calvo Reference Sanz-Serna and Calvo1994, Section 13.1). Clearly, increasing the number of stages does not lead to integrators with more complicated implementation; regardless of the value of $s$ , numerical integrations only consist of a sequence of flows of the split systems.

It is not necessary to add that to each of the integrators just described there corresponds a second integrator found by swapping the roles of the systems $\text{A}$ and $\text{B}$ , just as (3.15) corresponds to (3.13).

So far our attention has focused on palindromic formulas as these will be important later in connection with the property of reversibility (Theorem 4.2). It is also possible to consider splittings of the form

$$\begin{eqnarray}\unicode[STIX]{x1D713}_{h}=\unicode[STIX]{x1D711}_{b_{r}h}^{(\text{B})}\circ \unicode[STIX]{x1D711}_{a_{r}h}^{(\text{A})}\circ \cdots \circ \unicode[STIX]{x1D711}_{b_{2}h}^{(\text{B})}\circ \unicode[STIX]{x1D711}_{a_{2}h}^{(\text{A})}\circ \unicode[STIX]{x1D711}_{b_{1}h}^{(\text{B})}\circ \unicode[STIX]{x1D711}_{a_{1}h}^{(\text{A})},\end{eqnarray}$$

which we denote as

(3.21) $$\begin{eqnarray}(b_{r},a_{r},\ldots ,b_{2},a_{2},b_{1},a_{1})\end{eqnarray}$$

( $\sum _{i}a_{i}=\sum _{i}b_{i}=1$ ). The palindromic formulas in (3.19)–(3.20) are particular instances of this general format because it is always possible to set $a_{1}=0$ in (3.21).

Remark 3.5. When a splitting integrator of the general form (3.21) is applied to the splitting (3.9)–(3.10) of the system (2.2), the result coincides with the application of a so-called symplectic, explicit Runge–Kutta–Nyström (RKN) method. The properties (order, stability, etc.) of such an integrator may therefore be studied either by using techniques pertaining to splitting integrators (as is done in this paper) or by employing an RKN approach. The second methodology was favoured in the early years of geometric integration (Sanz-Serna Reference Sanz-Serna1992, Sanz-Serna and Calvo Reference Sanz-Serna and Calvo1994).

3.4 Fixed $h$ stability

If two numerical schemes are candidates to integrate an initial value problem for a given system (2.1), then the scheme that leads to smaller global errors for a given computational cost may seem more desirable (but, in the context of geometric integration, the geometric properties may play a role when choosing the integrator: see Section 4 below). Even though global errors may be bounded as in (3.6), in practice it is almost always impossible, for the problem at hand, to estimate realistically the error constant implied in the $O(h^{\unicode[STIX]{x1D708}})$ notation in the bound. For this reason, the literature on numerical integrators has traditionally resorted to well-chosen model problems where both the numerical and true solutions, $x_{n}$ and $x(t_{n})$ , may be written down in closed form; the performance of the various integrators on the model problem may then be investigated analytically and is taken as an indication of their performance when applied to realistic problems. Note that the actual numerical solution of the model problem cannot be expected to be of real practical interest, since in that problem the true solution is available analytically.

3.4.1 One degree of freedom

For our purposes, the model problem of choice is the harmonic oscillator (2.3). For Runge–Kutta methods, for the splitting algorithms in Section 3.3 (and in fact for all one-step integrators of practical interest) a time-step $(q_{n+1},p_{n+1})=\unicode[STIX]{x1D713}_{h}(q_{n},p_{n})$ may be expressed as

(3.22) $$\begin{eqnarray}\left[\begin{array}{@{}c@{}}q_{n+1}\\ p_{n+1}\end{array}\right]=\tilde{M}_{h}\left[\begin{array}{@{}c@{}}q_{n}\\ p_{n}\end{array}\right]\!,\quad \tilde{M}_{h}=\left[\begin{array}{@{}cc@{}}A_{h} & B_{h}\\ C_{h} & D_{h}\end{array}\right]\end{eqnarray}$$

for suitable method-dependent coefficients $A_{h}$ , $B_{h}$ , $C_{h}$ , $D_{h}$ . The evolution over $n$ time-steps is then given by

(3.23) $$\begin{eqnarray}\left[\begin{array}{@{}c@{}}q_{n}\\ p_{n}\end{array}\right]=\tilde{M}_{h}^{n}\left[\begin{array}{@{}c@{}}q_{0}\\ p_{0}\end{array}\right]\!,\end{eqnarray}$$

an expression to be compared with (2.4).

If a given $h>0$ is such that $|\tilde{M}_{h}^{n}|\rightarrow \infty$ as $n\rightarrow \infty$ , the magnitude of numerical solution $(q_{n},p_{n})$ will grow unboundedly, while the true solution of course remains bounded as $t\rightarrow \infty$ . Necessarily, global errors will be large for large $n$ . Then the integrator is said to be unstable for that particular choice of $h$ .

Example 3.6. For Euler’s rule (3.2) we find

$$\begin{eqnarray}A_{h}=D_{h}=1,\quad B_{h}=-C_{h}=h.\end{eqnarray}$$

The eigenvalues of $\tilde{M}_{h}$ are $1\pm \text{i}h$ with modulus $(1+h^{2})^{1/2}$ and, therefore, for any fixed $h>0$ , the powers $\tilde{M}_{h}^{n}$ grow exponentially as $n$ increases. Due to this numerical instability, Euler’s rule is completely unsuitable for integrating the harmonic oscillator, and this rules it out as a method to integrate more complicated oscillatory problems.

More precisely, in the step $n\rightarrow n+1$ , the radius $r=(q^{2}+p^{2})^{1/2}$ , which remains constant for the true solution, grows like

(3.24) $$\begin{eqnarray}r_{n+1}=(1+h^{2})^{1/2}r_{n}=\biggl(1+\displaystyle \frac{h^{2}}{2}+O(h^{4})\biggr)r_{n}\end{eqnarray}$$

for the Euler solution, so that $r_{n}=(1+h^{2})^{n/2}r_{0}$ . Note that, taking limits as $n\rightarrow \infty$ , $h\rightarrow 0$ , with fixed $nh$ , we find $r_{n}\rightarrow r_{0}=r(0)=r(nh)$ , as it corresponds to a convergent method.

Example 3.7. The midpoint rule (3.3) has

$$\begin{eqnarray}A_{h}=D_{h}=\displaystyle \frac{1-h^{2}/4}{1+h^{2}/4},\quad B_{h}=-C_{h}=\displaystyle \frac{h}{1+h^{2}/4}.\end{eqnarray}$$

The characteristic equation of $\tilde{M}_{h}$ is

(3.25) $$\begin{eqnarray}\unicode[STIX]{x1D706}^{2}-2A_{h}\unicode[STIX]{x1D706}+1,\end{eqnarray}$$

and therefore the product of the eigenvalues is 1. For $h\neq 0$ , $|A_{h}|<1$ and the matrix has a pair of complex conjugate eigenvalues of unit modulus. Then the powers $\tilde{M}_{h}^{n}$ remain bounded as $n$ increases and the method is stable for any $h$ .

Example 3.8. Let us next consider Strang’s splitting (3.13), applied with the splitting (3.9)–(3.10), which yields the velocity Verlet algorithm. This has

$$\begin{eqnarray}A_{h}=D_{h}=1-h^{2}/2,\quad B_{h}=h,\quad C_{h}=-h+h^{3}/4;\end{eqnarray}$$

the characteristic equation is again of the form (3.25). For $h>2$ , $A_{h}<-1$ and the eigenvalues are real and distinct, so that one of them has modulus ${>}1$ , and therefore the powers $\tilde{M}_{h}^{n}$ grow exponentially. For $h<2$ the eigenvalues are complex of unit modulus and the powers $\tilde{M}_{h}^{n}$ remain bounded. For $h=2$ , $\tilde{M}_{h}$ is a non-trivial Jordan block whose powers grow linearly (weak instability). To summarize, the integrator is unstable for $h\geq 2$ and has the stability interval $0<h<2$ . Integrations with $h>2$ lead to extremely large global errors, as we shall see below.

Example 3.9. For the alternative Strang formula (3.15), which yields the position Verlet algorithm, the coefficients are

$$\begin{eqnarray}A_{h}=D_{h}=1-h^{2}/2,\quad B_{h}=h-h^{3}/4,\quad C_{h}=-h.\end{eqnarray}$$

The discussion is almost identical to the one in the preceding example. The stability interval is also $0<h<2$ , as one may have guessed from the equal role that $q$ and $-p$ play in the harmonic oscillator (2.3).

Table 3.1.

Velocity Verlet integration of the harmonic oscillator. Relative errors after one or ten oscillation periods.

Example 3.10. In the situation of Example 3.8, are values of $h$ below the upper limit 2 satisfactory? The answer of course depends on the accuracy required. For the initial condition $q=1$ , $p=0$ , and different stable values of $h$ , Table 3.1 gives the relative error in the Euclidean norm

$$\begin{eqnarray}\displaystyle \frac{|(q_{n}-q(t_{n}),p_{n}-p(t_{n}))|}{|(q(t_{n}),p(t_{n}))|}\end{eqnarray}$$

at the final integration time, when the integration is carried out over an interval of length either one oscillation period (second column) or ten oscillation periods (third column). Columns are consistent with the order of convergence as in (3.6). The last rows reveal that the error increases linearly with $t$ (this will not be true when integrating other differential systems). In the first row the error grows more slowly than $t$ : the numerical solution $(q_{n},p_{n})$ remains close to the unit circle for all values of $n$ and therefore errors cannot be substantially larger than the diameter of the circle. From the table we see that if we are interested in errors below $10\%$ over an interval of length equal to ten oscillation periods, then $h$ has to be taken below $2\unicode[STIX]{x1D70B}/32\approx 0.20$ , that is, well below the upper end, $h=2$ , of the stability interval. To provide an indication of the effect of using unstable values of $h$ , we mention that with $h=\unicode[STIX]{x1D70B}$ the error after one period is ${\approx}46.4$ and after ten periods ${\approx}4.68\times 10^{17}$ .

Before we move to more complicated models, let us make two observations, valid for Runge–Kutta, splitting integrators and any other method of practical interest.

Remark 3.11. Because the problem is linear and rotationally invariant, the magnitude of the relative errors is independent of the initial condition $(q_{0},p_{0})\neq (0,0)$ .

Remark 3.12. Replacing the model (2.3) with the apparently more general system

(3.26) $$\begin{eqnarray}\displaystyle \frac{\text{d}q}{\text{d}t}=\unicode[STIX]{x1D714}p,\quad \displaystyle \frac{\text{d}p}{\text{d}t}=-\unicode[STIX]{x1D714}q,\quad \unicode[STIX]{x1D714}>0,\end{eqnarray}$$

with oscillation period $2\unicode[STIX]{x1D70B}/\unicode[STIX]{x1D714}$ , does not really change things: in view of Remark 3.1, integrating (3.26) with step size $h$ is equivalent to integrating (2.3) with step size $h/\unicode[STIX]{x1D714}$ . Because the length of the integration interval and the step size in Table 3.1 are given in terms of the oscillation period, the results displayed are valid for (3.26) for any choice of the value of $\unicode[STIX]{x1D714}$ . Regardless of the initial condition, if we are interested in relative errors below $10\%$ over an interval of length $20\unicode[STIX]{x1D70B}/\unicode[STIX]{x1D714}$ (equal to ten oscillation periods), then $h$ has to be taken below $2\unicode[STIX]{x1D70B}/(32\unicode[STIX]{x1D714})$ . Of course, for (3.26), stability requires that $h<2/\unicode[STIX]{x1D714}$ .

3.4.2 Several degrees of freedom; stability restrictions on $h$

Let us now move to the model with $d$ degrees of freedom,

$$\begin{eqnarray}\displaystyle \frac{\text{d}}{\text{d}t}q=M^{-1}p,\quad \displaystyle \frac{\text{d}}{\text{d}t}p=-Kq,\end{eqnarray}$$

where $M$ and $K$ are $d\times d$ , symmetric positive-definite matrices. This is the Hamiltonian system with

(3.27) $$\begin{eqnarray}H=\displaystyle \frac{1}{2}p^{T}M^{-1}p+\displaystyle \frac{1}{2}q^{T}Kq.\end{eqnarray}$$

In mechanics, $M$ and $K$ are the mass and stiffness matrices respectively.

This model may transformed into $d$ uncoupled one-degree-of-freedom oscillators. In fact, factor $M=LL^{T}$ (which may be done in infinitely many ways) and diagonalize the symmetric positive-definite matrix $L^{-1}KL^{-T}$ as $U^{T}L^{-1}KL^{-T}U=\unicode[STIX]{x1D6FA}^{2}$ , with $U$ orthogonal and $\unicode[STIX]{x1D6FA}$ diagonal with diagonal entries $\unicode[STIX]{x1D714}_{i}$ , $i=1,\ldots ,d$ . A simple computation yields the next result.

Proposition 3.13. With the notation above, the (non-canonical) change of dependent variables, $q=L^{-T}U\bar{q}$ , $p=LU\unicode[STIX]{x1D6FA}\bar{p}$ , decouples the system into a collection of $d$ harmonic oscillators:

$$\begin{eqnarray}\displaystyle \frac{\text{d}\bar{q}^{i}}{\text{d}t}=\unicode[STIX]{x1D714}_{i}\bar{p}^{i},\quad \displaystyle \frac{\text{d}\bar{p}^{i}}{\text{d}t}=-\unicode[STIX]{x1D714}_{i}\bar{q}^{i},\quad i=1,\ldots ,d\end{eqnarray}$$

(superscripts denote components).

Now, for all integrators of practical interest, decoupling and numerical integration commute: carrying out the integration in the old variables $(q,p)$ yields the same result as successively (i) changing variables in the system, (ii) integrating each of the uncoupled oscillators, and (iii) translating the result to the old variables. Therefore stability may be analysed under the assumption that the integration is performed in the uncoupled version.

Example 3.14. Consider a particular case with $d=2$ , $M=L=I$ , $\unicode[STIX]{x1D714}_{1}=1$ and $\unicode[STIX]{x1D714}_{2}=100$ and the initial condition

$$\begin{eqnarray}\bar{q}^{1}(0)=1,\quad \bar{p}^{1}(0)=0,\quad \bar{q}^{2}(0)=0.01,\quad \bar{p}^{2}(0)=0.\end{eqnarray}$$

We wish to integrate with the velocity Verlet algorithm over $0\leq t\leq 20\unicode[STIX]{x1D70B}$ (ten periods of the slower oscillation) and aim at absolute errors of magnitude ${\approx}0.1$ in the Euclidean norm in the $\mathbb{R}^{4}$ space of the variables $(q^{1},q^{2},p^{1},p^{2})$ or, equivalently, because here $q=U\bar{q}$ with $U$ orthogonal, in the Euclidean norm of the variables $(\bar{q}^{1},\bar{q}^{2},\bar{p}^{1},\bar{p}^{2})$ . On accuracy grounds it would be sufficient to take $h\approx 0.2$ : from Example 3.10 we know this is enough to integrate the first uncoupled oscillator with the desired accuracy, and the second oscillator should not contribute significantly to the error, due to the smallness of $\bar{q}^{2}(t)$ and $\bar{p}^{2}(t)$ for all $t$ . However, unless we take $h\unicode[STIX]{x1D714}_{2}<2$ , i.e.  $h<0.02$ , the errors in $\bar{q}^{2}(t)$ and $\bar{p}^{2}(t)$ will grow exponentially due to instability. In this example, and in many situations arising in practice, to avoid instabilities, the value of $h$ has to be chosen much smaller than accuracy would require. As a result the computational effort to span the time interval of interest would be much higher than might be expected on accuracy grounds. These situations are called stiff. To deal with stiffness one may resort to suitable implicit integrators such as the midpoint rule, or explicit integrators with large stability intervals.

Remark 3.15. It may seem that the initial condition in the preceding example is somewhat contrived, as the sizes of two components of $\bar{q}$ are so unbalanced. This is not so: it is easily checked that, in terms of the original energy in (3.27), the contributions to $H$ of the $\unicode[STIX]{x1D714}_{1}$ and $\unicode[STIX]{x1D714}_{2}$ oscillations are both equal to $1/2$ . Those oscillations are called normal modes in mechanics.

The present discussion is continued in Section 4.4.

We emphasize that, while the preceding material refers to the quadratic Hamiltonian (3.27), one may expect that it has some relevance for Hamiltonians close to that model. However, there are cases where the behaviour of a numerical integrator departs considerably from its behaviour when applied to (3.27). This point is discussed in Section 9.2.

4.1 Hamiltonian problems

Splitting integrators are symplectic in the following sense. Assume that the system (3.7) is Hamiltonian and that it is split in such a way that both split systems (3.8) are also Hamiltonian. Then the splitting integrator mapping $\unicode[STIX]{x1D713}_{h}$ is symplectic, as a composition (Proposition 2.3) of flows that are individually symplectic (Theorem 2.5). This ensures that the numerical solution shares the specific properties of Hamiltonian flows that derive from symplecticness. Note that for $\unicode[STIX]{x1D713}_{h}$ to be symplectic it is not enough that the system being integrated is Hamiltonian; if the split vector fields $f^{(\text{A})}$ and $f^{(\text{B})}$ are not Hamiltonian themselves, then $\unicode[STIX]{x1D713}_{h}$ cannot be expected to be symplectic. Of course, all possible splittings of a Hamiltonian vector field as a sum of Hamiltonian vector fields, $f=f^{(\text{A})}+f^{(\text{B})}$ , may be obtained by splitting the Hamiltonian $H=H^{(\text{A})}+H^{(\text{B})}$ in an arbitrary way and then setting $f^{(\text{A})}=J^{-1}\unicode[STIX]{x1D6FB}H^{(\text{A})}$ , $f^{(\text{B})}=J^{-1}\unicode[STIX]{x1D6FB}H^{(\text{B})}$ . We sum up as follows.

Note that the $n$ -fold composition $\unicode[STIX]{x1D713}_{h}^{n}$ that advances the numerical solution over $n$ time-steps is then also symplectic (Proposition 2.3).

Some implicit Runge–Kutta methods, including the midpoint rule and Gauss methods, are also symplectic: $\unicode[STIX]{x1D713}_{h}$ is a symplectic map whenever the system (2.1) being integrated is Hamiltonian as in (2.8). No explicit Runge–Kutta method is symplectic. In particular, Euler’s rule is not symplectic; according to Example 3.6 the Euler discretization of the harmonic oscillator increases area, as $\tilde{M}_{h}$ has determinant $1+h^{2}$ . This increase in area is related to the estimate (3.24) that shows that the Euler solution spirals outward, a non-Hamiltonian behaviour.

In the particular case of the Hamiltonian in (2.10), taking ${\mathcal{T}}$ and ${\mathcal{U}}$ to play the roles of $H^{(\text{A})}$ and $H^{(\text{B})}$ respectively leads to the splitting of the differential system given in (3.9)–(3.10) with $F=-\unicode[STIX]{x1D6FB}{\mathcal{U}}$ .

It would also be desirable to have integrators that preserved energy when applied to the Hamiltonian system (2.8), that is, $H\circ \unicode[STIX]{x1D713}_{h}=H$ . Unfortunately, for realistic problems such a requirement is incompatible with $\unicode[STIX]{x1D713}_{h}$ being symplectic (Sanz-Serna and Calvo Reference Sanz-Serna and Calvo1994, Section 10.3). It is then standard practice to insist on symplecticness and sacrifice conservation of energy. There are several reasons for this. Symplecticness plays a key role in the Hamiltonian formalism (see Theorem 2.5). In addition, while, as we have seen, it is not difficult to find symplectic formulas, standard classes of integrators do not include energy-preserving schemes unless the energy is assumed to have particular forms. Finally, as we shall discuss below, symplectic schemes have small energy errors even when the integration interval is very long.

4.2 Reversible problems

For reversible systems we have, from Proposition 2.9 and Theorem 2.10, the following theorem.

The midpoint rule and Gauss Runge–Kutta formulas also generate reversible mappings $\unicode[STIX]{x1D713}_{h}$ whenever they are applied to a reversible system. No explicit Runge–Kutta method is reversible.

The $n$ -fold composition $\unicode[STIX]{x1D713}_{h}^{n}$ that advances the numerical solution over $n$ time-steps is then also reversible (Proposition 2.9). However, note that if variable time-steps were taken, then the mapping $\unicode[STIX]{x1D713}_{h_{n}}\circ \cdots \circ \unicode[STIX]{x1D713}_{h_{1}}$ , which advances the solution from $t_{0}$ to $t_{n+1}$ , would not be reversible. This is one of the reasons for not considering variable time-steps here; see Remark 3.4.

The use of reversible integrators (with constant step sizes) ensures that the numerical solution inherits relevant geometric properties of the true solution of the differential system (Cano and Sanz-Serna Reference Cano and Sanz-Serna1997, Cano and Sanz-Serna Reference Cano and Sanz-Serna1998).

4.3 Modified equations

Modified equations are rather old (see references in Griffiths and Sanz-Serna Reference Griffiths and Sanz-Serna1986); however, their use as a means to analyse numerical integrators has only become popular in the past twenty years, since the emergence of geometric integration (Sanz-Serna Reference Sanz-Serna and Marsden1996).

4.3.1 Motivation

Let us first look at some examples.

Example 4.3. Suppose that the system (2.1) is solved with Euler’s rule (3.2) and let $\unicode[STIX]{x1D713}_{h}$ denote the corresponding map $x+hf(x)$ . Consider the new system, parametrized by $h$ ,

(4.1) $$\begin{eqnarray}\displaystyle \frac{\text{d}}{\text{d}t}x=\tilde{f}_{h}(x),\quad \tilde{f}_{h}(x)=f(x)-\displaystyle \frac{h}{2}f^{\prime }(x)f(x),\end{eqnarray}$$

with flow $\tilde{\unicode[STIX]{x1D711}}_{t}$ (for simplicity, the dependence of this flow on the parameter $h$ is not incorporated in the notation). By proceeding as in the derivation of (3.5), we find that the Taylor expansion of $\tilde{\unicode[STIX]{x1D711}}_{h}$ in powers of $h$ is

$$\begin{eqnarray}\displaystyle \tilde{\unicode[STIX]{x1D711}}_{h}(x) & = & \displaystyle x+h\tilde{f}_{h}(x)+\displaystyle \frac{h^{2}}{2}\tilde{f}_{h}^{\prime }(x)\tilde{f}_{h}(x)+O(h^{3})\nonumber\\ \displaystyle & = & \displaystyle x+h\biggl(f(x)-\displaystyle \frac{h}{2}f^{\prime }(x)f(x)\biggr)+\displaystyle \frac{h^{2}}{2}f_{h}^{\prime }(x)f_{h}(x)+O(h^{3})\nonumber\\ \displaystyle & = & \displaystyle x+hf(x)+O(h^{3}).\nonumber\end{eqnarray}$$

Thus, the Euler mapping $\unicode[STIX]{x1D713}_{h}(x)$ , which differs from the flow of the system (2.1) being solved in $O(h^{2})$ terms (first-order consistency), differs from the flow $\tilde{\unicode[STIX]{x1D711}}_{h}$ of the so-called modified or shadow system (4.1) in $O(h^{3})$ terms (second-order consistency). According to Theorem 3.2 (see Griffiths and Sanz-Serna Reference Griffiths and Sanz-Serna1986 for details), over bounded time intervals, the Euler solution differs from the corresponding solution of the modified system in $O(h^{2})$ terms. Therefore we may expect that the properties of the Euler solutions for (2.1) will resemble the properties of the solutions of (4.1) more than they resemble the properties of solutions of (2.1) itself. When studying the properties of Euler’s rule, working with (4.1) rather than with $\unicode[STIX]{x1D713}_{h}$ should be advantageous, because differential equations are simpler to analyse than maps. An illustration is provided next.

Example 4.4. Let us particularize the preceding example to the case of the harmonic oscillator (2.3). The modified system (4.1) is found to be

$$\begin{eqnarray}\displaystyle \frac{\text{d}q}{\text{d}t}=p+\displaystyle \frac{h}{2}q,\quad \displaystyle \frac{\text{d}p}{\text{d}t}=-q+\displaystyle \frac{h}{2}p,\end{eqnarray}$$

and from here a simple computation yields for the radius $r=(q^{2}+p^{2})^{1/2}$

$$\begin{eqnarray}\displaystyle \frac{\text{d}r}{\text{d}t}=\displaystyle \frac{h}{2}r,\end{eqnarray}$$

so that, over a time interval of length $h$ , $r$ is multiplied by the factor $\exp (h^{2}/2)=1+h^{2}/2+O(h^{3})$ . This is precisely what we found in (3.24) for the Euler solution of the harmonic oscillator (but the $O(h^{3})$ remainder here will not coincide with the one in (3.24), because there is an $O(h^{3})$ difference between $\unicode[STIX]{x1D713}_{h}$ and $\tilde{\unicode[STIX]{x1D711}}_{h}$ ).

Example 4.5. In (4.1), $\tilde{f}_{h}(x)$ is a first-degree polynomial in $h$ . If $\tilde{f}_{h}(x)$ is chosen to be quadratic in $h$ , that is, $\tilde{f}_{h}=f-(h/2)f^{\prime }f+h^{2}f_{2}$ , it is possible to determine $f_{2}$ so as to achieve $\unicode[STIX]{x1D713}_{h}-\tilde{\unicode[STIX]{x1D711}}_{h}=O(h^{4})$ for the Euler map $\unicode[STIX]{x1D713}_{h}(x)=x+hf(x)$ . Similarly, taking $\tilde{f}_{h}(x)$ as a suitable chosen polynomial in $h$ with degree $\unicode[STIX]{x1D707}=3,4,\ldots ,$ it is possible to achieve $\unicode[STIX]{x1D713}_{h}-\tilde{\unicode[STIX]{x1D711}}_{h}=O(h^{\unicode[STIX]{x1D707}+2})$ .

4.3.2 Definition

Given a (consistent) integrator $\unicode[STIX]{x1D713}_{h}$ for system (2.1), there exists a (unique) formal series in powers of $h$ ,

(4.2) $$\begin{eqnarray}\tilde{f}_{h}^{\infty }(x)=f(x)+hg^{[1]}(x)+h^{2}g^{[2]}(x)+\cdots\end{eqnarray}$$

(the $g^{[\unicode[STIX]{x1D707}]}$ map the space $\mathbb{R}^{D}$ into itself), with the property that, for each $\unicode[STIX]{x1D707}=0,1,\ldots ,$ the flow $\tilde{\unicode[STIX]{x1D711}}_{h}^{[\unicode[STIX]{x1D707}]}$ of the modified system

(4.3) $$\begin{eqnarray}\displaystyle \frac{\text{d}}{\text{d}t}x=\tilde{f}_{h}^{[\unicode[STIX]{x1D707}]}(x),\quad \tilde{f}_{h}^{[\unicode[STIX]{x1D707}]}(x)=f(x)+hg^{[1]}(x)+\cdots +h^{\unicode[STIX]{x1D707}}g^{[\unicode[STIX]{x1D707}]}(x),\end{eqnarray}$$

satisfies

(4.4) $$\begin{eqnarray}\unicode[STIX]{x1D713}_{h}(x)-\tilde{\unicode[STIX]{x1D711}}_{h}^{[\unicode[STIX]{x1D707}]}(x)=O(h^{\unicode[STIX]{x1D707}+2}).\end{eqnarray}$$

Furthermore, for a symmetric method, the odd-numbered $g^{[\unicode[STIX]{x1D707}]}(x)$ are identically zero:

(4.5) $$\begin{eqnarray}g^{[1]}(x)=g^{[3]}(x)=g^{[5]}(x)=\cdots =0.\end{eqnarray}$$

As $\unicode[STIX]{x1D707}$ increases, the solutions of the modified system (4.3) provide better and better approximations to $\unicode[STIX]{x1D713}_{h}$ . For some vector fields $f$ and some integrators (see an example in Proposition 4.9 below), it may happen that the formal series (4.2) converges for each $x$ (so that $\tilde{f}_{h}^{\infty }$ is a well-defined mapping $\mathbb{R}^{D}:\rightarrow \mathbb{R}^{D}$ ) and that furthermore the flow $\tilde{\unicode[STIX]{x1D711}}_{h}^{[\infty ]}$ of

(4.6) $$\begin{eqnarray}\displaystyle \frac{\text{d}}{\text{d}t}x=\tilde{f}_{h}^{\infty }(x)\end{eqnarray}$$

coincides with $\unicode[STIX]{x1D713}_{h}$ . In those cases, (4.6) provides an exact modified system to study $\unicode[STIX]{x1D713}_{h}$ . However, such situations are exceptional, because, in general, discrete dynamical systems (such as the one generated by $\unicode[STIX]{x1D713}_{h}$ ) possess features that cannot appear in flows of autonomous differential systems. It is the lack of convergence of (4.2) that makes it necessary to consider the truncations in (4.3) which do not exactly reproduce $\unicode[STIX]{x1D713}_{h}$ but approximate it with an error as in (4.4). It is possible, by increasing $\unicode[STIX]{x1D707}$ as $h$ becomes smaller and smaller, to render the discrepancy $\unicode[STIX]{x1D713}_{h}(x)-\tilde{\unicode[STIX]{x1D711}}_{h}^{[\unicode[STIX]{x1D707}(h)]}(x)$ exponentially small with respect to $h$ , as first proved by Neishtadt (Sanz-Serna and Calvo Reference Sanz-Serna and Calvo1994, Section 10.1).

4.3.3 Finding explicitly the modified equations

For splitting integrators, the terms of the series (4.2) may be found explicitly by using the Baker–Campbell–Hausdorff formula. The next theorem provides a summary of some key points; a more detailed description is given in Section 9.1.

Theorem 4.6. Assume that system (3.7) is integrated by means of a (consistent) splitting algorithm of the general format (3.21). The series (4.2) is of the form

$$\begin{eqnarray}\displaystyle & & \displaystyle (f^{(\text{A})}+f^{(\text{B})})+hC_{1,1}[f^{(\text{A})},f^{(\text{B})}]\nonumber\\ \displaystyle & & \displaystyle \qquad +h^{2}(C_{2,1}[f^{(\text{A})},[f^{(\text{A})},f^{(\text{B})}]]+C_{2,2}[f^{(\text{B})},[f^{(\text{A})},f^{(\text{B})}]])+\cdots \nonumber\end{eqnarray}$$

where, for each $\unicode[STIX]{x1D707}=1,2\ldots ,$ the coefficient of $h^{\unicode[STIX]{x1D707}}$ is a linear combination of linearly independent iterated commutators involving $\unicode[STIX]{x1D707}+1$ fields $f^{(\text{A})}$ , $f^{(\text{B})}$ . The coefficients $C_{i,j}$ that appear in the linear combinations are known polynomials on the coefficients $a_{k}$ and $b_{\ell }$ that appear in the formula (3.21).

If the splitting is palindromic then all the coefficients $C_{2i+1,j}$ corresponding to odd powers of $h$ vanish.

For Runge–Kutta methods it is also possible to give expressions for the terms in the series (4.2): see Calvo, Murua and Sanz-Serna (Reference Calvo, Murua, Sanz-Serna, Kloeden and Palmer1994) or Sanz-Serna (Reference Sanz-Serna and Marsden1996).

4.3.4 Modified equations and order

If in (4.2) the functions $g^{[1]},\ldots ,g^{[\unicode[STIX]{x1D707}-1]},\unicode[STIX]{x1D707}\geq 2$ , vanish so that

$$\begin{eqnarray}\tilde{f}_{h}^{[\unicode[STIX]{x1D707}]}(x)-f(x)=h^{\unicode[STIX]{x1D707}}g^{[\unicode[STIX]{x1D707}]}(x)+O(h^{\unicode[STIX]{x1D707}+1}),\end{eqnarray}$$

then the flow $\tilde{\unicode[STIX]{x1D711}}_{h}^{[\unicode[STIX]{x1D707}]}$ of (4.3) and the flow $\unicode[STIX]{x1D711}_{h}$ of (2.1) will satisfy

$$\begin{eqnarray}\tilde{\unicode[STIX]{x1D711}}_{h}^{[\unicode[STIX]{x1D707}]}(x)-\unicode[STIX]{x1D711}_{h}(x)=h^{\unicode[STIX]{x1D707}+1}g^{[\unicode[STIX]{x1D707}]}(x)+O(h^{\unicode[STIX]{x1D707}+2}),\end{eqnarray}$$

which in view of (4.4) implies that

$$\begin{eqnarray}\unicode[STIX]{x1D713}_{h}(x)-\unicode[STIX]{x1D711}_{h}(x)=h^{\unicode[STIX]{x1D707}+1}g^{[\unicode[STIX]{x1D707}]}(x)+O(h^{\unicode[STIX]{x1D707}+2}),\end{eqnarray}$$

that is, the integrator $\unicode[STIX]{x1D713}_{h}(x)$ has order $\unicode[STIX]{x1D707}$ (or higher). The converse is also true, because the argument may be reversed: order $\geq \unicode[STIX]{x1D707}$ implies that $g^{[1]},\ldots ,g^{[\unicode[STIX]{x1D707}-1]}$ must vanish. From (4.5), the order of a symmetric method must be an even integer. To sum up, we have the following result, which makes it possible to use the series (4.2) rather than the mapping $\unicode[STIX]{x1D713}_{h}$ to study the order of a given integrator.

Theorem 4.7. A (consistent) integrator has order $\geq \unicode[STIX]{x1D708}$ , $\unicode[STIX]{x1D708}\geq 2$ , if and only if the functions $g^{[1]},\ldots ,g^{[\unicode[STIX]{x1D708}-1]}$ appearing in (4.2) are identically zero.

When the order is exactly $\unicode[STIX]{x1D708}$ , $\unicode[STIX]{x1D708}=1,2,\ldots ,$ the leading term of the truncation error is $h^{\unicode[STIX]{x1D708}+1}g^{[\unicode[STIX]{x1D708}]}(x)$ .

A (consistent) symmetric integrator has even order.

This theorem, in tandem with Theorem 4.6, provides the standard way to write down the order conditions for splitting methods, that is, the relations on the coefficients $a_{k}$ , $b_{\ell }$ that are necessary and sufficient for an integrator to have order $\unicode[STIX]{x1D708}$ . Specifically, order $\geq \unicode[STIX]{x1D708}$ is equivalent to $C_{i,j}=0$ whenever $i<\unicode[STIX]{x1D708}$ . Here are the simplest illustrations.

• ∙ The Lie–Trotter formula (3.11) has $\unicode[STIX]{x1D708}=1$ , and, according to Theorems 4.6 and 4.7, the leading term of the local error is a constant multiple of $h^{2}[f^{(\text{A})},f^{(\text{B})}]$ , which matches our earlier finding in (3.12).

• ∙ The condition $C_{1,1}=0$ is necessary and sufficient for an integrator to have order $\unicode[STIX]{x1D708}\geq 2$ . It is automatically satisfied for palindromic formulas.

• ∙ Strang’s formula (3.13) has order exactly 2. The coefficient of $h^{3}$ in the expansion of the local error is a combination of iterated commutators with three fields. This was already found in (3.14).

• ∙ To have order $\unicode[STIX]{x1D708}\geq 3$ we have to impose the order conditions $C_{1,1}=0$ , $C_{2,1}=0$ , $C_{2,2}=0$ .

4.3.5 Modified equations and geometric integrators

The geometric properties of integrators have a clear impact on their modified systems.

Let us begin with the Hamiltonian case. If the split fields in (3.8) are Hamiltonian with Hamiltonian functions $H^{(\text{A})}$ and $H^{(\text{B})}$ , then by invoking Theorem 2.13 it is feasible to work with iterated Poisson brackets of $H^{(\text{A})}$ and $H^{(\text{B})}$ rather than with the iterated Lie brackets of $f^{(\text{A})}$ and $f^{(\text{B})}$ . Then the expansion in Theorem 4.6 is replaced by

$$\begin{eqnarray}\displaystyle & & \displaystyle (H^{(\text{A})}+H^{(\text{B})})-hC_{1,1}\{H^{(\text{A})},H^{(\text{B})}\}\nonumber\\ \displaystyle & & \displaystyle \qquad +h^{2}(C_{2,1}\{H^{(\text{A})},\{H^{(\text{A})},H^{(\text{B})}\}\}+C_{2,2}\{H^{(\text{B})},\{H^{(\text{A})},H^{(\text{B})}\}\})+\cdots \,.\nonumber\end{eqnarray}$$

For each $\unicode[STIX]{x1D707}$ , the modified system (4.3) is a Hamiltonian system. This is a reflection of the fact that, according to Theorem 4.1, splitting integrators give rise to mappings $\unicode[STIX]{x1D713}_{h}$ that are symplectic. In the class of Runge–Kutta methods it is also true that symplectic integrators have modified systems that are Hamiltonian.

Speaking informally, we may say that all integrators change the system being integrated into a modified system; in non-symplectic methods the (perhaps small) change is such that the modified system is no longer Hamiltonian. Symplectic methods are those that change Hamiltonian systems into Hamiltonian systems. This heuristic description cannot be made entirely rigorous because, as pointed out above, the exact modified system (4.6) only exists in a formal sense due to the lack of convergence of the series in (4.2). The existence of a Hamiltonian modified system is at the basis of many favourable properties of symplectic integrators.

Remark 4.8. The considerations above only make sense if the step size $h$ is held constant along the integration interval. Since the modified system changes with $h$ , variable step size implementations of symplectic integrators do not have well-defined modified systems and in fact their behaviour is closer to that of non-symplectic integrators than to that of symplectic integrators used with constant step sizes (Calvo and Sanz-Serna Reference Calvo and Sanz-Serna1993).

For the reversible case, consider the situation of Theorem 4.2. From Proposition 2.14 all iterated commutators involving an odd number of fields are themselves $S$ -reversible. Since, according to Theorem 4.6, for palindromic splitting integrators the iterated commutators with an even number of fields enter the expansion (4.2) with null coefficients, then the modified systems (4.3) are reversible.

4.4 Geometric integrators and the harmonic model problem

We take again the integration of the harmonic oscillator as a model problem (see (3.22)–(3.23)); now our interest is in studying in detail the behaviour of geometric integrators.

We focus on (consistent) integrators that are both symplectic and reversible. In terms of the matrix $\tilde{M}_{h}$ , the first of these properties corresponds to $A_{h}D_{h}-B_{h}C_{h}=1$ and, when this condition holds, reversibility is equivalent to $A_{h}=D_{h}$ . Our treatment follows that of Blanes et al. (Reference Blanes, Casas and Sanz-Serna2014).

The characteristic polynomial of $\tilde{M}_{h}$ is of the form (3.25) and there are four possibilities, the first two corresponding to unstable simulations and the last two corresponding to stable simulations.

1. (1) $h$ is such that $|A_{h}|>1$ . In that case $\tilde{M}_{h}$ has spectral radius ${>}1$ and therefore the powers $\tilde{M}_{h}^{n}$ grow exponentially with $n$ .

2. (2) $A_{h}=\pm 1$ and $|B_{h}|+|C_{h}|>0$ . The powers $\tilde{M}_{h}^{n}$ grow linearly with $n$ .

3. (3) $A_{h}=\pm 1$ , $B_{h}=C_{h}=0$ , that is, $\tilde{M}_{h}=\pm I$ , $\tilde{M}_{h}^{n}=(\pm I)^{n}$ .

4. (4) $h$ is such that $|A_{h}|<1$ . In that case, $\tilde{M}_{h}$ has complex conjugate eigenvalues of unit modulus and the powers $\tilde{M}_{h}^{n}$ , $n=0,1,\ldots ,$ remain bounded.

Comparing (3.22) with the result of setting $t=h$ in (2.4), we see that, by consistency, $B_{h}=h+O(h^{2})$ and $C_{h}=-h+O(h^{2})$ , and therefore $A_{h}=(1+B_{h}C_{h})^{1/2}=1-h^{2}/2+O(h^{3})$ . Thus, for $h>0$ sufficiently small, $A_{h}<1$ and the integration will be stable. The stability interval of the integrator is the longest interval $(0,h_{\max })$ such that integrations with $h\in (0,h_{\max })$ are stable. For reasons discussed in Example 3.14, methods with long stability intervals are often to be favoured. From Example 3.7, the midpoint rule has stability interval $(0,\infty )$ . Explicit integrators have stability intervals of finite length.

For each $h$ such that $|A_{h}|\leq 1$ , it is expedient to introduce $\unicode[STIX]{x1D703}_{h}\in \mathbb{R}$ such that $A_{h}=D_{h}=\cos \unicode[STIX]{x1D703}_{h}$ . For $|A_{h}|<1$ , we have $\sin \unicode[STIX]{x1D703}_{h}\neq 0$ and we may define

(4.7) $$\begin{eqnarray}\unicode[STIX]{x1D712}_{h}=B_{h}/\sin \unicode[STIX]{x1D703}_{h}.\end{eqnarray}$$

In terms of $\unicode[STIX]{x1D703}_{h}$ and $\unicode[STIX]{x1D712}_{h}$ , the matrices in (3.22) and (3.23) are then

(4.8) $$\begin{eqnarray}\tilde{M}_{h}=\left[\begin{array}{@{}cc@{}}\cos \unicode[STIX]{x1D703}_{h} & \unicode[STIX]{x1D712}_{h}\sin \unicode[STIX]{x1D703}_{h}\\ -\unicode[STIX]{x1D712}_{h}^{-1}\sin \unicode[STIX]{x1D703}_{h} & \cos \unicode[STIX]{x1D703}_{h}\end{array}\right]\end{eqnarray}$$

and

(4.9) $$\begin{eqnarray}\tilde{M}_{h}^{n}=\left[\begin{array}{@{}cc@{}}\cos (n\unicode[STIX]{x1D703}_{h}) & \unicode[STIX]{x1D712}_{h}\sin (n\unicode[STIX]{x1D703}_{h})\\ -\unicode[STIX]{x1D712}_{h}^{-1}\sin (n\unicode[STIX]{x1D703}_{h}) & \cos (n\unicode[STIX]{x1D703}_{h})\end{array}\right]\!.\end{eqnarray}$$

For a value of $h$ in the (stable) case $A_{h}=\pm 1$ , $B_{h}=C_{h}=0$ , we have $\sin \unicode[STIX]{x1D703}_{h}=0$ , so (4.7) does not make sense. However, the matrix $\tilde{M}_{h}$ is of the form (4.8) for any choice of $\unicode[STIX]{x1D712}_{h}$ . (Typically, for such a value of $h$ , one may avoid the indeterminacy in the value of $\unicode[STIX]{x1D712}_{h}$ by taking limits as $\unicode[STIX]{x1D716}\rightarrow 0$ in $\unicode[STIX]{x1D712}_{h+\unicode[STIX]{x1D716}}=B_{h+\unicode[STIX]{x1D716}}/\sin \unicode[STIX]{x1D703}_{h+\unicode[STIX]{x1D716}}$ .)

For a method of order $\unicode[STIX]{x1D708}$ , $\unicode[STIX]{x1D712}_{h}=1+O(h^{\unicode[STIX]{x1D708}})$ , $\unicode[STIX]{x1D703}_{h}=h+O(h^{\unicode[STIX]{x1D708}+1})$ as $h\rightarrow 0$ . By comparing the numerical $\tilde{M}_{h}^{n}$ in (4.9) with the true $M_{nh}$ in (2.4), one sees that a method with $\unicode[STIX]{x1D703}_{h}=h$ would have no phase error: the angular frequency of the rotation of the numerical solution would coincide with the true angular frequency of the harmonic oscillator. More generally, the difference $\unicode[STIX]{x1D703}_{h}-h$ governs the phase error. According to (4.9), this phase error grows linearly with $n$ (recall Table 3.1). On the other hand, a method with $\unicode[STIX]{x1D712}_{h}=1$ would have no energy error: the numerical solution would remain on the correct level curve of the Hamiltonian, i.e. on the circle $p^{2}+q^{2}=p_{0}^{2}+q_{0}^{2}$ . The discrepancy between $\unicode[STIX]{x1D712}_{h}$ and $1$ governs the energy errors. In (4.9) we see that these are bounded as $n$ grows.

The preceding considerations may alternatively be understood by considering the modified Hamiltonian given in the next result.

Proposition 4.9. Consider the application to the harmonic oscillator (2.3) of a (consistent) reversible volume-preserving integrator (3.22), and assume that the step size $h$ is stable, so that $\tilde{M}_{h}$ may be written in the form (4.8). Then $\unicode[STIX]{x1D713}_{h}$ exactly coincides with the $h$ flow of the modified Hamiltonian

$$\begin{eqnarray}\tilde{H}_{h}=\displaystyle \frac{\unicode[STIX]{x1D703}_{h}}{2h}\biggl(\unicode[STIX]{x1D712}_{h}p^{2}+\displaystyle \frac{1}{\unicode[STIX]{x1D712}_{h}}q^{2}\biggr).\end{eqnarray}$$

Figure 4.1.

Modified Hamiltonian for (a) velocity and (b) position Verlet. The dots are the first five points along a discrete orbit of the velocity/position Verlet integrator initiated at $(q,p)=(1,0)$ and the grey line provides the associated level set of the modified Hamiltonian. For comparison, the black contour line shows the corresponding level set of the true Hamiltonian.

In particular, numerical trajectories are contained in ellipses,

(4.10) $$\begin{eqnarray}\unicode[STIX]{x1D712}_{h}p^{2}+\displaystyle \frac{1}{\unicode[STIX]{x1D712}_{h}}q^{2}=\unicode[STIX]{x1D712}_{h}p_{0}^{2}+\displaystyle \frac{1}{\unicode[STIX]{x1D712}_{h}}q_{0}^{2},\end{eqnarray}$$

rather than in circles (Figure 4.1).

Remark 4.10. A comparison of a given integrator (3.22) with (2.4) shows that

(4.11) $$\begin{eqnarray}\left[\begin{array}{@{}c@{}}q_{n+1}\\ p_{n+1}\end{array}\right]=\tilde{M}_{h}\left[\begin{array}{@{}c@{}}q_{n}\\ p_{n}\end{array}\right]\!,\quad \tilde{M}_{h}=\left[\begin{array}{@{}cc@{}}A_{h} & -C_{h}\\ -B_{h} & D_{h}\end{array}\right]\end{eqnarray}$$

is a second integrator of the same order of accuracy. We might think that (4.11) arises from (3.22) by changing the roles of the variables $q$ and $p$ . In the particular case of splitting integrators, (4.11) arises from (3.22) after swapping the roles of the split systems $\text{A}$ and $\text{B}$ . The integrators (3.22) and (4.11) share the same interval of stability and the same $\unicode[STIX]{x1D703}_{h}$ . The function $\unicode[STIX]{x1D712}_{h}$ of (4.11) is obtained by changing the sign of the reciprocal of the function $\unicode[STIX]{x1D712}_{h}$ of (3.22). The important function $\unicode[STIX]{x1D70C}(h)$ to be introduced in Proposition 6.4 is also the same for (3.22) and (4.11). The velocity Verlet algorithm and the position Verlet algorithm provide an example of this kind of pair of integrators (see Examples 3.8 and 3.9).

The results we have just presented may be extended to quadratic Hamiltonians with $d$ degrees of freedom (3.27): it is sufficient to use diagonalization as in Proposition 3.13. In particular, for stability we require that the stability interval of the integrator contains all products $h\unicode[STIX]{x1D714}_{i}$ , where the frequencies $\unicode[STIX]{x1D714}_{i}$ are the square roots of the eigenvalues of $L^{-1}KL^{-T}$ , with $LL^{T}=M$ .

4.5 Optimal stability of Strang’s method

Let us fix an integer $N$ and consider consistent palindromic splitting integrators (3.19)–(3.20) with $s=N$ ; these use $N$ evaluations of $F$ per step when applied to problems of the form (3.9)–(3.10). The corresponding coefficient $A_{h}$ in (3.22) is a polynomial of degree $N$ in the variable $z=h^{2}$ (for obvious reasons $A_{h}$ is often called the stability polynomial of the integrator). We pointed out above that consistency imposes the relation $A_{h}=1-z/2+O(h^{2})$ . Our aim is to identify, among the class just described, the polynomial $A_{h}(z)$ that satisfies $|A_{h}(z)|<1$ for $0<h<h_{\text{max}}$ with $h_{\text{max}}$ as large as possible. As we shall show presently, the $A_{h}(z)$ sought corresponds to the integrator

(4.12) $$\begin{eqnarray}\unicode[STIX]{x1D713}_{h}=\overbrace{\unicode[STIX]{x1D713}_{h/N}^{V}\circ \cdots \circ \unicode[STIX]{x1D713}_{h/N}^{V}}^{N\text{ times}},\end{eqnarray}$$

where $\unicode[STIX]{x1D713}_{h}^{V}$ is the mapping associated with the Strang/Verlet formula (3.13).Footnote ⁷ Note that to carry out a step of length $h$ with the method in (4.12), one just has to take $N$ consecutive steps of length $h/N$ of standard velocity Verlet. In other words, subject to stability, if one wishes to take as long a step as possible with a budget of $N$ evaluations of the force per step, the best choice is to concatenate $N$ steps of Strang/Verlet. Footnote ⁸

To see the optimality of $\unicode[STIX]{x1D713}_{h}$ in (4.12), we first note that, after expressing $\unicode[STIX]{x1D713}_{h/N}^{V}$ in terms of the $\text{A}$ and $\text{B}$ flows and merging consecutive $\text{B}$ flows, the mapping (4.12) corresponds indeed to a palindromic splitting with $N$ stages. Then, from Example 3.8 we know that Verlet is stable for $0<h<2$ and this implies that (4.12) is stable for $0<h/N<2$ (for (4.12) the powers of $\tilde{M}_{h}$ are powers of the Verlet matrix $\tilde{M}_{h/N}^{V}$ ). In this way (4.12) has stability interval $(0,2N)$ and we shall prove next that this is the longest possible. From Example ?? 3.8, Verlet with step size $h/N$ has $A_{h/N}^{V}=1-h^{2}/(2N^{2})$ , which, in view of (4.8)–(4.9), implies that for (4.12) the coefficient $A_{h}$ has the expression

$$\begin{eqnarray}\cos \unicode[STIX]{x1D703}_{h}=\cos (N\unicode[STIX]{x1D703}_{h/N}^{V}).\end{eqnarray}$$

Recalling the definition $T_{N}(\cos \unicode[STIX]{x1D6FC})=\cos (N\unicode[STIX]{x1D6FC})$ of the Chebyshev polynomial $T_{N}$ , we observe that for (4.12)

$$\begin{eqnarray}A_{h}(z)=T_{N}\biggl(1-\displaystyle \frac{z}{2N^{2}}\biggr).\end{eqnarray}$$

Well-known properties of $T_{N}$ imply that no other polynomial $A_{h}(z)$ of degree $\leq N$ with $A_{h}=1-z/2+O(h^{2})$ has modulus $\leq 1$ in the interval $-1<1-z/(2N^{2})<1$ , that is, when $0<h<2N$ .

5.1 Simple Monte Carlo methods

Given a probability distribution $\unicode[STIX]{x1D707}$ in $\mathbb{R}^{D}$ (the target distribution) and a function $F:\mathbb{R}^{D}\rightarrow \mathbb{R}$ , the problem addressed by the algorithms considered here is to numerically estimate the following $D$ -dimensional integral with respect to $\unicode[STIX]{x1D707}$ :

(5.1) $$\begin{eqnarray}\unicode[STIX]{x1D707}(F)=\int _{\mathbb{R}^{D}}F(x)\unicode[STIX]{x1D707}(\text{d}x).\end{eqnarray}$$

In general, $\unicode[STIX]{x1D707}(F)$ cannot be determined analytically. Moreover, since the dimension $D$ might not be small, conventional numerical quadrature is likely not to be practical or even feasible.

The simple Monte Carlo method approximately computes $\unicode[STIX]{x1D707}(F)$ in (5.1) by generating $N$ independent and identically distributed (i.i.d.) samples $X_{1},\ldots ,X_{N}$ from $\unicode[STIX]{x1D707}$ , evaluating the function $F$ at these samples and using the estimator

(5.2) $$\begin{eqnarray}\bar{F}_{N}=\displaystyle \frac{1}{N}\mathop{\sum }_{i=1}^{N}F(X_{i}).\end{eqnarray}$$

Assuming that $\unicode[STIX]{x1D707}(F)<\infty$ , the law of large numbers states that

$$\begin{eqnarray}\lim _{N\rightarrow \infty }\bar{F}_{N}=\unicode[STIX]{x1D707}(F)\quad \text{as }N\rightarrow \infty ,\end{eqnarray}$$

almost surely. If, in addition, the standard deviation $\unicode[STIX]{x1D70E}_{0}(F)$ of the random variable $F(X)$ , $X\sim \unicode[STIX]{x1D707}$ , defined by

(5.3) $$\begin{eqnarray}\unicode[STIX]{x1D70E}_{0}(F)^{2}=\int _{\mathbb{ R}^{D}}(F(x)-\unicode[STIX]{x1D707}(F))^{2}\unicode[STIX]{x1D707}(\text{d}x)\end{eqnarray}$$

is finite, the central limit theorem ensures the following distributional limit:

$$\begin{eqnarray}N^{1/2}(\bar{F}_{N}-\unicode[STIX]{x1D707}(F))\overset{d}{\rightarrow }{\mathcal{N}}(0,\unicode[STIX]{x1D70E}_{0}(F)^{2})\quad \text{as }N\rightarrow \infty .\end{eqnarray}$$

Loosely speaking, this may be interpreted as stating that the distribution of $\bar{F}_{N}$ is approximately ${\mathcal{N}}(0,\unicode[STIX]{x1D70E}_{0}(F)^{2}/N)$ . Hence, the standard deviation of the error $\bar{F}_{N}-\unicode[STIX]{x1D707}(F)$ decreases like the inverse square root of the number of samples. Often this standard deviation is referred to as Monte Carlo error. Thus, to halve the Monte Carlo error the number of i.i.d. samples needs to be quadrupled.

In most cases of practical interest, we cannot directly generate i.i.d. samples $X_{i}$ from $\unicode[STIX]{x1D707}$ and must therefore resort to Markov chain Monte Carlo methods.

5.2 Markov chain Monte Carlo methods

Recall that a Markov chain with state space $\mathbb{R}^{D}$ is a sequence of random $D$ -vectors $\{X_{i}\}_{i\in \mathbb{N}}$ that satisfies the Markov property

$$\begin{eqnarray}\mathbb{P}(X_{i+1}\in A\mid X_{1},\ldots ,X_{i})=\mathbb{P}(X_{i+1}\in A\mid X_{i})\end{eqnarray}$$

for all measurable sets $A$ . In other words, given the past history of the chain $X_{1},\ldots ,X_{i}$ , the only information required to update the state of the chain is the current state $X_{i}$ . Here our interest is restricted to (time-) homogeneous chains, that is, to cases where $\mathbb{P}(X_{i+1}\in A\mid X_{i})$ is independent of $i$ .

Typically one constructs a homogeneous Markov chain in terms of its transition probabilities $\unicode[STIX]{x1D6F1}_{x}$ , $x\in \mathbb{R}^{D}$ . These are the probabilities

$$\begin{eqnarray}\unicode[STIX]{x1D6F1}_{x}(A)=\mathbb{P}(X_{i+1}\in A\mid X_{i}=x),\end{eqnarray}$$

with $i\in \mathbb{N}$ and $A$ measurable. Often $\unicode[STIX]{x1D6F1}_{x}$ may be computed as

$$\begin{eqnarray}\unicode[STIX]{x1D6F1}_{x}(A)=\int _{A}\unicode[STIX]{x1D6F1}_{x}(\text{d}x^{\prime }),\end{eqnarray}$$

for a suitable kernel $\unicode[STIX]{x1D6F1}_{x}(\text{d}x^{\prime })$ . Clearly the chain is determined once the transition probabilities and the distribution of $X_{1}$ are known. In practice, the term ‘chain’ is used in a wide sense to refer to the transition probabilities without specifying the distribution of $X_{1}$ .

A probability distribution $\unicode[STIX]{x1D708}$ is an invariant or stationary distribution of a Markov chain with transition probabilities $\unicode[STIX]{x1D6F1}_{x}$ if

$$\begin{eqnarray}\int _{\mathbb{R}^{D}}\unicode[STIX]{x1D6F1}_{x}(A)\unicode[STIX]{x1D708}(\text{d}x)=\int _{A}\unicode[STIX]{x1D708}(\text{d}x)\end{eqnarray}$$

holds for all measurable sets $A$ . We also say that $\unicode[STIX]{x1D6F1}_{x}$ preserves $\unicode[STIX]{x1D708}$ . In the situations we are interested in, a Markov chain will have a unique invariant distribution. If $\unicode[STIX]{x1D708}$ is the invariant distribution of the chain and in addition $X_{1}\sim \unicode[STIX]{x1D708}$ , one says that the chain is at stationarity.

Markov chain Monte Carlo (MCMC) methods generate a Markov chain $\{X_{i}\}_{i\in \mathbb{N}}$ that has the target $\unicode[STIX]{x1D707}$ as an invariant distribution and estimate $\unicode[STIX]{x1D707}(F)$ by the average (5.2). By analogy to the simple i.i.d. situation described above, one would like to have MCMC methods that meet two basic requirements.

• ∙ For each $F:\mathbb{R}^{D}\rightarrow \mathbb{R}$ such that $\unicode[STIX]{x1D707}(F)<\infty$ ,

  (5.4) $$\begin{eqnarray}\bar{F}_{N}\overset{\text{a.s.}}{\rightarrow }\unicode[STIX]{x1D707}(F)\quad \text{as }N\rightarrow \infty .\end{eqnarray}$$
  This is the MCMC analogue of the law of large numbers.

• ∙ For each $F:\mathbb{R}^{D}\rightarrow \mathbb{R}$ such that $\unicode[STIX]{x1D707}(F)<\infty$ and $\unicode[STIX]{x1D70E}_{0}(F)<\infty$ ,

  (5.5) $$\begin{eqnarray}N^{1/2}(\bar{F}_{N}-\unicode[STIX]{x1D707}(F))\overset{d}{\rightarrow }{\mathcal{N}}(0,\unicode[STIX]{x1D70E}(F)^{2})\quad \text{as }N\rightarrow \infty ,\end{eqnarray}$$
  for some $\unicode[STIX]{x1D70E}(F)$ . This is the MCMC analogue of the central limit theorem.

For each fixed function $F$ and Markov chain $\{X_{i}\}_{i\in \mathbb{N}}$ , the constant $\unicode[STIX]{x1D70E}(F)^{2}$ appearing in (5.5) is called the asymptotic variance of the MCMC estimator $\bar{F}_{N}$ . A straightforward calculation shows that this asymptotic variance satisfies

$$\begin{eqnarray}\unicode[STIX]{x1D70E}(F)^{2}=\unicode[STIX]{x1D70E}_{0}(F)^{2}+2\mathop{\sum }_{i>1}\operatorname{cov}_{\unicode[STIX]{x1D707}}(F(X_{1}),F(X_{i})),\end{eqnarray}$$

where $\unicode[STIX]{x1D70E}_{0}(F)$ is defined in (5.3) and the covariances are computed assuming that the chain is at stationarity. If the $X_{i}$ were independent, all the covariances $\operatorname{cov}_{\unicode[STIX]{x1D707}}(F(X_{1}),F(X_{i}))$ would vanish and we would recover the standard central limit theorem. Since, in most interesting cases, the $X_{i}$ in the Markov chain are not mutually independent, $\unicode[STIX]{x1D70E}(F)$ is often larger than $\unicode[STIX]{x1D70E}_{0}(F)$ . Generally speaking it is desirable to have low values of $\operatorname{cov}_{\unicode[STIX]{x1D707}}(X_{1},X_{i})$ so that $\unicode[STIX]{x1D70E}(F)$ is not far away from $\unicode[STIX]{x1D70E}_{0}(F)$ for each $F$ .

In practice, the inputs that the user has to supply to an MCMC algorithm include the following.

• ∙ A sample of the initial state $X_{1}$ . Ideally, this sample should be taken in a domain of state space of high probability. Otherwise the chain may need many steps to start generating useful samples. A discussion of this issue is beyond the scope of this paper.

• ∙ A possibly non-normalized density function $\unicode[STIX]{x1D70C}$ of the target $\unicode[STIX]{x1D707}$ , that is, the probability density function is $Z^{-1}\unicode[STIX]{x1D70C}(x)$ , where $Z=\int _{\mathbb{R}^{D}}\unicode[STIX]{x1D70C}(x)\,\text{d}x$ is not assumed to be $1$ ; the value of $Z$ is not required to run the algorithms.

5.3 Metropolis method for reversible maps

The replacement of the i.i.d. variables that simple Monte Carlo uses in the estimator (5.2) with variables of a Markov chain is of interest because it is not difficult to construct a chain that has a given target $\unicode[STIX]{x1D707}$ as an invariant distribution. The key to this construction is the Metropolis–Hastings accept/reject mechanism, which turns a given proposal chain (for which $\unicode[STIX]{x1D707}$ is not invariant) into a Metropolized chain, which leaves $\unicode[STIX]{x1D707}$ invariant. The simplest Metropolis rule was introduced in 1953 in a landmark paper (Metropolis et al. Reference Metropolis, Rosenbluth, Rosenbluth, Teller and Teller1953); later, Hastings provided an important generalization (Hastings Reference Hastings1970).

A review of the Metropolis–Hastings rule is not required for our purposes. However, we shall present a trimmed-down variant of Metropolis–Hastings that we will use to define HMC. This variant works in the special case where the target $\unicode[STIX]{x1D707}$ is invariant with respect to a linear involution, takes as input a reversible deterministic map, and manufactures a Markov chain that preserves $\unicode[STIX]{x1D707}$ . While the chain that we construct is not expected to satisfy a law of large numbers and therefore has no practical merit, Proposition 5.2 will be used later to analyse HMC methods.

The technique, patterned after Fang, Sanz-Serna and Skeel (Reference Fang, Sanz-Serna and Skeel2014), requires a non-normalized density function $\unicode[STIX]{x1D70C}:\mathbb{R}^{D}\rightarrow \mathbb{R}$ of the target, and, as pointed out above, a map $\unicode[STIX]{x1D6F7}:\mathbb{R}^{D}\rightarrow \mathbb{R}^{D}$ that is reversible with respect to a linear involution $S$ that preserves probability, that is, $\unicode[STIX]{x1D70C}\circ S=\unicode[STIX]{x1D70C}$ . The accept/reject mechanism is based on the acceptance probability $\unicode[STIX]{x1D6FC}:\mathbb{R}^{D}\rightarrow [0,1]$ defined as

(5.6) $$\begin{eqnarray}\unicode[STIX]{x1D6FC}(x)=\min \biggl\{1,\displaystyle \frac{\unicode[STIX]{x1D70C}(\unicode[STIX]{x1D6F7}(x))}{\unicode[STIX]{x1D70C}(x)}|\text{det}\,\unicode[STIX]{x1D6F7}^{\prime }(x)|\biggr\}\end{eqnarray}$$

( $\unicode[STIX]{x1D6F7}^{\prime }$ is the Jacobian matrix of $\unicode[STIX]{x1D6F7}$ ).

We consider the following algorithm.

Step 2 contains the accept/reject mechanism. In the case of acceptance the updated state coincides with the state proposed $\tilde{X}_{1}$ from Step 1; in the case of rejection the updated state is $S(X_{0})$ . Note that for rejection, conventional Metropolis mechanisms set the updated state of the chain to be $X_{0}$ .

Proof. The transition kernel of $\{X_{i}\}_{i\in \mathbb{N}}$ is given by

$$\begin{eqnarray}\unicode[STIX]{x1D6F1}_{x}(\text{d}x^{\prime })=\unicode[STIX]{x1D6FC}(x)\unicode[STIX]{x1D6FF}(x^{\prime }-\unicode[STIX]{x1D6F7}(x))\,\text{d}x^{\prime }+(1-\unicode[STIX]{x1D6FC}(x))\unicode[STIX]{x1D6FF}(x^{\prime }-S(x))\,\text{d}x^{\prime },\end{eqnarray}$$

where $\unicode[STIX]{x1D6FF}(\cdot )$ is the Dirac-delta function. Hence, for any measurable set $A$ , if $1_{A}(x)$ denotes the corresponding indicator function,

$$\begin{eqnarray}\displaystyle & & \displaystyle \int _{\mathbb{R}^{D}}\unicode[STIX]{x1D6F1}_{x}(A)\unicode[STIX]{x1D70C}(x)\,\text{d}x=\int _{A}\unicode[STIX]{x1D70C}(x)\,\text{d}x\nonumber\\ \displaystyle & & \displaystyle \qquad +\int _{\mathbb{R}^{D}}\unicode[STIX]{x1D6FC}(x)\unicode[STIX]{x1D70C}(x)1_{A}(\unicode[STIX]{x1D6F7}(x))\,\text{d}x-\int _{\mathbb{R}^{D}}\unicode[STIX]{x1D6FC}(x)\unicode[STIX]{x1D70C}(x)1_{A}(S(x))\,\text{d}x.\nonumber\end{eqnarray}$$

By change of variables in the third integral,

$$\begin{eqnarray}\displaystyle & & \displaystyle \int _{\mathbb{R}^{D}}\unicode[STIX]{x1D6F1}_{x}(A)\unicode[STIX]{x1D70C}(x)\,\text{d}x=\int _{A}\unicode[STIX]{x1D70C}(x)\,\text{d}x\nonumber\\ \displaystyle & & \displaystyle \qquad +\int _{\mathbb{R}^{D}}\unicode[STIX]{x1D6FC}(x)\unicode[STIX]{x1D70C}(x)1_{A}(\unicode[STIX]{x1D6F7}(x))\,\text{d}x\nonumber\\ \displaystyle & & \displaystyle \qquad -\int _{\mathbb{R}^{D}}\unicode[STIX]{x1D6FC}(S(\unicode[STIX]{x1D6F7}(x)))|\text{det}\,\unicode[STIX]{x1D6F7}^{\prime }(x)|\unicode[STIX]{x1D70C}(\unicode[STIX]{x1D6F7}(x))1_{A}(\unicode[STIX]{x1D6F7}(x))\,\text{d}x,\nonumber\end{eqnarray}$$

where we used the hypothesis that $\unicode[STIX]{x1D70C}\circ S=\unicode[STIX]{x1D70C}$ . The last two terms on the right-hand side of this equation cancel because

$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FC}(S(\unicode[STIX]{x1D6F7}(x))) & = & \displaystyle \min \biggl\{1,\displaystyle \frac{\unicode[STIX]{x1D70C}(x)}{\unicode[STIX]{x1D70C}(\unicode[STIX]{x1D6F7}(x))}|\text{det}\,\unicode[STIX]{x1D6F7}^{\prime }(x)|^{-1}\biggr\}\nonumber\\ \displaystyle & = & \displaystyle \unicode[STIX]{x1D6FC}(x)\unicode[STIX]{x1D70C}(x)\displaystyle \frac{1}{\unicode[STIX]{x1D70C}(\unicode[STIX]{x1D6F7}(x))}|\text{det}\,\unicode[STIX]{x1D6F7}^{\prime }(x)|^{-1},\nonumber\end{eqnarray}$$

which follows from the reversibility of $\unicode[STIX]{x1D6F7}$ , the hypothesis $\unicode[STIX]{x1D70C}\circ S=\unicode[STIX]{x1D70C}$ and Proposition 2.8.☐

5.4 The HMC method: basic idea

We consider a target distribution $\unicode[STIX]{x1D6F1}$ in $\mathbb{R}^{d}$ . If ${\mathcal{U}}:\mathbb{R}^{d}\rightarrow \mathbb{R}$ denotes the negative logarithm of the (not necessarily normalized) probability density function of the target, then

$$\begin{eqnarray}\unicode[STIX]{x1D6F1}(\text{d}q)=Z_{q}^{-1}\exp (-{\mathcal{U}}(q))\,\text{d}q,\quad \text{where }Z_{q}=\int _{\mathbb{R}^{d}}\exp (-{\mathcal{U}}(q))\,\text{d}q.\end{eqnarray}$$

The Monte Carlo algorithms studied here use ${\mathcal{U}}$ but do not require knowledge of the normalization factor $Z_{q}$ . In HMC, regardless of the application in mind, ${\mathcal{U}}$ is seen as the potential energy of a mechanical system with coordinates $q$ . Then auxiliary momenta $p\in \mathbb{R}^{d}$ and a quadratic kinetic energy function ${\mathcal{T}}(p)=(1/2)p^{T}M^{-1}p$ are introduced as in (2.10) ( $M$ is a positive-definite symmetric matrix chosen by the user).Footnote ⁹ The total energy of this fictitious mechanical system is $H={\mathcal{T}}+{\mathcal{U}}$ and the equations of motion are given in (2.9).

Example 5.3. In the particular case where the target $\unicode[STIX]{x1D707}$ is Gaussian with non-normalized density $\exp (-(1/2)q^{T}Kq)$ , $q\in \mathbb{R}^{d}$ , we have

$$\begin{eqnarray}{\mathcal{U}}(q)=(1/2)q^{T}Kq,\end{eqnarray}$$

which leads to the Hamiltonian (3.27) discussed above. For a univariate standard normal target, $H=(1/2)(p^{2}+q^{2})$ and Hamilton’s equations reduce to those of harmonic oscillator (2.3).

The Boltzmann–Gibbs distribution in $\mathbb{R}^{2d}$ corresponding to $H$ was discussed earlier in connection with Theorem 2.7. This distribution is defined as

(5.7) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6F1}_{BG}(\text{d}q,\text{d}p) & = & \displaystyle (2\unicode[STIX]{x1D70B})^{-d/2}|\text{det}\,M|^{-1/2}\exp \biggl(-\displaystyle \frac{1}{2}p^{T}M^{-1}p\biggr)\nonumber\\ \displaystyle & & \displaystyle \quad \times \,Z_{q}^{-1}\exp (-{\mathcal{U}}(q))\,\text{d}q\,\text{d}p\end{eqnarray}$$

(for simplicity the inverse temperature is taken here to be $\unicode[STIX]{x1D6FD}=1$ ). Clearly the target $\unicode[STIX]{x1D6F1}$ is the $q$ -marginal of $\unicode[STIX]{x1D6F1}_{BG}$ . The $p$ -marginal is Gaussian with zero mean and covariance matrix $M$ ; therefore samples from this marginal are easily available (and will be put to use in the algorithms below). A key fact for our purposes: Hamilton’s equations of motion (2.9) preserve $\unicode[STIX]{x1D6F1}_{BG}$ (Theorem 2.7).

HMC generates (correlated) samples $(q_{i},p_{i})\in \mathbb{R}^{2d}$ by means of a Markov chain that leaves $\unicode[STIX]{x1D6F1}_{BG}$ invariant; the corresponding marginal $q_{i}\in \mathbb{R}^{d}$ chain then leaves invariant the target distribution $\unicode[STIX]{x1D6F1}$ . The basic idea of HMC is encapsulated in the following algorithm (the duration $\unicode[STIX]{x1D706}>0$ is a – deterministic – parameter, whose value is specified by the user).

Note that $p_{0}$ plays no role, since the initial condition starts from $\unicode[STIX]{x1D709}_{0}$ . Step 1 is referred to as momentum refreshment or momentum randomization.

It is easy to see that this algorithm succeeds in preserving the distribution $\unicode[STIX]{x1D6F1}_{BG}$ .

The transition kernel of this chain is given by

$$\begin{eqnarray}\unicode[STIX]{x1D6F1}_{(q,p)}(\text{d}q^{\prime },\text{d}p^{\prime })=\mathbb{E}\{\unicode[STIX]{x1D6FF}((q^{\prime },p^{\prime })-\unicode[STIX]{x1D711}_{\unicode[STIX]{x1D706}}(q,\unicode[STIX]{x1D709}))\}\,\text{d}q^{\prime }\,\text{d}p^{\prime },\end{eqnarray}$$

where the expected value is over $\unicode[STIX]{x1D709}\sim {\mathcal{N}}(0,M)$ .

The most appealing feature of the algorithm is that if $\unicode[STIX]{x1D706}$ is sufficiently large, we may hope that the Markov transitions $i\rightarrow i+1$ produce values $q_{i+1}$ far away from $q_{i}$ , thus reducing the correlations in the chain and facilitating the exploration of the target distribution.

5.5 Numerical HMC

Algorithm 5.4 cannot be used in practice because in the cases of interest the exact solution flow $\unicode[STIX]{x1D711}_{\unicode[STIX]{x1D706}}$ of Hamilton’s equations is not available. It is then necessary to resort to numerical approximations to $\unicode[STIX]{x1D711}_{\unicode[STIX]{x1D706}}$ , but, as pointed out in Section 4, numerical methods cannot preserve volume in phase space and energy and therefore do not preserve exactly the Boltzmann–Gibbs distribution. To correct the bias introduced by the time discretization error, the numerical solution is Metropolized using Algorithm 5.1. However, this requires the numerical integrator to be reversible.

Let $\unicode[STIX]{x1D6F9}_{\unicode[STIX]{x1D706}}$ denote a numerical approximation of $\unicode[STIX]{x1D711}_{\unicode[STIX]{x1D706}}$ (more precisely, if the step size is $h$ and $n=\lfloor \unicode[STIX]{x1D706}/h\rfloor$ is the number of steps required to integrate up to $t=\unicode[STIX]{x1D706}$ , then $\unicode[STIX]{x1D6F9}_{\unicode[STIX]{x1D706}}=\unicode[STIX]{x1D713}_{h}^{n}$ ). In order to use Algorithm 5.1 with $\unicode[STIX]{x1D6F9}_{\unicode[STIX]{x1D706}}$ playing the role of $\unicode[STIX]{x1D6F7}$ and the momentum flip involution (2.13) playing the role of $S$ , we first note (Proposition 2.11) that the momentum flip involution preserves $H$ in (2.10) and, as a consequence, it preserves the Boltzmann–Gibbs distribution. In addition Theorem 2.10 ensures that the Hamiltonian flow is reversible with respect to this involution; it then makes sense (Theorem 4.2) to assume that the integrator chosen is such that $\unicode[STIX]{x1D6F9}_{\unicode[STIX]{x1D706}}$ is also reversible. The acceptance probability in (5.6) now reads

(5.8) $$\begin{eqnarray}\unicode[STIX]{x1D6FC}(q,p)=\min \{1,\text{e}^{-\unicode[STIX]{x1D6E5}H(q,p)}|\text{det}\,\unicode[STIX]{x1D6F9}_{\unicode[STIX]{x1D706}}^{\prime }(q,p)|\},\end{eqnarray}$$

where

$$\begin{eqnarray}\unicode[STIX]{x1D6E5}H(q,p)=H(\unicode[STIX]{x1D6F9}_{\unicode[STIX]{x1D706}}(q,p))-H(q,p)\end{eqnarray}$$

is the energy error (recall that if the integrator were exact $H(\unicode[STIX]{x1D6F9}_{\unicode[STIX]{x1D706}}(q,p))$ would coincide with $H(q,p)$ by conservation of energy, Theorem 2.6).

The transition kernel of the chain is given by

$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6F1}_{(q,p)}(\text{d}q^{\prime },\text{d}p^{\prime }) & = & \displaystyle \mathbb{E}\{\unicode[STIX]{x1D6FC}(q,\unicode[STIX]{x1D709})\unicode[STIX]{x1D6FF}((q^{\prime },p^{\prime })-\unicode[STIX]{x1D6F9}_{\unicode[STIX]{x1D706}}(q,\unicode[STIX]{x1D709}))\}\,\text{d}q^{\prime }\,\text{d}p^{\prime }\nonumber\\ \displaystyle & & \displaystyle \quad +\,\mathbb{E}\{(1-\unicode[STIX]{x1D6FC}(q,\unicode[STIX]{x1D709}))\unicode[STIX]{x1D6FF}((q^{\prime },p^{\prime })-S(q,\unicode[STIX]{x1D709}))\}\,\text{d}q^{\prime }\,\text{d}p^{\prime }.\nonumber\end{eqnarray}$$

In practice, the acceptance probability (5.8) may not be readily available due to the need to compute $\det \unicode[STIX]{x1D6F9}_{\unicode[STIX]{x1D706}}^{\prime }(q,p)$ . If the numerical approximation $\unicode[STIX]{x1D6F9}_{\unicode[STIX]{x1D706}}$ , in addition to being assumed reversible, is also volume-preserving (as it would be for splitting integrators according to Theorem 4.1), then the determinant drops from the formula, and then the acceptance probability

(5.9) $$\begin{eqnarray}\unicode[STIX]{x1D6FC}(q,p)=\min \{1,\text{e}^{-\unicode[STIX]{x1D6E5}H(q,p)}\}\end{eqnarray}$$

becomes easily computable. Variants where preservation of volume does not take place are studied by Fang et al. (Reference Fang, Sanz-Serna and Skeel2014).

5.6 Exact randomized HMC

The Hamiltonian flow in Step 2 of Algorithm 5.4 is what, in principle, enables HMC to make large moves in state space that reduce correlations in the Markov chain $\{q_{i}\}_{i\in \mathbb{N}}$ . Roughly speaking, one may hope that by increasing the duration $\unicode[STIX]{x1D706}$ , $q_{1}$ moves away from $q_{0}$ , thus reducing correlation. However, simple examples show that this outcome is far from assured.

Indeed, for the univariate standard normal target distribution in Example 5.3, the Hamiltonian flow is a rotation in the $(q,p)$ -plane with period $2\unicode[STIX]{x1D70B}$ . It is easy to see that if $q_{0}$ is taken from the target distribution, as $\unicode[STIX]{x1D706}$ increases from $0$ to $\unicode[STIX]{x1D70B}/2$ , the correlation between $q_{1}$ and $q_{0}$ decreases and for $\unicode[STIX]{x1D706}=\unicode[STIX]{x1D70B}/2$ , $q_{1}$ and $q_{0}$ are independent. However, increasing $\unicode[STIX]{x1D706}$ beyond $\unicode[STIX]{x1D70B}/2$ will cause an increase in the correlation, and for $\unicode[STIX]{x1D706}=\unicode[STIX]{x1D70B}$ , $q_{1}=-q_{0}$ and the chain is not ergodic. For general distributions, it is likely that a small $\unicode[STIX]{x1D706}$ will lead to a highly correlated chain, while choosing $\unicode[STIX]{x1D706}$ too large may cause the Hamiltonian trajectory to make a U-turn and fold back on itself, thus increasing correlation (Homan and Gelman Reference Homan and Gelman2014). Generally speaking, the performance of HMC may be very sensitive to changes in $\unicode[STIX]{x1D706}$ , as first noted by Mackenzie (Reference Mackenzie1989). In order to increase the robustness of the algorithm, Mackenzie suggested randomly varying $\unicode[STIX]{x1D706}$ from one Markov transition to the next, and for that purpose he used a uniform distribution in the interval $[\unicode[STIX]{x1D706}_{\text{min}},\unicode[STIX]{x1D706}_{\text{max}}]$ .

Recently Bou-Rabee and Sanz-Serna (Reference Bou-Rabee and Sanz-Serna2017) have studied an algorithm where the lengths of the time intervals of integration of the Hamiltonian dynamics at the different transitions of the Markov chain are independent and identically distributed exponential random variables with mean $\unicode[STIX]{x1D706}$ ; these durations are of course taken to be independent of the state of the chain. The algorithm is then as follows.

Analogous to Theorem 5.5, the probability distribution $\unicode[STIX]{x1D6F1}_{\operatorname{BG}}$ in (5.7) is an invariant distribution of the Markov chain $\{(q_{i},p_{i})\}_{i\in \mathbb{N}}$ defined by iterating Algorithm 5.10. Analytical results and numerical experiments (Bou-Rabee and Sanz-Serna Reference Bou-Rabee and Sanz-Serna2017, §4–5) show that the dependence of the performance of Algorithm 5.10 on the mean duration parameter $\unicode[STIX]{x1D706}$ is simpler than the dependence of the performance of Algorithm 5.4 on its constant duration parameter.

5.7 Numerical randomized HMC

Unfortunately, the complex dependence of correlation on the duration parameter $\unicode[STIX]{x1D706}$ of Algorithm 5.4 is not removed by time discretization and is therefore inherited by Algorithm 5.6. For instance, for the univariate standard normal target, it is easy to check that if $\unicode[STIX]{x1D706}$ is close to an integer multiple of $\unicode[STIX]{x1D70B}$ and $h>0$ is suitably chosen, then a Verlet numerical integration will result, for each $q_{0}$ , in $q_{1}=-q_{0}$ (a move that will be accepted by the Metropolis–Hastings step).

To improve its performance, Algorithm 5.6 is typically operated with values of $h$ that are randomized (Neal Reference Neal and Brooks2011). Since, due to stability restrictions, explicit integrators cannot be used with arbitrarily large values of the time-step, $h$ is typically chosen from a uniform distribution in an (often narrow) interval $(\unicode[STIX]{x1D6E5}t_{\text{min}},\unicode[STIX]{x1D6E5}t_{\text{max}})$ . The number of time-steps $N$ in each integration leg is kept constant, and therefore the length of the integration intervals is random with a uniform distribution in $(N\unicode[STIX]{x1D6E5}t_{\text{min}},N\unicode[STIX]{x1D6E5}t_{\text{max}})$ . Even after such a randomization, the fact remains that increasing the duration parameter will increase the computational cost of each integration leg and may impair the quality of the sampling.

Several modifications of Algorithm 5.10 that use numerical integration are suggested in Bou-Rabee and Sanz-Serna (Reference Bou-Rabee and Sanz-Serna2017). The most obvious of them approximates the Hamiltonian flow in Step 3 by a volume-preserving reversible integrator (such as Verlet) operated with a fixed step size $h$ and with the number of integration steps $m$ at the different transitions of the Markov chain being independent and identically distributed geometric random variables with mean $\unicode[STIX]{x1D706}/h$ . These random numbers are of course taken to be independent of the state of the chain. As in Algorithm 5.6, one needs an accept/reject mechanism to remove the bias due to the energy errors introduced by this integrator.

Analogous to Theorem 5.7, the probability distribution $\unicode[STIX]{x1D6F1}_{\operatorname{BG}}$ in (5.7) is an invariant distribution of the Markov chain $\{(q_{i},p_{i})\}_{i\in \mathbb{N}}$ defined by iterating Algorithm 5.11.

In the remainder of this paper our attention is focused on Algorithm 5.6. Experiments based on Algorithm 5.11 are reported in Section 8.2 and an example that illustrates the effects of randomization of $h$ and $\unicode[STIX]{x1D706}$ is presented in Section 9.3.

6.1 Mean energy error

Step 2 in Algorithm 5.6 requires a volume-preserving reversible integrator if (5.9) is to be used. As we shall discuss now, those geometric properties have a direct impact on the mean energy error. We begin with an auxiliary result (Beskos et al. Reference Beskos, Pillai, Roberts, Sanz-Serna and Stuart2013, Lemma 3.3) that holds for any volume-preserving reversible map.

Proposition 6.1. Let $\unicode[STIX]{x1D6F9}:\mathbb{R}^{2d}\rightarrow \mathbb{R}^{2d}$ be a volume-preserving bijection that is reversible with respect to the momentum flip involution (2.13), and set $\unicode[STIX]{x1D6E5}(q,p)=H(\unicode[STIX]{x1D6F9}(q,p))-H(q,p)$ . If $g$ is an odd real function of a real variable, then

$$\begin{eqnarray}\int _{\mathbb{R}^{2d}}g(\unicode[STIX]{x1D6E5}(q,p))\,\text{e}^{-H(q,p)}\,\text{d}q\,\text{d}p=-\int _{\mathbb{ R}^{2d}}g(\unicode[STIX]{x1D6E5}(q,p))\,\text{e}^{-H(\unicode[STIX]{x1D6F9}(q,p))}\,\text{d}q\,\text{d}p,\end{eqnarray}$$

provided that one of the integrals exists. If $g$ is even, then

$$\begin{eqnarray}\int _{\mathbb{R}^{2d}}g(\unicode[STIX]{x1D6E5}(q,p))\,\text{e}^{-H(q,p)}\,\text{d}q\,\text{d}p=\int _{\mathbb{ R}^{2d}}g(\unicode[STIX]{x1D6E5}(q,p))\,\text{e}^{-H(\unicode[STIX]{x1D6F9}(q,p))}\,\text{d}q\,\text{d}p,\end{eqnarray}$$

provided that one of the integrals exists.

Proof. By a change of variables under the map $(q,p)\mapsto S(\unicode[STIX]{x1D6F9}(q,p))$ ,

$$\begin{eqnarray}\displaystyle \int _{\mathbb{R}^{2d}}g(\unicode[STIX]{x1D6E5}(q,p))\,\text{e}^{-H(q,p)}\,\text{d}q\,\text{d}p & = & \displaystyle \int _{\mathbb{R}^{2d}}g(\unicode[STIX]{x1D6E5}(S(\unicode[STIX]{x1D6F9}(q,p))))\,\text{e}^{-H(S(\unicode[STIX]{x1D6F9}(q,p)))}\,\text{d}q\,\text{d}p\nonumber\\ \displaystyle & = & \displaystyle \int _{\mathbb{R}^{2d}}g(-\unicode[STIX]{x1D6E5}(q,p))\,\text{e}^{-H(\unicode[STIX]{x1D6F9}(q,p))}\,\text{d}q\,\text{d}p,\nonumber\end{eqnarray}$$

where in the first step we used that both $S$ and $\unicode[STIX]{x1D6F9}$ preserve volume and in the second we took into account that $S$ leaves $H$ invariant and that the reversibility of $\unicode[STIX]{x1D6F9}$ implies

(6.1) $$\begin{eqnarray}\unicode[STIX]{x1D6E5}(q,p)=-\unicode[STIX]{x1D6E5}(S(\unicode[STIX]{x1D6F9}(q,p))).\end{eqnarray}$$

☐

Figure 6.1.

A univariate target with probability modes at $q=\pm 1$ , leading to a double-well potential ${\mathcal{U}}$ . The continuous lines are contours of constant $H$ . The symmetry of the contours with respect to the axis $p=0$ is a consequence of the reversibility of the Hamiltonian flow. The solutions of Hamilton’s equations move from left to right when $p>0$ and from right to left when $p<0$ , so when contours are reflected over the horizontal axis the arrow of time is reversed. If a reversible $\unicode[STIX]{x1D6F9}$ maps $(q_{0},p_{0})$ into $(q^{\ast },p^{\ast })$ , it has to map $(q^{\ast },-p^{\ast })$ into $(q_{0},p_{0})$ , so as to preserve the symmetry of the figure. The transition $(q_{0},p_{0})\mapsto (q^{\ast },p^{\ast })$ has an increase in energy and $(q^{\ast },-p^{\ast })\mapsto (q_{0},p_{0})$ decreases energy in exactly the same amount; this is the content of formula (6.1).

We now use the proposition to bound the average of $\unicode[STIX]{x1D6E5}$ .

Theorem 6.2. Let $\unicode[STIX]{x1D6F9}:\mathbb{R}^{2d}\rightarrow \mathbb{R}^{2d}$ be a bijection that is volume-preserving and reversible with respect to the momentum flip involution (2.13). If the integral

(6.2) $$\begin{eqnarray}m_{\unicode[STIX]{x1D6E5}}=\int _{\mathbb{R}^{2d}}\unicode[STIX]{x1D6E5}(q,p)\,\text{e}^{-H(q,p)}\,\text{d}q\,\text{d}p,\end{eqnarray}$$

with $\unicode[STIX]{x1D6E5}(q,p)=H(\unicode[STIX]{x1D6F9}(q,p))-H(q,p)$ , exists, then

(6.3) $$\begin{eqnarray}0\leq m_{\unicode[STIX]{x1D6E5}}\leq \int _{\mathbb{R}^{2d}}\unicode[STIX]{x1D6E5}(q,p)^{2}\,\text{e}^{-H(q,p)}\,\text{d}q\,\text{d}p.\end{eqnarray}$$

Furthermore the first inequality is strict except in the trivial case where $\unicode[STIX]{x1D6E5}(q,p)$ vanishes for each $(q,p)$ .