Preface for the special issue “Geometric numerical integration, twenty-five years later”

Although for many years the community of numerical analysts has been implicitly assumed that numerical integrators applied to differential equations should preserve as many properties of the system as possible (and thus follow the celebrated motto The Purpose of Computing is Insight, Not Numbers [7]), it was during the 1950s in the field of accelerator physics where this necessity was perceived as more pressing. Thus, in the words of an early pioneer, if one wishes to examine solutions to differential equations, adoption of a ‘Hamiltonian’ or ‘canonical’ integration algorithm would be reassuring [10]. In fact, as pointed out by Forest in his enlightening review [5] referring to several particle accelerator physicists active during that period, they knew that ignoring the symplectic condition of the Poincaré map could lead to spurious effects. This was indeed the main motivation for one of the first symplectic integrators designed as such by R. De Vogelaere in connection with the design of a particle accelerator in the Midwestern United States: De Vogelaere thought that if the numerical integrator possessed the same qualitative properties as the Hamiltonian system, then the effect of discretization error would be no worse than the effect of errors in the model [18]. The schemes designed by De Vogelaere in his unpublished report [2] in 1956 were of first and second order, and it was another accelerator physicist, R. Ruth, who presented in 1983 a third-order symplectic method which turned out to be a splitting method [15]. This paper can be considered as the actual starting point in the systematic exploration of symplectic integrators along several parallel avenues: (i) the use of generating functions in the context of Hamiltonian mechanics to produce appropriate canonical transformations approximating the exact flow in each integration step [1,3]; (ii) the conditions that standard integrators such as the Runge–Kutta methods have to satisfy to be symplectic [9,16,19]; (iii) the design of explicit symplectic methods of order 4 and higher for Hamiltonian systems that can be split into two pieces which can be solved exactly when considered as independent systems [4,14], with the help of the Lie formalism. This notion of splitting was then further clarified and put to use by Suzuki [20] and Yoshida [23] to construct methods of arbitrarily high order.

The early 1990s saw a dramatic increase in the interest and applications of symplectic integrators in several fields, and in particular in dynamical astronomy. Thus, numerical integrations of the planetary equations over very large time intervals carried out by low order symplectic integrators revealed the existence of chaotic phenomena in the Solar System [22]. On the other hand, another symplectic integrator, the Verlet method, has for many years been the method of choice for simulating problems in molecular dynamics [21]. As a matter of fact, neither the accuracy nor the linear stability of these schemes help to explain their success in practice.

It was in this context that a new paradigm in the numerical integration of differential equations was deemed necessary, far beyond the classical consistency/stability approach. It was J.M. Sanz-Serna who first formulated these ideas in his contribution to the conference on the State of the Art in Numerical Analysis organized by the Institute of Mathematics and its Applications at the University of York (UK) in 1996. In fact, ‘Geometric Integration’ is the title of the subsequent paper in the proceedings [17] and the concept that encapsulates this new paradigm.

What, then, is geometric (numerical) integration? Perhaps the simplest definition can be found in the following sentence by McLachlan and Quispel: ‘Geometric integration’ is the term used to describe numerical methods for computing the solution of differential equations, while preserving one or more physical/mathematical properties of the system exactly (i.e. up to round-off error) [13]. Thus, rather than taking primarily into account prerequisites such as consistency and stability, the aim is to reproduce the qualitative features of the solution of the differential equation being discretized, in particular its geometric properties, such as the symplectic character (for Hamiltonian systems), the phase-space volume (for divergence-free vector fields), time-reversal symmetries, first integrals of motion (energy, linear and angular momentum), Casimirs, Lyapunov functions, etc. In these structure-preserving methods one tries to incorporate as many of these properties as possible and, as a result, they exhibit an improved qualitative behaviour. Not only that, it turns out that it typically also allows for a significantly more accurate integration for long-time intervals than with general-purpose methods.

Perhaps the first property one should ask of any integrator when applied to a given dynamical system is preservation of the resulting phase space. This is satisfied, in particular, by any Runge–Kutta method when the phase space is ${\mathbb{R}}^n$, but more sophisticated numerical schemes are required when the system evolves on a more complicated manifold such a Lie group. The time dependent Schrödinger equation describing the evolution of quantum systems is a case in point, and numerical integrators based on the much honoured Magnus expansion [11] are particular instances of such Lie-group methods [8].

The success of symplectic integrators for simulating Hamiltonian systems in classical mechanics led in a natural way to other lines of research. In particular, given the close connection between Hamiltonian and Lagrangian mechanics and the fact that the corresponding equations of motion are derived from Hamilton's principle of critical action, it is reasonable to discretize this variational principle and obtain a map which, by construction, approximates the exact trajectory corresponding to the solution of the (continuous) Euler–Lagrange equations. This, in a nutshell, is the basis of the variational methods [12], which have been widely used for the simulation of mechanical systems.

The year 2021 therefore marks the 25th anniversary of the formal creation of ‘Geometric Numerical Integration’. Along this period, it has been an active and interdisciplinary area of research, and is nowadays the subject of intensive development. Just as an illustration, typing ‘Geometric Numerical Integration’ into Google brings up more than 24 million results, whereas the most authoritative reference in the field [6] counts two editions and more than 5700 citations in Google Scholar. At this stage, it seems natural, in our opinion, to review how structure-preserving methods and techniques provide more accurate descriptions when simulating physical systems of different nature but nevertheless possessing some qualitative properties which must be preserved by discretization.

This Special Issue contains two review articles and five research papers on geometric integration and related issues. Abell and Dragt review the key role played by structure-preserving techniques in accelerator physics, beyond the origina motivation for designing symplectic integrators. Oteo and Ros, for their part, provide a lively portrait of the great mathematician Wilhelm Magnus, the expansion that bears his name and the ulterior development of a well known class of Lie-group methods. The application of Lie-group methods to mechanical systems is the subject of the research paper by Celledoni, Cokaj, Leone, Murari and Owren, whereas Tran and Leok focus on the formulation of variational integrators in the setting of Lagrangian and Hamiltonian partial differential equations.

Problems in celestial mechanics and dynamical astronomy, and in particular the gravitational N-body problem, have traditionally been a strong motivation for designing new types of geometric numerical integrators. The analysis carried out in the paper by Antoñana, Chartier and Murua on the application of majorant series in the N-body problem constitutes the first step towards a new generation of numerical schemes in this setting.

The importance of preservation of properties in the field of stochastic differential equations has been only recently recognized. Cohen and Vilmart address this topic for the class of stochastic Poisson systems.

Last but not least, Böhle, Kuehn and Thalhammer investigate how to efficiently integrate the Kuramoto model for coupled oscillators.