Fast Switch and Spline Function Inversion Algorithm with Multistep Optimization and k-Vector Search for Solving Kepler’s Equation in Celestial Mechanics

Abstract

Obtaining the inverse of a nonlinear monotonic function f(x) over a given interval is a common problem in pure and applied mathematics, the most famous example being Kepler’s description of orbital motion in the two-body approximation. In traditional numerical approaches, this problem is reduced to solving the nonlinear equation f(x)−y=0 in each point y of the co-domain. However, modern applications of orbital mechanics for Kepler’s equation, especially in many-body problems, require highly optimized numerical performance. Ongoing efforts continually attempt to improve such performance. Recently, we introduced a novel method for computing the inverse of a one-dimensional function, called the fast switch and spline inversion (FSSI) algorithm. It works by obtaining an accurate interpolation of the inverse function f−1(y) over an entire interval with a very small generation time. Here, we describe two significant improvements with respect to the performance of the original algorithm. First, the indices of the intervals for building the spline are obtained by k-vector search combined with bisection, thereby making the generation time even smaller. Second, in the case of Kepler’s equation, a multistep method for the optimized calculation of the breakpoints of the spline polynomial was designed and implemented in Cython. We demonstrate results that accurately solve Kepler’s equation for any value of the eccentricity e∈[0,1−ϵ], with ϵ=2.22×10−16, which is the limiting error in double precision. Even with modest current hardware, the CPU generation time for obtaining the solution with high accuracy in a large number of points of the co-domain can be kept to around a few nanoseconds per point.

Bivariate Infinite Series Solution of Kepler’s Equations

“…The resulting polynomials will form a 2-D spline [23], generalizing the 1-D spline that has been proposed in Refs. [24,25] for solving KE for every M when e is fixed. This bivariate spline may be used for accelerating computations involving the repetitive solution of Kepler's equation for several different values of e and M [23], as for exoplanet search [27,28] or for the implementation of Enke's method [1].…”

“…with coefficients c k,q depending on the choice of the base values (e c , M c ). These solutions converge locally around (e c , M c ), and can be used to devise 2-D spline algorithms for the numerical computation of the eccentric anomaly E for every (e, M) [23], generalizing the 1-D spline methods that have been described recently [24,25]. Since they do not require the evaluation of transcendental functions in the generation procedure, splines based on polynomial expansions, such as the 1-D cubic spline of Refs.…”

“…Since they do not require the evaluation of transcendental functions in the generation procedure, splines based on polynomial expansions, such as the 1-D cubic spline of Refs. [24,25], or the 2-D quintic spline of Ref. [23], which is based on the solutions presented here, are more convenient for numerical computations than expansions in terms of trigonometric functions.…”

Comment on "An efficient code to solve the Kepler equation. Elliptic case"

“…Other methods, such as inverse series [20,3,21], or splines [23,22], will also be affected by similar limits since they use divisions by f ′ to compute the coefficients of their expansions. Moreover, CORDIC-like methods for solving the KE have been proposed and shown to be affected by a limiting accuracy ∼ 10 −15…”

“…[17]. We also argue that this limiting accuracy is expected to affect also other methods that use derivatives, such as inverse series [20,3,21], or splines [23,22], or divisions by differences of values of f for points that are close, as occurs in the Secant method or Inverse Quadratic Interpolation appearing in Brent's scheme [2]. Moreover, CORDIC-like methods for solving the KE have also a similar limiting accuracy [25,26].…”