Institute of Fundamental Sciences, Massey University, Palmerston North

For most physical systems, from the very large (galaxies) to the very small (molecules) to the very complicated (the weather), it is not too difficult to write down the equations of motion. These usually take the form of ordinary or partial differential equations. The problem lies in solving these equations: usually it is impossible. In fact, with the new insights brought by chaos theory since the 1960s, we can see that it would be hopeless to expect most equations to have a simple solution; and since the evolution of most systems is exquisitely sensitive to tiny perturbations or errors (the ‘butterfly effect’), one wonders if it is even worth trying.

Therefore in practice, most equations of motion are solved numerically (‘integrated’) on computers. The development of methods for doing this goes back more than 100 years (at first doing the calculations by hand or with mechanical calculators), and huge strides have been made even in recent decades. There are now some very robust and sophisticated computer packages that will try and do a good job with any equation you throw at them. Of course, they do make some errors, and the errors inevitably build up. Again, because of results in chaos theory, this was believed to be unavoidable.

Meanwhile, far away from the world of numerical analysis textbooks, many applied scientists were very happy using methods of their own devising, and getting excellent results. One of the simplest of these is called the leapfrog method. For the equations of motion

dq/dt = p (1a)

dp/dt = f(q) (1b)

where q is a vector representing the positions of the particles in the system (stars, molecules, or whatever), p their momenta (velocities), and f(q) the forces acting on the particles, the leapfrog method looks like this: from the state at step n, namely (q_n, p_n), we calculate the new state as follows:

q_n+1 = q_n + h p_n (2a)

p_n+1 = p_n + h f(q_n+1)(2b)

Here the number h is called the time step, whose value depends a lot on the problem. It might be 10^-15 seconds for molecular dynamics or 100 years for galactic evolution. It turns out that this method gives excellent results for very long simulations precisely because of the q _n+1 in equation (2b): we first update the q’s using the current p, and then update the p’s using the new value of q. (This alternating process gives the leapfrog method its name.) This in turn is related to the structure of the original equations (1). These are not any old equations but come from a law of physics, Newton’s second law force = mass x acceleration. It was gradually realized during the 19^th century that this law embodies extra, hidden, properties called symplectic conservation, in addition to more well-known conservation laws like conservation of energy and angular momentum. By the late 1960s the study of this structure gave rise to a new branch of geometry, symplectic geometry, which eventually, by the late 1980s, began to influence numerical analysts who developed special methods that obey the law of symplectic conservation. The leapfrog method is the simplest method of this class, and still one of the best. At last, after being in popular use for several decades, the surprisingly good performance of the leapfrog method was explained. It’s been possible to show that, although you can’t stop the quantitative errors building up during the simulation, by obeying the law of symplectic conservation, you won’t commit any qualitative errors and so the overall results will be far more reliable.

So the key turned out to be to study the equations of motion (1) in depth and to discover their internal mathematical structure. It is that structure, rather than, say, any particular numbers appearing in the equations, that determines the overall behaviour of the system. What kinds of internal structure might equations of motion have? Some of them were already well-known: they can have conserved quantities like energy, and they can have symmetries like left-right, forwards-backwards, or rotational symmetry. During the 1990s many people began studying properties like these and how to design cheap computer methods which inherit them, and the field of geometric integration began to take shape.

One of the most popular ways of designing geometric integrators is by "splitting". One takes the system and breaks it into a sum of simpler pieces, each piece so simple that it can be solved exactly with no error. This is hardly an original idea: Descartes in his Discourse on Method took as an "inviolable precept" to "divide each of the difficulties which I encountered into as many parts as possible, and as might be required for an easier solution." (And even Henry Ford is supposed to have said that "Nothing is particularly hard if you divide it into small jobs." ) In the case of equation (1), the leapfrog method works by splitting into (1a) and (1b) and solving each part exactly. The only error is due to the interaction between these two parts. This simple idea has turned out to have great power and generality and has been a driving force in our research.

By the late 1990s many structural properties of differential equations were known and many methods found which could preserve at least some of them. An early obstacle had appeared with a theorem that no method can preserve energy and symplecticity simultaneously. That is, there is a law of physics which actual systems obey effortlessly but which we know it is impossible to mimic on a computer. However, a general theory of what all possible such geometric properties are, how they arise in practice, and which of them can be preserved in practical algorithms was (and is) still lacking.

One way forward, pursued in research at Massey, is suggested by looking at the algebraic structure of the systems of the form (1). Given two such systems, there is a way of "multiplying" them to get a third system, and this system has the same structure as the original two. One says the systems form a Lie algebra, introduced by the Norwegian mathematician Sophus Lie in the 1870s and which became one of the fundamental structures used by physics in the 20^th century. We can then go hunting to see what all possible Lie algebras are, and, armed with this list, see if we can find them turning up in applications. This approach is proving very successful as (just as with the example above) the structures are not always immediately obvious unless you know what you are looking for. We are also coming up with some new structures which have not been considered in general dynamical systems theory. Some of these occur in algorithms used in weather forecasting, in the equations of mathematical physics like the Maxwell equations, and in systems with special forms of friction. So there is an interplay between the pure mathematics of classifying Lie algebras, the numerical methods, the applications, and the study of dynamical systems and chaos. Geometric integration is only just beginning.