A MULTIRATE STÖRMER ALGORITHM FOR CLOSE ENCOUNTERS^*

Much of our understanding of the evolution of the solar system comes from long N-body simulations of planets and test particles. If the stepsize for these simulations is chosen so the local truncation error is below machine precision, the error in the numerical solution is that due to round-off. When the integration method used for such simulations is implemented in a standard way, the round-off error will typically be systematic, see for example the test results for implicit Gauss Runge–Kutta methods in Hairer et al. (2008). The round-off error will grow as t for conserved quantities where t is time, and as t² for other dynamical variables such as the position.

Brouwer (1937) showed that if the round-off error is stochastic, the power laws become t^1/2 and t^3/2, respectively. At least four integration schemes that achieve these power laws have been developed.

The first was the order 13 Störmer method of Grazier et al. (2005b). The test results in Grazier et al. (2005a) for simulations of the Sun and the four Jovian planets, Jupiter, Saturn, Uranus, and Neptune, showed the energy and phase error grew as approximately t^1/2 and t^3/2, respectively, when the stepsize was 1/1024 (4.1 days) of Jupiter's orbital period. Laskar et al. (2004) presented a symplectic method of the order of O(h⁸) + O(h⁴²) where h is the time step and is a typical planetary mass in solar masses. They performed a simulation of the Sun, the eight planets, Pluto, and the Moon over 500 million years using a stepsize of 1.83 days and found that the error in the energy was well approximated by a Gaussian distribution. Hairer et al. (2008) investigated the propagation of round-off error for implicit Gauss Runge–Kutta methods when fixed-point iteration is used to solve for the stages. They showed the order 12 method could be implemented so that the error in the energy grow approximately as t^1/2 when performing simulations of the Sun, the Jovian planets, and Pluto with a stepsize of 167 days (each step requires considerably more function evaluations than the order 13 Störmer method). Hayes (2008) showed that the Taylor method of Jorba & Zou (2005) obeys Brouwer's law when the method is used with a sufficiently small local error tolerance.

Grazier et al. (1999) used the Störmer method of Grazier et al. (2005b) to perform simulations of the Sun, the Jovian planets and 100,000 test particles over 10⁹ years. As was common at the time for simulations of this type, a test particle was removed if it came close to a planet. This removal mechanism ensured the simulations could be completed in a reasonable time but meant no insight was gained about the orbital evolution of test particles after a close encounter.

We present, analyze, and test a multirate algorithm based on the Störmer method of Grazier et al. (2005b) that accurately integrates close encounters between test particles and planets, and that retains the good error propagation of the Störmer method in Grazier et al. (2005b) when a test particle is not in a close encounter. The algorithm uses stepsizes H and h_i, i = 1, ..., N_p, where N_p is the number of planets, h_i ≪ H, and H/h_i is an integer. We refer to H and h_i as the full and reduced stepsizes, respectively. The full stepsize is used to integrate the orbits of the Sun and planets at all times, and the orbital motion of the test particles when they are not undergoing a close encounter. The reduced stepsize h_i is used to integrate the orbital motion of test particles during a close encounter with the ith planet. The position of the Sun and all planets during a close encounter is calculated using Hermite interpolation.

The use of multirate algorithms for N-body simulations is not new. Aarseth (1971) presented an order four scheme that employed different timesteps for each body. This scheme has been used extensively, see for example Weibel et al. (1990), and Schindler & Müller (1993). More recently, Saha & Tremaine (1994), and Emel'yanenko (2007) presented multirate symplectic methods. The method of Saha & Tremaine exhibits short-term errors that are O(Ω²h²) where is a typical planetary mass in solar masses and Ω is a typical orbital frequency. When a special starting procedure is used, the long-term errors are O(²Ω³h²). The method of Emel'yanenko is based on the leapfrog scheme and is second order. Our multirate algorithm differs significantly from the above algorithms because it is intended for accurate simulations.

When designing our algorithm we took advantage of two features of simulations of planets and test particles. The orbital motion of each test particle is independent of all other test particles. This enables us to integrate the motion of each test particle independently of the other test particles. The second feature is that only a small percentage of test particles undergo close encounters at a given time. We performed simulations to specifically measure this percentage and found it was never greater than one percent for the type of simulations we are interested in. This small percentage means that having h_i ≪ H will typically increase the CPU time by just a small amount compared with that when the stepsize H is used for all bodies. This leads to a large reduction in CPU time for the simulation compared to the existing method, which uses the same stepsize for all bodies because the stepsize must be noticeably smaller than H to ensure the close encounters are integrated accurately.

We define the equations of motion in Section 2 and give a brief overview of the implementation of the Störmer methods of Grazier et al. (2005b). We follow this in Section 3 with a description of how we integrate close encounters. An important part of our algorithm is an accurate interpolation scheme for the position of the Sun and planets. If the interpolation is inaccurate, the integration of close encounters will be inaccurate. We analyze our interpolation scheme in Section 4 and discuss other aspects of the close encounter algorithm in Section 5. We present a summary of our numerical testing in Section 6 and end with a discussion in Section 7.

Throughout the remainder of this paper we refer to the Sun and planets as the massive bodies, and use the name CloseEncounter to denote our algorithm for handling close encounters.

We assume the simulations are performed in three-dimensional Cartesian coordinates with the origin at the center of mass of the massive bodies. The description below is readily changed if the origin is centered on the Sun. Let R₀ be the position of the Sun at time t and R_i, i = 1, ..., N_p, the position of the ith planet at time t. A subscript of 0 and not 1 is used for the Sun to simplify the description of CloseEncounter. The equations of motion for the massive bodies are

$$\begin{equation} \ddot{{\mathbf {R}}}_i = \sum _{l=0,l \ne i}^{N_p} \frac{\mu _{{ l}} ({\mathbf {R}}_{{ l}} - {\mathbf {R}}_i)}{\Vert {\mathbf {R}}_{{ l}}-{\mathbf {R}}_i\Vert _2^3}, \quad i=0,\ldots,N_p, \end{equation} \tag{ 1 }$$

where μ_l is G times the mass of the lth body, G being the gravitational constant.

Let N_t be the number of test particles at time t, and r_j, j = 1, ..., N_t, the position of the jth test particle at time t. The equations of motion for the jth test particle are

$$\begin{equation} \ddot{{\mathbf {r}}}_j = \sum _{l=0}^{N_p} \frac{\mu _{{ l}} ({\mathbf {R}}_{{ l}} - {\mathbf {r}}_{{ j}})}{\Vert {\mathbf {R}}_{{ l}} -{\mathbf {r}}_{{ j}}\Vert _2^3},\quad j = 1,\ldots ,N_t. \end{equation} \tag{ 2 }$$

We observe that the right side of Equation (2) does not depend on r_l, l ≠ j.

Both Equations (1) and (2) have the form

$$\begin{equation} \ddot{{\mathbf {x}}}= \mathbf {f}(t,{\mathbf {x}}). \end{equation} \tag{ 3 }$$

The update formula for the kth order Störmer method applied to Equation (3) is

$$\begin{equation} {\mathbf {x}}_{n+1} - 2 {\mathbf {x}}_n + {\mathbf {x}}_{n-1} = h^2 \sum _{l=0}^{k-1} \alpha _{{ l}} \mathbf {f}_{n-l}, \quad n \ge k-1, \end{equation} \tag{ 4 }$$

where x_l is the numerical approximation to x(t) at t = lh ≡ t_l, l = 0, 1, ..., and f_l = f(t_l, x_l). The starting values x_l, l = 0, ..., k − 1, are assumed sufficiently accurate. The calculation of these values is discussed in a later section.

The coefficients α_l in Equation (4) for high-order Störmer methods are large in magnitude and alternate in sign. These two properties induce substantial growth in the round-off error on long N-body simulations of the solar system; see for example Quinlan (1994). This growth is greatly reduced if method (4) is written in the backward difference form

$$\begin{equation} {\mathbf {x}}_{n+1} - 2 {\mathbf {x}}_n + {\mathbf {x}}_{n-1} = h^2 \sum _{l=0}^{k-1} \beta _l \nabla ^l \mathbf {f}_n, \end{equation} \tag{ 5 }$$

where ∇^lf_n is the lth backward difference for the Störmer method on the step from t_n to t_{n + 1}. The backward differences are defined by ∇⁰f_n = f_n and ∇^{l + 1}f_n = ∇^lf_n − ∇^lf_{n − 1}, l = 0, 1, ..., and are related to terms in a Taylor Series expansion. The coefficients β_l are all positive and decrease monotonically with l. The analysis and test results in Grazier et al. (2005b) show that these two properties of the β_l lead to large reductions in the round-off error.

The round-off error can be further reduced by writing Equation (5) in the summed form, see for example p. 472 of Hairer et al. (1993),

$$\begin{equation} {\mathbf {x}}_{n+1} = {\mathbf {x}}_n + h {\mathbf {v}}_n, \quad {\mathbf {v}}_n = {\mathbf {v}}_{n-1} + h \sum _{l=0}^{k-1} \beta _l \nabla ^l \mathbf {f}_n, \end{equation} \tag{ 6 }$$

and by forming the sum in Equation (6) from the highest to lowest term. This final implementation with h chosen so the truncation error is below machine precision is used by Grazier et al. (2005b) for the order 13 Störmer method and achieves Brouwer's Law.

Our integration of the orbital motion of the bodies from t = t_n to t = t_{n + 1} consists of two parts. The massive bodies are integrated using the full stepsize H. The test particles are then integrated one particle at a time from t_n to t_{n + 1}. If a particle is not undergoing a close encounter at t = t_n, it is integrated using the full stepsize H and the same algorithm as for the massive bodies. Otherwise, CloseEncounter is used.

Table 1 gives pseudo-code for CloseEncounter applied to the ith planet and jth particle. When the integration of the jth particle from t_n to t_{n + 1} is required, the 3-vectors d_{k − 2} and d_{k − 1} are first calculated. The components of d_{k − 2} are the three components for the jth particle in the backward difference ∇^{k − 2}f_n; the components of d_{k − 1} are the corresponding components in the backward difference ∇^{k − 1}f_n. Next, the scalar D

$$\begin{equation*} D \equiv \max \lbrace \Vert \mathbf {d}_{k-2}\Vert _2, \Vert \mathbf {d}_{k-1}\Vert _2\rbrace , \end{equation*}$$

is formed, where || · ||₂ denotes the unweighted L₂ norm.

The quantity D is a measure of how rapidly the orbit of the test particle is changing with t. The closer a particle is to a planet, the faster the change and the larger D will be. We test if D is greater than a threshold value d_thres. If it is, the particle is deemed to be undergoing a close encounter with the planet and CloseEncounter is invoked. If not, the full stepsize H is used for the integration step from t_n to t_{n + 1}.

We use max {||d_{k − 2}||₂, ||d_{k − 1}||₂} and not just ||d_{k − 1}||₂ in the definition of D to reduce the likelihood D is zero fortuitously. If D is zero, the full stepsize H will be used on the next step which will often be the right choice.

If the particle was integrated from t_{n − 1} to t_n using a stepsize of H, CloseEncounter first uses the one-step integrator OneStep, discussed in Section 5, to calculate accurate values for the position of the jth particle at t = t_n − lh_i, l = 1, ..., k − 1. These positions along with that at t = t_n form the starting values for the integration of the close encounter using the Störmer method. OneStep integrates just the equations of motion for the particle; the position of the massive bodies for these calculations are found by interpolation, as described later. The calculation of the particle's position at t = t_n − lh_i, l = 1, ..., k − 1, is not required if a stepsize of h_i was used when integrating the particle from t_{n − 1} to t_n i.e., if the particle is already undergoing a close encounter.

Once the starting values are known, the main loop of CloseEncounter is executed. On each trip through the loop, the integrator Störmer first advances the particle h_i in time by integrating the equations of motion for the particle. Störmer is the same integrator as that used for the massive bodies. As with OneStep, interpolation is used get the position of the massive bodies. The distance of the particle from the center of the planet at the end of the step is then calculated. If this distance is no larger than the planet's radius R_i, the particle is removed from the simulation. If not, $({\mathbf {r}}_j-{\mathbf {R}}_i) \cdot (\dot{{\mathbf {r}}}_j-\dot{{\mathbf {R}}}_i)$ is calculated and its sign compared with the sign on the previous step. If the sign has changed, the particle has changed from inbound to outbound relative to the planet. Kepler's two-body problem is used to estimate the minimum distance of the particle from the planet on the step just taken. If this distance is no larger than R_i, the particle is removed from the simulation. The Keplerian approximation is very accurate because h_i is small and there are no other massive bodies nearby.

The integration using the reduced stepsize h_i is continued until either the particle is removed or t_{n + 1} is reached. We do not check if the particle leaves the close encounter partway across the interval [t_n, t_{n + 1}]. This decision means up to M_i − 1 steps with D ⩽ d_thres could be taken. We decided the increase in CPU time was not large enough to warrant complicating CloseEncounter by permitting a switch from a stepsize of h_i to H partway across [t_n, t_{n + 1}].

If upon reaching t_{n + 1}, the particle has left the close encounter, the integration of the particle will be continued with Störmer and a stepsize H. This continuation and the calculation of D necessitates saving the position of the particle at t = t_{n + 1} − lH, l = 1, ..., k − 1.

As noted in the Introduction, an important part of the algorithm described above is the interpolation scheme used to approximate the position of the massive bodies on the interval [t_n, t_{n + 1}]. Our requirement was that the interpolation scheme use as little CPU time as possible subject to the interpolation error being insignificant compared with the accumulated error in the solution. This requirement on the interpolation error is met for all t, including small t, by having the error below machine precision.

We first implemented and tested a high order difference scheme, an obvious choice given the Störmer method employs differences. Our testing in quadruple precision showed that the small orbital eccentricity of the planets meant the truncation error for the difference scheme was typically six orders of magnitude smaller than machine precision for double precision, a result we could have anticipated from the analysis in Grazier et al. (2005b). This meant that the difference scheme provided considerably more accuracy than required. This is especially so for large t when the accumulated error in the solution is noticeably larger than unit round-off. We therefore decided against using the difference scheme and investigated Hermite interpolation because it is of lower order and has a smaller overhead.

Let w denote one component of the position of a massive body. The cubic Hermite polynomial P_{3, n} for the time step from t = t_n to t = t_{n + 1} interpolates the numerical approximations w_n, $\dot{\omega }_n$, w_{n + 1}, $\dot{\omega }_{n+1}$ to w and w' at t_n and t_{n + 1}, respectively. P_{3, n} can be written in several forms. Given the way our Störmer method is implemented, the form that uses the least CPU time is

$$\begin{equation} P_{3,n} = d_0 w_n + d_1 H \dot{\omega }_n + d_2 w_{n+1} + d_3 H \dot{\omega }_{n+1}, \end{equation} \tag{ 7 }$$

where d₀ = (τ − 1)²(2τ + 1), d₁ = (τ − 1)²τ, d₂ = τ²(3 − 2τ), d₃ = τ²(τ − 1), and τ = (t − t_n)/H.

The one-step quintic Hermite polynomial $P_{5,n}^h$ for the time step from t_n to t_{n + 1} interpolates w_n, $\dot{\omega }_n$ $\ddot{\omega }_n$, w_{n + 1}, $\dot{\omega }_{n+1}$, and $\ddot{\omega }_{n+1}$. As with P_{3, n}, P_{5, n} can be written in several forms. The form analogous to that in Equation (7) is

$$\begin{eqnarray} P_{5,n} &=& d_0 w_n + d_1 H \dot{\omega }_n + d_2 H^2 \ddot{\omega }_n + d_3 w_{n+1}\nonumber\\ &&+ d_4 H \dot{\omega }_{n+1} + d_5 H^2 \ddot{\omega }_{n+1}, \end{eqnarray} \tag{ 8 }$$

where d₀ = (1 − τ)³(6τ² + 3τ + 1), d₁ = (1 − τ)³τ(3τ + 1), d₂ = (1 − τ)³τ²/2, d₃ = τ³(6τ² − 15τ + 10), d₄ = τ³(1 − τ)(3τ − 4), d₅ = τ³(τ − 1)²/2.

Since the massive bodies do not collide, the components of their position are smooth functions of t. This means that the standard theorem for the interpolation error in Hermite interpolation, see for example Theorem 2.1.5.9 on p. 57 of Stoer & Bulirsch (2002), can be applied. The interpolation error E₃(t) and E₅(t) for cubic and one-step quintic Hermite interpolation can be written as

$$\begin{eqnarray} E_3(t) &=& \frac{(t-t_n)^2(t-t_{n+1})^2w^{(4)}(\xi _3)}{4!},\nonumber\\ E_5(t) &=& \frac{(t-t_n)^3(t-t_{n+1})^3w^{(6)}(\xi _5)}{6!}, \end{eqnarray} \tag{ 9 }$$

where ξ₃, ξ₅ ∈ [t_n, t_{n + 1}]. From this we have the upper bounds on the interval [t_n, t_{n + 1}] of

$$\begin{eqnarray} &&|E_3(t)| \le \frac{H^4}{2^4 4!} \max |w^{(4)}(u)|,\nonumber\\ &&|E_5(t)| \le \frac{H^6}{2^6 6!} \max |w^{(6)}(u)|, \end{eqnarray} \tag{ 10 }$$

where the maxima are taken over the interval [t_n, t_{n + 1}] and we have used H = t_{n + 1} − t_n.

We now investigate the above bounds for the planets (the conservation laws for the center of mass imply the bounds for the Sun are closely related to those for the planets). The stepsize t_{n + 1} − t_n is small compared with the period of a planet. This means that the orbital motion of each planet on the interval [t_n, t_{n + 1}] is well approximated by Kepler's two-body problem. For example, if in a simulation of the Sun and Jovian planets, the mass of Saturn, Uranus, and Neptune is set to zero, the error in Jupiter's position over 4.1 days is approximately 10⁻⁸ astronomical units. A good estimate of |E₃(t)| and |E₅(t)| can then be obtained by applying the Hermite interpolation schemes to Kepler's problem. This problem can be stated as

$$\begin{eqnarray} &&\ddot{\omega }_1(t) = -\frac{\omega _1}{(\omega _1^2+\omega _2^2)^{3/2}}, \quad \ddot{\omega }_2(t) = -\frac{\omega _2}{(\omega _1^2+\omega _2^2)^{3/2}},\nonumber \\ &&\omega (0) = \left[1-e,0\right]^T, \quad \dot{\omega }(0) =[0,\sqrt{(1+e)(1-e)^{-1}}]^T,\nonumber\\ \end{eqnarray} \tag{ 11 }$$

where e is the eccentricity of the orbit.

The analytical solution to Equation (11) is $\omega (t) = [\cos (u) - e,\sqrt{(}1-e^2)\sin (u)]^T$ where u is the solution of Kepler's equation t = u − esin (u). The expressions for w⁽⁴⁾ and w⁽⁶⁾ in Equation (9) are readily obtained from the analytical solution for Kepler's problem using a symbolic manipulation package. For x we have

$$\begin{eqnarray} w^{(4)} &=& \frac{8\,e{c}^{2}-12\,{e}^{3}{c}^{2}+c+6\,{e}^{2}{c}^{3}-11\,{e}^{2}c+15\,{e}^{3}-7\,e}{\left(1-ec\right)^7}, \nonumber\\ w^{(6)} &=& {\frac{P^x_{5,5}}{\left(1-ec \right) ^{11}}}, \end{eqnarray} \tag{ 12 }$$

where c ≡ cos (u), and $P^x_{5,5}(c,e)$ is the quintic bi-variate polynomial given in the Appendix. For y we have

$$\begin{eqnarray} w^{(4)} &=& {\frac{(6\,{e}^{2}{c}^{2}-15\,{e}^{2}+8\,ec+1) \sqrt{1-e^2} \sin (u ) }{(1-ec) ^{7}}}, \nonumber\\ w^{(6)} &=& {\frac{P^y_{4,4}\sqrt{1-e^2}\sin \left(u\right)}{\left(1-ec \right) ^{11}}}, \end{eqnarray} \tag{ 13 }$$

where $P^y_{4,4}(c,e)$ is the quartic bi-variate polynomial given in the Appendix.

We attempted to use analytical methods to obtain the global maximum of |w⁽⁴⁾| and |w⁽⁶⁾| as a function of e. We did not obtain an elegant proof and resorted to finding the maximum numerically for e = (i/1000)e_max, i = 0, ..., 1000. The maximum for each e was found by evaluating |w⁽⁴⁾| and |w⁽⁶⁾| at 10,001 uniformly spaced values of u on [0, 2π] and taking the maximum of these values.

An important question is what value to use for the maximum eccentricity e_max. Since our algorithm is intended for simulations with the Sun and Jovian planets as the massive bodies, we used the eccentricities of the Jovian planets to help us choose e_max. N-body simulations, see for example Ito & Tanikawa (2002), show the eccentricities of the Jovian planets oscillate slowly with time. The eccentricity for Jupiter oscillates between approximately 0.025 and 0.065, and that for the remaining Jovian planets between approximately 0.00 and 0.09 (Saturn), 0.00 and 0.08 (Uranus), and 0.00 and 0.025 (Neptune). These eccentricities are bounded above by 0.1 and we used e_max = 0.1.

A stepsize of 1/1024 Jupiter's period is equivalent to a stepsize of 2π/1024 when modeling Jupiter's orbital motion using Kepler's two-body problem. If we take e = 0.065 and H = 2π/1024,

$$\begin{equation} |E_{3}(t)| \le 6.17 \times 10^{-12}\quad {and}\quad |E_{5}(t)| \le 5.45 \times 10^{-18}\quad \end{equation} \tag{ 14 }$$

for x and

$$\begin{equation} |E_3(t)| \le 4.90 \times 10^{-12}\quad {and}\quad |E_5(t)| \le 4.48 \times 10^{-18}\quad \end{equation} \tag{ 15 }$$

for y, where the numerical values in the upper bounds are given to three significant figures.

We calculated the above upper bounds for Saturn, Uranus, and Neptune using the maximum eccentricities listed above and the appropriate stepsize. The upper bounds for the three planets were all less than the upper bounds for Jupiter in Equations (14) and (15).

To assess how tight the upper bounds were, we first found very accurate solutions to Kepler's two-body problem at t_i = 2iπ/1024, i = 0, ..., 1024. We then formed the cubic and quintic Hermite polynomials on the intervals [t_i, t_{i + 1}], i = 0, ..., 1023, and used numerical techniques to estimate the maximum of the interpolation error on [0, 2π]. We found this estimate agreed with the upper bounds in Equations (14) and (15) to three significant figures.

For Kepler's two-body problem with e small, $\sqrt{\vphantom{A^A}\smash{\hbox{${\omega _1^2+\omega _2^2}$}}} \approx 1$. This means the upper bounds (14) and (15) can then be taken as relative errors for ω₁ and ω₂. This implies that the interpolation error in CloseEncounter will be below machine precision for quintic Hermite interpolation but not for cubic Hermite interpolation, suggesting quintic Hermite interpolation should be used. Our numerical tests quickly established that cubic Hermite interpolation was not sufficiently accurate for small t.

During a close encounter, the local truncation error of a Störmer method will vary approximately as a high inverse power of the distance of the test particle from the planet. This implies the largest local truncation error for a given planet will occur when the particle is just above the surface of the planet. This observation leads to the following way of estimating h_i. Choose it so that the Störmer method accurately integrates a low eccentricity Keplerian orbit with semi-minor axis equal to the radius of the planet. We used the analytical techniques in Grazier et al. (2005b) and found the error was below machine precision for eccentricities of zero and 1/10 if the stepsize was respectively 1/128 and 1/256 of the orbital period. Table 2 lists the resulting h_i and M_i for each of the Jovian planets. The entries are given to three significant figures. We return to estimating a suitable value for h_i in the next section.

The integrator OneStep needs to calculate numerical solutions to machine precision or close to this precision. OneStep need not be overly efficient because it is used just to initialize the Störmer method at the start of an integration or when the stepsize for a particle is reduced from H to h_i. Suitable integration methods for OneStep include polynomial extrapolation, see for example Hairer et al. (1993), and high-order explicit Runge–Kutta Nyström methods; see for example Dormand et al. (1987). We chose the order 12 formula of Dormand et al. (1987) for the testing described in the next section because the formula was more accurate at limiting precision than polynomial extrapolation.

We performed extensive numerical testing of CloseEncounter. We established the accuracy and efficiency of the new method, and investigated the dependence of its performance on H, M_i, d_thres and the interpolation scheme. We also compared the new method with the existing method (the method where the same small stepsize is used for all bodies). Because the orbital motion of a test particle is independent of other particles we used just one particle in our testing.

Here we present a summary of the test results for two contrasting and demanding problems of an asteroid, modeled as a test particle, having close encounters with Jupiter. The problems are denoted by AST1 and AST2. Figure 1 gives the graph of the separation between Jupiter and the asteroid for AST1 and AST2. The initial conditions for the two problems together with the μ are listed in the Appendix. Quintic Hermite interpolation was used for all test results presented here.

Zoom In Zoom Out Reset image size

The asteroid for AST1 has close encounters with Jupiter at 2316, 2998, 3999, 4850, 5610, and 6970 days, where the times are given to the nearest day and the last close encounter is mild. The closest the asteroid gets to the center of Jupiter is 76.4 R_J where R_J is the radius of Jupiter. This occurs during the first close encounter. AST1 is a worthy test problem because it shows how the error in the position of the asteroid can be rapidly magnified by multiple close encounters.

The asteroid for AST2 makes a single close encounter at 1926 days. The asteroid's distance from the center of Jupiter at this time is 1.4 R_J, making the close encounter markedly more difficult numerically than the individual close encounters of AST1.

We estimate the error in the position of the asteroid at time t as r₁(t) − r_ref(t) where r_ref(t) is the position taken from a reference solution, and denote ||r₁(t) − r_ref(t)||₂ by E(t). We calculated a reference solution for AST1 and AST2 by solving them in quadruple precision using the integrator ODEX2 of Hairer et al. (1993) with tolerances of 10⁻²ⁱ, i = 7, ..., 13, and using the results of these integrations to obtain a solution accurate to 16 significant figures. Our choice of units for distance and time means the error in the velocity of the asteroid is smaller than that for the position.

We begin with the results for AST1. The left half of Figure 2 gives the graph of E(10000) as a function H for the existing method. The graph has two distinct parts. For large H, E(10000) varies approximately as a smooth power law, indicating that the truncation error is dominating the round-off error. For small H, E(10000) oscillates rapidly with H about a mean value that increases slowly with decreasing H, indicating the round-off error is dominating the truncation error.

Zoom In Zoom Out Reset image size

The stepsize at the knee between the two parts is optimal in that E will not be smaller for smaller H except fortuitously if the stepsize corresponds to a trough of an oscillation in E. We estimated this optimal stepsize by using linear least squares to fit a power law to both parts in the graph and finding H at their intersection. This gave a stepsize of 0.64 days.

The right side of Figure 2 gives graphs of E for CloseEncounter as a function of M_i. The graphs are for d_thres equal to 2.0 × 10⁻¹⁶, 1.6 × 10⁻¹⁶, 1.3 × 10⁻¹⁶, and 1.0 × 10⁻¹⁶. The full stepsize H is 10000/1020 = 9.80 for all graphs and is optimal in that smaller values lead to little reduction in E.

We observe that as M_i increases from a small value, E decreases rapidly with M_i to approximately 10⁻¹⁰ when M_i = 15, and then undergoes oscillations with M_i. We redid the simulations in quadruple precision and the oscillations disappeared, showing they are caused by round-off error.

For each value of d_thres used in the right of Figure 2, CloseEncounter was invoked for six sub-intervals of t on [0, 10000]. The end-points of the sub-intervals varied with d_thres but the variation was insignificant for our purposes. The sub-intervals bracketed the six close encounters well. The first five sub-intervals were of similar width, and the sixth sub-interval was half this width. The width is smaller because the sixth close encounter is milder than the first five.

The optimal full stepsize H of 9.80 days for the new method is larger than the stepsize of 4.1 days used by Grazier et al. (1999) in their simulations. The main reason for this difference is that Grazier et al. (1999) did not optimize the stepsize in the same way we did. A stepsize of 9.80 days means the bounds |E₅(t)| in Section 4 for the interpolation error are a little larger than the unit round-off, indicating the optimal stepsize for the new scheme also depends on the integration error of the Störmer method when t is not small.

Figure 3 gives the graphs of E for the existing and new method applied to AST2. As for AST1, we used a full stepsize 9.80 days for the new method. The numerical severity of the close encounter means the optimal stepsize for the existing method is over two orders of magnitude smaller than that for AST1, and the optimal M_i for the new method is over two orders of magnitude larger.

Zoom In Zoom Out Reset image size

Our least squares fit of the power laws to the graph in the left of Figure 3 gave an estimated optimal stepsize of 0.00177, which is approximately 360 times smaller than for AST1. The estimated error for this stepsize is 3.6 × 10⁻⁸. We observe from the right of Figure 3 that the optimal M_i is approximately 6000, in reasonable agreement with the values for Jupiter in Table 2, and that the typical error for the new method is smaller than that for the existing method when M_i is greater than approximately 5000. The oscillations in E for the two methods means that this difference in the error is of marginal significance for our work.

The above results for AST1 and AST2 highlight one potential difficulty with the existing and new methods. The optimal value of the main parameter, H for the existing method and M_i for the new method, is problem dependent. We performed numerical tests with the existing method and found a stepsize of 0.00172 days minimized the maximum of the error for AST1 and AST2. We then performed tests with the new method with a full stepsize of 9.80 days and found M_i = 6250 minimized the maximum of the error for AST1 and AST2. The error for both methods was approximately 7 × 10⁻⁸.

The above results show that the minimum possible error for the existing and new methods applied to AST1 and AST2 are comparable. To establish the relative efficiency of the two methods we measured the CPU time. In all measurements the full stepsize for the new method was as in our previous tests (9.80 days).

Table 3 gives the ratio of the CPU time for the existing and new methods. The CPU times used to calculate the ratios are an average over five runs; the sample standard deviation over five runs was no more than three percent. The column in the table headed "Individual" is when H for the existing method and M_i for the new method is chosen to minimize the error on AST1 and AST2 separately. The new method is five times faster than the existing method on AST1 and 31 times faster on AST2. The main contributions to the six-fold increase in going from AST1 to AST2 are

• 1.  
  the stepsize used by the existing method to integrate all six bodies in AST2 was approximately 360 times smaller than that used to integrate AST1;
• 2.  
  the stepsize used by the new method to integrate five of the six bodies in AST2 was the same as that for the five bodies in AST1; and
• 3.  
  the stepsize used by the new method to integrate the sixth body in AST2 was 400 times smaller than that used to integrate the body in AST1, similar to the 360 factor for the existing method.

The cost of performing the interpolation and the number of steps by the new method with the reduced stepsize also contributed to the increase.

The column in Table 3 headed "Combined" gives the ratio of the CPU time when H for the existing method and M_i for the new method is chosen to minimize the maximum error for AST1 and AST2. The new method is 72 times faster than the existing method on AST1 and 20 times faster on AST2.

We began using a multirate Störmer algorithm for simulations of the outer solar system several years ago. We improved the scheme as we gained experience with it and we believe the algorithm as described here is near optimal for an order 13 Störmer integrator that does not vary the stepsize from step to step.

Our numerical comparisons show that the integration error of the new method and the existing method (where the same stepsize is used for all bodies) are comparable if the stepsizes are chosen to minimize the error. As part of our extensive numerical testing, we found that the new method satisfied Brouwer's Law for the intervals of time when a test particle was not undergoing a close encounter. This was expected because outside close encounters the new existing is essentially the existing method which does obey Brouwer's Law.

The speed-up of the new method relative to the existing method when the stepsizes for the two methods are chosen to minimize the error depends on the number of test particles, the severity of the close encounters, and the percentage of time spent in close encounters. The speed-up for the test results presented in the previous section varied from 5 to 72.

Our analysis and test results showed that quintic Hermite interpolation was more accurate than cubic Hermite interpolation. This negative result about cubic Hermite interpolation does not mean it should be rejected. The error bounds are sufficiently small that cubic Hermite interpolation could be used for large time when the accumulated error in the position of the massive bodies far exceeds the bounds. This observation suggests the hybrid scheme of using quintic Hermite interpolation in the first part of a simulation and cubic Hermite interpolation for the remainder of the simulation.

In the description of our algorithm, we assumed as often will be the case that all massive bodies other than the Sun were planets. Our algorithm extends without modification to simulations involving non-planetary massive bodies provided the massive bodies can be integrated throughout the simulation with the full stepsize and their position found sufficiently accurately using quintic Hermite interpolation. Our algorithm also extends to simulations where test particles have close encounters with the Sun. We have not included this because few such close encounters occur for the simulations we perform. If a test particle has a close encounter with the Sun, we save the information about its orbit and remove the test particle. At the end of the simulation, we perform a separate integration of the removed test particles if their number is significant.

The authors thank the referee for the detailed report. The referee's comments about the new method relative to the existing method were especially invaluable. The authors wish to acknowledge the contribution of the NeSI high-performance computing facilities and the staff at the Centre for eResearch at the University of Auckland. New Zealand's national facilities are provided by the New Zealand eScience Infrastructure (NeSI) and funded jointly by NeSI's collaborator institutions and through the Ministry of Business, Innovation and Employment's Infrastructure programme (http://www.nesi.org.nz).

Our analysis of the interpolation scheme in Section 4 involved the quintic bi-variate polynomial $P^x_{5,5}(c,e)$ in w⁽⁶⁾ for x and the quartic bi-variate polynomial $P^y_{4,4}(c,e)$ in w⁽⁶⁾ for y. We found these polynomials using the symbolic manipulation package MAPLE. $P^x_{5,5}(c,e)$ is

$$\begin{equation} P^x_{5,5}(c,e) = \sum _{i=0}^5 \alpha _i c^i \end{equation} \tag{ A1 }$$

where α₀ = 945e⁵ − 630e³ + 31e, α₁ = −945e⁴ + 352e² − 1, α₂ = −1260e⁵ + 1138e³ − 52e, α₃ = 954e⁴ − 328e², α₄ = 360e⁵ − 444e³, and α₅ = −120e⁴. $P^x_{4,4}(c,e)$ is

$$\begin{equation} P^x_{4,4}(c,e) = \sum _{i=0}^4 \alpha _i c^i \end{equation} \tag{ A2 }$$

where α₀ = 945e⁴ − 210e⁴ + 1, α₁ = −840e³ + 52e, α₂ = −840e⁴ + 328e², α₃ = 444e³, and α₄ = 120e⁴.

Table 4 lists the initial position and velocity of the Sun, the Jovian planets and the asteroid for problem AST1. The position and velocity for the Jovian planets and asteroid were supplied by H. F. Levison (2003, private communication). The position and velocity of the Sun were calculated using the conservation laws for the center of mass. The position is in astronomical units and the velocity is in astronomical units per day.

The initial conditions for problem AST2 were the same as those for AST1 except the position and velocity of the asteroid were [ − 3.99, 0.400, 0.240]^T and [ − 0.001812, −0.0085, 0.0001]^T, respectively.

Our μ values for the massive bodies were

$$\begin{eqnarray*} \begin{array}{lrlr} \rm Sun & {\tt 0.2959139769527998E-03} \nonumber\\ {\rm Jupiter} & {\tt 0.2825345909524226E-06} \\ {\rm Saturn} & {\tt 0.8459715185680657E-07}\nonumber\\ {\rm Uranus} & {\tt 0.1292024916781969E-07} \\ {\rm Neptune} & {\tt 0.1524358900784276E-07} \end{array} \end{eqnarray*}$$

The μ for the Sun includes that for Mercury, Venus, Earth, and Mars. The μ values are known more accurately than G and M separately.