IT is now well known that numerical methods such as the ordinary Runge--Kutta methods are not ideal for integrating Hamiltonian systems, because Hamiltonian systems are not generic in the set of all dynamical systems, in the sense that they are not structurally stable against non-Hamiltonian perturbations. The numerical approximation to a Hamiltonian system obtained from an ordinary numerical method does introduce a non-Hamiltonian perturbation. This means that a Hamiltonian system integrated using an ordinary numerical method will become a dissipative (non-Hamiltonian) system, with completely different long-term behaviour, since dissipative systems have attractors and Hamiltonian systems do not.
This problem has led to the introduction of methods of symplectic integration for Hamiltonian systems, which do preserve the features of the Hamiltonian structure by arranging that each step of the integration be a canonical or symplectic transformation [Menyuk1984,Feng1986,Sanz-Serna & Vadillo1987,Itoh & Abe1988,Lasagni1988,Sanz-Serna1988,Channell & Scovel1990,Forest & Ruth1990,MacKay1990,Yoshida1990,Auerbach & Friedmann1991,Candy & Rozmus1991,Feng & Qin1991,Miller1991,Marsden et al. 1991,Sanz-Serna & Abia1991,Maclachlan & Atela1992].
A symplectic transformation satisfies
where M is the Jacobian of the map for the integration step, and J is the matrix
with I being the identity matrix. Preservation of the symplectic form is equivalent to preservation of the Poisson bracket operation, and Louiville's theorem is a consequence of it.
Many different symplectic algorithms have been developed and discussed, and many of them are Runge--Kutta methods [Lasagni1988,Sanz-Serna1988,Channell & Scovel1990,Forest & Ruth1990,Yoshida1990,Candy & Rozmus1991,Sanz-Serna & Abia1991,Maclachlan & Atela1992]. Explicit symplectic Runge--Kutta methods have been introduced for separable Hamiltonians of the form . Fourth-order explicit symplectic Runge--Kutta methods for this case are discussed by Candy & Rozmus candy, Forest & Ruth forest, and Maclachlan & Atela maclachlan, and sixth-order and eighth-order methods are developed by Yoshida yoshida. In the special separable case where , and we have a Hamiltonian of potential form, even more accurate methods of fourth-order and fifth-order have been developed by Maclachlan & Atela maclachlan. No explicit symplectic Runge--Kutta methods exist for general Hamiltonians which are not separable. Lasagni lasagni and Sanz-Serna sanz both discovered that the implicit Gauss--Legendre Runge--Kutta methods are symplectic. Maclachlan & Atela maclachlan find these Gauss--Legendre Runge--Kutta methods to be optimal for general Hamiltonians. Thus symplectic integration proves to be a situation where implicit Runge--Kutta methods find a use, despite the computational penalty involved in implementing them compared to explicit methods.
A positive experience with practical use of these methods in a problem from cosmology has been reported by Santillan Iturres et al. santillan. They have used the methods described by Channell & Scovel channell to integrate a rather complex Hamiltonian, discovering a structure (suspected to be there from nonnumerical arguments) which nonsymplectic methods were unable to reveal.
Although symplectic methods of integration are undoubtedly to be preferred in dealing with Hamiltonian systems, it should not be supposed that they solve all the difficulties of integrating them; they are not perfect. Channell & Scovel channell give examples of local structure introduced by the discretization. For another example, integration of an integrable Hamiltonian system, where the solution of Newton's equations is reducible to the solution of a set of simultaneous equations, followed by integration over single variables, and trajectories lie on invariant tori, will cause a nonintegrable perturbation to the system. For a small perturbation, however, such as we should get from a good symplectic integrator, the KAM theorem tells us that most of the invariant tori will survive. Nevertheless, the dynamical behaviour of the symplectic map is qualitatively different to that of the original system, since in addition to invariant tori, the symplectic map will possess island chains surrounded by stochastic layers. Thus the numerical method perturbing the nongeneric integrable system restores genericity.
There is a more important reason why care is needed in integrating Hamiltonian systems, even with symplectic maps, and that is the lack of energy conservation in the map. It would seem to be an obvious goal for a Hamiltonian integration method both to preserve the symplectic structure and to conserve the energy, but it has been shown that this is in general impossible, because the symplectic map with step length h would then have to be the exact time-h map of the original Hamiltonian. Thus a symplectic map which only approximates a Hamiltonian cannot conserve energy [Zhong & Marsden1988,MacKay1990,Marsden et al. 1991]. Algorithms have been given which are energy conserving at the expense of not being symplectic, but for most applications retaining the Hamiltonian structure is more important than energy conservation. Marsden et al. marsden mention an example where using an energy-conserving algorithm to integrate the equations of motion of a rod which can both rotate and vibrate leads to the absurd conclusion that rotation will virtually cease almost immediately in favour of vibration.
In fact, the symplectic map with step length h is the exact time-h map of a time-dependent Hamiltonian with period h, and is near to the time-h map of the original Hamiltonian , so that the quantity , which is measuring energy conservation, is a good guide to the accuracy of the method. In many cases, the lack of energy conservation is not too much of a problem, because if the system is close to being integrable, and has less than two degrees of freedom, there will be invariant tori in the symplectic map which the orbits cannot cross, and so the energy can only undergo bounded oscillations. This is in contrast to integrating the same system with a nonsymplectic method, where there would be no bound on the energy, which could then increase without limit. This is a major advantage of symplectic methods. However, consider a system which has two degrees of freedom. The phase space of the symplectic map is extended compared to that of the original system, so that an N-degree-of-freedom system becomes an -degree-of-freedom map. (The extra degree of freedom comes from t and .) Now in the case where N=2, the original system, if it were near integrable, would have two-dimensional invariant tori acting as boundaries to motion in the three-dimensional energy shell, but in the map the extra degree of freedom would mean that the three-dimensional invariant tori here would no longer be boundaries to motion in the five-dimensional energy shell, so Arnold diffusion would occur. This is a major qualitative difference between the original system and the numerical approximation. It has been shown to occur for two coupled pendulums by Maclachlan & Atela maclachlan, and proves that symplectic methods should not blindly be relied upon to provide predictions of long-time behaviour for Hamiltonian systems.
There is a further point about symplectic maps that affects all numerical methods using floating-point arithmetic, and that is round-off error. Round-off error is a particular problem for Hamiltonian systems, because it introduces non-Hamiltonian perturbations despite the use of symplectic integrators. The fact that symplectic methods do produce behaviour that looks Hamiltonian shows that the non-Hamiltonian perturbations are much smaller than those introduced by nonsymplectic methods. However, it is shown by Earn & Tremaine earn that round-off error does adversely affect the long-term behaviour of Hamiltonian maps like the standard map, by introducing dissipation. To iterate the map they instead use integer arithmetic with Hamiltonian maps on a lattice that they construct to be better and better approximations to the original map as the lattice spacing is decreased. They show that these lattice maps are superior to floating-point maps for Hamiltonian systems. Possibly a combination of the techniques of symplectic methods and lattice maps may lead to the numerical integration of Hamiltonian systems being possible without any non-Hamiltonian perturbations.