Table of contents

Please see PDF of preface.

Geometric integration is the numerical integration of a differential equation, while preserving one or more of its 'geometric' properties exactly, i.e. to within round-off error. Many of these geometric properties are of crucial importance in physical applications: preservation of energy, momentum, angular momentum, phase-space volume, symmetries, time-reversal symmetry, symplectic structure and dissipation are examples. In this paper we present a survey of geometric numerical integration methods for ordinary differential equations. Our aim has been to make the review of use for both the novice and the more experienced practitioner interested in the new developments and directions of the past decade. To this end, the reader who is interested in reading up on detailed technicalities will be provided with numerous signposts to the relevant literature.

The paper provides an introduction and survey of conservative discretization methods for Hamiltonian partial differential equations. The emphasis is on variational, symplectic and multi-symplectic methods. The derivation of methods as well as some of their fundamental geometric properties are discussed. Basic principles are illustrated by means of examples from wave and fluid dynamics.

This paper is a very personal view of the field of geometric integration in accelerator physics—a field where often work of the highest quality is buried in lost technical notes or even not published; one has only to think of Simon van der Meer Nobel prize work on stochastic cooling—unpublished in any refereed journal. So I reconstructed the relevant history of geometrical integration in accelerator physics as much as I could by talking to collaborators and using my own understanding of the field. The reader should not be too surprised if this account is somewhere between history, science and perhaps even fiction.

In this paper, we review and extend recent research on averaging integrators for multiple time-scale simulation such as are needed for physical N-body problems including molecular dynamics, materials modelling and celestial mechanics. A number of methods have been proposed for direct numerical integration of multiscale problems with special structure, such as the mollified impulse method (Garcia-Archilla, Sanz-Serna and Skeel 1999 SIAM J. Sci. Comput.20 930–63) and the reversible averaging method (Leimkuhler and Reich 2001 J. Comput. Phys.171 95–114). Features of problems of interest, such as thermostatted coarse-grained molecular dynamics, require extension of the standard framework. At the same time, in some applications the computation of averages plays a crucial role, but the available methods have deficiencies in this regard. We demonstrate that a new approach based on the introduction of shadow variables, which mirror physical variables, has promised for broadening the usefulness of multiscale methods and enhancing accuracy of or simplifying computation of averages. The shadow variables must be computed from an auxiliary equation. While a geometric integrator in the extended space is possible, in practice we observe enhanced long-term energy behaviour only through use of a variant of the method which controls drift of the shadow variables using dissipation and sacrifices the formal geometric properties such as time-reversibility and volume preservation in the enlarged phase space, stabilizing the corresponding properties in the physical variables. The method is applied to a gravitational three-body problem as well as a partially thermostatted model problem for a dilute gas of diatomic molecules.

A large number of splitting methods for autonomous separable systems exist in the literature which have been designed for many different structures of the vector field. However, the performance of most of these methods is diminished and their orders of accuracy are frequently reduced when applied to non-autonomous problems. Based on the formal solution obtained from the Magnus series expansion, we show how to modify a standard splitting method for autonomous problems to treat non-autonomous systems with similar or better efficiency. We illustrate this technique to build new fourth- and sixth-order schemes whose performance is then illustrated on several numerical examples.

This paper constructs and analyses an adaptive moving mesh scheme for the numerical simulation of singular PDEs in one or more spatial dimensions. The scheme is based on computing a Legendre transformation from a regular to a spatially non-uniform mesh via the solution of a relaxed form of the Monge–Ampère equation. The method is shown to preserve the inherent scaling properties of the PDE and to identify natural computational coordinates. Numerical examples are presented in one and two dimensions which demonstrate the effectiveness of this approach.

In this paper we develop and analyse new explicit Magnus expansions for the nonlinear equation Y' = A(t, Y)Y defined in a matrix Lie group. In particular, integration methods up to order four are presented in terms of integrals which can be either evaluated exactly or replaced by conveniently adapted quadrature rules. The structure of the algorithm allows us to change the step size and even the order during the integration process, thus improving its efficiency. Several examples are considered, including isospectral flows and highly oscillatory nonlinear differential equations.

If the three moments of inertia are distinct, the solution to the Euler equations for the free rigid body is given in terms of Jacobi elliptic functions. Using the arithmetic–geometric mean algorithm (Abramowitz and Stegun 1992 Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables (New York: Dover)), these functions can be calculated efficiently and accurately. Compared to standard numerical ODE and Lie–Poisson solvers, the overall approach yields a faster and more accurate numerical solution to the Euler equations. This approach is designed for mass asymmetric rigid bodies. In the case of symmetric bodies, the exact solution is available in terms of trigonometric functions, see Dullweber et al (1997 J. Chem. Phys.107 5840–51), Reich (1996 Fields Inst. Commun.10 181–91) and Benettin et al (2001 SIAM J. Sci. Comp.23 1189–203) for details. In this paper, we consider the case of asymmetric rigid bodies subject to external forces. We consider a strategy similar to the symplectic splitting method proposed in Reich (1996 Fields Inst. Commun.10 181–91) and Dullweber et al (1997 J. Chem. Phys.107 5840–51). The method proposed here is time-symmetric. We decompose the vector field of our problem into a free rigid body (FRB) problem and another completely integrable vector field. The FRB problem consists of the Euler equations and a differential equation for the 3 × 3 orientation matrix. The Euler equations are integrated exactly while the matrix equation is approximated using a truncated Magnus series. In our experiments, we observe that the overall numerical solution benefits greatly from the very accurate solution of the Euler equations. We apply the method to the heavy top and the simulation of artificial satellite attitude dynamics.

In this paper we discuss the conservation of wave action under numerical discretization by variational and multisymplectic methods. Both the abstract wave action conservation defined with respect to a smooth, periodic, one-parameter ensemble of flow realizations and the specific wave action based on an approximated and averaged Lagrangian are addressed in the numerical context. It is found that the discrete variational formulation gives rise in a natural way not only to the discrete wave action conservation law, but also to a generalization of the numerical dispersion relation to the case of variable coefficients. Indeed a fully discrete analogue of the modulation equations arises. On the other hand, the multisymplectic framework gives easy access to the conservation law for the general class of multisymplectic Runge–Kutta methods. A numerical experiment confirms conservation of wave action to machine precision and suggests that the solution of the discrete modulation equations approximates the numerical solution to order on intervals of .

In this paper, we analyse a family of exponential integrators for second-order differential equations in which high-frequency oscillations in the solution are generated by a linear part. Conditions are given which guarantee that the integrators allow second-order error bounds independent of the product of the step size with the frequencies. Our convergence analysis generalizes known results on the mollified impulse method by García-Archilla, Sanz-Serna and Skeel (1998, SIAM J. Sci. Comput.30 930–63) and on Gautschi-type exponential integrators (Hairer E, Lubich Ch and Wanner G 2002 Geometric Numerical Integration (Berlin: Springer), Hochbruck M and Lubich Ch 1999 Numer. Math.83 403–26).

The main contribution of this paper is to present a canonical choice of a Hamiltonian theory corresponding to the theory of discrete Lagrangian mechanics. We make use of Lagrange duality and follow a path parallel to that used for construction of the Pontryagin principle in optimal control theory. We use duality results regarding sensitivity and separability to show the relationship between generating functions and symplectic integrators. We also discuss connections to optimal control theory and numerical algorithms.

This paper develops the theory of Abelian Routh reduction for discrete mechanical systems and applies it to the variational integration of mechanical systems with Abelian symmetry. The reduction of variational Runge–Kutta discretizations is considered, as well as the extent to which symmetry reduction and discretization commute. These reduced methods allow the direct simulation of dynamical features such as relative equilibria and relative periodic orbits that can be obscured or difficult to identify in the unreduced dynamics. The methods are demonstrated for the dynamics of an Earth orbiting satellite with a non-spherical J₂ correction, as well as the double spherical pendulum. The J₂ problem is interesting because in the unreduced picture, geometric phases inherent in the model and those due to numerical discretization can be hard to distinguish, but this issue does not appear in the reduced algorithm, where one can directly observe interesting dynamical structures in the reduced phase space (the cotangent bundle of shape space), in which the geometric phases have been removed. The main feature of the double spherical pendulum example is that it has a non-trivial magnetic term in its reduced symplectic form. Our method is still efficient as it can directly handle the essential non-canonical nature of the symplectic structure. In contrast, a traditional symplectic method for canonical systems could require repeated coordinate changes if one is evoking Darboux' theorem to transform the symplectic structure into canonical form, thereby incurring additional computational cost. Our method allows one to design reduced symplectic integrators in a natural way, despite the non-canonical nature of the symplectic structure.

The theory of modified equations (MEs) for discretizations of ODEs is reconsidered. Obstructions to convergence of series expansions of MEs are pinpointed and alternative approaches are presented which provide more accurate descriptions of numerical approximations through MEs. We emphasize how structural assumptions on the ODE can be used to improve estimates. Then we give arguments for a slightly alternative approach based on time-dependent MEs which avoids the asymptotic nature traditionally associated with MEs. Some applications of the theory are also provided.

In this paper we review some group theoretic techniques applied to discretizations of PDEs. Inspired by the recent years active research in Lie group- and exponential-time integrators for differential equations, we will in the first part of the paper present algorithms for computing matrix exponentials based on Fourier transforms on finite groups. As an example, we consider spherically symmetric PDEs, where the discretization preserves the 120 symmetries of the icosahedral group. This motivates the study of spectral element discretizations based on triangular subdivisions. In the second part of the paper, we introduce novel applications of multivariate non-separable Chebyshev polynomials in the construction of spectral element bases on triangular and simplicial sub-domains. These generalized Chebyshev polynomials are intimately connected to the theory of root systems and Weyl groups (used in the classification of semi-simple Lie algebras), and these polynomials share most of the remarkable properties of the classical Chebyshev polynomials, such as near-optimal Lebesgue constants for the interpolation error, the existence of FFT-based algorithms for computing interpolants and pseudo-spectral differentiation and existence of Gaussian integration rules. The two parts of the paper can be read independently.

We derive order conditions for commutator-free Lie group integrators. For certain problems, these schemes can be good alternatives to the Runge–Kutta–Munthe-Kaas schemes, especially when applied to stiff problems or to homogeneous manifolds with large isotropy groups. The order conditions correspond to certain subsets of the set of ordered rooted trees. We discuss ways to select these subsets and their combinatorial properties. We also suggest how the reuse of flow calculations can be included in order to reduce the computational cost. In the case that at most two flow calculations are admitted in each stage, the order conditions simplify substantially. We derive families of fourth-order schemes which effectively use only five flow calculations per step.

In this paper, we study local structures of volume-preserving maps and source-free vector fields, which are defined in the Euclidean n-space Rⁿ with n ⩾ 3. First, we prove that any volume-preserving map, defined in some neighbourhood of the origin, can be represented as a composition of n − 1 essentially two-dimensional area-preserving maps. This result can be viewed as an analogue of the following known fact (Feng and Shang 1995 Volume-preserving algorithms for source-free dynamical systems Numer. Math.71 451–63): any source-free vector field on Rⁿ can be represented as a sum of n − 1 essentially two-dimensional Hamiltonian vector fields. Then, we present a local representation of source-free vector fields under volume-preserving coordinate changes. Finally, we construct a Lie algebra of skew-symmetric tensor potentials of second order associated with source-free vector fields. The Lie algebra turns out to be isomorphic to the Lie algebra of source-free vector fields.

Quantum analysis is reformulated to clarify its essence, namely the invariance of quantum derivative for any choice of definitions of the differential df(A) satisfying the Leibniz rule. This formulation with use of the inner derivation δ_A is convenient to study quantum corrections in contrast to the Feynman operator calculus. The present analysis can also be used to find a general scheme of constructing exponential product formulae of higher order. General recursive schemes are also reviewed with an emphasis to standard symmetric splitting formulae. Multiple integral representations of q-derivatives are derived using such general integral formulae of quantum derivatives as are expressed by hyperoperators. A simple explanation of the connection between quantum derivatives and q-derivative is also given.

The constant-pressure, constant-temperature (NPT) molecular dynamics approach is re-examined from the viewpoint of deriving a new measure-preserving reversible geometric integrator for the equations of motion. The underlying concepts of non-Hamiltonian phase-space analysis, measure-preserving integrators and the symplectic property for Hamiltonian systems are briefly reviewed. In addition, current measure-preserving schemes for the constant-volume, constant-temperature ensemble are also reviewed. A new geometric integrator for the NPT method is presented, is shown to preserve the correct phase-space volume element and is demonstrated to perform well in realistic examples. Finally, a multiple time-step version of the integrator is presented for treating systems with motion on several time scales.