High order explicit Runge-Kutta pairs for ephemerides of the Solar System and the Moon

Numerically integrated ephemerides of the Solar System and the Moon require very accurate integrations of systems of second order ordinary di erential equations. We present a new family of 8-9 explicit Runge-Kutta pairs and assess the performance of two new 8-9 pairs on the equations used to create the ephemeris DE102. Part of this work is the introduction of these equations as a test problem for integrators of initial value ordinary di erential equations.


Introduction
Ephemerides for the solar system and the Moon can be obtained by n umerically integrating a system of ordinary di erential equations which model all signi cant gravitational attractions between the bodies.To take full advantage of the accuracy of modern astronomical observations and to distinguish between competing analytical theories for the motion of the planets and the Moon, the global error in the integrations must be very small.Another characteristic of the integrations is that they often span a large interval of astronomical time, necessitating many i n tegration steps.
The requirement of a very small global error and the need to take many steps often means the integrations must be performed using extended precision such as 80-bit arithmetic or quadruple precision.
The ordinary di erential equations for ephemerides are non-sti and hence explicit Runge-Kutta (ERK) pairs are suitable methods for performing the integrations.Of the many ERK pairs available, the 13-stage 7-8 pair of Prince and Dormand 5] has proven to be as e cient a s a n y other on many problems when using double precision arithmetic, except possibly at low accuracy requirements.
In particular, the pair is noticeably more e cient than 8-9 pairs.We investigate whether this result holds for numerically integrated ephemerides.
We beginour investigation with a search of an existing family of 16-stage 8-9 pairs for a near optimal pair.Next we derive a family of 17-stage 8-9 pairs and search this for a near optimal pair.We then compare the new pairs and the 7-8 pair of Prince and Dormand on the equations of DE102.
2 Order nine pairs 2.1 De nitions Consider the initial value problem y 0 = f (x y) y(x 0 ) = y 0 (1)   where 0 d=dx, f : R R n !R n and the solution y(x) is su ciently di erentiable.
The 8-9 ERK pairs we investigate have s-stages and generate an order nine approximation y i and an order eight a p p r o ximation b y i to y(x i ), i = 1 2 : : : according to where h = x i ; x i;1 and f j = f (x i;1 + hc j y i;1 + h j;1 X k=1 a j k f k ) j = 1 : : : s (c 1 = 0 ) : We refer to c j j = 1 : : : s , as the abscissae and a ij j = 1 : : : i ; 1 i = 2 : : : s as the interior weights.To retain the one step nature of the methods, we restrict the abscissae to the interval 0 1].The e ciency of 8-9 pairs can beassessed using the size of the principal error coecients.These can bewritten as e 10 ( j ) = ( j ) 10! ( where T 10 is the set of rooted trees of order ten, (t j ) and (t j ) are positive i n tegers and k ( j ) are functions of the interior weights and abscissae.We use two measures of the size of the principal error coe cients (4) 2.2 Sixteen stages Verner 7] derived a family of 16-stage 8-9 pairs with c 2 , c 5 , c 9 , c 10 , c 11 , c 13 , c 14 and a 11 6 as free parameters (To simplify what follows, we have interchanged the coe cients for the fourteen and sixteenth stages, this can be done without changing the properties of the pairs.)The order nine formula in the pairs uses the rst fteen stages and the order eight formula uses all sixteen stages.The coe cients b j b b j , j = 2 : : : 7, b 16 , b b 14 and b b 15 are identically zero.
Verner presented the coe cients of a pair from this family which had c 2 = 1=12, c 5 = (2 + 2 p 6)=15, c 9 = 1=2, c 10 = 1=3, c 11 = 1=4, c 12 = 4=3, c 13 = 5=6, c 14 = 1=6 and a 11 6 = 0 .The pair has E 2 10 = 6 :1 10 ;5 and E 1 10 = 3 :1 10 ;5 and has beenused when comparing the numerical performance of 8-9 pairs with pairs of other orders.However the pair was intended as an illustration of the derivation and not as one with a near optimal value of E 2 10 or E 1 10 .To assess in a problem independent way if the 8-9 family of Verner contains more e cient pairs, and if so, how much more e cient, we performed a minimisation of E 2 10 over the free parameters, subject to the constraint that the coe cients be no larger than M in magnitude.This constraint is commonly used when selecting a pair from a family and is intended to prevent the selection of a pair with poor round-o error properties.Although no one value of M is used, it is often 20 or 30 and we chose 30.
A slightly smaller value of E 2 10 is possible if a smaller grid size is used, but since the numberof derivative evaluations varies approximately as the ninth root of E 2 10 , the gain in e ciency is small.A signi cantly smaller value of E 2 10 , approximately twice as small, is possible if the abscissae are not constrained to the interval 0 1], but this choice means the pair is no longer a o n e step method.
A comparison of E 2 10 for the new pair and the one of Verner suggests the new pair will beapproximately 70 percent more e cient at small stepsizes, raising the possibility o f i t being competitive with pairs of other orders.

Seventeen stages
The work of Sharp and Smart 6] for 4-5 and 5-6 ERK pairs shows a gain in e ciency is possible if an extra stage is used to form the pair.The extra stage means more free parameters are available, permitting a smaller value of E 2 10 , but this is at the expense of increasing (by one) the numberof function evaluations required to take a step.
To investigate if a gain in e ciency was possible for 8-9 pairs, we derived a family of 17-stage 8-9 pairs.The family has six more free parameters (three abscissae, three interior weights) than the 16-stage 8-9 pairs.
A comparison of E 2 10 for the new 16-stage and 17-stage pairs suggests, after scaling by the number of stages, the 17-stage pair will be 18 percent more e cient at small stepsizes.

Generalised
The families of 8-9 pairs described in the previous sub-section are readily generalised to include either one or two extra free parameters.
One generalisation is to replace b b j , j = 1 : : : s by the convex linear combination b j + ( 1 ; ) b b j .This substitution is equivalent to making one of the previously identically zero b b a free parameter.The local error estimate for the pair is changed, but since b j , j = 1 : : : s ; 1 remain the same, the error coe cients of the order nine formulae and hence E 2 10 (and E 1 10 ) are unchanged.The second generalisation is based on a transformation obtained by V erner 8] for two families of 8-stage 5-6 (ERK) pairs.Verner showed the family of Prince and Dormand 5] which has c 2 , c 3 , c 5 , c 6 , b 8 and b b 7 as free parameters can be obtained from the class in Verner 7] which has c 2 , c 3 , c 5 and c 6 as free parameters using a simple transformation on the last two rows of interior weights.
This transformation generalises (Verner,  The introduction of these two free parameters changes the local error estimate and the principal error coe cients of the order nine formula.However, as is the case for the 5-6 pairs in 7], the change in E 2 10 and E 1 10 is small for near optimal pairs.

DE102
Newhall, Standish and Williams 4] presented DE 102, an ephemeris of the solar system and the Moon, obtained by integrating a system of 33 second order ordinary di erential equations of the form y 00 = f (t y y 0 ): (8) The system (8) consists of equations of motion for the nine planets, the Moon together with three equations for the lunar physical librations.The motion of the Sun is found from the de nition of the solar system barycentre.The equations contain contributions from point-mass interactions, gure e ects for Earth and the Moon, Earth tides and perturbations from the ve asteroids (1) Ceres, (2) Pallas, (4) Vesta, (7) Iris and (324) Bamberga.
The calculations required for one evaluation of the second derivative for (8) are described in Figure 1.A fuller description is given in 4] and by inference in the program DE118i.ARC o f Moshier, available on the internet.

Numerical experiments
We conducted numerical tests of the two new 8-9 pairs and the 7-8 pair of Prince and Dormand on the DE102 equations described in the previous section.The results are illustrated below.The pairs are denoted by PD78 (Prince and Dormand 7-8), P16 (new 16-stage) and P17 (new 17-stage).
A computer which performed quadruple precision in hardware was unavailable and hence we used the the multiprecision Fortran90 package MPFUN90 of Bailey 1], with the precision level set at 35 digits, approximately that of quadruple precision.The multi-1.Initialise a) Calculate the heliocentric position and velocity for the asteroids and transform to approximate barycentric values.These values are corrected once the correct position of the Sun is known.b) Calculate the distance between the bodies.The distances involving the Sun or asteroids are estimates only.These distances are corrected once the correct position of the Sun is known.c) Use xed-point iteration to nd the correct position and velocity of the Sun and asteroids, then correct the distances involving the Sun or asteroids.d) Calculate the cube of the distances between all bodies.The coe cients of the 7-8 pair as speci ed in 5] are accurate to approximately 18 digits.We recalculated the coe cients in 100 digit arithmetic, using the values of the free parameters in 5], and used these coe cients, rounded to 35 digits.The global error in a numerical solution was obtained by calculating a more accurate solution and taking the di erence between the two solutions.

Point-mass acceleration
The rst example is for an integration interval of 20 and local error tolerances of 10 i , i = ;14 : : : ;22. Figure 2 contains the log-log graph of the number of derivative evaluations against the norm of the end-point global error (the points have been joined for clarity).Pair P17 is more e cient than P16 suggesting the e ciency is improved by adding a stage.The gain in e ciency varies from 15 to 20 percent, in good agreement with that predicted using E 2 10 .The pairs P16 and P17 are more e cient than PD8 for global errors smaller than (approximately) 10 ;16 , and 10 ;18:5 respectively.In addition and as can beexpected from the order of the pairs, the e ciency of 8-9 pairs relative to the 7-8 pair increases as the global error decreases.For example, P17 is 16 percent more e cient for a global error of 10 ;20 and 29 percent more e cient for a global error of 10 ;22 .The second example is for an integration interval of 50 using the same local error tolerances as in the rst example.The results are given in Figure 3. P16 was excluded because our test results such as those in Figure 2 showed P16 was less e cient than P17 for the local error tolerances we were using.
The e ciency of P17 relative to PD78 as a function of the global error is similar to that for the rst example, except for a minor di erence at the larger global errors.The global errors are larger than in the rst example, a result which is consistent with a larger interval of integration.

Discussion
The main aim of our work was to investigate if 8-9 explicit Runge-Kutta pairs were more e cient than lower order pairs, principally 7-8 pairs, for numerically integrated ephemerides.We derived a new family of 8-9 pairs, obtained near optimal 8-9 pairs from this family and an existing one, and compared the performance of these pairs and the 7-8 pair of Prince and Dormand on the equations of the ephemerides DE102.
Our testing showed the 8-9 pairs were usually more e cient than the 7-8 pairs.The gain in e ciency was not large, but given the amount of CPU time required to produce ephemerides, the gain is signi cant.Our testing also showed that near optimal 17-stage 8-9 can bemore e cient than near optimal 16-stage 8-9 pairs.As part of this work we introduced the equations of DE102 as a test problem for integrators of initial value ordinary di erential equations.This problem, in addition to being a realistic one, has several interesting numerical aspects.For example, the position and velocity of the Sun is found by solving a system of nonlinear (algebraic) equations.The question arises as to the most e cient way of solving these equations.
fje 10 ( j )jg : private communication) to other families of pairs, including the 8-9 pairs in this paper.This means b 16 and b b 15 , previously zero in the 16-stage 8-9 pairs, and b 17 and b b 16 , previously zero in the 17-stage 8-9 pairs, can befree parameters.
a) Calculate the Newtonian acceleration of all bodies.b) Calculate the post-Newtonian acceleration of the planets and the Moon.

3 .
Figure of the Moon a) Form the rotation matrix for the transformation from space to body coordinates.b) Calculate the e ect of the point-mass Earth on the lunar gure and add this to the lunar acceleration.c) Calculate the torque on the Moon due to the gravitational interaction between the lunar gure and the external point-mass Earth.d) The acceleration from b) induces an acceleration in the Earth -add this to the Earth's acceleration.e) Calculate the e ect of the point-mass Sun on the lunar gure and add this to the lunar acceleration.f )Calculate the torque on the Moon due to the gravitational interaction between the lunar gure and the external point-mass Sun.g) Calculate the acceleration of the libration angles.4. Figure of the Earth a) Calculate the e ect of the point-mass Moon on the Earth's gure and add this to the Earth's acceleration.b) The acceleration from a) induces an acceleration in the Moon -add this to the lunar acceleration.c) Calculate the contribution to the acceleration of the Moon and the Earth due to the Earth tides.d) Calculate the e ect of the point-mass Sun on the Earth's gure and add this to the Earth's acceleration.The accelerations in this section are adjusted for the precession and nutation of the equinox and obliquity of the ecliptic.

Figure 1 :
Figure 1: A summary of the calculations required for one evaluation of the second derivative in the mathematical model of DE 102.

Figure 2 :
Figure 2: A log-log graph of the number of derivative e v aluations against the norm of the end-point global error for DE102 with a integration interval of 20.Prince and Dormand 7-8 pair -dashed line, new 16-stage pair -dotted line, new 17-stage pair -solid line.
integration interval of 50 Base 10 logarithm of the norm of the end−point global error Base 10 logarithm of the number of derivative evaluations 13−stage 17−stage

Figure 3 :
Figure 3: A log-log graph of the number of derivative e v aluations against the norm of the end-point global error for DE102 with a integration interval of 50.Prince and Dormand 7-8 pair -dashed line, new 17-stage pair -solid line.