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.