Efficient orbit integration by scaling for Kepler energy consistency

Extending Nacozy's idea of manifold correction by using the concept of the integral invariant relation, we propose a new approach to numerically integrate quasi-Keplerian orbits. The method integrates the time evolution of the Kepler energy and the usual equation of motion simultaneously. Then it directly adjusts the integrated position and velocity by a spatial scale transformation in order to satisfy the Kepler energy relation rigorously at every integration step. The scale factor is determined by solving an associated cubic equation precisely with the help of Newton's method. In treating multiple bodies, the Kepler energies are integrated for each body and the scale factors are adjusted separately. The implementation of the new method is simple, the added cost of computation is low, and its applicability is wide. Numerical experiments show that the scaling reduces the integration error drastically. In the case of pure Keplerian orbits, the truncation error grows linearly with respect to time, and the round-off error grows more slowly than that. When perturbations exist, a component that grows with the second or a higher power of time appears in the truncation error, but its magnitude is reduced significantly as compared with the case without scaling. The rate of decrease varies roughly as the 5/4 to 5/2 power of the strength of the perturbing acceleration, where the power index depends on the type of perturbation. The method seems to suppress the accumulation of roundoff errors in the perturbed cases, although the details remain to be investigated. The new approach provides a fast and high-precision device with which to simulate the orbital motions of major and minor planets, natural and artificial satellites, comets, and space vehicles at a negligible increase in computational cost.