Numerical integration of relativistic equations of motion for Earth satellites (Q1011852)

This paper is concerned with the numerical investigation of the motion of an artificial satellite around the Earth taking into account the relativistic effects due to the Earth oblateness and the gravitational action of the Sun according to the model proposed by \textit{V. A. Brumberg} [Celest. Mech. Dyn. Astron. 88, No.~2, 209--225 (2004; Zbl 1131.70311)]. Denoting by \( \mathbf{q} = \mathbf{q} (t) \) the position of the satellite in a geocentric frame, the differential equations have the form \( \mathbf{q}'' (t) = \nabla_{\mathbf{q}} V(t,\mathbf{q}) + \mathbf{f},\) where the first term in the right hand side includes the forces derivable from a potential force, and \( \mathbf{f} \) the remaining relativistic forces. Here the author writes the above equation in the Lagrangian form \( ( L_{\mathbf{q}'} )' - L_{\mathbf{q}} + \mathbf{f} = 0 \) with \( L = (1/2) \| \mathbf{q}' \|^2 + V \), and an auxiliary energy equation \( E' = - \partial_t L + \mathbf{f}\mathbf{q}'\). For the long-time numerical integration of these equations, the author proposes a partitioned Runge-Kutta symplectic formula derived by \textit{L. Jay} [SIAM J. Numer. Anal. 33, No.~1, 368--387 (1996; Zbl 0903.65066)], which combines the four-stage Lobatto IIIA and IIIB formulas for positions and velocities. The main advantage of this method, as shown in several numerical experiments with energy equations, is that due to its geometric properties it provides an accurate numerical integration in contrast with standard non-symplectic integrators of high order. Finally, since the Lobatto IIIA and IIIB methods are fully implicit, a starting algorithm proposed by \textit{I. Higueras} and \textit{T. Roldán} [Appl. Numer. Math. 56, No.~1, 1--18 (2006; Zbl 1083.65067)] is used to improve the computational efficiency of the new method.