next up previous
Next: Miscellaneous comments concerning dynamical Up: Contents Previous: Exponential divergence in N-body

Integrators

Unless otherwise stated, the term ``system'' in this section will mean an autonomous Hamiltonian system. The discussion of symplectic integration for Hamiltonian systems in this section draws heavily from Sanz-Serna [66], which is a very readable introduction to the subject. Another seminal paper is Channel and Scovel [13].

The question of what kind of integrator to use to integrate N-body systems is one that has had plenty of attention. No discussion of stellar system N-body integrators would be complete without mention of Aarseth's classic integration scheme for collisional systems [1], which includes regularization. Makino [55] derived formulae for optimal order and timestep criteria for Aarseth-type integrators. For collisionless N-body systems, and indeed many non-astronomical Hamiltonian systems, symplectic integrators have enjoyed widespread popularity. Recent years have also seen a growing interest in exactly conserving integrators.

The popularity of symplectic integrators stems from the fact that the time-t flow of any Hamiltonian system is a symplectic or canonical mapping. One consequence of this is that the mapping is volume preserving in phase space. This is a consequence of Liouville's Theorem (see, for example, [56]), which states that the phase-space density of trajectories inside an ensemble remains constant with time. A symplectic integrator also preserves phase space volume to within the machine precision. As a result, symplectic integrators tend to preserve qualitative properties of phase space trajectories: trajectories do not cross, and although energy is not exactly conserved, energy fluctuations are bounded. In fact, it can be shown that the phase space of a linear Hamiltonian problem simply undergoes slight stretches and/or contractions along various axes, so for example a circular trajectory in 2 dimensions becomes an ellipse. Furthermore, if the Hamiltonian function is sufficiently smooth, then the discrete solution produced by a symplectic integrator lies on the exact solution of a problem that has a slightly perturbed Hamiltonian. In general, even for non-smooth Hamiltonians, the discrete solution lies exponentially close in tex2html_wrap_inline2154 to an exact solution to a nearby Hamiltonian problem, where tex2html_wrap_inline2154 is the local error [66].

Non-symplectic integrators generally turn a conservative system into a dissipative one, while a symplectic integrator, on average over time, preserves the conservative nature of the system. Earn and Tremaine [23] point out that roundoff error actually destroys the exact symplecticness of a time- and/or space-continuous system. They demonstrate how one can build a precisely symplectic discrete map over a lattice discretization of space using integer arithmetic, so that if the modeled problem is Hamiltonian, then the discrete model is also exactly Hamiltonian, but with a slightly different Hamiltonian function. This builds confidence in the qualitative nature of the solution being similar to that of the continuous system.

Since symplectic integrators preserve the topological structure of trajectories in phase space, they are probably more reliable for long-term integrations than non-symplectic integrators, since a non-symplectic integrator causes the system to become dissipative, which in turn can badly distort the topological structure of the phase flow over long periods. However, symplectic integrators are often implicit, and require more function evaluations and smaller timesteps to produce the same computational accuracy as compared to standard non-symplectic integrators [66]. Thus, for short-term integrations where accuracy is paramount, a standard non-symplectic integrator may be more efficient.

Another drawback of symplectic integration is that any integration scheme which is symplectic for constant timestep does not remain symplectic if the timestep is changed based solely and explicitly on the state vector [71]. Recall that the time- tex2html_wrap_inline2158 map produced by a symplectic integrator lies exactly (or exponentially close to) the true time- tex2html_wrap_inline2158 flow of a slightly perturbed Hamiltonian problem. However, if tex2html_wrap_inline2158 changes, then the problem we are close to also changes. Thus, if we vary the timestep, it appears that there is no single Hamiltonian problem our symplecticly-generated solution remains close to.

Nonetheless, several researchers have begun to experiment with variable timestep symplectic integrators. Hut, Makino and McMillan [40] introduce an implicit method of symmetrization which is intended to preserve approximate time symmetry. For the Leapfrog integrator, their method reduces to

equation187

where h(y) is a function used to compute the timestep, based solely on the state vector. They tested their symmetrized leapfrog on a Kepler problem with eccentricity e=0.9, using 1 implicit iteration per timestep and tex2html_wrap_inline2168 timesteps per revolution. The errors in energy, semi-major axis, eccentricity, time of perihelion passage and longitude of pericenter all have moderately sized periodic errors, but secular drift is nil or linear with a small slope, compared to the non-symmetrized error of linear or quadratic, respectively. They also demonstrate how a 4th-order generalization of Leapfrog can use their time-symmetrization technique, and use it to perform an impressive integration of an orbit with eccentricity tex2html_wrap_inline2170 . With 3 iterations per timestep of the implicit time-stepping symmetrization equation, they were able to integrate this trajectory using tex2html_wrap_inline2172 timesteps per orbit, with a drift in the above quantities of only one part in tex2html_wrap_inline2174 , over tex2html_wrap_inline2168 revolutions. This seems very impressive. They do not mention any theoretical results on the symplecticness of their method, however.

Skeel and Biesiadecki [72] introduce another method for using different timesteps while preserving the symplectic nature of the integrator. They show that, if it is possible to split the Hamiltonian into independent additive terms, then each term may be integrated with a different, but constant, timestep, and the integration remains exactly symplectic. Each timestep must be an integer multiple of the smallest timestep. For example, if there are only two step sizes and the large timestep is a multiple M of the smaller, then the force at timestep n is computed as

equation193

They show that an artificial splitting of the Hamiltonian into ``quickly'' and ``slowly'' varying components can also utilize their method. They test their method on the Kepler problem, setting the ``quickly'' varying component to be the half of the orbit that is closest to the central mass, and the slowly varying component to be the other half of the orbit. Their results are similar to the second-order results of Hut et al.  [40]. It would be interesting to see if a 4th order version of their method would also give comparable results to Hut et al.  [40]. The error of this method is analyzed more fully in Littell, Skeel, and Zhang [54], where they come to the unsurprising conclusion that the smallest stepsize times the highest frequency of the ``quickly'' changing components must be small. In other words, their method still needs an adequate sample of the highest frequency changes in order to integrate them. They prove the also unsurprising result that, if tex2html_wrap_inline2182 is the maximum frequency in the system, integrated with timestep h, and tex2html_wrap_inline2186 is the maximum frequency of one of the ``slower'' components, then the most efficient timstep for that component is tex2html_wrap_inline2188 .

In addition to symplectic integration, there is an active and quickly blooming body of literature on exactly conserving algorithms. In general, there does not exist an integrator that is both symplectic and exactly energy conserving for Hamiltonian systems [69, 71]. However, since energy is not the only quantity conserved in most systems, there is still the possibility of building a symplectic integrator that conserves other quantities, although there are limits to what the integrator can conserve in comparison to the real system [25]. Furthermore, if the system is not Hamiltonian, then a symplectic integrator provides no added benefit, whereas a conserving integrator may still be of great value. One must be careful, however, how one enforces energy conservation. It is plausible that a method which naïvely enforces global energy conservation may introduce spurious mechanisms of energy transport within the system [69].

Shadwick, Bowman, and Morrison [69] introduce a non-traditional integration scheme which is explicit, and exactly conserving of all quantities one is willing to program it to conserve. The drawback is that the method is not general: one needs to derive a complicated set of explicit formulae for each new integrator-problem pair. One applies the integration method to one's problem on paper, and then derives exactly how the conserved quantities of interest are not conserved. Then one simply adds terms to the integration routine that compensate for the non-conservation. If the standard version of the integrator is explicit, then so is the modified version. Although it is a very tedious process, it seems to have impressive results. The authors present 3 examples: the ``three-wave'' problem, the Lotka-Volterra predator-prey problem, and the Kepler problem, each demonstrated with both Euler's method and a 2nd order predictor-corrector. In the Kepler problem they force conservation of energy, angular momentum, and a quantity known as the Runge-Lenz vector. They note that non-conservation of phase-space volume in Hamiltonian systems can be used as a measure of error, since one can no longer use conservation of other quantities as a measure of error. The authors, however, do not clearly state the limitations of their method. For example, it is not clear precisely why one couldn't use their method to force conservation of phase space volume as well as energy, which we know to be impossible [25, 71].

We will close our discussion of integrators with a remark by Sanz-Serna [66, p.278,] (d is the dimension of configuration space):

Conservation of energy restricts the dynamics of the numerical solution by forcing the computed points to be on the correct (2d-1)-dimensional manifold H= constant, but otherwise poses no restriction to the dynamics: within the manifold the points are free to move anywhere and only motions orthogonal to the manifold are forbidden. When d is large this is clearly a rather weak restriction. On the other hand, symplecticness restricts the dynamics in a more global way: all directions in phase space are taken into account.
Athough I think he makes an extremely interesting point, I note that his remark does not address the possibility that Shadwick et al.  [69] have demonstrated: that energy is not the only quantity that can be conserved. Thus, it is this author's opinion that much interesting work remains in the area of integration of conservative systems.


next up previous
Next: Miscellaneous comments concerning dynamical Up: Contents Previous: Exponential divergence in N-body

Wayne Hayes
Fri Dec 27 17:41:39 EST 1996

Access count (updated once a day) since 1 Jan 1997: 16452