Symplectic Integrators: T + V Revisited and Round-Off Reduced

TL;DR

This paper revisits T+V symplectic integrators, enhancing their efficiency and accuracy with higher-order algorithms, variable step sizes, and a novel method to reduce round-off errors, making them competitive with MVS methods for long-term simulations.

Contribution

The paper introduces new 6th-order T+V algorithms with force gradients and symplectic correctors, and a simple modification to significantly reduce round-off errors in T+V integrators.

Findings

6th-order T+V algorithms are competitive with MVS in some cases

Variable step sizes improve T+V efficiency

Round-off errors can be reduced by several orders of magnitude

Abstract

Symplectic integrators separate a problem into parts that can be solved in isolation, alternately advancing these sub-problems to approximate the evolution of the complete system. Problems with a single, dominant mass can use mixed-variable symplectic (MVS) integrators that separate the problem into Keplerian motion of satellites about the primary, and satellite-satellite interactions. Here, we examine T+V algorithms where the problem is separated into kinetic T and potential energy V terms. T+V integrators are typically less efficient than MVS algorithms. This difference is reduced by using different step sizes for primary-satellite and satellite-satellite interactions. The T+V method is improved further using 4th and 6th-order algorithms that include force gradients and symplectic correctors. We describe three 6th-order algorithms, containing 2 or 3 force evaluations per step, that are competitive with MVS in some cases. Round-off errors for T+V integrators can be reduced by several orders of magnitude, at almost no computational cost, using a simple modification that keeps track of accumulated changes in the coordinates and momenta. This makes T+V algorithms desirable for long-term, high-accuracy calculations.

Equations

$$\frac{dq}{dt}$$

$$q(\tau )=e^{\tau F}q(0)=\left(1+\tau F+\frac{{\tau }^2}{2}{F}^2+\cdots \right)q(0)$$

$$q(\tau )=e^{\tau (A+B)}q(0)$$

$$e^{\tau (A+B)}$$

$$e^{\tau A}{e}^{\tau B}=\exp \left\{\tau (A+B)+\frac{{\tau }^2}{2}[A,B]+\cdots \right\}$$

$$e^A{e}^B$$

$$S_2=e^{\tau B/2}{e}^{\tau A}{e}^{\tau B/2}=\exp \left\{\tau (A+B)+O({\tau }^3)+\cdots \right\}$$

$$S_4$$

$$S_{4G}$$

$$e^C{e}^K{e}^{-C}=\exp \left(K+[C,K]+\frac{1}{2}[C,C,K]+\frac{1}{6}[C,C,C,K]+\cdots \right)$$

$$S_{4C}$$

$$\exp \left(\frac{{\tau }^2}{12}[A,B]\right)\simeq$$

$$\exp \left(\frac{\tau A}{4}\right)\exp \left(\frac{\tau B}{6}\right)\exp \left(-\frac{\tau A}{4}\right)\exp \left(-\frac{\tau B}{6}\right)\exp \left(-\frac{\tau A}{4}\right)$$

$$\exp \left(-\frac{\tau B}{6}\right)\exp \left(\frac{\tau A}{4}\right)\exp \left(\frac{\tau B}{6}\right)$$

$$S_{6A}$$

$$a$$

$$g$$

$$h$$

$$j$$

$$k$$

$$l$$

$$S_{6B}$$

$$30a^4-90a^3+78a^2-26a+3=0$$

$$b$$

$$g$$

$$k$$

$$l$$

$$h$$

$$S_{6C}$$

$$15a^5-60a^4+90a^3-60a^2+18a-2=0$$

$$b$$

$$g$$

$$j$$

$$l$$

$$h$$

$$e^C$$

$$\exp (X_i)=\exp ({\alpha }_i\tau A)\exp ({\beta }_i\tau B)$$

$$j$$

Full text

††pubyear: 2018††pagerange: Symplectic Integrators: T + V Revisited and Round-Off Reduced–Symplectic Integrators: T + V Revisited and Round-Off Reduced

1 Introduction

Symplectic integrators are often favored for integrating the orbits of planetary and satellite systems. These integrators have two important advantages over most other algorithms: (i) they show good long-term energy conservation properties, and (ii) they are computationally efficient for problems that involve a dominant central mass.

The usual strategy when devising a symplectic integrator is to divide the problem of interest into 2 or more parts that can each be solved easily in isolation. The algorithm then advances the parts alternately in a series of sub-steps that combine to approximate the evolution of the whole system. One possibility, that can be applied to any non-dissipative N-body problem, is to divide the problem into parts involving the kinetic energy $T$ and potential energy $V$ respectively (Gladman et al., 1991). This kind of algorithm is often referred to as a “T+V” integrator as a result. It is straightforward to advance each of these parts separately using Hamilton’s equations. Either the momenta stay fixed while the coordinates change, or vice versa.

T+V integrators are commonly used in some fields of physics, such as classical and quantum mechanical molecular dynamics problems (Bandrauk & Shen, 1993; Forbert & Chin, 2001; Omelyan et al., 2002a). However, their use for studying planetary and satellite systems in celestial mechanics has largely been superseded by “mixed-variable” symplectic (MVS) integrators (Wisdom & Holman, 1991). MVS algorithms take advantage of the presence of a dominant central body in the system, and split the problem into (i) Keplerian motion about the central body, and (ii) direct and indirect perturbations due to interactions between the less massive objects. This separation makes it possible to use substantially longer steps for the same level of accuracy as a T+V integrator.

One way to improve the efficiency of T+V algorithms is to use different step sizes for motion around the central body, and for other interactions. Since interactions between satellites are typically much weaker than the force from the central body, the satellite interactions can be advanced using a longer time step without compromising the accuracy of the algorithm. We will explore this possibility in this paper, and show that T+V integrators can be competitive with MVS algorithms in some circumstances, and that a simple modification makes them much more resistant to round-off errors.

The rest of this paper is organized as follows. Section 2 describes how symplectic integrators can be constructed, and describes some examples that are accurate to second, fourth, and sixth-order in the step size. In Section 3, we apply these integrators to the 2-body Kepler problem, while Section 4 looks at systems containing more than two bodies. In Section 5, we show how T+V integrators can be modified to greatly reduce round-off errors. Section 6 contains a summary.

2 Devising Symplectic Integrators

Consider a system of $N$ bodies moving in 3 spatial dimensions described by a Hamiltonian $H_{F}$. Using Hamilton’s equations, the evolution of any quantity $q$ can be expressed as

[TABLE]

where $t$ is the time, and ${\bf x}_{i}$ and ${\bf p}_{i}$ are coordinates and momenta of body $i$. Here $F$ is an operator that encapsulates the evolution of the system and depends on $H_{F}$.

The value of $q$ after one time step $\tau$ can be expressed as

[TABLE]

For most $N$-body systems, $F$ is too complicated to allow an exact solution. However, as we noted in the introduction, we can usually separate the problem into 2 or more parts that are easy to solve. For example if we separate the Hamiltonian so that $H_{F}=H_{A}+H_{B}$, with corresponding operators $A$ and $B$, then $q(\tau)$ is given by

[TABLE]

where

[TABLE]

We note that the operators $A$ and $B$ do not commute in general, so that $AB\neq BA$.

A simple integrator consists of just two steps in which the system is advanced for one time step under each of the two sub-problems separately. Using the operators $A$ and $B$, we can describe this integrator as

[TABLE]

where the square brackets indicate a commutator defined by $[A,B]=AB-BA$.

The algorithm described by Eqn. 5 differs from the true system, described by Eqn. LABEL:eq_full_system by a factor $O(\tau^{2})$, so the integrator is accurate to first order in the step size.

2.1 Second and 4th-Order Integrators

We can devise higher-order integrators using the Baker-Campbell-Hausdorff formula (Yoshida, 1990), which gives the product of exponential operators in terms of a series of commutators:

[TABLE]

where we use Yoshida’s compact commutator notation $[A,B,C]\equiv[A,[B,C]]$ and $[A,B,C,D]\equiv[A,[B,[C,D]]]$ etc.

The well-known second-order leapfrog integrator consists of 3 sub-steps:

[TABLE]

Using this algorithm, the sub-system $H_{B}$ is advanced for half a time step, then the sub-system $H_{A}$ is advanced for a full step, followed by another half step for sub-system $H_{B}$.

Leapfrog is an example of a time-symmetric integrator where the algorithm is unchanged if the sequence of sub-steps is reversed. Symmetric algorithms have the advantage that they contain no error terms with even powers of the timestep $\tau$ (Yoshida, 1990), and we will only consider symmetric integrators from now on.

Forest & Ruth (1990) described a 4th-order algorithm that contains 7 sub-steps:

[TABLE]

where $a=1/(4-2^{4/3})=0.6756\ldots$. Although this integrator is accurate to 4th-order in the step size, it does not perform as well as one might expect since the sub-steps are large (and some travel backwards in time), so the coefficients of the leading error terms are large (Chambers, 2003).

2.2 Force Gradients

One way to reduce the size of the sub-steps is to use terms with “force gradients” (Omelyan et al., 2002b). Many $N$-body problems contain only quadratic momentum terms. If we split the Hamiltonian into two parts such that one of them $H_{A}$ contains all of the momentum terms, and the other $H_{B}$ depends only on coordinates, then the operator $[B,B,A]$ will depend on coordinates only. A sub-step consisting of the $[B,B,A]$ operator can be advanced easily according to Hamilton’s equations since the coordinates remain fixed for the duration of that sub-step.

A simple 4th-order integrator that includes a force-gradient term is

[TABLE]

Note that the $[B,B,A]$ force gradient term can be advanced at the same time as the $2\tau B/3$ term since both these sub-steps only alter the momenta. We will examine the form of $[B,B,A]$ in more detail in Sections 3 and 4.

2.3 Symplectic Correctors

Another strategy for developing integrators is to use a “symplectic corrector” (Wisdom et al., 1996). Correctors are a sequence of sub-steps applied before and after an integration step such that the terms at the end of one step exactly cancel those applied at the start of the next step. In practice, this means the corrector or its inverse needs only be applied at the start of an integration and immediately before output is required. Correctors typically consist of a complicated series of sub-steps, and can be computationally expensive. However, they are worth the cost if output is not needed very often.

Wisdom et al. (1996) provide a useful formula for the effect of a corrector $C$ on an unmodified integrator or “kernel” $K$:

[TABLE]

Using this formula, we can obtain the following 4th-order integrator that includes a force gradient and a corrector:

[TABLE]

In general, we can only get an approximate expression for $[A,B]$ accurate to some order in the step size. The following corrector is accurate to $O(\tau^{3})$, which is sufficient to make the whole integrator 4th-order accurate:

[TABLE]

The inverse corrector reverses the order of these sub-steps as well as their signs.

2.4 Sixth-Order Algorithms

Sixth order symplectic integrators typically require many substeps. For example, Yoshida (1990) gives three examples that each contain 15 substeps. The number of sub-steps can be reduced substantially when force gradients and symplectic correctors are included. In this subsection, we will examine three sixth-order algorithms, one that consists of 5 sub-steps, and two that contain 7 sub-steps. All of these use ordinary force-gradient terms as well as a higher-order derivative of the force gradient that is proportional to the operator $[B,B,A,A,B]$. For systems with quadratic momenta, this operator consists only of coordinates, and can be advanced in the same way as $[B,B,A]$.

The following sixth-order algorithm contains 5 sub-steps plus a corrector:

[TABLE]

where

[TABLE]

This integrator is accurate to sixth order in the step size provided that the Hamiltonian for the system contains only quadratic momenta, in which case $[B,B,B,A]=[A,B,B,B,A]=[B,B,B,B,A]=0$. Note that unlike the previous algorithms, the kernel of this integrator begins with a sub-step involving operator $A$ rather than $B$. There is no equivalent 6th-order algorithm beginning with $B$ that has sub-steps with real coefficients.

Two other 6th-order algorithms will be useful later due to the special form of the corrector in each case. These additional constraints on the correctors mean that each integrator kernel requires 7 sub-steps instead of 5.

The first 7-step integrator has the following form

[TABLE]

where $a$ is the smaller real root of

[TABLE]

which gives $a\simeq 0.577,953,138,043,435\cdots$, and

[TABLE]

Another 7-step integrator has the following form (note the different location of the $g$ and $h$ terms):

[TABLE]

where $a$ is the real root of

[TABLE]

which gives $a\simeq 0.567,040,718,865,478\cdots$, and

[TABLE]

For each of the 6th-order algorithms, we need a corrector composed of alternating $A$ and $B$ operators that is accurate to $O(\tau^{5})$. These correctors contain only even powers of $\tau$. We note that a corrector of the following form can be adapted to give any corrector with even powers of $\tau$ up to 4th-order by choosing appropriate values of the free parameters $\alpha_{1,2}$ and $\beta_{1,2}$:

[TABLE]

where

[TABLE]

These correctors contain 32 terms and are expensive to calculate, but the resulting integrators remain efficient as long as output is not required too often. (It is possible that the desired correctors could be constructed with fewer terms, but we do not explore this here.)

The values of the $\alpha$ and $\beta$ are related to the coefficients $j$, $k$ and $l$ of the $[A,B]$, $[A,A,A,B]$ and $[A,B,B,A]$ terms in the corrector by the following equations:

[TABLE]

Note that the sign of the $l$ term can be changed by reversing the order of $A$ and $B$ in Eqn. 22

3 The Kepler Problem

In this section, we apply the integration algorithms described above to the Kepler problem in which a test particle orbits a point mass $M$. The Hamiltonian in this case is

[TABLE]

where $r$ is the distance between the objects.

We separate the Hamiltonian into kinetic and potential energy terms:

[TABLE]

With this separation, $H_{A}$ can be advanced by keeping the momenta fixed and linearly increasing the positions at constant velocity. Similarly, $H_{B}$ can be advanced by keeping the coordinates fixed and applying a constant acceleration to the velocities. The force gradient operators used by some of the higher-order integrators are particularly simple in this case:

[TABLE]

and it is straightforward to advance the system under either of these operators.

Figure 1 shows the performance of some of the T+V integrators described above when applied to the Kepler problem. The figure shows the rms energy error as a function of computer time for integrations lasting 100,000 orbital periods for an orbit with an eccentricity of 0.1. From the slope of the curves, we see that the 2nd, 4th and 6th-order algorithms are behaving as expected. For example, decreasing the step size by a factor of 2 doubles the computation time while roughly reducing the error by factors of 4, 16 and 64 respectively. The error for the 6th-order integrator levels off at very small step sizes due to round-off error.

4 The $N$-Body Problem

In this section, we consider the evolution of a system containing several bodies orbiting a dominant central mass. We will work with democratic heliocentric coordinates (Duncan et al., 1998), which consist of coordinates ${\bf X}$ with respect to the central body, and momenta ${\bf P}$ with respect to the center of mass. Using these coordinates, the Hamiltonian can be split into the following parts:

[TABLE]

where $N$ is the number of satellites orbiting the central body which has index 0.

Note that $H_{I}$ is typically much smaller than $H_{A}$ and $H_{B}$ for small satellite-to-primary mass ratios. We can make use of this difference to produce a modified leapfrog algorithm with the following steps

• •

Apply a corrector $C_{I}$.

• •

Advance $H_{I}$ for $\tau/2$

• •

Do the following $M$ times:

• –

Advance $H_{B}$ for $\tau/(2M)$

• –

Advance $H_{A}$ for $\tau/M$

• –

Advance $H_{B}$ for $\tau/(2M)$

• •

Advance $H_{I}$ for $\tau/2$

• •

Apply an inverse corrector $-C_{I}$.

where $M$ is an integer and $\tau$ is the step size. When $M=1$, we have the usual leapfrog algorithm. More efficient algorithms will use $M>1$. The speed-up can be significant when $N$ is larger than a few, since advancing $H_{I}$ requires $O(N^{2})$ operations compared with $O(N)$ for $H_{A}$ and $H_{B}$.

Note that we have included a corrector $C_{I}$ to eliminate the leading error term that contains a single factor of $I$, where

[TABLE]

We do not need to include $B$ here since $[B,I]=0$. This corrector doesn’t do much to improve the performance of leapfrog, but it will improve some of the higher-order algorithms discussed below, so we include it here for consistency.

We can use a similar procedure for the conventional 4th-order integrator $S_{4}$. The 4th-order gradient integrator $S_{4G}$ requires minor modification since the gradient is more complicated than for the Kepler problem of the previous section. The 4th-order gradient algorithm has the following steps:

• •

Apply a corrector $C_{I}$.

• •

Advance $H_{I}$ for $\tau/2$

• •

Do the following $M$ times:

• –

Advance $H_{B}$ for $\tau/(6M)$

• –

Advance $H_{A}$ for $\tau/(2M)$

• –

Advance $H_{B}$ for $2\tau/(3M)$ and $[B,B,A]$ for $-\tau^{3}/(72M)$

• –

Advance $H_{A}$ for $\tau/(2M)$

• –

Advance $H_{B}$ for $\tau/(6M)$

• •

Advance $H_{I}$ for $\tau/2$

• •

Apply an inverse corrector $-C_{I}$.

where the gradient operator is now given by

[TABLE]

where

[TABLE]

Figure 2 shows the performance of the T+V integrators when integrating the Sun and the 8 planets of the Solar System. The integrations last for 100,000 years. The step size for the direct terms $H_{I}$ is 1.8 days, while cases for multiple values of $M$ are shown. The figure also shows integrations using a 2nd-order leapfrog MVS integrator, and a 2nd-order MVS leapfrog with a corrector that is included with the Mercury $N$-body integration package (Chambers, 1999, 2010). For the MVS integrators, we advance the Keplerian orbits using the accurate and efficient routine described by Rein & Tamayo (2015).

The classical second and fourth order integrators $S_{2}$ and $S_{4}$ perform poorly compared to the MVS integrator. However, the other fourth order algorithms and the sixth order integrator $S_{\rm 6A}$ are more efficient (requiring less CPU time for a given accuracy) than the standard MVS algorithm. MVS with a corrector is more efficient than all the T+V integrators, although it is only slightly better than $S_{\rm 6A}$. Thus, the usual great speed disadvantage of T+V compared to MVS integrators can be substantially reduced by using these 4th and 6th order algorithms.

The accuracy of the second and fourth order integrators varies with the step size roughly as expected. Doubling the step size roughly doubles the integration cost, while reducing the error by factors of about 4 and 16 for the second and fourth order cases respectively. However, $S_{\rm 6A}$ doesn’t perform as well as expected. In fact, the slope of the curve in Figure 2 shows that $S_{\rm 6A}$ mostly behaves as a 2nd order integrator rather than 6th order. In the following subsections we explore the reason for this behavior and show how it can be fixed.

4.1 Improved Sixth-Order Integrators

We can get a sense of the problem with the 6th-order algorithm $S_{\rm 6A}$ by examining the simpler 4th-order algorithm, $S_{\rm 4C}$, described by Eqn. 11, with the inclusion of direct terms that are advanced by the operator $I$. This algorithm also uses a corrector. Consider the simplest case in which the step size for $I$ is the same as for $A$ and $B$. The integrator kernel $K$ is then

[TABLE]

We then want to apply a corrector $C$ outside the kernel. The most general corrector available at second order in the step size has the form

[TABLE]

where $j$ and $k$ are constants. Note that $[B,I]=0$ since both the corresponding pieces of the Hamiltonian depend only on the coordinates. Using this general corrector, a complete step of the integrator is

[TABLE]

The $[A,A,B]$ and $[B,B,A]$ terms are both eliminated by choosing $j=1/12$ as in the original S4C integrator given in Eqn. 11. The $[I,I,A]$ term can be neglected since it contains 2 factors of the small quantity $I$. However, it is impossible to eliminate both the remaining $O(\tau^{3})$ terms with a single value of $k$. This means that, in principle, the integrator will behave as second order rather than 4th-order.

The underlying cause of the problem is that the $[A,B]$ term in the corrector should really be applied in between the sub-steps involving $I$ in Eqn. 31 rather than outside them. However, this would entail applying the corrector and its inverse at the same frequency as the $I$ sub-steps, which would be inefficient since the corrector is typically expensive. The results shown in Figure 2 suggest that the error incurred by moving the corrector outside the $I$ steps is not a problem for the 4th-order integrator $S_{4C}$, presumably because the uncorrected terms contain a factor of $I$ which is small. However, the problem becomes obvious for the 6th-order algorithm, since multiple terms at $O(\tau^{3})$ and $O(\tau^{5})$ are not properly eliminated when the corrector is moved.

One way to overcome this problem is to use an integrator kernel that doesn’t require an $[A,B]$ term in its corrector. For example, if the integrator kernel without the $I$ steps has the following form

[TABLE]

then it will have the following form when the $I$ steps are added

[TABLE]

so that both the $[A,A,I]$ and $[B,A,I]$ terms can be eliminated by a corrector of the form $(\tau^{2}/12)[A,I]$. Some terms at $O(\tau^{5})$ will remain uncorrected regardless of what other terms are in the corrector, but these uncorrected terms will contain at least one factor of $I$, and should be small enough not to degrade the integrator’s performance.

The integrator $S_{\rm 6B}$, described by Eqn. 15 in Section 2.4 has a corrector with the necessary properties, and we now examine the performance of this integrator. Figure 3 shows the performance of $S_{\rm 6B}$ compared to $S_{\rm 6A}$ for an integration of the Sun and 8 planets of the Solar System. The integrations last for 100,000 years. The step size for the direct planet-planet terms, represented by $H_{I}$, is held constant at 1.8 days, while several values of the step size are considered for the other parts of the Hamiltonian.

For large step sizes, the energy error of $S_{\rm 6B}$ varies as roughly as the sixth power of the step size, so the algorithm is indeed behaving as a sixth-order integrator, as desired. The computational cost is generally smaller than for $S_{\rm 6A}$ for a given energy error. As we saw earlier, $S_{\rm 6A}$ mostly behaves as a second-order algorithm. For small step sizes, the accuracy of $S_{\rm 6B}$ saturates at about 1 part in $10^{-13}$. Further reductions in the step size do not improve the accuracy. This implies that the dominant source of error at this point is the direct planet-planet terms of $H_{I}$. Reducing the error further would require reducing the step size for these terms as well as the other parts of the Hamiltonian.

4.2 Simple, Exact Correctors

Another strategy for improving the 6th-order integrator $S_{\rm 6A}$ is to return the corrector to its correct position and find a way to implement the corrector more efficiently. This means that the corrector can be applied at every integration step as it should be. It turns out that this can be done, at least for the Kepler problem where the Hamiltonian equivalent to the operator $\tau^{2}[A,B]$ is integrable and can be solved efficiently.

For the Kepler problem, we have

[TABLE]

and the Hamiltonian equivalent to $\tau^{2}[A,B]$ is

[TABLE]

Advancing the system under $H_{AB}$ using Hamilton’s equation, we get

[TABLE]

Noting that the evolution of the coordinates does not depend on the momenta, we can solve these equations analytically to get

[TABLE]

where the subscript 0 indicates the initial values, and we have used the fact that $H_{AB}$ is a constant to solve the momentum equations.

A corrector with the form $\tau^{4}[A,B,B,A]$ can be advanced analytically in a similar way, but it is not obvious how to do the same for $\tau^{4}[A,A,A,B]$. Therefore, we will use integrator $S_{\rm 6C}$, described by Eqn. 18 in Section 2.4, which has a corrector that only contains terms that we can advance analytically.

To use this integrator, we once again adopt democratic heliocentric coordinates, but separate out the indirect momentum terms, so that the Hamiltonian has the following parts:

[TABLE]

where $N$ is the number of satellites orbiting the central body which has index 0. Note that $H_{A}$ and $H_{B}$ correspond to $N$ separate Kepler problems, so we can use the analytic corrector derived above.

One step of the integrator looks like this:

• •

Apply a corrector $C_{I}$

• •

Advance $H_{I}$ for $\tau/2$

• •

Advance $H_{S}$ for $\tau/2$

• •

Apply an analytic corrector proportional to $\tau^{2}[A,B]$

• •

Apply an analytic corrector proportional to $\tau^{4}[A,B,B,A]$

• •

Do the following $M$ times

• –

Apply the 7-step kernel of integrator $S_{\rm 6C}$ applied to $H_{A}$ and $H_{B}$ only.

• •

Apply an analytic inverse corrector proportional to $\tau^{4}[A,B,B,A]$

• •

Apply an analytic inverse corrector proportional to $\tau^{2}[A,B]$

• •

Advance $H_{S}$ for $\tau/2$

• •

Advance $H_{I}$ for $\tau/2$

• •

Apply an inverse corrector $-C_{I}$

where we have included the corrector $C_{I}$ to eliminate the leading error term involving $I$. This needs only be applied at the start of the integration and when output is required. We note that the $[A,B]$ and $[A,B,B,A]$ correctors can be applied one after the other without appreciably affecting the accuracy.

Figure 4 shows the performance of this algorithm compared to $S_{6A}$ for the same problem shown in Figures 2 and 3. In this case, the efficiency of the new 6th-order algorithm $S_{\rm 6C}$ is only a modest improvement over that of $S_{\rm 6A}$, presumably because applying the analytic corrector is still somewhat expensive. The algorithm does appear to behave as a 6th-order integrator for large step sizes. However, the error quickly saturates to a constant level at smaller step sizes due to the error incurred by the terms in $H_{I}$ and $H_{S}$, which are integrated with the same step size for all the cases in Figure 4. For this problem, at least, it appears that the method used by integrator $S_{\rm 6B}$ is more efficient than applying an analytic corrector as in $S_{\rm 6C}$.

5 Round-Off Error

For small step sizes, the main source of error in $N$-body integrations is “round-off” error caused by the limited precision of the computer rather than the accuracy of the integration algorithm. For T+V integrators, we can greatly reduce round-off error, at very little computational cost, using a simple procedure.

Consider a step in which the coordinates ${\bf X}$ are modified by an amount ${\bf dX}$. Typically, ${\bf dX}$ will be much smaller in magnitude than ${\bf X}$, especially for small step sizes. The information stored in the least significant digits of ${\bf dX}$ is lost when it is added to ${\bf X}$, and this error accumulates stochastically (or worse, depending on the computer) over the course of many steps. The same argument applies to changes in the momenta ${\bf P}$.

We can save much of this lost information by keeping track of the accumulated changes ${\bf dX}$ and ${\bf dP}$, updating these at each sub-step, and noting exactly how much of this information is transferred to ${\bf X}$ (or ${\bf P}$) when the quantities are updated. This is actually trivial to achieve in practice using the following procedure:

• •

At the start of an integration, set ${\bf dX}={\bf dP}=0$.

• •

Every time a new ${\bf dX}$ or ${\bf dP}$ is calculated, add it to the existing value of ${\bf dX}$ or ${\bf dP}$.

• •

When updating ${\bf X}$ (or ${\bf P}$), follow these steps:

• –

Store the coordinate values before the update ${\bf X}_{0}={\bf X}$.

• –

Update the coordinates: ${\bf X}={\bf X}_{0}+{\bf dX}$.

• –

Modify the changes: ${\bf dX}\rightarrow{\bf dX}+({\bf X}_{0}-{\bf X})$.

Note that the inclusion of the parentheses in the last step is essential for the procedure to work.

Following these steps often reduces round-off error by 2–3 decimal orders of magnitude, and requires minimal extra computational cost. Figure 5 shows an example using a million-year integration of the Sun and the 8 planets of the solar system. The figure shows the absolute energy error versus time for the sixth-order T+V integrator $S_{\rm 6B}$, using the same step size for $I$ as the other parts of the Hamiltonian. The result is compared to MVS leapfrog with a symplectic corrector. The step size for both algorithms is 0.23 days.

The energy error increases over time for the MVS algorithm due to a combination of round-off error and errors incurred by the routine that advances the Kepler problem. After 1 million years, the error is roughly 1 part in $10^{11}$, and it is likely to increase further for longer integrations. By contrast, the energy error for the T+V algorithm remains very small, less than 1 part in $10^{14}$ throughout the integration. The error at the end of the simulation is only slightly larger than that after only 100 years. This suggests that T+V integrators may be preferred to MVS algorithms for long-term integrations for which a high degree of accuracy is required.

6 Summary

Symplectic integrators separate the Hamiltonian for an $N$-body system into two or more parts that can be solved easily in isolation. The evolution of the complete system is approximated by combining multiple sub-steps that alternately advance one of the sub-systems. For $N$-body systems with a dominant central body, a common strategy is to separate the problem into (i) Keplerian orbits about the central body, and (ii) interactions between the satellites (Wisdom & Holman, 1991). These algorithms are called mixed-variable symplectic (MVS) integrators.

In this paper, we re-examine another class of symplectic integrators in which the Hamiltonian for the system is split into terms involving the kinetic energy $T$ and the potential energy $V$ respectively (Gladman et al., 1991).

The main conclusions of this study are

Classical second and fourth-order T+V integrators require substantially more computer time than MVS algorithms for the same level of accuracy. 2. 2.

The speed of T+V integrators can be improved by using different step sizes for strong and weak forces associated with primary-satellite and satellite-satellite terms respectively. 3. 3.

More efficient fourth and sixth-order T+V algorithms can be developed using force gradients and symplectic correctors. We describe 3 new sixth-order algorithms that require either 2 or 3 force evaluations per step, plus force gradients and their derivatives. 4. 4.

The fourth and sixth-order integrators are often more efficient than a second-order MVS (leapfrog) integrator, and are competitive with a second-order MVS algorithm that includes a symplectic corrector. 5. 5.

High-order T+V algorithms like these may be especially favorable compared to second-order MVS for systems containing many planets (which reduces the fractional cost of the Keplerian motion and indirect terms). 6. 6.

Round-off errors for T+V integrators can be reduced greatly, at little extra computational cost, using a simple modification that keeps track of the accumulated changes in the coordinates and momenta.

Acknowledgements

I would like to thank an anonymous referee for helpful comments on this paper.