ISSN 0012-2661, Differential Equations, 2017, Vol. 53, No. 7, pp. 923–934.
Pleiades Publishing, Ltd., 2017.
Original Russian Text
G.G. Elenin, T.G. Elenina, 2017, published in Differentsial’nye Uravneniya, 2017, Vol. 53, No. 7, pp. 950–961.
Adaptive Symplectic Conservative Numerical Methods
for the Kepler Problem
G. G. Elenin
and T. G. Elenina
Lomonosov Moscow State University, Moscow, 119991 Russia
Received February 2, 2017
Abstract— We suggest and substantiate a uniﬁed form of a family of adaptive conservative
numerical methods for the Kepler problem. The family contains methods of the second, fourth,
and sixth approximation order as well as an exact method. The methods preserve all the global
properties of the exact solution of the problem. The variable time step is chosen automatically
depending on the properties of the solution.
Mathematical models of motion of a set of interacting point masses are widely used in inter-
disciplinary studies. It suﬃces to mention ﬁelds of basic research such as celestial mechanics and
classical molecular dynamics [1–5]; mathematical models describe the motion of celestial bodies
and spacecraft in the former and the atomic and molecular motion in the latter. Modern theoreti-
cal studies of the evolution of such systems are based on the analysis of a numerical solution of the
Cauchy problem for the Hamiltonian equations.
Numerical methods commonly used for solving the Cauchy problem for systems of ordinary
diﬀerential equations distort the most important global properties of the exact solution for Hamil-
tonian systems; this manifests itself most clearly on large time intervals. It often turns out that,
starting from some instant of time, the approximate solution essentially describes the drawbacks of
the numerical method rather than the phenomenon to be studied. That is why it is important to
develop numerical methods that have a suﬃciently high approximation order and, in the framework
of exact arithmetics, preserve as many global properties of exact solutions of Hamiltonian systems
as possible. The global properties include the symplecticity of the map taking the initial state to
the current state, reversibility in time, conservation of the phase volume and the standard integrals
of motion (the total momentum, the total angular momentum, and the total energy), and also
additional conservation laws due to the symmetry of the potential energy of interaction of point
This problem of numerical mathematics has been paid considerable attention in the last three
decades [6–11]. Nevertheless, no methods preserving all global properties of solutions of the Cauchy
problem for Hamiltonian systems have been developed yet.
The Kepler problem is the simplest nontrivial example of problems of this class. Its solution de-
scribes the evolution of the state of two gravitationally interacting point masses in the plane of their
motion. The total momentum, the total angular momentum, the total energy, and the components
of the Laplace–Runge–Lenz vector are conserved on the solutions of the problem [6, p. 26]. Note
that the Verlet method , which is widely used in problems of molecular dynamics, conserves nei-
ther the total energy nor the components of the Laplace–Runge–Lenz vector in the Kepler problem.
The discrete gradient method  is not symplectic and does not conserve the components of the
Laplace–Runge–Lenz vector. As a result, the orbits of motion given by these methods qualitatively
diﬀer from the exact orbit.
In the recent years, there has been increasing interest in the development of numerical methods
for the Kepler problem [14–17]. In particular, new methods have been developed based on the idea