DESIGN OF A MECHANICAL N-BODY INTEGRATOR

BEN LEIMKUHLER

Abstract. This document describes a gravitational N-body integration algorithm incorporating the following features: (1) it is exactly (neglecting rounding errors) time-reversible, (2) angular and linear momentum are conserved, (3) smooth switching functions are used to split potential terms into local and weak parts so that weaker long-range forces are evaluated relatively rarely and close interactions are identified, (4) close approaches between bodies are resolved accurately, using an efficient integration method, (5) the stepsize varies automatically based on an appropriate Sundman time reparameterization. Although this method is formally second order, the most intensive computations (the close approach dynamics) are executed at higher order, thus improving the overall accuracy of the scheme.

1. Introduction

The N-body problem of celestial mechanics is described by a Hamiltonian

ℋ (p, q) = T (p) + f(q),

where T : ℝ^{3N} → ℝ and f : ℝ^{3N} → ℝ are smooth kinetic and potential energy functions defined in terms of the individual momenta p_{1}, p_{2}, …, p_{N} and positions q_{1}, q_{2}, …, q_{N} of the bodies by

T (p) = T (p_{1}, p_{2}, …, p_{N}) := ∑_{i=1}^{N} \frac{|p_{i}|^{2}}{2m_{i}} ,

and

f(q) = f(q_{1}, q_{2}, …, q_{N}) := − ∑_{ij} \frac{Gm_{i}m_{j}}{r_{ij}} , r_{ij} = |q_{i} − q_{j}|,

The idea for this article grew out of discussions with Martin Lo, Jet Propulsion Laboratory, Pasadena, CA. This work was performed while the author was a member at the Mathematical Sciences Research Institute, Berkeley, CA, in the fall of 1998. Research at MSRI is suppored in part by NSF grant DMS-9701755. The author was also supported in part by NSF Grant #DMS 9627330.

where G is the gravitational constant and m_{i} represents the mass of the ith body.

The N-body problem is the seminal problem of dynamical systems, and study of its solutions, ergodic properties, and topological structure continues to generate a great deal of mathematical and physical interest (see [14] and other recent references therein). The numerical simulation problem also remains of terrific importance. This problem is solved routinely in studies of steller and planetary dynamics [29, 21, 1], and related problems arise in quasiclassical studies of atomic systems [22, 5, 25].

While many numerical schemes for the N-body problem have been developed over the years [7, 1], these codes may exhibit deficiencies in very long time integration or in scattering studies. One approach to improving the qualitative behavior of numerical simulation methods is to incorporate some of the many geometric properties of the phase flow, such as time-reversal symmetry, symplectic structure, and integrals such as angular momentum. Methods which preserve symmetries and invariants are sometimes referred to as geometric or mechanical integrators [23, 16, 12]. Symplectic and symmetric algorithms for smooth N-body trajectories (without close approaches) have been successfully used for long-term simulations of the solar system [29].

Close approaches introduce both theoretical and computational challenges for the Coulombic N-body problem. In the approach described here, we assume that there are potentially a large number of bodies, but at any given time only a few bodies are engaged in very close approaches. None of the geometric algorithms can be expected to preserve structure or provide substantial efficiency improvements under frequent discontinuous switchings (such as may be used for close encounters). This problem was discovered and explained by Calvo and Sanz-Serna [3] in the context of variable stepsize symplectic integration: frequent changes in the timestepping map lead to a deteriation of the numerical stability of geometric methods. For variable stepsizes, it turns out to be more practical to use a Sundman or Poincaré-type time-transformation to rescale the vector field or Hamiltonian in such a way that (i) geometric structure is preserved, and (ii) the new system can be solved with a fixed timestep. Articles by Stoffer [26], Hut et al [9], and the current author with collaborators [10, 8] developed the use of time-reversible variable stepsizes for N-body problems. More recently several articles have indicated how regularization methods can be incorporated in a reversible or symplectic framework [19, 13, 11].

Another approach to variable stepsize integration, based on a hierarchical splitting, was developed in [24], and we exploit this idea in a somewhat different way below, to define a smooth switching of the close approaching bodies. We then separate the weaker interactions by splitting, and identify potential close approaches using Verlet lists [28], partitioning the bodies into small groups. The multiple timestepping framework of molecular dynamics (r-RESPA [27]) allows some control of the amount of computation performed between evaluations of the long-range forces. In order to correct for the increasing sensitivity of the dynamics of closeapproaching bodies, we incorporate a Sundman time transformation using the reversible-adaptive framework developed previously by the author[13].

For problems with only two-body close approaches, much of the computational work will be in the resolution of these small systems. The switches somewhat complicate the dynamics of the two-body pairs (they are no longer pure Kepler problems). We describe an efficient technique for recovering their dynamics, based on a reduced coordinate set, a regularization, and an efficient implementation of a higher-order implicit Gauss-Legendre integrator.

When three or more body interactions in close approach must be taken into consideration (in the localized subsystem), a system of relative variables can be introduced and a regularization developed in terms of twobody pairs. In the standard approach based on Kustaanheimo-Stiefel transformation, the coupling of variables in the regularized Hamiltonian makes high-accuracy geometric integrators very costly. We suggest instead an approach based on relative coordinates, splitting, and two-body pair integration.

Our integrator could easily be adapted to the Coulombic problems arising in quasiclassical simulations of atomic systems[22, 5, 25, 13], to systems with multiple fixed bodies, or to systems subject to applied fields or arbitrary perturbing potentials, using additional splittings.

This article describes the design of the integrator. A companion article will address the implementation issues associated to this method.

2. Smooth Switches, Splitting, Multiple-Timestepping and Reversible Adaptive Time-Stepping

We use the term smooth switch to describe a real function χ which smoothly passes from one to zero in some finite subinterval of ℝ_{+}:

χ ∈ C^{k}(ℝ_{+}, [0, 1]), χ(r) = 1, r < r_{−}, χ(r) = 0, r ≥ r_{+}, χ′(r) ≤ 0.

Piecewise polynomial switches of any desired smoothness can be constructed as follows. First define a polynomial

p(x) = (−1)^{k+1} \frac{(x^{2} − h^{2})^{k+1}}{h^{2k+2}} .

This function is easily seen to satisfy the conditions p(0) = 1, p^{(l)}(h) = p^{(l)}(−h) = 0, l = 0, 1, …, k. Normalizing the integral of this polynomial on the interval [−h, h] and subtracting from unity results in a new polynomial

q(x) = 1 − \frac{∫_{−h}^{x} p(x)dx}{∫_{−h}^{h} p(x)dx}.

Now for given r_{−} and r_{+} we set h = (r_{+} − r_{−})/2, r̄ = (r_{+} + r_{−})/2, then the function

χ(r) = \{ \begin{matrix}1 & r ≤ r_{−} \\ q(r − r̄), & r_{−} ≤ r < r_{+} \\ 0, & r ≥ r_{+} \end{matrix}

is a C^{k} switch (See Fig. 1). The fact that we can choose χ to be piecewise polynomial will be found to have some positive ramifications for the efficiency of our code.

Figure 1. A smooth switch.

The purpose of the smooth switch is to limit the high-accuracy resolution of trajectories to close approaching bodies, while reducing the number of force evaluations for weak interactions. We decompose each pair interaction in the gravitational potential using smooth switches:^{1}

f = \underbrace{− ∑_{ij} χ(r_{ij})\frac{Gm_{i}m_{j}}{r_{ij}}}_{f_{loc}} + \underbrace{− ∑_{ij}(1 − χ(r_{ij}))\frac{Gm_{i}m_{j}}{r_{ij}}}_{f_{weak}} .

In a two-body problem, there are three distinct regimes in the dynamics of f_{loc}. For separations r ≤ r_{−}, the dynamics are pure Keplerian. For r ≥ r_{+}, the bodies are in free motion. For r ∈ (r_{−}, r_{+}), the motion is determined by the (somewhat artificial) dynamics of the switch, integrable dynamics in a rapidly decaying force field.^{2}

To design a one-step method, we utilize a Hamiltonian splitting of the form ℋ = ℋ _{loc} + f_{weak}, where

ℋ _{loc} := T + f_{loc}.

Let Φ_{t,ℋ} represent the flow map on Hamiltonian ℋ , parameterized by t. We construct a splitting method using multiple time-stepping [27] with m inner substeps:

Φ_{Δt,ℋ} ≈ Φ̂_{Δt,ℋ} := Φ_{Δt,f_{weak}} [Φ̂_{\frac{1}{m} Δt,ℋ _{loc}}]^{m} .

(In most cases, we would take m = 1, but the incorporation of multiple timestepping adds some flexibility.) Here Φ̂_{Δt,ℋ _{loc}} is a time-reversible approximation method for solving the localized gravitational problem, as described in the next section.

The method Φ̂_{Δt,ℋ} is a first order fixed stepsize integrator, and it is not time-reversible. We use this method as the basis for constructing a second order variable stepsize reversible integrator according to the prescription of [8]:

1For simplicity, we here suppose that the parameters r_{−}, r_{+} of each switch are all identical. In many cases, it may be more practical to choose the parameters of the switch dependent on the masses of the particles of the pair. In the related classical atomic model, the switch might also depend on the product of the charges of the particle pair. This topic will be addressed in the companion article mentioned in the introduction.

2Some peculiar effects might be observed if the switched fields were integrated for any length of time. These effects will be eliminated in the limit of a small splitting timestep.

( \begin{matrix}q & n+1/2 \\ p & n+1/2 \end{matrix} ) = Φ̂_{\frac{1}{2} Δt_{n},ℋ} ( \begin{matrix}q & n \\ p & n \end{matrix} ) ,
\frac{1}{Δt_{n}} + \frac{1}{Δt_{n+1}} = \frac{2}{g_{n+1/2}Δt},
( \begin{matrix}q & n+1 \\ p & n+1 \end{matrix} ) = Φ̂_{\frac{1}{2} Δt_{n+1},ℋ}^{∗} ( \begin{matrix}q & n+1/2 \\ p & n+1/2 \end{matrix} ) .

where

g_{n+1/2} = g(q^{n+1/2}, p^{n+1/2})

and g(q, p) is the time-reparameterization function, a smooth, positive, scalar-valued function which is invariant under p → −p (see below). The adjoint method needed above is easily computed since Φ̂_{Δt,ℋ} is the composition of symmetric maps:

Φ̂_{Δt,ℋ}^{∗} = [Φ̂_{\frac{1}{m} Δt,ℋ _{loc}}]^{m} Φ_{Δt,f_{weak}}.

Although the method is formally second order (symmetric methods have even order [6]), we expect to utilize a higher-order integrator for the local interaction dynamics, as described in the next and following sections. Viewed as a second order method, the scheme is therefore expected to have a relatively small leading error constant compared to a method such as Störmer-Verlet which does nothing special to integrate close approaching bodies. If desired, a third order splitting method could be used in place of Φ̂_{Δt,ℋ} , resulting in a fourth order method overall (see [17] for a discussion of how to construct such splittings).

We anticipate that the considerations for the choice of time-transformation will be similar to those discussed in [2, 13]. In the latter reference, it was suggested to use a control of the form

g_{min} = \frac{1}{1 + r̄^{−3/2}}

where r̄ represents the smallest pair separation. In some cases, it may be desirable to have a smooth control, in which case we could use

g_{sm} = \frac{1}{1 + (∑_{i<j} r_{ij}^{−3m/2})^{1/m}} ,

where m is a positive integer [2]).

3. Integration of ℋ _{loc}.

We now describe the integrator for ℋ _{loc}, a key component of the overall scheme. We will assume throughout this section that the solution is desired on the time interval [0, Δt], that initial positions q^{0} and momenta p^{0} are provided and that new values q^{1} and p^{1} are to be computed.

Because of the use of switching functions, a relatively small number of the interaction terms are active in ℋ _{loc} at any given timestep; most of the particles simply drift linearly with fixed momenta. We would like to introduce the concept of a local interaction graph G based the positive potential interactions among the bodies. This graph will be used to limit the local computation (in the relatively expensive regularization variables) to minimal subgroups of bodies. However, there is some question about how to define G. Clearly the vertices of G are the body indices, but what should we take for the edges?

Given vertices i and j, it is evidently not sufficient to add the edge īj only when ∥q_{i}^{0} − q_{j}^{0}∥ ≤ r_{+}, since, during the timestep, the two bodies may move closer to each other and thus have a close interaction by the end of the timestep. Moreover, the one-step methods we propose for integrating the few-body problems will evaluate the forces at an intermediate point in the subinterval. Still it is not sufficient to simply link i and j based on the initial, final and intermediate steps, since it is necessary to know these interactions at the start of the timestep (after all, the whole purpose of constructing the graph is to restrict the computation to local groups!). To define the interaction graph, we therefore need a very efficient method to predict, based only on the initial positions and velocities, which particles have the potential to have a local interaction during the current timestep. The same technique used to construct Verlet lists[28] in molecular dynamics can be used for this purpose.

To identify potential close approaches, we proceed as follows for each of the i bodies (being careful to work with squared norms to avoid computation of extra square roots). We iterate over the bodies with index j > i. Define

Δq = q_{i}^{0} − q_{j}^{0}, Δp = p_{i}^{0} − p_{j}^{0},

We first compute σ^{0} = |Δq|^{2} Clearly if σ^{0} < r_{+}^{2}, we should add edge īj to G, but because our method is heuristic (based only on linear trajectories), we include a relaxation factor β > 1 and check instead σ^{0} < βr_{+}^{2}. Next, we compute the time of close approach of the linear trajectories through (q_{i}^{0}, p_{i}^{0}) and (q_{i}^{1}, p_{i}^{1}):

t^{∗} = −\frac{Δq · Δp}{|Δp|^{2}}.

If t^{∗} ∈ [0, Δt], then we compute the squared separation between these two trajectories at t^{∗}, σ^{∗} = |Δq + t^{∗}Δp|^{2}, and check if σ^{∗} < βr_{+}^{2}. Finally, we compute the separation at the right endpoint of the time interval, σ^{1} = |Δq + ΔtΔp|^{2}, adding the edge if this is less than βr_{+}^{2}.

After updating the interaction graph, we subdivide the indices 1, 2, …, N into L disjoint connected components Ω_{1}, Ω_{2}, … Ω_{L}. These components will generally consist of one or several indices. The cost of resolving the close approaches rises very rapidly with the number of bodies, while the attainable accuracy diminishes, but close approaches of three or more bodies are rare. In the next two sections, we describe the integration of the few-body problems, first for only two-body approaches, then for the more general situation.

4. Fast Integrator for the Switched Problem: Two-Body Case

The computation of the close encounters, is typically the most costly part of the calculation, so substantial algorithmic effort is warranted in the interest of overall efficiency. In this section we describe a technique for two-body approaches based on a reducing coordinate transformation, regularization, and the Gauss-Legendre family of higher-order implicit Runge-Kutta methods. In many cases, for example in tracking the motion of a body on a tour of the solar system, this is all that will be needed to integrate the close approaches.

Any Hamiltonian system of the form

H_{2bdy} = \frac{1}{2m_{1}} |p_{1}|^{2} + \frac{1}{2m_{2}} |p_{2}|^{2} + ϕ(|q_{1} − q_{2}|)

is integrable. For Kepler’s problem, an elegant solution was known at least to Gauss, and is detailed in [15]. A number of papers have considered efficient numerical methods for solving the associated nonlinear equation (see e.g. [19, 20]). A method avoiding all transcendental functions was given in [13]. Solutions for special cases of the two-body problem with various power potentials are discussed in [4], however for arbitrary ϕ, the computation of the associated quadrature and nonlinear equations becomes complicated by potential singularities and a sign change of the integrand, and an involved procedure is then needed to perform the computations, including a special quadrature algorithm and case-dependent changes of variables.

For the general case, it proves to be economical and robust to use coordinate transformations to reduce the dynamics of H_{2bdy} to a regularized one-degree of freedom problem which is then solved by an efficient numerical timestepping method. The original variables can be recovered by solving a linear system in tandem with the reduced problem, and the result is a fast 2-body solver.

We begin by changing to relative variables,

B ( \begin{matrix}q & 1 \\ q & 2 \end{matrix} ) = ( \begin{matrix}q & cm \\ δ \end{matrix} ) ,  B = ( 1  −1^{\frac{m_{1}}{m_{1}+m_{2}}  \frac{m_{2}}{m_{1}+m_{2}}} )

and

( \begin{matrix}p & 1 \\ p & 2 \end{matrix} ) = B^{T} ( \begin{matrix}p & cm \\ π \end{matrix} ) ,

taking the Hamiltonian to

H_{2bdy} = \frac{1}{2M} |p_{cm}|^{2} + \frac{1}{2m̄}|π|^{2} + ϕ(|δ|),

with M = m_{1} + m_{2}, 1/m̄ = 1/m_{1} + 1/m_{2}. We now concentrate on the decoupled reduced problem with

H_{red} = \frac{1}{2m̄}|π|^{2} + ϕ(|δ|),

The differential equations on this Hamiltonian are

(1)  \frac{d}{dt}δ = π/m̄,
(2)  \frac{d}{dt}π = −\frac{ϕ′(|δ|)}{|δ|}δ.

Now introduce new variables ξ = |δ|^{2}, η = δ · π. Equations for these variables are readily derived through differentiation, the use of (1)-(2), and the constant of motion H_{red} = E_{0}:

(3)  \frac{d}{dt}ξ = \frac{2}{m̄}η
(4)  \frac{d}{dt}η = −ξ^{1/2}ϕ′(ξ^{1/2}) + 2E_{0} − 2ϕ(ξ^{1/2}).

Again, this is an integrable system (it is Hamiltonian with a single degree-of-freedom). To solve it, we could proceed to define a potential ϕ̃(ξ) = ∫ (−ξ^{1/2}ϕ′(ξ^{1/2}) + 2E_{0} − 2ϕ(ξ^{1/2}))dξ, then recover the solution by piecing together solutions of the separable scalar ODEs

\frac{d}{dt}ξ = ±\sqrt{m̄(Ẽ_{0} − ϕ̃(ξ))},

where Ẽ_{0} is the energy in the associated Hamiltonian. Numerically speaking, this is rather challenging. Not only do we have to design an accurate numerical quadrature and a switching based on sign of ξ̇, but the solution is then only obtained in implicit form, defined by a certain nonlinear equation. The entire problem is further complicated by the singularity in the integrand. Despite the problems, a viable method for the general case could undoubtedly be derived based on quadrature, but we are skeptical that it would be as efficient and robust as the alternative we describe in the sequel.

Before proceeding, we must address the regularity of (3)-(4). Does the problem admit solutions which make arbitrarily close approaches to the singularity? The key to answering this question is to use the fact that the switch becomes one at sufficiently small interbody separation. Thus we need only answer the question for ϕ(ξ) = −γ/ξ, γ a positive constant. The Hamiltonian corresponding to (3)-(4) becomes in the close approach

H̃ = \frac{1}{m̄}η^{2} − 2γξ^{1/2} − 2E_{0}ξ

In the cartesian variable, it follows by a short calculation that

H̃ = (cos^{2} θ − 1)|π|^{2}|δ|^{2} = −|π × δ|^{2},

where θ is the angle between π and δ, thus the energy in the reduced problem is the negative square of the angular momentum of the 2-body system. Except along a zero angular momentum orbit H̃ < 0, and it follows that ξ is bounded strictly from zero as well. However, even if the total angular momentum were known to be nonzero, there is nothing to prevent the “local” angular momentum from approaching zero, so arbitrarily close approaches are possible.

The first step in solving (3)-(4) is thus to perform a regularization, which we implement at the level of the differential equations. We make a coordinate transformation replacing ξ^{1/2} by a new variable ζ, then we perform a Sundman time transformation, dt/dτ = ζ. The result of these changes of variables is a new system of the form

(5)  \frac{dζ}{dτ} = \frac{1}{m̄}η
(6)  \frac{dη}{dτ} = −ζ^{2}ϕ′(ζ) + 2E_{0}ζ − 2ζϕ(ζ),
(7)  \frac{dt}{dτ} = ζ

This system has the following properties: (i) the singularity is removed, so integration can in theory be performed even through head-on collisions, and (ii) the right hand side of (6) is piecewise polynomial, thus very efficient to compute. The solution of the regularized differential equations requires an accurate and stable geometric integration method such as the efficiently implemented implicit Gauss-Legendre method described below.

In tandem with the numerical solution of (3)-(4), the cartesian variable should be recovered from the same numerical method applied to the system

\frac{d}{dt}δ = π/m̄,
\frac{d}{dt}π = −\frac{ϕ′(ξ^{1/2})}{ξ^{1/2}}δ.

This computation can also be accelerated, as described below.

4.1. Accelerated RK Method for the Regularized Reduced Problem. We next describe an efficient implementation of the implicit Gauss-Legendre Runge-Kutta scheme of s-stages (order 2s) for the equations of motion (5)-(7) and the corresponding system for δ and π. An alternative would of course be to solve the equations (5)-(7) accurately in quadrature, including a nonlinear implicit equation for the fictive time interval; for the regularized equations, the quadrature approach is greatly improved. We favor the use of an RK method since it avoids (1) switching of the integrand and (2) square roots.

Because this is a such a critical component of the overall method, we discuss the computation in detail. After rewriting all equations in terms of the fictive time variable, the equations become:

(8)  \frac{d}{dτ} ζ = \frac{1}{m̄}η,
(9)  \frac{d}{dτ} η = g(ζ),
(10)  \frac{d}{dτ} t = ζ,
(11)  \frac{d}{dτ} δ = \frac{ζ}{m̄}π,
(12)  \frac{d}{dτ} π = α(ζ)δ.

where α and g are scalar functions. The solution of (8)-(10) decouples from (11)-(12). A curious twist is that the timestep Δt is known in advance, but the corresponding fictive timestep Δτ must be computed. In order to recover the maximal order of accuracy, the functions α and g must be C^{2s+1}, thus we need to use a C^{2s+1} switch.

Let the coefficients of the s-stage Gauss-Legendre RK method be described in the standard Butcher notation (see [6]) by a matrix A = (a_{ij})_{i,j=1}^{s}, and vectors b = (b_{1}, .., b_{s})^{T} and c = (c_{1}, .., c_{s})^{T}. Applied to u̇ = f(t, u), the method computes u_{n+1} from u_{n} in a timestep of size Δt according to the formulas

(13)  U_{i} = u_{n} + Δt ∑_{j=1}^{s} a_{ij}f(t_{n} + c_{j}Δt, U_{j}), i = 1, …, s

(14) u_{n+1} = u_{n} + Δt ∑_{j=1}^{s} b_{j}f(t_{n} + c_{j}Δt, U_{j}),

After application of (13)-(14) to the scalar equations (8)-(9) and simplification, there results a system of s equations for the values G_{i} of g at the stages,

G_{i} = g(Z_{i})
Z_{i} = ζ_{n} + \frac{Δτc_{i}}{m̄} η_{n} + \frac{Δτ^{2}}{m̄} ∑_{j=1}^{s} â_{ij}G_{j}, i = 1, …, s,

where the coefficients â_{ij} are those of the square of A, i.e. (â_{ij}) = (a_{ij})^{2}. Along with this system for Z_{i}, we must solve the equation for the timestep Δτ,

Δt = Δτ ∑_{j=1}^{s} b_{j}Z_{j}.

The Jacobian matrix of the resulting s + 1-dimensional nonlinear system for Z_{1}, …, Z_{s} and Δτ can be written in the block form

J = [ \begin{matrix}Ĵ & v \\ u^{T} & w \end{matrix} ] .

where elements of the s × s matrix Ĵ are

J_{ij} = δ_{ij} − \frac{Δτ^{2}}{m̄} â_{ij}g′(Z_{j})

where δ_{ij} is one if i = j and otherwise zero, v has elements

v_{i} = −\frac{c_{i}}{m̄}η_{n} − \frac{2Δτ}{m̄} ∑_{j=1}^{s} â_{ij}g(Z_{j}),

u = −Δτb, and w = − ∑_{j=1}^{s} b_{j}Z_{j}.

The nonlinear equations should be comparably easy to resolve in the regularized variables, but having a good predictor is very important for the ultimate success of the method. Previous step values along the numerical solution (or internal stages at the previous step of a Runge-Kutta method) are sometimes used for prediction, but we prefer to keep all computations localized to within the current timestep due to the potential changes in the local interaction graph from step to step. One possibility is to use a high-order Taylor series approximation (obtained through differentiation of the differential equation). This is especially viable because the function g is assumed to be C^{2s+1} and piecewise polynomial.

Ideally, the Jacobian need only be calculated once at the beginning of the step (from the predicted values), factored, and held fixed during a modified Newton iteration. In order to preserve the geometric invariants, it is essential that the equations be iterated to convergence (to rounding error or a small multiple thereof); this condition probably precludes the use of anything but a modified Newton iteration. In some cases, a damped or global Newton iteration may be needed to ensure convergence.

Our next task is to solve (11)-(12), for which only the internal (stage) variables of the previous calculation are actually needed. Observe that this system decouples into three linear single degree-of-freedom systems. Applying (13)-(14) to (11)-(12) therefore results in three decoupled linear systems of dimension s, each having the same coefficient matrix.

The cost of this procedure is moderate. For example, if s = 4, which yields an eighth order inner integrator, each pair integration will require the calculation and factoring of 5 × 5 and 4 × 4 matrices, and back-solves for a handful of right hand sides.

5. Higher-Body Close Encounters

We next consider numerical integrators for the case that a component of the local interaction graph can include three or more bodies. The equations in this case are no longer integrable, so the only possibility is to derive an appropriate numerical discretization scheme. In general, our problem is to integrate numerically an d-body problem with Hamiltonian of the form

\frac{1}{2} ∑_{i=1}^{d} \frac{|p_{i}|^{2}}{m_{i}} + ∑_{j=1}^{d−1} ∑_{k=j+1}^{d} ϕ(|q_{i} − q_{j}|)

Due to our use of splitting and time-transformation, we can restrict attention to the case where the number of bodies d is small.

As is the usual practice, we first introduce relative coordinates for the pair separations,

δ_{1} = q_{1} − q_{2}, δ_{2} = q_{1} − q_{3}, …, δ_{d−1} = q_{1} − q_{d},
δ_{d} = q_{2} − q_{3}, …, δ_{2d−3} = q_{2} − q_{d},
…, δ_{\frac{1}{2} (d^{2}−d)} = q_{d−1} − q_{d}.

Together with the center of mass, these variables, though degenerate (D = (d^{2} − d)/2 ≥ d, d ≥ 3), are adequate to describe any configuration. Canonical equations can be developed in terms of the relative variables, by introducing suitable momenta π_{i} canonically conjugate to the δ_{i}. Neglecting the center of mass motion, the Hamiltonian becomes

H_{rel} = \frac{1}{2} ∑_{i=1}^{D} \frac{|π_{i}|^{2}}{2m̄_{i}} + ∑_{j≠k} σ_{jk}π_{j} · π_{k} + ∑_{i=1}^{D} ϕ(|δ_{i}|),

where m̄_{i} = 1/(1/m_{α} + 1/m_{β}), where m_{α} and m_{β} are the masses of the bodies of the ith pair, and σ_{jk} is either zero (if the pairs on which δ_{j} and δ_{k} are based involve no common body) or else ± the reciprocal of the common mass of the jth and kth pairs of bodies, depending on the order of the differences.

Now there are several ways to proceed. A traditional course, described concisely in [18], is to apply the Kustaanheimo-Stiefel regularizing transformation to all the variables. The kinetic energy metric is then position dependent and rather complicated. The nonseparable nature of the Hamiltonian precludes the use of (known classes of) explicit geometric integrators, leaving us essentially with implicie methods such as the Gauss-Legendre Runge-Kutta methods or Lobatto IIIA-B partitioned Runge-Kutta pairs[6], but in this case, there are no simple reductions of the numbers of variables such as we found in the previous section. For example, a fourth-order method for the three-body problem would require the solution of a nonlinear system of dimension 2 × 24 = 48 at each step. This is clearly unacceptable in long-term simulations, for which many such nonlinear systems may need to be solved.

Instead, we will look for a splitting of H_{rel} that requires only the solution of two-body problems at each step. The numerical method of the last section can then be used to solve the two-body problems.

One such splitting would break H_{rel} into two parts,

H_{rel} = H_{1} + H_{2},

where

H_{1} = ∑_{i=1}^{D} \frac{|π_{i}|^{2}}{2m̄_{i}} + ∑_{i=1}^{D} ϕ(|δ_{i}|),

and

H_{2} = ∑_{j≠k} σ_{jk}π_{j} · π_{k}.

H_{1} consists of D decoupled two-body problems, whereas H_{2} is only P dependent. Hence both terms are integrable.

An alternative splitting method would divide H into D parts,

H_{rel} = G_{1} + G_{2} + … G_{M}

where

G_{i} = \frac{|π_{i}|^{2}}{2m̄_{i}} + ϕ(|δ_{i}|) + \frac{1}{2} ∑_{k≠i} σ_{ik}π_{i} · π_{k}.

Each of the D terms can be seen to be integrable, since only the ith momentum vector changes during the integration of G_{i}. In the case of a 3-body problem, this splitting has three terms. For higher order, the number would rapidly increase.

Both splittings are easily seen to conserve angular momentum, since each of the terms has this feature and the discretization is constructed by concatenation of the exact flows (or conserving approximations of the exact flows) on each term.

To obtain higher-order methods, we suggest to use the concatenation technique [30, 17, 23]. For the alternative decomposition, second and fourth order splitting methods are also easily derived.

References

[1] Aarseth, S.J., Direct methods for N-body simulation, in Multiple Time Scales, eds. J.U. Brackbill and B.I. Cohen, Academic Press, NY, 1985.

[2] Bond, S., and Leimkuhler, B., Time-transformations for reversible variable stepsize integration, Numerical Algorithms, to appear.

[3] Calvo, M.P. and Sanz-Serna, J.M., The development of variable stepsize symplectic integrators with applications to the two-body problem, SIAM J. Sci. Comput. 14, 936, 1993.

[4] Goldstein, H., Classical Mechanics, 2nd Ed., Addison Wesley, New York, 1980.

[5] Gu, Y., and Yuan, J.-M., Chaotic scattering of electrons with He^{+}, Phys. Rev. A 47, R2442, 1993.

[6] Hairer, E., Norsett, S.P., and Wanner, G., Solving Ordinary Differential Equations, I, 2nd Ed., Springer, 1994.

[7] Heggie, D.C., The N-body problem in stellar dynamics, in Long-Term Dynamical Behavior of Natural and Artificial N-Body Systems, Ed. A.E. Roy, NATO ASI Series, Kluwer, Dordnecht, 1988.

[8] Holder, T., Leimkuhler, B., and Reich, S., Explicit variable stepsize and timereversible integration, preprint, 1998.

[9] Hut, P., Makino, J., McMillan, S., Building a better leapfrog, The Astrophysical Journal 443, L93-L96, 1995.

[10] Huang, W., and Leimkuhler, B., The Adaptive Verlet Method, SIAM J. Sci. Comput. 18, 239, 1997.

[11] Hut,P., Funato, Y., Kokubo, E., Makino, J., McMillan, S., Time symmetrization meta-algorithms, in Computational Astrophysics, Proc. of the 12th ’Kingston meeting’ on Theoretical Astrophysics, Eds. D.A. Clarke and M.J. West, ASP Conf. Series. to appear.

[12] Budd, C.J. and Iserles, A., Geometric integration: numerical solution of differential equations on manifolds, Numerical Analysis Reports, DAMTP 1998/NA10, University of Cambridge, 1998.

[13] Leimkuhler, B., Reversible adaptive regularization: perturbed Kepler motion and classical atomic trajectories, Trans. Roy Soc., to appear.

[14] Klein, M., and Knauf, A., Classical planar scattering by Coulombic potentials, Lecture Notes in Physics m13, Springer-Verlag, Berlin, 1992.

[15] Landau, L.D., and Lifshitz, E.,M., Mechanics (Course in Theoretical Physics, Vol. 1), 3rd Ed., Pergamon Press, Oxford, 1976.

[16] Marsden, J., Lectures on Mechanics, Cambridge University Press, Cambridge, [17] McLachlan, R., On the numerical integration of ordinary differential equations

by symmetric composition methods, SIAM J. Sci. Comp., 16, 151, 1995.

[18] Mikkola, S., A practical and regular formulation of the N-body equations, Mon. Not. R. Astr. Soc., f̄ 215, 171–177, 1985.

[19] Mikkola, S., A cubic approximation for Kepler’s equation. Celestial Mechanics and Dynamical Astronomy 40, 303-312, 1987.

[20] Nijerhuis, A., Solving Kepler’s equation with high efficiency and accuracy, Celestial Mechanics and Dynamical Astronomy 51, 319-330, 1991.

[21] S.J. Peale, Some unsolved problems in evolutionary dynamics, Celestial Mechanics and Dynamical Astronomy 46, 253-275, 1989.

[22] Rzazewski, K., Lewenstein, M., and Salieres, P., Multielectron stabilization of atoms in a laser field: classical perspective, Phys. Rev. A 49, 1196, 1994.

[23] Sanz-Serna, J.M. and Calvo, M.P., Numerical Hamiltonian Problems, Chapman and Hall, 1994.

[24] Skeel, R.D. and Biesidecki, J., Symplectic integration with variable stepsize, Annals of Numerical Mathematics, 1, 191, 1994.

[25] Steckel, J., and Jaffe, C., The bifurcations of the Langmuir orbit in the twoelectron atom, Hamiltonian systems with three or more degrees of freedom, NATO ASI Series, Kluwer, Dordnecht, 1998.

[26] Stoffer, D., Variable steps for reversible integration methods, Computing 55, 1, 1995.

[27] Tuckerman, M. and Berne, B.J., J. Chem. Phys. 97, 1990 (1992).

[28] Verlet, L., Computer experiments on classical fluids I. Thermodynamical properties of Lennard-Jones molecules, Physical Review 159, 98 (1967).

[29] Wisdom, J. and Holman, M., Symplectic maps for the N-body problem, The Astronomical Journal 102, 1528, 1991.

[30] Yoshida, H., Construction of higher order symplectic integrators, Phys. Lett. A 150, 1990.

Department of Mathematics, 405 Snow Hall, University of Kansas, Lawrence, KS 66045 USA

  E-mail address: leimkuhl@math.ukans.edu