Comparison of numerical methods for the integration of natural satellite systems

We present a new implementation of the recurrent power series (RPS) method which we have developed for the integration of the system of N satellites orbiting a point-mass planet. This implementation is proved to be mon efficient than previously developed implementations of the same method. Furthermore, its comparison with two of the most popular numerical integration methods: the 10th-order Gauss-Jackson backward difference method and the Runge-Kutta-Nystrom RKN12(10)17M shows that the RPS method is more than one order of magnitude better in accuracy than the other two. Various test problems with one up to four satellites are used, with initial conditions obtained from ephemerides of the saturnian satellite system. For each of the three methods we find the values of the user-specified parameters (such as the method's step-size (h) or tolerance (TOL)) that minimize the global error in the satellites' coordinates while keeping the computer time within reasonable limits. While the optimal values of the step-sizes for the methods GJ and RKN are all very small (less than T/100), the ones that are suitable for the RPS method are within the range: T/13 < h < T/6 (T being the period of the innermost satellite of the problem). Comparing the results obtained by the three methods for these step-sizes and for the various test problems we observe the superiority of the RPS method over GJ in terms of accuracy and over RKN both in accuracy and in speed.