Bivariate Infinite Series Solution of Kepler's Equations
aa r X i v : . [ a s t r o - ph . I M ] F e b B IVAR IATE I NFINITE S ER IES S OLUTION OF K EPLER ’ S E QUATIONS
Daniele Tommasini
Applied Physics Department, School of Aeronautic and Space EngineeringUniversidade de Vigo. As Lagoas s/n, Ourense, 32004 Spain [email protected]
February 23, 2021 A BSTRACT
A class of bivariate infinite series solutions of the elliptic and hyperbolic Kepler equations is de-scribed, adding to the handful of 1-D series that have been found throughout the centuries. Thisresult is based on an iterative procedure for the analytical computation of all the higher-order par-tial derivatives of the eccentric anomaly with respect to the eccentricity e and mean anomaly M in a given base point ( e c , M c ) of the ( e, M ) plane. Explicit examples of such bivariate infinite se-ries are provided, corresponding to different choices of ( e c , M c ) , and their convergence is studiednumerically. In particular, the polynomials that are obtained by truncating the infinite series up tothe fifth degree reach high levels of accuracy in significantly large regions of the parameter space ( e, M ) . Besides their theoretical interest, these series can be used for designing 2-D spline numeri-cal algorithms for efficiently solving Kepler’s equations for all values of the eccentricity and meananomaly. K eywords Elliptic Kepler equation · hyperbolic Kepler equation · Orbital Mechanics · Astrodynamics
In the Newtonian approximation, the time dependence of the relative position of two distant or spherically symmetricbodies that move in each other’s gravitational field can be written with explicit analytical formulas involving a finitenumber of terms only when the eccentricity, e , is equal to or , corresponding to circular and parabolic orbits,respectively [Roy(2005)]. For < e < and for e > , such evolution can be obtained by solving for E one of thefollowing two Kepler Equations (KEs) (see e.g. Chap. 4 of Ref. [Roy(2005)]), M = f ( e, E ) = (cid:26) E − e sin E, for e < e sinh E − E, for e > , (1)where M and E are measures of the epoch and the angular position called the mean and the eccentric anomaly,respectively (for convenience, the same symbols are used here for the elliptic and hyperbolic anomalies, even thoughthey are defined in different ways).A handful of infinite series solutions of Eqs. (1) have been found throughout the centuries (see Chap. 3 inRef. [Colwell(1993)]). The solution for < e < has been written as an expansion in powers of e [Lagrange(1771)],or as an expansion in the basis functions sin( nM ) with coefficients proportional to the values J n ( e ) of Bessel func-tions [Bessel(1805), Fernande(1994)]. Levi-Civita [Levi-Civita(1904a), Levi-Civita(1904b)] described a series inpowers of the combination z = e exp( √ − e )1+ √ − e . Finally, Stumpff found an infinite series expansion in powers of M [Stumpff(1968)]. PREPRINT - F
EBRUARY
23, 2021This article describes a class of solutions of KEs, Eqs. (1), in terms of bivariate infinite series in powers of both e and M , E = ∞ X k =0 ∞ X q =0 c k,q ( e − e c ) k ( M − M c ) q , (2)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 ) [Tommasini and Olivieri(2021)], generalizing the 1-D spline methods that have been describedrecently [Tommasini and Olivieri(2020a), Tommasini and Olivieri(2020b)]. Since they do not require the evaluationof transcendental functions in the generation procedure, splines based on polynomial expansions, such as the 1-Dcubic spline of Refs. [Tommasini and Olivieri(2020a), Tommasini and Olivieri(2020b)], or the 2-D quintic splineof Ref. [Tommasini and Olivieri(2021)], which is based on the solutions presented here, are more convenient fornumerical computations than expansions in terms of trigonometric functions. Let the unknown exact solution of Eq. (1) be E = g ( e, M ) . (3)If the analytical expression of the partial derivatives of g ( e, M ) were known, a bivariate Taylor expansion couldbe written for any choice of base values e c , E c , M c = f ( e c , E c ) , so that Eq. (2) would be demonstrated with thecoefficients given by c k,q = 1 k ! q ! (cid:20) ∂ k + q g∂e k ∂M q ( e c , M c ) (cid:21) . (4)In order to obtain such derivatives, we notice that the definitions in Eqs. (1) and (3) imply the identity, E = g ( e, f ( e, E )) . (5)In this expression, e and E are considered as the independent variables. Therefore, by taking the differential, weobtain, d E = ∂g∂e ( e, f ( e, E )) d e + ∂g∂M ( e, f ( e, E )) (cid:20) ∂f∂e ( e, E ) d e + ∂f∂E ( e, E ) d E (cid:21) . (6)Since e and E are independent, the coefficients of d E and d e must cancel separately. This condition can be used toobtain the partial derivatives of g . Taking also into account Eqs. (1) and (3), the cancellation of the coefficient of d E implies, ∂g∂M ( e, M ) = 1 ∂f∂E ( e, g ( e, M )) ≡ λ − eC , (7)where we have defined the parameter λ and the functions C such that λ = 1 and C = cos g ( e, M ) , for e < , or λ = − and C = cosh g ( e, M ) , for e > . As it could be expected, Eq. (7) coincides with the usual rule for thederivative of the inverse function when e is considered to be a fixed parameter [Stumpff(1968), Colwell(1993)]. Thecancellation of the coefficient of d e in Eq. (6) implies, ∂g∂e ( e, M ) = − ∂g∂M ( e, M ) ∂f∂e ( e, g ( e, M )) ≡ S − eC , (8)with S defined as S = sin g ( e, M ) , for e < , or S = sinh g ( e, M ) , for e > . In the case of the elliptic KE, the resultof Eq. (8) was used in Ref. [Fernande(1994)] to derive and expansion in the basis sin nM .Eqs. (7) and (8), taken together with Eqs. (1) and (3), can be used for the iterative computation of all the higher orderderivatives entering Eq. (2) for e = 1 . The calculations can be simplified by expressing all the derivatives in terms ofonly λ , S , and C , and using the following identities, which can be derived from the definitions of S and C and Eqs.(7) and (8), ∂S∂e ( e, M ) = C S − e C , ∂S∂M ( e, M ) = λ C − e C , (9) ∂C∂e ( e, M ) = − λ S − e C , ∂C∂M ( e, M ) = − S − e C . (10)2 PREPRINT - F
EBRUARY
23, 2021The second order derivatives can then be obtained by applying the operators ∂∂e and ∂∂M and the rules of Eqs. (9) and(10) to the first order derivatives given in Eqs. (7) and (8). The result is, ∂ g∂e ( e, M ) = 2 C S (1 − e C ) − λ e S (1 − e C ) , (11) ∂ g∂M ( e, M ) = − λ e S (1 − e C ) , (12) ∂ g∂e∂M ( e, M ) = λ C (1 − e C ) − e S (1 − e C ) . (13)Similarly, the third order derivatives can be obtained by applying the operators ∂∂e and ∂∂M and the rules of Eqs. (9)and (10) to the second order derivatives, Eqs. (11), (12), and (13). The result is, ∂ g∂e ( e, M ) = 6 C S − λ S (1 − e C ) − λ e C S (1 − e C ) + 3 e S (1 − e C ) , (14) ∂ g∂M ( e, M ) = − e C (1 − e C ) + 3 λ e S (1 − e C ) , (15) ∂ g∂e ∂M ( e, M ) = 2 λ C − S (1 − e C ) − e C S (1 − e C ) + 3 λ e S (1 − e C ) , (16) ∂ g∂e∂M ( e, M ) = − λ S (1 − e C ) − λ e C S (1 − e C ) + 3 e S (1 − e C ) . (17)All the higher order derivatives can be obtained by iterating this procedure. These expressions are exact, but theydepend on the unknown function g through S and C . Nevertheless, taken together with Eq. (3), they can be usedto compute 2-D Taylor series solutions of KEs. This can be done by choosing a pair of base values, e c and E c ,corresponding to M c = f ( e c , E c ) . The values of the coefficients entering Eqs. (4) and (2) can then be computed bysubstituting λ = 1 , S = sin E c , C = cos E c , for e c < , or λ = − , S = sinh E c , C = cosh E c , for e c > , in theexpressions for the derivatives of g , and by defining the zeroth order term ∂ g∂e M ( e c , M c ) = g ( e c , M c ) = E c . Thisprocedure can be used to build a class of bivariate infinite series solutions of the elliptic and hyperbolic KEs, one forany given choice of base values. Three explicit examples will be given in Sec. 3.The determination of the radius of convergence for the univariate series has been a formidable mathematical problem(see Chap. 6 of Ref. [Colwell(1993)]). Although the rigorous generalization of such analyses to the bivariate series ofEqs. (2) and (4) surpasses the scope of this work, a numerical criterion for estimating the region of convergence willstill be given in Sec. 3. In this section, three examples of bivariate infinite series solutions of KEs are given. They have been obtained fromEqs. (2) and (4) by applying the methods discussed in Sec. 2 for the computation of the derivatives of g , for threedifferent choices of the base values ( e c , M c ) . All the non-vanishing terms up to fifth order are shown explicitly. Since g ( e, − M ) = − g ( e, M ) , it is sufficient to solve KEs only for positive values of M . Moreover, for e < the M domaincan be reduced to the interval ≤ M ≤ π , and then the solution for every M can be obtained by using the periodicityof f and g .In all cases, it is convenient to define approximate solutions S n obtained by truncating the infinite series of Eq. (2)keeping only the terms with k + q ≤ n , so that S n ( e, M ) = n X k =0 n − k X q =0 c k,q ( e − e c ) k ( M − M c ) q , (18)with coefficients given by Eq. (4). The errors E n of the approximate solutions S n can then be evaluated in a selfconsistent way, E n ( e, M ) = | S n ( e, M ) − S n ( e, f ( e, S n ( e, M ))) | . (19)3 PREPRINT - F
EBRUARY
23, 2021From a practical point of view, the convergence of the infinite series for certain values of ( e, M ) means that E n ( e, M ) should tend to decrease for increasing n . This idea is used for obtaining an estimate of the region of convergence inthe ( e, M ) parameter space by comparing the average errors for lower and higher values of n with the following ruleof thumb, E ( e, M ) + E ( e, M ) + E ( e, M ) >
32 [ E ( e, M ) + E ( e, M )] . (20) e c = 0 , M c = 0 Choosing e c = 0 , E c = 0 , so that M c = 0 rad, λ c = 1 , S c = 0 , C c = 1 , the series of Eqs. (2) and (4) becomes, E = M + eM + e M + e M − e M + e M − e M + · · · (21)Fig. 1 shows the variation of the errors E n ( e, M ) , Eq. (19), computed along the line M = πe , when the bivariatepolynomials S n of Eqs. (18) are obtained by truncating Eq. (21) to degree n . Along this line, the error E is at thelevel of arithmetic double precision ( ǫ double = 2 . × − ) for e = Mπ . . . For almost all values of e . . the error decreases as n increases, as expected for a convergent series. Below this value, inversions in the order of the E n s only occur around e ∼ . and e ∼ . , when the second order and the first order approximations, respectively,become smaller than higher order ones. For e & . , the inversions become more evident and cannot be justified anylonger with oscillatory terms in f and g . These considerations hint at the conclusion that in this specific direction thelimit for convergence should be somewhere in between e ∼ . and e ∼ . . A more precise, though still approximate,criterion of convergence may be obtained with the rule of thumb of Eq. (20). In the case of the M = π e line, this rulegives the values e < . (or M < . rad), in agreement with the conclusion obtained from the qualitative analysisof Fig. 1.Fig. 2 shows the contour levels in the ( e, M ) plane of the error E affecting the fifth degree polynomial approximation, S , as given by Eq. (21). The limit e < . for the convergence of Lagrange’s univariate series (see page26 of Ref. [Colwell(1993)]) is represented by the vertical dashed magenta line. The continuous black curve marksthe boundary of the region of convergence of the bivariate series of Eq. (21), as estimated with Eq. (20). The error E is kept below ∼ − rad for e . . and M . π/ , and is reduced to the level ∼ − rad for e ∼ . and M ∼ π/ . Moreover, the fifth order approximation reaches machine precision ǫ double in an entire neighborhood ofsize ∆ e ∼ × − , ∆ M ∼ × − rad around the point ( e c , M c ) . e c = , M c = π − Choosing e c = , E c = π , so that M c = π − , λ c = 1 , S c = 1 , C c = 0 , and defining δ = e − e c = e − and ∆ = M − M c = M − π − , the series of Eqs. (2) and (4) becomes, E = π δ − δ + 2 δ ∆ + ∆ − δ − δ ∆ − δ ∆ + ∆ δ + 244 δ ∆ + 222 δ ∆ + 52 δ ∆ −
192 ++ 37 δ − δ ∆ − δ ∆ − δ ∆ − δ ∆ + 9∆
384 + · · · (22)The errors E n ( e, M ) shown in Fig. 3 for five polynomials S n , obtained by truncating Eq. (22) to degree n , have beencomputed along the two lines connecting the opposite points (0 , and (1 , π ) of the ( e, M ) domain with the mid point ( e c , M c ) = ( , π − ) . For most values of e . . the errors tend to decrease as n increases, though the inversionsbecome more evident as e approaches . In this case, the rule of thumb of Eq. (20) gives the values e < . (or M < . rad) along this line.Fig. 4 shows the contour levels in the ( e, M ) plane of the error E affecting the fifth degree polynomial approximation, S , as given by Eq. (22). The continuous black curve marks the boundary of the region of convergence, as estimatedwith Eq. (20). It can be seen that the Taylor series based in the mid point ( , π − ) converges for most values of ( e, M ) ,although even in this case there are two small regions in which it is expected to diverge. Moreover, the fifth degreepolynomial reaches machine precision ǫ double in an elongated neighborhood of the point ( , π − ) along a diagonal linecrossing the entire e domain from ( e = 0 , M = 1 . rad ) to ( e ≃ , M = 0 . rad ) , with transverse size (along the M direction) ranging from ∼ × − rad close to the endpoints, to ∼ − rad around the center ( e c , M c ) .4 PREPRINT - F
EBRUARY
23, 2021 −3 −2 −1 e , M/π −15 −12 −9 −6 −3 e rr o r ( r a d ) Figure 1: Errors E n ( e, M ) affecting the approximate polynomial solutions S n ( e, M ) of KE, as a function of theeccentricity e (in logarithmic scale). The S n are obtained by truncating the infinite series of Eq. (21) up to degree n ,for n = 1 , · · · , , and the plotted errors have been computed along the line M = π e of the ( e, M ) plane. e c = 2 , M c = 0 Choosing e c = 2 , E c = 0 , so that M c = 0 , λ c = − , S c = 0 , C c = 1 , defining δ = e − e c = e − , the series of Eqs.(2) and (4) becomes, E = M − M δ + M δ − M − M δ + 7 M δ M δ − M δ M
60 + · · · , (23)where now E and M indicate the (dimensionless) hyperbolic anomalies.Fig. 5 shows the variation of the errors E n ( e, M ) computed along the line M = e − , for < M ≤ , when thebivariate polynomials S n of Eqs. (18) are obtained by truncating Eq. (23) to degree n . Along this line, the error E reaches machine precision ǫ double for e − M . . . It can be seen that the inversions in the order of the errors E n become more and more important for e − & . , in agreement with the estimate e < . (corresponding to M < . ) obtained by applying the condition of convergence Eq. (20) to the direction of the M = e − line.For the hyperbolic motion, the values of e and M can vary in infinite ranges, < e < ∞ , < M < ∞ (due tothe symmetry for M → − M ). In Fig. 6, the contour levels in the ( e, M ) plane of the error E for the solution (23)have been drawn in the region e ≤ , M ≤ . This region has been chosen in such a way that the plot contains thecontinuous black curve marking the boundary of the region of convergence, as estimated with Eq. (20). Moreover,the fifth degree polynomial reaches machine precision ǫ double in an entire neighborhood of size ∆ e ∼ × − , ∆ M ∼ × − , around the point ( e c , M c ) . 5 PREPRINT - F
EBRUARY
23, 2021 −3 −2 −1 e −3 −2 −1 M ( r a d ) −19 −16 −13 −10 −7 −4 −1 Figure 2: Contour levels of the error E ( e, M ) affecting the fifth degree polynomial approximation, Eq. (21), as afunction of the eccentricity e and the mean anomaly M (both in logarithmic scales). The continuous black curvemarks the boundary of the region of convergence for the bivariate series of Eq. (21), as estimated with Eq. (20). Thevertical dashed magenta line represents the limit of the region of convergence for Lagrange’s univariate series. I have described an analytical procedure for the exact computation of all the higher-order partial derivatives of theelliptic and hyperbolic eccentric anomalies with respect to both the eccentricity e and the mean anomaly M . Althoughsuch derivatives depend implicitly on the solution of KE, they can be computed explicitly by choosing a couple ofbase values e c and E c for the eccentricity and the eccentric anomaly, so that the corresponding value M c of the meananomaly can be obtained without solving KE. For any such choice of ( e c , M c ) , an infinite Taylor series expansion inboth M and e can then be written, which is expected to converge in a suitable neighborhood of ( e c , M c ) . A rule ofthumb for estimating the actual size of the region of convergence has also been given.Three explicit examples of such series have then been provided, two for the elliptic and one for the hyperbolic KE.Each of them, for fixed base point, turns out to converge in large parts of the ( e, M ) plane. For ( e, M ) close to ( e c , M c ) within a range ∆ e ∼ ∆ M/π = O (10 − ) , the polynomial obtained by truncating the infinite series up tothe fifth degree reaches an accuracy at the level of machine double precision. More far away from ( e c , M c ) , but stillwithin the region of convergence, higher order terms should be introduced to maintain such an accuracy.Since these new solutions converge locally around ( e c , M c ) , a suitable set of them, centered around dif-ferent ( e c , M c ) and truncated up to a certain degree, can be used to design an algorithm for the numeri-cal computation of the function E ( e, M ) for every value of ( e, M ) . The resulting polynomials will forma 2-D spline [Tommasini and Olivieri(2021)], generalizing the 1-D spline that has been proposed in Refs.[Tommasini and Olivieri(2020a), Tommasini and Olivieri(2020b)] for solving KE for every M when e is fixed. This6 PREPRINT - F
EBRUARY
23, 2021 e , a M + b −14 −11 −8 −5 −2 e rr o r ( r a d ) Figure 3: Errors E n ( e, M ) for the approximate polynomial solutions S n ( e, M ) of KE, obtained by truncating theinfinite series of Eq. (22) up to degree n , for n = 1 , · · · , . The errors have been computed along the two linesconnecting the points (0 , , ( , π − ) , and (1 , π ) of the ( e, M ) plane, so that a = π − and b = 0 for e < . , and a = b = π +1 for e > . . (Notice that here only the vertical axis is logarithmic, while the horizontal axis is chosen tobe linear.)bivariate spline can be used for accelerating computations involving the repetitive solution of Kepler’s equation for sev-eral different values of e and M [Tommasini and Olivieri(2021)], as for exoplanet search [Makarov and Veras(2019),Eastman et al.(2019) ] or for the implementation of Enke’s method [Roy(2005)]. I thank David N. Olivieri for discussions. This work was supported by grants 67I2-INTERREG, from Axencia Galegade Innovación, Xunta de Galicia, and FIS2017-83762-P from Ministerio de Economia, Industria y Competitividad,Spain.
References [Roy(2005)] Roy, A.E.
Orbital Motion , 4 ed.; Institute of Physics Publishing, Bristol and Philadelphia, 2005.[Colwell(1993)] Colwell, P.
Solving Kepler’s Equation Over Three Centuries ; Willmann-Bell Inc., Richmond, VA,1993.[Lagrange(1771)] Lagrange, J. Sur le problème de Kepler.
Memoires de l’Academie Royale des Sciences (Berlin) , , 204–233. 7 PREPRINT - F
EBRUARY
23, 2021 e M ( r a d ) −15 −12 −9 −6 −3 Figure 4: Contour levels of the error E affecting the fifth degree polynomial approximation of Eq. (22), as a functionof the eccentricity e and the mean anomaly M . The continuous black curve marks the boundary of the region ofconvergence, as estimated with Eq. (20). (Notice that here the axes for e and M are linear.)[Bessel(1805)] Bessel, F. Ueber die Berechnung der wahren Anomalie in einer von der Parabel nicht sehr verschei-denen Bahn. Monatliche Correspondenz zur Beforderung der Erd-und Himmels-Kunde herausgegeben vonFreiherrn von Zach , , 197–207.[Fernande(1994)] Fernande, S.D.S. Extension of the solution of kepler’s equation to high eccentricities. CelestialMechanics and Dynamical Astronomy , , 297–308.[Levi-Civita(1904a)] Levi-Civita, T. Sopra la equazione di Kepler. Astronomische Nachrichten , , 313–314.[Levi-Civita(1904b)] Levi-Civita, T. Sopra la equazione di Kepler. Rendiconti della Reale Accademia dei Lincei,Classe di scienze fisiche, matematiche e naturali E , , 260–268.[Stumpff(1968)] Stumpff, K. On the application of Lie-series to the problems of celestial mechanics. NationalAeronautics and Space Administration , Technical Note D-4460 .[Tommasini and Olivieri(2021)] Tommasini, D.; Olivieri, D.N. Univariate and Bivariate Quintic Spline Algorithmsfor Fast and Accurate Solution of the Elliptic Kepler Equation for All Values of the Eccentricity and MeanAnomaly (in preparation), 2021.[Tommasini and Olivieri(2020a)] Tommasini, D.; Olivieri, D.N. Fast Switch and Spline Scheme for Accurate Inver-sion of Nonlinear Functions: The New First Choice Solution to Kepler’s Equation.
Applied Mathematics andComputation , , 124677. .[Tommasini and Olivieri(2020b)] Tommasini, D.; Olivieri, D.N. Fast Switch and Spline Function Inversion Algo-rithm with Multistep Optimization and k-Vector Search for Solving Kepler’s Equation in Celestial Mechanics. Mathematics , , 2017. . 8 PREPRINT - F
EBRUARY
23, 2021 −3 −2 −1 e − 2 , M −16 −13 −10 −7 −4 −1 e rr o r Figure 5: Errors E n ( e, M ) for the approximate polynomial solutions S n ( e, M ) of the hyperbolic KE, obtained bytruncating the infinite series of Eq. (23), as functions of e − (in logarithmic scale). The errors have been computedalong the line M = e − .[Makarov and Veras(2019)] Makarov, V.V.; Veras, D. Chaotic Rotation and Evolution of Asteroids and Small Planetsin High-eccentricity Orbits around White Dwarfs. The Astrophysical Journal , , 1–7. .[Eastman et al.(2019)] Eastman, J.D.; Rodriguez, J.E.; Agol, E.; Stassun, K.G.; Beatty, T.G.; Vanderburg, A.; Gaudi,S.; Collins, K.A.; Luger, R. EXOFASTv2: A public, generalized, publication-quality exoplanet modelingcode, 2019. 9
PREPRINT - F
EBRUARY