In Chapter 2, we introduced the concept of a numerical integrator as a mapping which approximates the flow-map of a given system of differential equations. We have also seen a few instances of how such integrators behave, demonstrating concepts such as convergence and order of accuracy. We observed that the typical picture is a locally accurate approximation that gradually drifts further from the true trajectory (see Fig. 2.3, Fig. 2.5 and the left panel of Fig. 2.7); the rate of drift can be reduced by reducing the stepsize (and thereby also increasing the amount of computational work), but the qualitative picture does not change in any significant way.

What stands out as remarkable, therefore, is the behavior, illustrated in the right panel of Fig. 2.7, of the Euler-B method, which retains bounded trajectories when applied to the harmonic oscillator. In Chapter 2, we provided an explanation for this in the form of a linear stability analysis showing that certain methods, including Störmer–Verlet and Euler-B, have eigenvalues on the unit circle when applied to the harmonic oscillator (or any other oscillatory linear system), if the stepsize is below some threshold value. The Euler-B and Störmer-Verlet methods (among others) possess a strong asymptotic stability property for linear systems.

It is interesting to note that a related long-term stability property extends to nonlinear models. If we apply, for example, the Störmer–Verlet methods to the Lennard-Jones oscillator, we obtain the results illustrated in Fig. 4.1 (compare with Fig. 2.3 and Fig. 2.5).

In Chapter 4, we introduced several first- and second-order symplectic integration methods for Hamiltonian systems. In this chapter, we will discuss the construction of “higher-order” symplectic methods (with order greater than two), focusing in particular on those types of schemes that have been found to be most useful for practical computations. In traditional practice, higher-order integrators are employed for solving problems with relatively smooth solutions, such as gravitational simulations (solar system simulations, satellite trajectories). They are also traditionally used for many types of computations when very high accuracy (for example near the rounding error of the computer) is desired.

As we have seen in Chapter 2, the appropriateness of a given numerical method for a given computational task is a complicated issue. In some cases, the principles of geometric integration are in contradiction with the demand for high accuracy. If the purpose of simulation is to reconstruct, as exactly as possible, a particular trajectory segment, it may not matter what sort of qualitative features the integrator possesses: the efficiency of the integration method in terms of solution error per unit work is instead of paramount importance. Since the development of symplectic integrators adds a number of additional constraints on the design of the method, such schemes typically sacrifice something in efficiency compared with their nonsymplectic counterparts at similar accuracy, for example requiring an extra force evaluation or two at each timestep. Thus the problem of correctly determining the precise entry point and time instant that a space probe arrives at the Martian atmosphere is a task best handled by a standard integration method, for example, a high-order multistep integrator (for example, Diva) or explicit Runge–Kutta method (for example, RKSUITE).

In this chapter, we discuss formulation issues and symplectic integration methods for simulating the motion of a rigid body. Rigid bodies arise frequently in engineering, chemistry, and physics. For example, they occur in molecular simulation when the flexibility of small polyatomic units such as the water molecule, or CH4 is ignored. Cogwheels, space vehicles, and the planets are some other objects that are commonly modeled by rigid bodies.

Even in the absence of external applied forces, any rigid body more complicated than a uniform sphere will exhibit complicated motion, as defined by the moments of inertia of the body. A hint of the potential complexity of the motion is provided by the classic illustration using a hardbound book, which typically has three unequal moments of inertia I1 < I2 < I3 with I1 corresponding to an axis drawn along the binding, I2 to an axis across the cover, and I3 to an axis through the pages of the book (see Fig. 8.1). As the book is tossed up and spinning around each of the axes, the following dynamics are observed: around the first and third axes, the motion combines a stable periodic rotation with the rising and falling motion due to gravity, whereas the rotation with respect to the middle axis is much more complicated. (See, for example, for more explanation.) (It helps to place a rubber band around the book's cover to keep it closed while conducting experiments.)

Developing a method to simulate general rigid body motions, especially for long-term integration, proves an interesting and challenging task. The first issue we must confront is the selection of a set of coordinates that describe body orientation and spatial position.

About geometric integration

This book is about simulating dynamical systems, especially conservative systems such as arise in celestial mechanics and molecular models. We think of the integrator as the beating heart of any dynamical simulation, the scheme which replaces a differential equation in continuous time by a difference equation defining approximate snapshots of the solution at discrete timesteps. As computers grow in power, approximate solutions are computed over ever-longer time intervals, and the integrator may be iterated many millions or even billions of times; in such cases, the qualitative properties of the integrator itself can become critical to the success of a simulation. Geometric integrators are methods that exactly (i.e. up to rounding errors) conserve qualitative properties associated to the solutions of the dynamical system under study.

The increase in the use of simulation in applications has mirrored rising interest in the theory of dynamical systems. Many of the recent developments in mathematics have followed from the appreciation of the fundamentally chaotic nature of physical systems, a consequence of nonlinearities present in even the simplest useful models. In a chaotic system the individual trajectories are by definition inherently unpredictable in the exact sense: solutions depend sensitively on the initial data. In some ways, this observation has limited the scope and usefulness of results obtainable from mathematical theory. Most of the common techniques rely on local approximation and perturbation expansions, methods best suited for understanding problems which are “almost linear,” while the new mathematics that would be needed to answer even the most basic questions regarding chaotic systems is still in its infancy.

We have seen in the previous chapter that integrators preserving symplectic structure and/or first integrals can often be constructed in a straightforward way. In this chapter, we consider the properties of those methods and the implications for long-term simulations.

The traditional approach of numerical analysis generally assumes that the purpose of simulation is the faithful reproduction of a particular solution or trajectory, but individual trajectories typically are not of primary interest in most modern, scientific research; rather, the scientist typically treats the trajectory as a particular realization of a fundamentally stochastic evolution modelling in some way the myriad undetermined perturbations present in a “real-world” environment. It was the important discovery of Lorenz that differential equations can exhibit a chaotic solution behavior that includes an essentially stochastic or “random” component. The scientist views the model being analyzed as representative of a class of nearby models based on parameters which are typically only empirically (and approximately) determined. Furthermore, exact initial conditions are also typically not available. Some classical examples of such a scenario are molecular dynamics and numerical weather prediction.

It is now apparent that most modern large-scale simulations are conducted with timesteps and time intervals such that the numerical solution cannot be thought of as close to any particular model solution. The purpose of wedding the development of integrators to the standard axiomatic principle of timestepping – that one is attempting to approximate a particular trajectory – is thus called into question. Although high accuracy often is not needed in nonlinear dynamics computations, we must recognize certain important constraints imposed by the laws of nature.