Geometric Numerical Integration (GNI) is the study and design of numerical methods for differential equations that preserve the underlying geometric structures of the exact flow: symplectic forms, invariants, symmetries, constraints, and manifold structure. Rather than focusing solely on minimizing local truncation error at a fixed time, GNI asks:
Can a numerical method approximate not only the solution values, but the qualitative, structural behaviour of the continuous dynamical system, especially over long times?
This viewpoint is synthesized in the monograph by Hairer, Lubich, and Wanner [HLW06], and is presented concisely in Blanes & Casas [BC26], on which this companion is based. [oai_citation:0‡Université Paris 13](https://www.math.univ-paris13.fr/~cuvelier/docs/Enseignements/Energetique/AnaNumII/12-13/TP1/G3/Geometric_Numerical_Integration_-_Hairer_Chap1.pdf?utm_source=chatgpt.com)
To motivate the subject, we begin from first principles: flows of ordinary differential equations, their geometric properties, and the limitations of traditional error-centric numerical analysis.
Consider an autonomous system of ordinary differential equations (ODEs)
\[ \dot{y}(t) = f(y(t)), \qquad y(t) \in \mathbb{R}^d, \quad f : \mathbb{R}^d \to \mathbb{R}^d. \]
Under standard assumptions (e.g. locally Lipschitz \(f\)), for any initial state \(y_0 \in \mathbb{R}^d\) and time \(t\) in some interval containing \(0\), there exists a unique solution \(y(t; y_0)\). The mapping
\[ \varphi^t : \mathbb{R}^d \to \mathbb{R}^d, \qquad \varphi^t(y_0) := y(t; y_0) \]
is called the flow of the ODE at time \(t\). For an autonomous system, the family \(\{\varphi^t\}_{t\in I}\) satisfies the flow (group) properties:
At this level, numerical methods are approximations of the one-step map \(\varphi^h\) for some timestep \(h\):
\[ y_{n+1} = \Phi_h(y_n) \approx \varphi^h(y_n). \]
The key idea in GNI is to treat \(\Phi_h\) itself as a discrete dynamical system and to study how well it reproduces the geometric and qualitative properties of \(\varphi^t\), especially for long-time integration.
The phase space \(\mathbb{R}^d\) is often endowed with additional structure:
The exact flow \(\varphi^t\) typically preserves some or all of these structures:
Geometric integrators are numerical methods \(\Phi_h\) that exactly preserve some of these properties for the discrete dynamics \(y_{n+1} = \Phi_h(y_n)\), regardless of step-size (within stability limits).
Classical numerical analysis for ODEs emphasizes:
For many applications this is sufficient. However, in problems such as:
the long-time qualitative behaviour is crucial: boundedness of energy, correct qualitative phase portraits, existence of quasi-periodic tori, etc. Standard high-order methods can reproduce short-time trajectories very accurately yet fail dramatically over long times by slowly but steadily destroying invariants.
Geometric numerical analysis shifts the emphasis:
The outcome is a theory that explains why relatively low-order geometric methods can dramatically outperform higher-order non-geometric methods in long-time simulations of structured systems.
A canonical example is the one-dimensional harmonic oscillator:
\[ \dot{q} = p, \qquad \dot{p} = -q, \]
with Hamiltonian
\[ H(q,p) = \frac{1}{2}(p^2 + q^2). \]
The exact flow preserves energy exactly: \(H(q(t),p(t)) = H(q_0,p_0)\) for all \(t\).
The explicit Euler method for this system is
\[ \begin{aligned} q_{n+1} &= q_n + h p_n,\\ p_{n+1} &= p_n - h q_n. \end{aligned} \]
This map is not symplectic, not volume-preserving, and the discrete energy \(H(q_n,p_n)\) typically grows (for small \(h\)), causing a spiralling-out numerical trajectory in phase space.
One variant of the symplectic Euler method (sometimes called “implicit in \(p\)”) is
\[ \begin{aligned} p_{n+1} &= p_n - h q_n,\\ q_{n+1} &= q_n + h p_{n+1}. \end{aligned} \]
This method is symplectic: its Jacobian matrix \(D\Phi_h(q,p)\) satisfies \(J^\top \Omega J = \Omega\) with the standard symplectic matrix \(\Omega = \begin{pmatrix}0 & 1\\ -1 & 0\end{pmatrix}\). Although it does not conserve the original energy exactly, backward error analysis shows that it nearly conserves a nearby (modified) Hamiltonian over very long times.
Qualitatively, instead of a spiralling trajectory, we obtain a nearly closed invariant curve, and the numerical energy oscillates within a small band.
Task marker: in later sections we can formalize the proof of symplecticity for symplectic Euler and perform a rigorous backward error analysis for the oscillator.
Many important classes of differential equations are characterized by particular geometric structures. Table 1.1 summarizes a few, along with typical geometric integrator families designed for them.
| Class of system | Geometric structure preserved by flow | Typical geometric integrators |
|---|---|---|
| Hamiltonian ODEs | Symplectic form, volume, (approx.) Hamiltonian | Symplectic Runge–Kutta, splitting/composition methods, variational integrators |
| Reversible systems | Reversing symmetry \(R\) with \(R^2 = \mathrm{id}\) | Symmetric compositions, time-reversible integrators |
| ODEs on manifolds / Lie groups | Manifold constraint, group structure | Lie-group integrators (RKMK, Crouch–Grossman, Munthe-Kaas methods) |
| Gradient systems | Lyapunov functional decreasing along flow | Energy-diminishing / contractive methods, discrete gradient methods |
| Volume-preserving systems | Volume form | Volume-preserving integrators (e.g. splitting methods on divergence-free fields) |
| Highly oscillatory Hamiltonian systems | Adiabatic invariants, actions | Modulated Fourier integrators, tailored splitting schemes |
The modern theory of geometric numerical integration is primarily developed in:
The present document follows the structure of [BC26] and aims to provide a more expansive, first-principles exposition that:
Task marker: Section 1.2 will formalize flows, invariants, and symmetries using differential geometric tools (Lie derivatives, pullbacks, pushforwards), and will set up the precise definition of symplectic and other geometric structures.
In Section 1.1, we introduced the notion of a flow of an ODE and motivated the idea of preserving geometric properties in numerical integration. We now formalize these concepts using basic differential geometric language. Our goal is to define:
Throughout, we work primarily on \( \mathbb{R}^d \) for simplicity, but the concepts extend naturally to smooth manifolds; see e.g. [AM78].
Let \(M\) be a smooth manifold (for most of this section you may take \(M = \mathbb{R}^d\)). A vector field on \(M\) is a smooth assignment that to each point \(x \in M\) associates a tangent vector \(f(x) \in T_x M\). We write:
\[ f : M \to TM, \qquad x \mapsto f(x) \in T_x M, \]
and in a global coordinate chart on \(\mathbb{R}^d\) we simply view \(f(x)\) as a function \(f : \mathbb{R}^d \to \mathbb{R}^d\).
A smooth curve \(\gamma: I \to M\) is called an integral curve of a vector field \(f\) if it satisfies the ODE:
\[ \dot{\gamma}(t) = f(\gamma(t)), \qquad t \in I \subset \mathbb{R}. \]
Under standard existence and uniqueness conditions (e.g. \(f\) locally Lipschitz), for every point \(x \in M\) there is a maximal interval \(I_x\) and a unique integral curve \(\gamma_x : I_x \to M\) with \(\gamma_x(0) = x\). The flow of \(f\) is then the mapping
\[ \varphi : D \subset \mathbb{R} \times M \to M, \qquad \varphi(t,x) := \gamma_x(t), \]
defined on an open domain \(D\) such that \((0,x) \in D\) for all \(x\). For each fixed \(t\) such that \(\{t\} \times M \subset D\), we write \(\varphi^t : M \to M\), \(\varphi^t(x) = \varphi(t,x)\).
The flow satisfies the group property:
\[ \varphi^0 = \mathrm{id}_M, \qquad \varphi^{t+s} = \varphi^t \circ \varphi^s, \]
whenever the compositions are defined. Conversely, a one-parameter group of diffeomorphisms \(\{\varphi^t\}_{t}\) with \(\varphi^0 = \mathrm{id}\) arises from a unique (time-independent) vector field \(f\) given by the infinitesimal generator:
\[ f(x) = \left.\frac{d}{dt}\right|_{t=0} \varphi^t(x). \]
An autonomous ODE on \(M\) has the form
\[ \dot{y} = f(y), \]
with \(f\) independent of time. Its flow satisfies the group property above. By contrast, a non-autonomous ODE has the form
\[ \dot{y} = f(t, y), \]
and its evolution operator \(\Phi^{t,s}\) mapping \(y(s)\) to \(y(t)\) generally satisfies a two-parameter semigroup property:
\[ \Phi^{t,s} \circ \Phi^{s,r} = \Phi^{t,r}, \qquad \Phi^{s,s} = \mathrm{id}. \]
A common trick, central in geometric integration for non-autonomous systems, is to convert the system into an autonomous one on an extended phase space \(M \times \mathbb{R}\); we return to this in Section 7.
Task marker: later sections will replace the ad-hoc RK2 integrator used for plotting with structure-preserving methods and compare phase portraits.
An important qualitative feature of many dynamical systems is the existence of quantities that do not change along solutions.
Definition 1.1 (First integral). Let \(f\) be a vector field on \(M\) with flow \(\varphi^t\). A smooth function \(I : M \to \mathbb{R}\) is called a first integral or invariant if
\[ I(\varphi^t(x)) = I(x) \quad \text{for all } x \in M \text{ and all } t \text{ in the domain of the flow}. \]
Equivalently, the derivative of \(I\) along integral curves vanishes:
\[ \frac{d}{dt} I(\gamma(t)) = 0 \quad \text{whenever} \quad \dot{\gamma}(t) = f(\gamma(t)). \]
In local coordinates, this is expressed using the gradient:
\[ \frac{d}{dt} I(y(t)) = \nabla I(y(t)) \cdot \dot{y}(t) = \nabla I(y(t)) \cdot f(y(t)) = 0. \]
The differential-geometric formulation uses the Lie derivative of a function with respect to a vector field:
\[ \mathcal{L}_f I = f[I] := dI(f) = \nabla I \cdot f. \]
Then \(I\) is a first integral if and only if
\[ \mathcal{L}_f I = 0. \]
For geometric integrators, the question is: can we construct numerical methods that preserve such invariants exactly? For quadratic invariants, special Runge–Kutta and projection methods can achieve this; see [HLW06, Ch. 3].
The behaviour of a dynamical system is often constrained by its symmetries. These symmetries are crucial in geometric integration because numerical methods can be designed to respect them.
Definition 1.2 (Symmetry of a vector field). Let \(G\) be a Lie group acting smoothly on \(M\) via \(\Phi : G \times M \to M\), \((g,x) \mapsto \Phi(g,x)\). An element \(g \in G\) is a symmetry of the vector field \(f\) if
\[ T_x \Phi(g, \cdot) \, f(x) = f(\Phi(g,x)) \qquad \text{for all } x \in M, \]
i.e. the pushforward of \(f\) under the group action equals \(f\) at the transformed point. This implies that the flow commutes with the group action:
\[ \Phi(g,\varphi^t(x)) = \varphi^t(\Phi(g,x)). \]
Symmetries give rise to conserved quantities via Noether’s theorem when the system derives from a variational principle; see [AM78].
A particularly important type of symmetry is time-reversal.
Definition 1.3 (Reversible system). Let \(R : M \to M\) be an involution, i.e. \(R^2 = \mathrm{id}_M\). An autonomous system \(\dot{y} = f(y)\) is called reversible with respect to \(R\) if
\[ f(Ry) = - T_y R \, f(y). \]
This implies the reversing symmetry relation for the flow:
\[ R \circ \varphi^t \circ R = \varphi^{-t}. \]
For example, Hamiltonian systems with standard symplectic structure are reversible under the involution \(R(q,p) = (q,-p)\), when the Hamiltonian satisfies appropriate conditions. Reversibility implies the existence of symmetric orbits and constrains the long-time behaviour.
Geometric integrators can be constructed to be self-adjoint (time-symmetric), which is the discrete analogue of reversibility and plays a key role in achieving good long-time behaviour; see [HLW06, Ch. 2, 3].
To formalize the preservation of geometric quantities by flows (and later by numerical methods), we briefly recall the notions of pullback, pushforward, and Lie derivative.
Let \(\psi : M \to N\) be a smooth map between manifolds. Its differential at \(x\) is the linear map
\[ T_x \psi : T_x M \to T_{\psi(x)} N, \]
called the pushforward. If \(f\) is a vector field on \(M\), we define the pushforward vector field \(\psi_* f\) on \(\psi(M)\) (when \(\psi\) is a diffeomorphism) by
\[ (\psi_* f)(\psi(x)) = T_x \psi \, f(x). \]
For a smooth function \(g: N \to \mathbb{R}\), the pullback by \(\psi\) is simply composition:
\[ (\psi^* g)(x) := g(\psi(x)). \]
More generally, if \(\alpha\) is a differential 1-form (or higher-rank form) on \(N\), the pullback \(\psi^*\alpha\) is a form on \(M\) defined so that it “undoes” the effect of the pushforward; see [AM78] for the standard construction.
Let \(f\) be a vector field with flow \(\varphi^t\). For a smooth function \(I : M \to \mathbb{R}\), the Lie derivative of \(I\) along \(f\) is defined by
\[ (\mathcal{L}_f I)(x) = \left.\frac{d}{dt}\right|_{t=0} I(\varphi^t(x)) = \nabla I(x)\cdot f(x). \]
For a differential form \(\alpha\), the Lie derivative is defined by
\[ (\mathcal{L}_f \alpha)_x = \left.\frac{d}{dt}\right|_{t=0} (\varphi^t)^* \alpha_{\varphi^t(x)}. \]
A form \(\alpha\) is preserved by the flow if and only if \((\mathcal{L}_f \alpha) = 0\), which in turn is equivalent to \((\varphi^t)^*\alpha = \alpha\) for all \(t\).
These notions are fundamental in defining symplectic and volume-preserving structures and in proving their invariance under flows and numerical methods.
Two particularly important geometric structures in GNI are:
A symplectic form on a manifold \(M\) is a closed, nondegenerate 2-form \(\omega\), that is:
On \(M = \mathbb{R}^{2m}\) with coordinates \(z=(q,p)\), the standard symplectic form is
\[ \omega = \sum_{i=1}^m dq_i \wedge dp_i, \]
or in matrix form, \(\omega(u,v) = u^\top \Omega v\) with the canonical symplectic matrix \(\Omega = \begin{pmatrix} 0 & I \\ -I & 0\end{pmatrix}\).
Given a smooth function \(H : M \to \mathbb{R}\) called the Hamiltonian, the associated Hamiltonian vector field \(X_H\) is uniquely defined by
\[ i_{X_H}\omega = dH, \]
i.e. \(\omega(X_H, \cdot) = dH\). In canonical coordinates this yields Hamilton’s equations:
\[ \dot{q} = \frac{\partial H}{\partial p}, \qquad \dot{p} = - \frac{\partial H}{\partial q}. \]
The flow \(\varphi^t\) of a Hamiltonian vector field is symplectic, meaning
\[ (\varphi^t)^* \omega = \omega, \]
or in coordinates \(\,(D\varphi^t(z))^\top \Omega\, D\varphi^t(z) = \Omega\). This symplecticity implies conservation of phase-space volume (Liouville’s theorem) and many other qualitative properties; see [HLW06, Ch. 1].
On \(M = \mathbb{R}^d\), the standard volume form is \(\mathrm{d}V = dx_1 \wedge \cdots \wedge dx_d\). A diffeomorphism \(\psi : M \to M\) is volume-preserving if
\[ \psi^*(\mathrm{d}V) = \mathrm{d}V \quad \Longleftrightarrow \quad \det D\psi(x) = 1 \quad \forall x. \]
For a vector field \(f\), the Lie derivative condition \(\mathcal{L}_f(\mathrm{d}V) = 0\) is equivalent to \(\nabla \cdot f = 0\): a divergence-free vector field generates a volume-preserving flow.
Hamiltonian vector fields are divergence-free with respect to the canonical volume form, hence Hamiltonian flows are volume-preserving. There exist other classes of volume-preserving flows that are not symplectic, and GNI has developed specialized integrators for them; see [HLW06, Ch. VI].
| Structure | Defining property | Flow condition |
|---|---|---|
| Symplectic | Closed, nondegenerate 2-form \(\omega\) | \((\varphi^t)^*\omega = \omega\) |
| Volume-preserving | Volume form \(\mathrm{d}V\) | \((\varphi^t)^*\mathrm{d}V = \mathrm{d}V\) |
| Reversible | Involution \(R : M \to M\) | \(R \circ \varphi^t \circ R = \varphi^{-t}\) |
| Invariant function | First integral \(I : M \to \mathbb{R}\) | \(I \circ \varphi^t = I\) |
Task marker: later chapters will introduce symplectic integrators whose discrete maps \(\Phi_h\) satisfy the same algebraic symplectic condition, and volume-preserving schemes with \(\det D\Phi_h = 1\).
In this section we have:
Geometric numerical integrators are, in essence, numerical maps \(\Phi_h\) that attempt to mimic these preservation properties at the discrete level. In the next sections, we focus specifically on Hamiltonian systems and symplectic geometry (Chapter 2) and then on symplectic integrators (Chapter 3), following and expanding the treatment in [BC26, HLW06].
The purpose of this section is to explain how flows of ODEs can be understood geometrically: as one-parameter groups of diffeomorphisms generated by vector fields, acting on functions, forms, and structures on phase space.
This viewpoint is essential in geometric numerical integration because it allows us to interpret a numerical method not merely as a recursion formula, but as a discrete flow map whose geometric properties can be studied, compared, and preserved. Our objective is to formalize the following ideas:
Let \(M\) be a smooth manifold and \(f\) a smooth vector field on \(M\). As previously defined, the flow of \(f\) is a smooth mapping
\[ \varphi : D \subset \mathbb{R} \times M \to M, \qquad (t,x) \mapsto \varphi(t,x), \]
where \(D\) is an open neighbourhood of \(\{0\}\times M\) and \(\varphi(0,x)=x\). For each fixed \(t\), define the time-\(t\) map \(\varphi^t(x):=\varphi(t,x)\).
When the system is autonomous, the collection \(\{\varphi^t\}_{t\in\mathbb{R}}\) satisfies the group property:
\[ \varphi^{t+s} = \varphi^t \circ \varphi^s, \qquad \varphi^0 = \mathrm{id}_M. \]
Thus each \(\varphi^t\) is a diffeomorphism and the flow is a one-parameter subgroup of \(\mathrm{Diff}(M)\), the infinite-dimensional Lie group of diffeomorphisms of \(M\). See Abraham–Marsden [AM78] for classical foundations.
Differentiating the flow map with respect to initial conditions yields the variational equation:
\[ \frac{d}{dt}\Big(D\varphi^t(x)\Big) = Df(\varphi^t(x)) \, D\varphi^t(x). \]
This governs the evolution of small perturbations. The matrix \(D\varphi^t(x)\) is the key object in deciding whether the flow preserves structures such as a symplectic form \(\omega\) or a volume form \(\mathrm{d}V\).
A fundamental concept in differential geometry is that a vector field acts as an infinitesimal derivation on functions, tensor fields, and forms. This is formalised via the Lie derivative.
For a smooth function \(g : M \to \mathbb{R}\), the Lie derivative of \(g\) along \(f\) is defined by
\[ (\mathcal{L}_f g)(x) = \left.\frac{d}{dt}\right|_{t=0} g(\varphi^t(x)) = \nabla g(x)\cdot f(x). \]
Thus \(f\) acts as a first-order differential operator. This is the basis for the concept of invariant functions, \(\mathcal{L}_f g = 0\), as covered in §1.2.
For a vector field \(X\) on \(M\), its Lie derivative with respect to \(f\) is the commutator:
\[ \mathcal{L}_f X = [f,X]. \]
The commutator captures algebraic structure of vector fields and is the starting point for the Baker–Campbell–Hausdorff series used in splitting methods.
A central conceptual tool is the formal identification of the flow \(\varphi^t\) of an autonomous vector field \(f\) with the exponential of the derivation \(f\cdot\nabla\). Let
\[ F := f\cdot\nabla = \sum_{i=1}^d f_i(x)\frac{\partial}{\partial x_i}. \]
The flow operator acting on functions satisfies the ODE
\[ \frac{d}{dt}(g \circ \varphi^t) = F(g\circ \varphi^t), \qquad (g\circ\varphi^0)=g. \]
Solving gives the formal exponential representation:
\[ g(\varphi^t(x)) = \big(e^{tF}g\big)(x), \]
where
\[ e^{tF} = \sum_{n=0}^\infty \frac{t^n}{n!} F^n. \]
The operator \(e^{tF}\) is thus the abstract “flow operator.” In geometric numerical integration, splitting methods approximate \(e^{h(A+B)}\) by compositions of \(e^{hA}\) and \(e^{hB}\).
The set of vector fields on a manifold forms a Lie algebra under the commutator \([f,g]\). The exponentials \(e^{tF}\) form local flows, and compositions of flows can be described using the Baker–Campbell–Hausdorff formula:
\[ e^{hA} e^{hB} = \exp\!\left( h(A+B) + \tfrac{h^2}{2}[A,B] + \tfrac{h^3}{12}[A,[A,B]] - \tfrac{h^3}{12}[B,[A,B]] + \cdots \right). \]
In our context, \(A\) and \(B\) are differential operators associated with vector fields (or Hamiltonian vector fields), and the BCH formula gives the exact vector field whose flow is approximated by a composition of subflows.
This is the algebraic backbone of splitting methods, composition methods, and Magnus integrators. See Blanes & Casas [BC26] or Hairer–Lubich–Wanner [HLW06].
For non-autonomous systems,
\[ \dot{y} = f(t,y), \]
the evolution operator \(\Phi^{t,s}\) (flow from time \(s\) to \(t\)) satisfies
\[ \frac{\partial}{\partial t}(\Phi^{t,s}g) = F(t)\,\Phi^{t,s}g, \qquad F(t) := f(t,\cdot)\cdot\nabla, \]
with \(\Phi^{s,s}=\mathrm{id}\). The formal solution is a time-ordered exponential:
\[ \Phi^{t,s} = \mathcal{T}\exp\!\left( \int_s^t F(\tau)\, d\tau \right), \]
where \(\mathcal{T}\) denotes chronological ordering:
\[ \mathcal{T}\exp\!\left(\int_s^t F(\tau)\,d\tau\right) = \mathrm{id} + \int_s^t F(\tau_1)d\tau_1 + \int_s^t\!\!\int_s^{\tau_1} F(\tau_1)F(\tau_2)\,d\tau_2\,d\tau_1 + \cdots \]
This formula is computationally useless in raw form, but it is the foundation for the Magnus expansion, which provides computable closed-form commutator expressions for accurate integrators of non-autonomous linear systems. See [BC26, Ch. 5].
The example below shows how the composition of two simple planar flows behaves differently from their combined flow, illustrating non-commutativity of vector fields and the emergence of commutator terms.
We consider vector fields
\[ f(x,y) = (1,\,0), \qquad g(x,y) = (0,\,x), \]
whose flows do not commute.
The separation of the red and blue points visualises the Lie bracket \([f,g]\neq 0\). In splitting methods of order > 1, these commutators appear explicitly in the error terms and are managed by higher-order compositions.
In this section we have:
These ideas form the algebraic-geometric backbone of modern structure-preserving algorithms. We now move to Chapter 2, where symplectic geometry and Hamiltonian flows are developed rigorously.