Significant Improvement in Planetary System Simulations from Statistical Averaging

Symplectic maps (Hairer et al. 2006) are well suited for studying the astrophysical N-body problem; they have been used to study the evolution of our solar system and exoplanetary systems (e.g., Wisdom & Holman 1991). In this note, we show that an application of a symplectic integrator, regardless of its time step or order, can give biased errors and trajectories. This bias can be mitigated by considering a statistical sample of trajectories.

Consider the problem of the Sun plus the four giant planets. The bodies are treated as point particles obeying Newtonian gravity. We integrate this system for 10,000 yr using the symplectic Wisdom–Holman map (WH) in canonical Democratic Heliocentric coordinates (Wisdom & Holman 1991; Duncan et al. 1998). If h is the timestep and is the ratio of the maximum planetary mass to the Solar mass, the error scales as ${ \mathcal O }(\epsilon {h}^{2})$. ³ As a function of phase space coordinates, the N-body Hamiltonian is,

$$\begin{eqnarray}\begin{array}{rcl}H & = & \displaystyle \sum _{i}\displaystyle \frac{{{\boldsymbol{P}}}_{i}^{2}}{2{m}_{i}}-\displaystyle \frac{{{Gm}}_{0}{m}_{i}}{{Q}_{i}}\\ & & +\displaystyle \frac{1}{2{m}_{0}}{\left(\displaystyle \sum _{i\ne 0}{{\boldsymbol{P}}}_{i}\right)}^{2}-\displaystyle \sum _{0\lt i\lt j}\displaystyle \frac{{{Gm}}_{i}{m}_{j}}{{Q}_{{ij}}},\end{array}\end{eqnarray} \tag{ 1 }$$

while the error in the Hamiltonian is (Hernandez & Dehnen 2017),

$$\begin{eqnarray}&&\begin{array}{l}{H}_{\mathrm{err}}=\displaystyle \frac{{h}^{2}}{24}\left(-\displaystyle \frac{1}{{m}_{0}}{\left(\displaystyle \sum _{i\ne 0}\displaystyle \frac{{{Gm}}_{0}{m}_{i}}{{Q}_{i}^{3}}{{\boldsymbol{Q}}}_{i}\right)}^{2}\right.\\ -\displaystyle \sum _{0\lt i\lt j}\displaystyle \frac{{{Gm}}_{i}{m}_{j}}{{Q}_{{ij}}^{5}}\left[{V}_{{ij}}^{2}{Q}_{{ij}}^{2}-3{\left({{\boldsymbol{V}}}_{{ij}}\cdot {{\boldsymbol{Q}}}_{{ij}}\right)}^{2}\right]\\ +\displaystyle \sum _{\displaystyle \genfrac{}{}{0em}{}{i,j\gt 0}{i\ne j}}\displaystyle \frac{{G}^{2}{m}_{0}{m}_{i}{m}_{j}}{{Q}_{{ij}}^{3}{Q}_{i}^{3}}{{\boldsymbol{Q}}}_{{ij}}\cdot {{\boldsymbol{Q}}}_{i}\\ +\displaystyle \sum _{i\ne 0}\left(\displaystyle \frac{G}{{Q}_{i}^{5}}\left[({{\boldsymbol{P}}}_{\odot }\cdot {{\boldsymbol{P}}}_{i}){Q}_{i}^{2}-3({{\boldsymbol{P}}}_{\odot }\cdot {{\boldsymbol{Q}}}_{i})({{\boldsymbol{P}}}_{i}\cdot {{\boldsymbol{Q}}}_{i})\right]\right.\\ +\displaystyle \frac{2{{Gm}}_{i}}{{m}_{0}{Q}_{i}^{5}}\left[{P}_{\odot }^{2}{Q}_{i}^{2}-3{\left({{\boldsymbol{P}}}_{\odot }\cdot {{\boldsymbol{Q}}}_{i}\right)}^{2}\right]\\ +\left.\left.\displaystyle \frac{2}{{m}_{i}}{\left(\displaystyle \sum _{j\ne 0,i}\displaystyle \frac{{{Gm}}_{i}{m}_{j}}{{Q}_{{ij}}^{3}}{{\boldsymbol{Q}}}_{{ij}}\right)}^{2}\right)\right)+{ \mathcal O }({h}^{4}).\end{array}\end{eqnarray} \tag{ 2 }$$

Here, G is the gravitational constant, Q _i and P _i are the ith canonical coordinate vectors, m_i is the ith mass, V _i  =  P _i /m_i , and P _⊙ = ∑_j≠0 P _j . V _ij  =  V _i  −  V _j and Q _ij  =  Q _i  −  Q _j . Some terms in Equation (2) scale as and are dominant while others scale as ².

We choose h = 1 yr and generated output every 10 yr. The initial conditions are at a time randomly chosen between 0 and 100 yr, t_start. To obtain the initial conditions, we run WH with 17th order symplectic correctors (Wisdom et al. 1996) and h = 0.005 yr from time 0 to t_start (we apply a corrector to the initial conditions and an inverse corrector on output). Integrating in this way from 0 to 100 yr gives energy error δ = ΔH/H = 1.9 × 10⁻¹³, close to the limit from rounding errors.

The top left panel in Figure 1 displays δ over the 10,000 yr. The error is not symmetric about 0, indicating a trajectory with biased errors. This is confirmed with the adjacent histogram of the errors. The mean δ, $\overline{\delta }$, is 8.2 × 10⁻⁷. This curve is described exactly by the Hamiltonian,

$$\begin{eqnarray}&&\tilde{H}=H+{H}_{\mathrm{err}}.\end{eqnarray} \tag{ 3 }$$

Note $\tilde{H}$, not H, is conserved by the mapping, as can be confirmed numerically (Hernandez & Dehnen 2017).

Zoom In Zoom Out Reset image size

Next, we perform an additional 99 experiments. These are identical to the first except that t_start is chosen randomly each time. The results of the combined 100 experiments are plotted in the bottom left panel of Figure 1. The errors are more symmetric about 0, resulting in $\overline{\delta }=1.5\times {10}^{-7}$, which is 5.4 times smaller. We find evidence that $\overline{\delta }$ is 0 as the number of initial conditions tends to ∞ . For 10,000 initial conditions, $\overline{\delta }=-7.5\times {10}^{-9}$, consistent with a standard error scaling as n^−1/2, and a population $\overline{\delta }=0$, where n is the number of initial conditions.

The adjacent plot shows a histogram that is approximately normally distributed around 0 with skewness of 0.042. This skewness was not improved by using the 10,000 initial conditions or by allowing t_start to take values between 0 and 1000 yr, to capture secular frequency timescales.

Our experiment demonstrates that the trajectories of a symplectic integrator depend on the initial phase, which is undesired. We can understand this by noting that the initial phases are found by solving H which, in turn, changes the value of $\tilde{H}$. The subsequent evolution of H from the mapping depends on the initial $\tilde{H}$. At leading order in , the energy error is a weighted sum of the semimajor axes errors, so a biased energy error implies at least one biased semimajor axis error. Biases in other action-like variables may also exist.

Our averaging procedure improved the distribution of errors. Such averaging procedures could have a significant impact, for instance, in improving calculation of secular frequencies in the solar system (Laskar et al. 2011).

Substantial effort has been invested into higher accuracy symplectic integrators (Yoshida 1990; Wisdom et al. 1996; Laskar & Robutel 2001; Blanes et al. 2013; Rein et al. 2019). We investigate also their biases, using WH with a 17th order symplectic corrector. The error scales as ${ \mathcal O }(\epsilon {h}^{18})$ + ${ \mathcal O }({\epsilon }^{2}{h}^{2})$. The integration time remains the same; output is generated every 10 yr, and the timestep is h = 0.1 yr. The right panels of Figure 1 show the results for one and 100 initial conditions, respectively. Again, histograms of the errors are also shown. The error is dominated by the ${ \mathcal O }({\epsilon }^{2}{h}^{2})$ term. Compared to the low accuracy experiments, the order in the error has changed by $\epsilon {\left(0.1\right)}^{2}\,=\,{10}^{-5}$, explaining the change in scale. This is a factor smaller than if we had just reduced the time step. Our goal is to reduce the error, which we have done by using a smaller step and symplectic correctors, but other ways of reducing error exist.

For one initial condition, the error is again not symmetric about 0, and $\overline{\delta }=1.2\times {10}^{-10}$. When the 100 initial conditions are studied, $\left|\overline{\delta }\right|$ decreases by a factor 5.0; $\overline{\delta }$ is now − 2.5 × 10⁻¹¹. The errors are approximately normally distributed around 0 with skewness of 0.022. An ensemble of initial conditions again improved bias for higher accuracy methods.

Using standard statistical tools to distinguish action-like variables between the low and high accuracy integrations is not helpful. We computed the probability distribution functions (PDFs) in eccentricity and semimajor axis of Jupiter for the 100 low accuracy and 100 high accuracy integrations. Performing a 2-sample Kolmogorov–Smirnov test finds the distributions to be statistically consistent at the 5% significance level. Small differences in the PDFs due to the accuracy of the integrator are not statistically significant. Statistical agreement between disparate integrators in solar system dynamics has also been shown by Hernandez et al. (2020).

High accuracy symplectic integrators have a number of disadvantages compared to conventional integrators. The high accuracy symplectic integrators,

• 1.  
  can at best only match the error behavior of conventional methods obeying Brouwer's Law (Hairer et al. 2008; Bartram & Wittig 2021; Hernandez & Holman 2021),
• 2.  
  do not readily accomodate adaptive stepping, which can make computation more efficient, and
• 3.  
  unlike conventional integrators, typically assume a conservative Hamiltonian system is being solved.

Our results lead us to conclude that we can achieve a high degree of precision and accuracy in measurement of some solar system quantities through ensembles of low accuracy integrations.