Continue the Series Solution to Keplers Equation Begun in Equation 267

Abstract

A new approach for solving Kepler equation for elliptical orbits is developed in this paper. This new approach takes advantage of the very good behaviour of the modified Newton–Raphson method when the initial seed is close to the looked for solution. To determine a good initial seed the eccentric anomaly domain [0, π] is discretized in several intervals and for each one of these intervals a fifth degree interpolating polynomial is introduced. The six coefficients of the polynomial are obtained by requiring six conditions at both ends of the corresponding interval. Thus the real function and the polynomial have equal values at both ends of the interval. Similarly relations are imposed for the two first derivatives. In the singular corner of the Kepler equation, M ≪ 1 and 1 −e ≪ 0 an asymptotic expansion is developed. In most of the cases, the seed generated leads to reach machine error accuracy with the modified Newton–Raphson method with no iterations or just one iteration. This approach improves the computational time compared with other methods currently in use.

1 INTRODUCTION

The basic goal of the Kepler equation (KE) is to provide the value of the eccentric anomaly E in terms of the mean anomaly M and the eccentricity e of the orbit:

\begin{equation*} M = E -e\, \sin E. \end{equation*}

Over the centuries the resolution of the KE has been studied for a wide variety of scientists such that in virtually every decade from 1650 to the present there have appeared papers devoted to the Kepler problem and its solution (Colwell1993). Although there exists a unique solution of the KE, there are many methods to describe or approximate it. Basically it depends on the solvers's motivations and the mathematical tools available in the epoch of study.

Before 1860 the work was motivated directly by celestial mechanics, but it showed results less sensitive to the needs of astronomers. After 1860 the development and refinement of what we now call the special functions of applied mathematics became the involvement of mathematicians in solving KE more outstanding. In the case of the astronomers, they were generally more interested in being able to solve KE quickly and repeatedly, but with only modest accuracy (Rambaut1906). From 1860 to sometime after 1890, the work on KE was mostly derivative and computational (Colwell1993). We can highlight the work developed by Cassini (1719) providing a non-analytic solution, the infinite series solution developed by Levi-Civita (1904) or the iterative methods based on the Newton's Principia (1686).

From mid-20th century, the classic and more recent methods focused on the development of optimal procedures where the accuracy of the algorithm and the computational effort become the main aspects to consider. Most of the algorithms are iterative and they depend on how the starting estimate for E is determined. Nevertheless, nowadays it is possible to find analytical algorithms as the method described in Pál (2009) to solve Kepler's problem. In what follows we briefly describe some of the approaches currently used to solve equation (2).

• The Serafin approach described in Serafin (1986). The bounds for E are chosen through several inequalities to successively divide the root interval and improve the lower and upper bounds. To do that, two linear methods based on elementary properties of convex functions are used.

• The Gooding and Odell approach described in Gooding & Odell (1988). The starter value for E is provided by solving a cubic equation. After that, a corrected value is given through the Halley's method which is used as input to approximate the evaluation of equation (2) and the derivative by their Taylor series. The final value is computed with the Newton–Raphson method. It is an iterative method which requires several transcendental function evaluations.

• The Nijenhuis approach described in Nijenhuis (1991). The starter for E is computed depending on four different regions of (M, e)-space. In three of these regions the starter value is defined with a Halley step based on a version of the KE in which the sine function is replaced by a simple approximation. In the remaining region, which contains the singular point, the starter value is given by the solution of Mikkola's cubic equation and refined through a Newton step applied to quintic equation. An mth order generalized Newton method is applied thereafter to obtain the final solution, where m depends on the accuracy required. It is an iterative method which requires two trigonometric evaluations.

• The Markley approach described in Markley (1995). The starter point for E is the root of a cubic equation, whose expression is derived through a Padé approximation for sinE, which is replaced in equation (2). After that, a fifth-order refinement of the solution is applied only once. It is a no iterative method which requires four transcendental function evaluations (a square root, a cube root and two trigonometric functions).

• The Fukushima approach described in Fukushima (1996). The starter value for E is a trivial upper bound of the solution: π. The approximate solution is determined by transforming the Newton–Raphson method from the continuous x-space to a discretized j-space, where x represents E and j the solution index. The corrected value for E is another approximation of the Newton–Raphson method by approximating the evaluation of equation (2) and the derivative by their truncated Taylor series around the approximate solution. It is an iterative method which does not require transcendental functions.

• The Mortari and Elipe approach described in Mortari & Elipe (2014). Two ranges of mean anomaly are defined to select where the solution should be determined. The lower bound for E is derived through two implicit functions, which are non-rational Bézier functions, linear or quadratic, depending on the derivatives of the initial bound values. The upper bound for E is estimated by applying the Newton–Raphson method with the lower bound as the starting value. After that, if the machine error accuracy is not reached, the lower and upper bounds define a new range of searching. It is an iterative method which requires several transcendental function evaluations (one quadratic iteration requires 19 additions, 20 products and one square root, while one linear iteration requires seven additions and five products).

In the last century, with the advent of calculators and computers, there is no impediment to achieving quick solutions of great accuracy. In particular, the symbolic manipulators like matlab, maple or mathematica are very powerful calculators which are easy to use and have a very intuitive syntax. Besides, they can be directed to provide the corresponding code in some of the standard programming languages such as c, c++ or fortran, which are specially adapted for numerical calculus. In this paper we develop a procedure to solve the KE for elliptic orbits taking advantage of the full potential of the symbolic manipulators. We focus the study on the determination of an appropriate seed to initialize the numerical method for solving the KE, considering the optimization already tested of the well-known Newton–Raphson method.

2 ELLIPTIC KEPLER EQUATION

The KE for the elliptical case,

\begin{equation} x = y-e\, \sin y, \end{equation}

(1)

determines a non-linear function y =y(x, e), where y =E is the unknown eccentric anomaly, M =x the mean anomaly which is known and e < 1 the eccentricity. For given values of e ∈ [0, 1[ and M ∈ [0, 2π] the function y =y(x, e) is univocally defined by equation (1) and this property can be deduced from the Banach fixed-point theorem.

The function g(y) is contractive if there is a non-negative real number k ∈ ]0, 1[ such that ∀y ₁, y ₂ ∈ [0, 2π]:

\begin{equation*} \vert g(y_{1})-g(y_{2})\vert \le k \vert y_{1}-y_{2}\vert . \end{equation*}

The smallest value for k is the Lipschitz constant of g.

If we consider the equation g(y) =x +e siny, the function g(y) is a contractive mapping for every value of x, since its derivative g΄(y) =e cosy satisfies the condition |g΄(y)| < 1 when e ∈ [0, 1]. In effect, applying the mean value theorem to function g(y), there exist an intermediate value y★ ∈ [y ₁, y ₂] which verifies

\begin{equation*} g(y_{2})-g(y_{1}) = g^{\prime }(y^*)(y_{2}-y_{1}). \end{equation*}

Taking into account that g΄(y) =e cosy, we have

\begin{equation*} \vert g(y_{2})-g(y_{1})\vert = \vert e \cos y^*\vert \,\vert y_{2}-y_{1}\vert \le e \vert y_{2}-y_{1}\vert , \end{equation*}

i.e. k, the Lipschitz constant of g, verifies k <e, proving the property.

Finally, applying the Banach fixed-point theorem for contractive functions, it can be concluded that there is a fixed point which verifies y =g(y) and this fixed point is the looked for solution of the KE. Starting from a given value y ₀, for example y ₀ =x, the sequence y _n + 1 =g(y _n ) turns out to be convergent and its limit is the fixed point.

However, this sequence is slowly convergent – especially when e is close to unity – and this way to solve KE is not effective, from a numerical point of view. Thus, the interest of this procedure is much more theoretical than practical.

The solution of the KE can be seen as a root finder problem of the function

\begin{equation} f(y) \equiv y-e\, \sin y-x=0 \end{equation}

(2)

for given values of e and x. Because of the transcendental character of the equation its solution is not trivial. Nevertheless, by introducing the function

\begin{equation*} F(y) = y - \frac{f(y)}{f^{\prime }(y)} = \frac{e(\sin y - y\cos y) +x}{1 - e\cos y}, \end{equation*}

where f΄(y) is the derivative of f(y), the solution of the KE reduces to the search of the fixed point of F(y); note that y =F(y)⇔f(y) = 0. Again, we can form a convergent sequence y _n + 1 =F(y _n ) whose limit is the root of the function f(y); the convergence of this new sequence is faster than the previous one (see Stumpf1999) but still is insufficient from a practical point of view.

There are several analytical solutions given by expansions in terms of Bessel functions of first kind like

\begin{eqnarray*} E &=& M + 2 \sum _{1}^{\infty } J_k(k\,e) \frac{\sin (k\,M)}{k}, \\ J_k(k\,e) &=& \frac{1}{\pi } \int _0^{\pi } \cos (k\,u - k\,e \sin u)\, \mathrm{d} u. \end{eqnarray*}

More analytical solutions can be found in references Battin (1999), Siewert & Burniston (1972) and Colwell (1993). However, these expansions are not competitive among the numerical methods of solution, as J. M. A. Danby comments in his book (Danby1992).

2.1 Domain of the equation

Initially, x and y range in the interval [0, 2π]; however, the transformation

\begin{equation} y = 2\pi - \eta , \qquad x=2\pi - \xi \end{equation}

(3)

convert equation (1) into

\begin{equation*} \xi = \eta - e\, \sin \eta . \end{equation*}

As a consequence, if x, y ∈ ]π, 2 π] ⇒ ξ, η ∈ [0, π] and the mapping (3) reduces the problem to the interval [0, π].

3 NEWTON–RAPHSON ALGORITHMS

The root searched is solution of equation (2) and we wish to establish a successive approximation method starting from a seed y ₀. In the classical Newton–Raphson (CNR) algorithm a sequence is formed whose elements become arbitrarily close to each other as the sequence progresses:

\begin{equation*} y_{n+1} = y_n - \frac{f(y_n)}{f^{\prime }(y_n)}. \end{equation*}

This Cauchy sequence converges to the solution and it is very well known that the convergence rate of this algorithm is quadratic; thus the number of correct figures is doubled, approximately, in each iteration of the sequence.

A better convergence rate is obtained with the modified Newton–Raphson (MNR) algorithm which we use in our code to solve the KE. If y _n in one of the terms of the sequence, the next term will be y _n + 1 =y _n  + Δy _n , where Δy _n should verify

\begin{equation*} f(y_{n+1})=f(y_{n}+\Delta y_{n})=0, \end{equation*}

with f given by equation (2). By approximating f(y _n + 1) by its second-order Taylor expansion about y _n ,

\begin{equation*} f(y_{n+1}) \approx f(y_{n})+f^{\prime }(y_{n})\Delta y_{n}+\frac{1}{2}f^{\prime \prime }(y_{n}) \Delta ^2 y_{n}, \end{equation*}

the condition f(y _n + 1) = 0 provides the value of Δy _n :

\begin{equation*} \Delta y_{n} = \frac{-f^{\prime }(y_{n}) \pm \sqrt{{f^{\prime }}^2(y_{n})-2f(y_{n})f^{\prime \prime }(y_{n})}}{f^{\prime \prime }(y_{n})}. \end{equation*}

Then, we multiply numerator and denominator by

\begin{equation*} -f^{\prime }(y_{n}) \mp \sqrt{{f^{\prime }}^{2}(y_{n})-2f(y_{n})f^{\prime \prime }(y_{n})} \end{equation*}

to obtain

\begin{equation} \Delta y_{n} = \frac{-2f(y_{n})}{f^{\prime }(y_{n}) \pm \sqrt{\vert {f^{\prime }}^{2}(y_{n})-2f(y_{n})f^{\prime \prime }(y_{n})\vert }}, \end{equation}

(4)

where we have to select the sign ( + ) when f΄(y _n ) is positive; if f΄(y _n ) < 0 the sign ( − ) must be taken.

The use of absolute value in equation (4) does not affect the convergence of the algorithm. It is introduced to avoid the algorithm fails during the series generation when the square root is a complex number. The first and second derivatives of equation (2) are given by the relations

\begin{eqnarray} f^{\prime }(y_{n})=1- e\, \cos y_{n}, \end{eqnarray}

(5)

\begin{eqnarray} f^{\prime \prime }(y_{n}) = e\, \sin y_{n}. \end{eqnarray}

(6)

Equation (4) summarizes the MNR algorithm. It should be noticed that the algorithm becomes the CNR method, in which the function is approximated by its first-order Taylor expansion, when f″(y _n ) = 0. In general, the convergence rate of the MNR algorithm is cubic.

The MNR algorithm is a particular case of the algorithm used in the Conway method. Bruce A. Conway, of the University of Illinois, applied a root-finding method of Laguerre (1834–1886) to the solution of KE (Conway1986). Although the method is intended for finding the roots of a polynomial, it works equally as well for KE. The algorithm associated with that method summarizes in the equation

\begin{equation*} \Delta y_{n}= \frac{-(1+p)f(y_{n})}{f^{\prime }(y_{n})\pm \sqrt{\vert p[p\,{f^{\prime }}^{2}(y_{n})-(1+p)f(y_{n})f^{\prime \prime }(y_{n})]\vert }}, \end{equation*}

with p =m − 1 and m is the degree of the polynomial. The convergence of the algorithm is almost independent of m; Conway opts for m = 5. As we can see, the MNR method is a particular case of the Conway method for m = 2.

The main difference between the CNR and the Conway or the MNR methods is associated with its reliability. The three methods need an initial seed from which start the generation of the convergent sequence. The CNR method is much more affected by the quality of the initial seed than the Conway or the MNR methods. If the quality of the initial seed is poor the CNR method could not converge to the solution.

Table1 shows the different sequences generated by the CNR and the MNR algorithms when trying to solve the KEs with values M = 2.5 and e = 0.8 and two different seeds: the bad one, y ₀ = 0.1, and another not too bad, y ₀ = 0.3. Only the first six terms of the sequence has been calculated.

Table 1.

The right answer is y = 2.78172231.

Sequences for M = 2.5 and e = 0.8
Seed y 0 = 0.1		Seed y 0 = 0.3
CNR	MNR	CNR	MNR
12.25640804	5.82975453	10.6355865	3.95106621
−29.74030805	1.97868144	3.70452359	2.86354381
0.74394754	2.76120459	2.73142300	2.78176216
6.32984624	2.78172169	2.78194606	2.78172231
−12.55060918	2.78172231	2.78172231	2.78172231
62.72807874	2.78172231	2.78172231	2.78172231

Sequences for M = 2.5 and e = 0.8
Seed y 0 = 0.1		Seed y 0 = 0.3
CNR	MNR	CNR	MNR
12.25640804	5.82975453	10.6355865	3.95106621
−29.74030805	1.97868144	3.70452359	2.86354381
0.74394754	2.76120459	2.73142300	2.78176216
6.32984624	2.78172169	2.78194606	2.78172231
−12.55060918	2.78172231	2.78172231	2.78172231
62.72807874	2.78172231	2.78172231	2.78172231

Table 1.

The right answer is y = 2.78172231.

Sequences for M = 2.5 and e = 0.8
Seed y 0 = 0.1		Seed y 0 = 0.3
CNR	MNR	CNR	MNR
12.25640804	5.82975453	10.6355865	3.95106621
−29.74030805	1.97868144	3.70452359	2.86354381
0.74394754	2.76120459	2.73142300	2.78176216
6.32984624	2.78172169	2.78194606	2.78172231
−12.55060918	2.78172231	2.78172231	2.78172231
62.72807874	2.78172231	2.78172231	2.78172231

Sequences for M = 2.5 and e = 0.8
Seed y 0 = 0.1		Seed y 0 = 0.3
CNR	MNR	CNR	MNR
12.25640804	5.82975453	10.6355865	3.95106621
−29.74030805	1.97868144	3.70452359	2.86354381
0.74394754	2.76120459	2.73142300	2.78176216
6.32984624	2.78172169	2.78194606	2.78172231
−12.55060918	2.78172231	2.78172231	2.78172231
62.72807874	2.78172231	2.78172231	2.78172231

For the bad seed the CNR does not converge to the real solution; on the contrary, the MNR always converges to the real answer even though we start from a very poor seed.

However, the behaviour of the algorithms is quite different if we use a very good seed as starting point. Table2 shows the first six terms of the sequences generated for the same values M = 2.5, e = 0.8 when the initial seed is very good: y ₀ = 2.8. Both algorithms converge quickly to the real answer.

Table 2.

The right answer is y = 2.781722308989884.

Sequences for M = 2.5 and e = 0.8
Seed y 0 = 2.8
CNR	MNR
2.781748270190563	2.781722746897827
2.781722309044170	2.781722308989884
2.781722308989884	2.781722308989884
2.781722308989884	2.781722308989884
2.781722308989884	2.781722308989884
2.781722308989884	2.781722308989884

Sequences for M = 2.5 and e = 0.8
Seed y 0 = 2.8
CNR	MNR
2.781748270190563	2.781722746897827
2.781722309044170	2.781722308989884
2.781722308989884	2.781722308989884
2.781722308989884	2.781722308989884
2.781722308989884	2.781722308989884
2.781722308989884	2.781722308989884

Table 2.

The right answer is y = 2.781722308989884.

Sequences for M = 2.5 and e = 0.8
Seed y 0 = 2.8
CNR	MNR
2.781748270190563	2.781722746897827
2.781722309044170	2.781722308989884
2.781722308989884	2.781722308989884
2.781722308989884	2.781722308989884
2.781722308989884	2.781722308989884
2.781722308989884	2.781722308989884

Sequences for M = 2.5 and e = 0.8
Seed y 0 = 2.8
CNR	MNR
2.781748270190563	2.781722746897827
2.781722309044170	2.781722308989884
2.781722308989884	2.781722308989884
2.781722308989884	2.781722308989884
2.781722308989884	2.781722308989884
2.781722308989884	2.781722308989884

The accuracy of the Newton–Raphson algorithms is jeopardized when M ≪ 1 and 1 −e ≪ 1; we call this region of the plane (e, M), the singular corner. In such a region y is also small, and the expansion of the equation for small values of y provides

\begin{equation*} y \approx \frac{x}{1-e}, \end{equation*}

i.e. the solution is bad determined from a numerical point of view since it is the ratio of two infinitesimals.

The central idea, and also the goal of this work, is to generate a very good initial seed in all possible cases, also in the singular corner. Then we will use one of the Newton–Raphson algorithms to find the right solution of the KE. Roughly, we exploit the very good characteristics of the Newton–Raphson algorithms (CNR, MNR or Conway) because they assure the convergence and provide the solution with accuracy and velocity if the quality of the initial seed is high.

In this paper, we try to demonstrate that it is much more cost-effective, from a computational point of view, to improve the initial seed and then use it to feed a Newton–Raphson algorithm than to look for new algorithms which would be able to beat the simplicity and performances of the Newton–Raphson algorithms.

Moreover, the solution of the KE is an issue with more than three centuries (see Colwell1993); along such a long period of time, several methods have been proposed most of them associated with mathematical tools or techniques of its epoch. At present, one important mathematical tool which characterizes this epoch is the symbolic manipulator (macsyma, maxima, derive, maple, mathematica, matlab, etc.). We try to exploit these tools to the full, in order to generate a c or fortran code robust, reliable, exact and fast.

4 THE SEED VALUE

Our approach starts discretizing the eccentric anomaly domain [0, π] in 12 intervals of 15° of longitude, where the corresponding M _i ≡ x _i is obtained from equation (1) (see Table3). This way, we defined 12 intervals [x _i , x _i + 1], i = 1, …, 12.

In each one of these intervals we introduce a fifth-degree polynomial p _i (x, e),i = 1, …, 12 to interpolate the eccentric anomaly:

\begin{eqnarray*} p_i(x,e)&=&a^{(i)}_{0}(e)+a^{(i)}_{1}(e)x+a^{(i)}_{2}(e)x^{2}+a^{(i)}_{3}x^{3}(e)\nonumber\\ &&+a^{(i)}_{4}(e)x^{4}+a^{(i)}_{5}(e)x^{5}. \end{eqnarray*}

Notice that the polynomials used in the method have coefficients that depend on the eccentricity (see AppendixA). In order to simplify the equation, we abuse notation by writing the coefficients only as a function of i, i.e. the polynomial they are related,

\begin{equation} p_i(x)=a^{(i)}_{0}+a^{(i)}_{1}x+a^{(i)}_{2}x^{2}+a^{(i)}_{3}x^{3}+a^{(i)}_{4}x^{4}+a^{(i)}_{5}x^{5}. \end{equation}

(7)

The polynomial p _i (x) is only valid inside the interval [x_i , x _{i + 1}],i = 1, …, 12, where we approximate the function by the polynomial: y(x) ≈ p _i (x).

In order to determine the six coefficients of p _i (x) the following six boundary conditions should be imposed

\begin{eqnarray*} p_i(x_{i})&=& y(x_{i}), \qquad \quad p_i(x_{i+1})=y(x_{i+1}), \\ p^{\prime }_i(x_{i})&=& y^{\prime }(x_{i}), \qquad \quad p^{\prime }_i(x_{i+1})=y^{\prime }(x_{i+1}), \\ p^{\prime \prime }_i(x_{i})&=& y^{\prime \prime }(x_{i}), \qquad \quad p^{\prime \prime }_i(x_{i+1})=y^{\prime \prime }(x_{i+1}), \end{eqnarray*}

where the two first derivatives of y are given by

\begin{eqnarray} y^{\prime } = \frac{\mathrm{d} y}{\mathrm{d} x} =\frac{1}{1-e\, \cos y}, \end{eqnarray}

(8)

\begin{eqnarray} y^{\prime \prime } = \frac{\mathrm{d}^2 y}{\mathrm{d} x^2} =\frac{-e\, \sin y}{(1- e\, \cos y)^{3}}, \end{eqnarray}

(9)

and they can be easily calculated at the ends of the interval. These conditions assure that the function y(x) and the polynomial p _i (x) share the same values and the same first and second derivatives at the ends of the interval [x _i , x _i + 1].

The polynomials p _i (x),i = 1, …, 12, can be generated easily with the help of a symbolic manipulator only once; moreover, the symbolic manipulator provides the fortran or c code of the polynomials that will be appropriately stored in one module to be used later. We use the version 18.02 of maple.

Consequently, given e and M, we determine the starter value y ₀ following the next procedure.

• First, given e the ends x _i , x _i + 1 of the intervals described in Table3 are calculated.

• The interval [x _i , x _i + 1] should be selected in such a way that M ∈ [x _i , x _i + 1]. Note that the corresponding values [y _i , y _i + 1] are known and the seed y ₀ satisfies: y ₀ ∈ [y _i , y _i + 1].

• Use the corresponding polynomial p _i (x) to obtain the seed to be used with the Newton–Raphson algorithm:

  \begin{equation} y_0 = p_i(M). \end{equation}

  (10)

Table 3.

Important values for the polynomial approach, where γ _is = sin ((i − 1)π/12) and γ _ic = cos ((i − 1)π/12).

Intervals defined in the solution
i	x i  =y i  −esiny i	|$y_{i}=(i-1)\frac{\pi }{12}$|	|$\frac{\mathrm{d} y}{\mathrm{d} x}\Big |_{x_i}$|	|$\frac{\mathrm{d}^2 y}{\mathrm{d} x^2}\Big |_{x_i}$|
1	0	0	|$\frac{1}{1-e}$|	0
2	|$\frac{\pi }{12}-\gamma _{2s}e$|	|$\frac{\pi }{12}$|	|$\frac{1}{1-\gamma _{2c}e}$|	|$\frac{-\gamma _{2s}e}{(1-\gamma _{2c}e)^3}$|
3	|$\frac{\pi }{6}-\frac{1}{2}e$|	|$\frac{\pi }{6}$|	|$\frac{1}{1-\gamma _{3c}e}$|	|$\frac{-e}{2(1-\gamma _{3c}e)^3}$|
4	|$\frac{\pi }{4}-\gamma _{4s}e$|	|$\frac{\pi }{4}$|	|$\frac{1}{1-\gamma _{4c}e}$|	|$\frac{-\gamma _{4s}e}{(1-\gamma _{4c}e)^3}$|
5	|$\frac{\pi }{3}-\gamma _{5s}e$|	|$\frac{\pi }{3}$|	|$\frac{2}{2-e}$|	|$\frac{-8\gamma _{5s}e}{(2-e)^3}$|
6	|$\frac{5\pi }{12}-\gamma _{6s}e$|	|$\frac{5\pi }{12}$|	|$\frac{1}{1-\gamma _{6c}e}$|	|$\frac{-\gamma _{6s}e}{(1-\gamma _{6c}e)^3}$|
7	|$\frac{\pi }{2}-e$|	|$\frac{\pi }{2}$|	1	−e
8	|$\frac{7\pi }{12}-\gamma _{6s}e$|	|$\frac{7\pi }{12}$|	|$\frac{1}{1+\gamma _{6c}e}$|	|$\frac{-\gamma _{6s}e}{(1+\gamma _{6c}e)^3}$|
9	|$\frac{2\pi }{3}-\gamma _{5s}e$|	|$\frac{2\pi }{3}$|	|$\frac{2}{2+e}$|	|$\frac{-8\gamma _{5s}e}{(2+e)^3}$|
10	|$\frac{3\pi }{4}-\gamma _{4s}e$|	|$\frac{3\pi }{4}$|	|$\frac{1}{1+\gamma _{4c}e}$|	|$\frac{-\gamma _{4s}e}{(1+\gamma _{4c}e)^3}$|
11	|$\frac{5\pi }{6}-\frac{1}{2}e$|	|$\frac{5\pi }{6}$|	|$\frac{1}{1+\gamma _{3c}e}$|	|$\frac{-e}{2(1+\gamma _{3c}e)^{3}}$|
12	|$\frac{11\pi }{12}-\gamma _{2s}e$|	|$\frac{11\pi }{12}$|	|$\frac{1}{1+\gamma _{2c}e}$|	|$\frac{-\gamma _{2s}e}{(1+\gamma _{2c}e)^3}$|
13	π	π	|$\frac{1}{1+e}$|	0

Intervals defined in the solution
i	x i  =y i  −esiny i	|$y_{i}=(i-1)\frac{\pi }{12}$|	|$\frac{\mathrm{d} y}{\mathrm{d} x}\Big |_{x_i}$|	|$\frac{\mathrm{d}^2 y}{\mathrm{d} x^2}\Big |_{x_i}$|
1	0	0	|$\frac{1}{1-e}$|	0
2	|$\frac{\pi }{12}-\gamma _{2s}e$|	|$\frac{\pi }{12}$|	|$\frac{1}{1-\gamma _{2c}e}$|	|$\frac{-\gamma _{2s}e}{(1-\gamma _{2c}e)^3}$|
3	|$\frac{\pi }{6}-\frac{1}{2}e$|	|$\frac{\pi }{6}$|	|$\frac{1}{1-\gamma _{3c}e}$|	|$\frac{-e}{2(1-\gamma _{3c}e)^3}$|
4	|$\frac{\pi }{4}-\gamma _{4s}e$|	|$\frac{\pi }{4}$|	|$\frac{1}{1-\gamma _{4c}e}$|	|$\frac{-\gamma _{4s}e}{(1-\gamma _{4c}e)^3}$|
5	|$\frac{\pi }{3}-\gamma _{5s}e$|	|$\frac{\pi }{3}$|	|$\frac{2}{2-e}$|	|$\frac{-8\gamma _{5s}e}{(2-e)^3}$|
6	|$\frac{5\pi }{12}-\gamma _{6s}e$|	|$\frac{5\pi }{12}$|	|$\frac{1}{1-\gamma _{6c}e}$|	|$\frac{-\gamma _{6s}e}{(1-\gamma _{6c}e)^3}$|
7	|$\frac{\pi }{2}-e$|	|$\frac{\pi }{2}$|	1	−e
8	|$\frac{7\pi }{12}-\gamma _{6s}e$|	|$\frac{7\pi }{12}$|	|$\frac{1}{1+\gamma _{6c}e}$|	|$\frac{-\gamma _{6s}e}{(1+\gamma _{6c}e)^3}$|
9	|$\frac{2\pi }{3}-\gamma _{5s}e$|	|$\frac{2\pi }{3}$|	|$\frac{2}{2+e}$|	|$\frac{-8\gamma _{5s}e}{(2+e)^3}$|
10	|$\frac{3\pi }{4}-\gamma _{4s}e$|	|$\frac{3\pi }{4}$|	|$\frac{1}{1+\gamma _{4c}e}$|	|$\frac{-\gamma _{4s}e}{(1+\gamma _{4c}e)^3}$|
11	|$\frac{5\pi }{6}-\frac{1}{2}e$|	|$\frac{5\pi }{6}$|	|$\frac{1}{1+\gamma _{3c}e}$|	|$\frac{-e}{2(1+\gamma _{3c}e)^{3}}$|
12	|$\frac{11\pi }{12}-\gamma _{2s}e$|	|$\frac{11\pi }{12}$|	|$\frac{1}{1+\gamma _{2c}e}$|	|$\frac{-\gamma _{2s}e}{(1+\gamma _{2c}e)^3}$|
13	π	π	|$\frac{1}{1+e}$|	0

Table 3.

Important values for the polynomial approach, where γ _is = sin ((i − 1)π/12) and γ _ic = cos ((i − 1)π/12).

Intervals defined in the solution
i	x i  =y i  −esiny i	|$y_{i}=(i-1)\frac{\pi }{12}$|	|$\frac{\mathrm{d} y}{\mathrm{d} x}\Big |_{x_i}$|	|$\frac{\mathrm{d}^2 y}{\mathrm{d} x^2}\Big |_{x_i}$|
1	0	0	|$\frac{1}{1-e}$|	0
2	|$\frac{\pi }{12}-\gamma _{2s}e$|	|$\frac{\pi }{12}$|	|$\frac{1}{1-\gamma _{2c}e}$|	|$\frac{-\gamma _{2s}e}{(1-\gamma _{2c}e)^3}$|
3	|$\frac{\pi }{6}-\frac{1}{2}e$|	|$\frac{\pi }{6}$|	|$\frac{1}{1-\gamma _{3c}e}$|	|$\frac{-e}{2(1-\gamma _{3c}e)^3}$|
4	|$\frac{\pi }{4}-\gamma _{4s}e$|	|$\frac{\pi }{4}$|	|$\frac{1}{1-\gamma _{4c}e}$|	|$\frac{-\gamma _{4s}e}{(1-\gamma _{4c}e)^3}$|
5	|$\frac{\pi }{3}-\gamma _{5s}e$|	|$\frac{\pi }{3}$|	|$\frac{2}{2-e}$|	|$\frac{-8\gamma _{5s}e}{(2-e)^3}$|
6	|$\frac{5\pi }{12}-\gamma _{6s}e$|	|$\frac{5\pi }{12}$|	|$\frac{1}{1-\gamma _{6c}e}$|	|$\frac{-\gamma _{6s}e}{(1-\gamma _{6c}e)^3}$|
7	|$\frac{\pi }{2}-e$|	|$\frac{\pi }{2}$|	1	−e
8	|$\frac{7\pi }{12}-\gamma _{6s}e$|	|$\frac{7\pi }{12}$|	|$\frac{1}{1+\gamma _{6c}e}$|	|$\frac{-\gamma _{6s}e}{(1+\gamma _{6c}e)^3}$|
9	|$\frac{2\pi }{3}-\gamma _{5s}e$|	|$\frac{2\pi }{3}$|	|$\frac{2}{2+e}$|	|$\frac{-8\gamma _{5s}e}{(2+e)^3}$|
10	|$\frac{3\pi }{4}-\gamma _{4s}e$|	|$\frac{3\pi }{4}$|	|$\frac{1}{1+\gamma _{4c}e}$|	|$\frac{-\gamma _{4s}e}{(1+\gamma _{4c}e)^3}$|
11	|$\frac{5\pi }{6}-\frac{1}{2}e$|	|$\frac{5\pi }{6}$|	|$\frac{1}{1+\gamma _{3c}e}$|	|$\frac{-e}{2(1+\gamma _{3c}e)^{3}}$|
12	|$\frac{11\pi }{12}-\gamma _{2s}e$|	|$\frac{11\pi }{12}$|	|$\frac{1}{1+\gamma _{2c}e}$|	|$\frac{-\gamma _{2s}e}{(1+\gamma _{2c}e)^3}$|
13	π	π	|$\frac{1}{1+e}$|	0

Intervals defined in the solution
i	x i  =y i  −esiny i	|$y_{i}=(i-1)\frac{\pi }{12}$|	|$\frac{\mathrm{d} y}{\mathrm{d} x}\Big |_{x_i}$|	|$\frac{\mathrm{d}^2 y}{\mathrm{d} x^2}\Big |_{x_i}$|
1	0	0	|$\frac{1}{1-e}$|	0
2	|$\frac{\pi }{12}-\gamma _{2s}e$|	|$\frac{\pi }{12}$|	|$\frac{1}{1-\gamma _{2c}e}$|	|$\frac{-\gamma _{2s}e}{(1-\gamma _{2c}e)^3}$|
3	|$\frac{\pi }{6}-\frac{1}{2}e$|	|$\frac{\pi }{6}$|	|$\frac{1}{1-\gamma _{3c}e}$|	|$\frac{-e}{2(1-\gamma _{3c}e)^3}$|
4	|$\frac{\pi }{4}-\gamma _{4s}e$|	|$\frac{\pi }{4}$|	|$\frac{1}{1-\gamma _{4c}e}$|	|$\frac{-\gamma _{4s}e}{(1-\gamma _{4c}e)^3}$|
5	|$\frac{\pi }{3}-\gamma _{5s}e$|	|$\frac{\pi }{3}$|	|$\frac{2}{2-e}$|	|$\frac{-8\gamma _{5s}e}{(2-e)^3}$|
6	|$\frac{5\pi }{12}-\gamma _{6s}e$|	|$\frac{5\pi }{12}$|	|$\frac{1}{1-\gamma _{6c}e}$|	|$\frac{-\gamma _{6s}e}{(1-\gamma _{6c}e)^3}$|
7	|$\frac{\pi }{2}-e$|	|$\frac{\pi }{2}$|	1	−e
8	|$\frac{7\pi }{12}-\gamma _{6s}e$|	|$\frac{7\pi }{12}$|	|$\frac{1}{1+\gamma _{6c}e}$|	|$\frac{-\gamma _{6s}e}{(1+\gamma _{6c}e)^3}$|
9	|$\frac{2\pi }{3}-\gamma _{5s}e$|	|$\frac{2\pi }{3}$|	|$\frac{2}{2+e}$|	|$\frac{-8\gamma _{5s}e}{(2+e)^3}$|
10	|$\frac{3\pi }{4}-\gamma _{4s}e$|	|$\frac{3\pi }{4}$|	|$\frac{1}{1+\gamma _{4c}e}$|	|$\frac{-\gamma _{4s}e}{(1+\gamma _{4c}e)^3}$|
11	|$\frac{5\pi }{6}-\frac{1}{2}e$|	|$\frac{5\pi }{6}$|	|$\frac{1}{1+\gamma _{3c}e}$|	|$\frac{-e}{2(1+\gamma _{3c}e)^{3}}$|
12	|$\frac{11\pi }{12}-\gamma _{2s}e$|	|$\frac{11\pi }{12}$|	|$\frac{1}{1+\gamma _{2c}e}$|	|$\frac{-\gamma _{2s}e}{(1+\gamma _{2c}e)^3}$|
13	π	π	|$\frac{1}{1+e}$|	0

Most of the times, this initial seed is a very good approximation to the solution of the KE; especially when M is close to the ends of the interval [x _i , x _i + 1] since this seed is practically the solution of the KE. However, it is well known the singular behaviour the KE in the singular corner. In this zone the approximation of y(x) by means of the polynomials p _i (x) it is not very effective and a special analysis should be performed.

5 ANALYSIS IN THE SINGULAR CORNER

To describe the solution close to the singular corner we introduce the value ε = 1 −e assuming that ε ≪ 1:

\begin{equation} x = y - (1 - {\epsilon })\, \sin y = y - \sin y + {\epsilon }\, \sin y. \end{equation}

(11)

Our goal is to describe numerically the exact solution y _v (x) of equation (11) with enough accuracy to be part of the seed used to start the Newton–Raphson convergent process. In order to do that, an asymptotic expansion in power of the small parameter ε will be obtained.

In the singular corner, when x is small, y is also small. However it is possible to distinguish several regimes inside this region.

• The inner region, where x ≈ O(ε²) and y ≈ O(ε). In this inner region the most important term on the right-hand side of equation (11) is εsiny.

• The intermediate region, where |$x \approx {\rm O}({\epsilon }^{\frac{3}{2}})$| and |$y \approx {\rm O}({\epsilon }^{\frac{1}{2}})$|⁠. In this intermediate region both terms on the right-hand side of equation (11) are of the same order: y − siny ≈ ε siny.

• The outer region, where |$x \approx {\rm O}({\epsilon }^{\frac{3}{4}})$| and |$y \approx {\rm O}({\epsilon }^{\frac{1}{4}})$|⁠. In this outer region the most important term on the right-hand side of equation (11) is y − siny.

In each one of these regions it is possible to obtain an asymptotic expansion of the exact solution y _v (x) of equation (11) in terms of some power of the small parameter ε. During the study we assessed the boundaries of the regions, obtaining the best accuracy when

• the inner region is defined, approximately, for (x, e) satisfying |$x < 0.001\,{\epsilon }^{\frac{3}{2}}$|⁠;

• the intermediate and outer region are defined in the same region for (x, e) satisfying |$x \ge 0.001\,{\epsilon }^{\frac{3}{2}}$|⁠.

As a consequence, only two regions in the singular corner are considered at the end: inner and intermediate-outer; these approximate solutions permit to generate a very good seed to feed the Newton–Raphson convergent process (see Fig. 1).

Figure 1.

The exact solution of the KE and the asymptotic solution – up to third order – obtained in the singular corner. The agreement is excellent. Here five values of ε have been considered: ε = 0.001, 0.051, 0.101, 0.151, 0.201.

5.1 Inner region

In the inner region we introduce the following change of variables:

\begin{equation} x = {\epsilon }^2\, \xi , \qquad y = {\epsilon }\, \eta , \qquad \mbox{with} \qquad \xi , \eta \approx {\rm O}(1). \end{equation}

(12)

Note that (ξ, η) are both positive amounts. We look for a solution of equation (11) given by

\begin{equation*} \xi = \eta + a_1 {\epsilon }+ a_2 {\epsilon }^2 + a_3 {\epsilon }^3 + \cdots , \end{equation*}

where the coefficients a_i , i = 1, 2, … can be obtained expanding appropriately the right-hand side of equation (11). The solution is

\begin{equation} \xi = \eta + \frac{1}{6}\eta ^3\, {\epsilon }- \frac{1}{6}\eta ^3\, {\epsilon }^2 - \frac{1}{120}\eta ^5\, {\epsilon }^3 + \frac{1}{120}\eta ^5\, {\epsilon }^4 + \cdots . \end{equation}

(13)

The next step is to invert this solution to obtain η as a function of ξ; to do that we introduce the expansion

\begin{equation*} \eta = \xi + b_1 {\epsilon }+ b_2 {\epsilon }^2 + b_3 {\epsilon }^3 + \cdots , \end{equation*}

where the coefficients b_i , i = 1, 2, … are functions of ξ that must be calculated. Introducing this expansion in equation (13) and requiring that the different coefficients of each order match it is possible to obtain the values of b_i , i = 1, 2, …. The result is

\begin{eqnarray} \eta \!=\! \xi \!-\! \frac{1}{6}\xi ^3\, {\epsilon }\!+\! \frac{1}{12} \xi ^3 \,(\xi ^2+2)\, {\epsilon }^2 \!-\!\frac{1}{360}\xi ^5\,(20\xi ^2+57)\, {\epsilon }^3 \!+\! \cdots . \nonumber\\ \end{eqnarray}

(14)

5.2 Intermediate-outer region

In the intermediate region we introduce the following change of variables:

\begin{equation} x = {\epsilon }^{\frac{3}{2}}\, \chi , \qquad y = {\epsilon }^{\frac{1}{2}}\, \sigma , \qquad \mbox{with} \qquad \chi , \sigma \approx {\rm O}(1). \end{equation}

(15)

Note that (χ, σ) are both positive amounts. We look for a solution of equation (11) given by

\begin{equation*} \chi = a_0 + a_1 {\epsilon }+ a_2 {\epsilon }^2 + a_3 {\epsilon }^3 + \cdots , \end{equation*}

where the coefficients a_i , i = 1, 2, … can be obtained expanding appropriately the right-hand side of equation (11). The solution is

\begin{eqnarray} \chi \!=\! \frac{1}{6}\,\sigma ^3 \!+\! \sigma \!-\! \frac{1}{120}\sigma ^3\,(\sigma ^2+20)\, {\epsilon }\!+\!\frac{1}{5040}\sigma ^5 (\sigma ^2+42)\, {\epsilon }^2 \!+\! \cdots . \nonumber\\ \end{eqnarray}

(16)

Let σ₀ > 0 be the real and positive solution of the cubic equation

\begin{eqnarray} f(\sigma ) = \frac{1}{6}\,\sigma ^3 \!+\! \sigma \!-\! \chi =0 \Leftrightarrow f(\sigma ) \!=\! \sigma ^3 \!+\! 6\,\sigma - 6\,\chi =0. \nonumber\\ \end{eqnarray}

(17)

This root always exist; in effect, the function f(σ) verifies

\begin{eqnarray*} f(0) &=& - 6\,\chi < 0, \\ f(\chi ) &=& \chi ^3 > 0, \end{eqnarray*}

therefore in the interval [0, χ] there is a root of f(σ) = 0 at least. But there is only one root, because the derivative

\begin{equation*} f^{\prime }(\sigma ) = 3\,\sigma ^2 + 6 >0 \end{equation*}

is always positive in that interval. That root is given by

\begin{equation} \sigma _0 = S-\frac{2}{S}, \end{equation}

(18)

where |$S=(\Lambda + 3\chi )^{\frac{1}{3}}$| and |$\Lambda = \sqrt{8+9\chi ^2}$|⁠. It should be noticed that expression (18) can present a rounding-error defect when |$S\approx \sqrt{2}$| (σ₀ ≪ 1). In order to minimize this error, we develop an alternative procedure to provide σ₀. We rewrite the cubic equation (17) as

\begin{equation*} \sigma (\sigma ^2 + 6) - 6\,\chi =0 \Leftrightarrow \sigma =\frac{6\chi }{\sigma ^2+6}, \end{equation*}

such that σ is replaced by equation (18) on the right-hand side of the equation, obtaining

\begin{equation*} \sigma _0=\frac{6\chi }{2+S^2+\frac{4}{S^2}} \end{equation*}

which provides a smaller rounding-error when χ ≪ 1, i.e. when |$\Lambda \approx 2\,\sqrt{2}$|⁠.

Once the value of σ₀ is known it is possible to invert the solution (16) to obtain σ as a function of χ; to do that we introduce the expansion

\begin{equation*} \sigma = \sigma _0 + b_1 {\epsilon }+ b_2 {\epsilon }^2 + b_3 {\epsilon }^3 + \cdots , \end{equation*}

where the coefficients b_i , i = 1, 2, … are functions of σ₀ that must be calculated. Introducing this expansion in equation (16) and requiring that the different coefficients of each order match it is possible to obtain the values of b_i , i = 1, 2, …. The result is

\begin{eqnarray} \sigma \!=\! \sigma _0 \!+\! \frac{\sigma _0^3\,(\sigma _0^2+20)}{60(\sigma _0^2+2)}\, {\epsilon }\!+\! \frac{\sigma _0^5\,(\sigma _0^6 \!+\! 25 \sigma _0^4 \!+\! 340 \sigma _0^2\!+\!840)}{1400(\sigma _0^2+2)^3 } {\epsilon }^2 + \cdots . \nonumber\\ \end{eqnarray}

(19)

6 ANALYSIS OF RESULTS

We performed an exhaustive numerical analysis of the algorithm that we propose in these pages. Two main aspects have been considered: (1) accuracy and (2) speed. There are additional considerations that must be taken into account like e.g. the reliability of the procedure, the ease of codification and so on. We start analysing the speed.

6.1 Speed of the process

The KE for the elliptical case has been solved with the method proposed in this paper by using the MNR algorithm. The KE has been solved ≈4 × 10⁶ times, using double and quadruple precision and a tolerance equal to ε_tol = 1.11 × 10⁻¹⁵ and 1.0 × 10⁻²⁴, respectively. It should be noticed that when working with tolerances close to the zero of the machine, any numerical procedure based on consecutive approximations could be affected by artificial numerical chaos. By using an appropriate tolerance it is possible to escape from such a numerical chaos. In each run we count the number of iterations needed to reach a solution with a residual ρ = |y −e siny −x| lower than the tolerance ε_tol. This number is taken as a measure of the speed of the procedure.

Fig.2 shows in the plane (e, M) the number of iterations needed to reach a residual ρ lower than the tolerance ε_tol; this number is always lower or equal unity. Notice that in the colour code of the figure, white colour corresponds to one iteration and black colour to zero iterations. The green colour (two iterations) and red colour (three iterations) do not appear in the figure. The averaged value obtained in our calculations is 0.9843 iterations. The reliability of this code is remarkable because after many millions of runs, the algorithm never enters in dead loops or other kind of non-convergent phases.

Figure 2.

Results with the sdg code (12 polynomials) using the modified Newton–Raphson (MNR) iteration scheme.

Two different calculations have been carried out: (1) using standard double precision and (2) using quadruple precision. In both cases the results can be summarized in Fig.2 because from the point of view of the number of iterations there are no differences in both calculations.

6.2 Accuracy analysis

In this section we focus the analysis on the sdg code with the iteration scheme provided by the MNR algorithm. In this code the iteration ends when the residual ρ = |y −esiny −x| is lower that the zero of the machine ε_tol = 2.22 × 10⁻¹⁶.

Let us consider the true solution y _v (x, e) for given values of e and the mean anomaly M =x. Let y _c (x, e) the numerical solution provided by the sdg code for these particular values. Obviously, due to truncation and round-off errors it is possible to write

\begin{equation*} y_c = y_v + \varepsilon _{{\rm abs}}, \end{equation*}

where ε_abs is the total absolute error associated with the numerical solution y _c (x, e). If the numerical solution y _c (x, e) is a good approximation of the true solution, the value of ε_abs is an infinitesimal: ε_abs ≪ 1. In such a case, this absolute error is closely related with the residual ρ. In effect, for the numerical solution y _c (x, e) the residual is given by

\begin{equation*} \rho = | y_v + \varepsilon _{{\rm abs}} - e\, \sin (y_v + \varepsilon _{{\rm abs}}) -x |. \end{equation*}

Taking into account that the true solution verifies

\begin{equation*} y_v - e\, \sin y_v - x = 0, \end{equation*}

the residual takes the form

\begin{equation*} \rho = | \varepsilon _{{\rm abs}} | |1 - e\,\cos y_v | + o(\varepsilon _{{\rm abs}}^2). \end{equation*}

As a consequence the total absolute error is given by

\begin{equation} | \varepsilon _{{\rm abs}} | = \frac{\rho }{|1 - e\, \cos y_v |}. \end{equation}

(20)

Basically, it depends on the residual ρ, just like the total relative error

\begin{equation} |\varepsilon _{{\rm rel}}| = \frac{|y_c - y_v|}{|y_v|} = \frac{\rho }{|y_v\, (1 - e\,\cos y_v) |}. \end{equation}

(21)

Let us fix the value of the eccentricity e; then we scan the whole interval M ∈ [0, π] calculating the residual ρ after zero iteration (the residual provided by the starting seed), after one iteration, two iterations and so on. Let ρ_max the maximum residual that we found in the scanning of the whole interval M ∈ [0, π] for each iteration number. These maximum residuals permit to calculate the maximum values of the absolute error ε_max for different number of iterations, and these absolute errors can be associated with the value of the eccentricity e.

If we plot the values of such maximum absolute errors versus the eccentricity e we have a very good idea of the accuracy involved in our calculations after zero, one, two, … iterations. Notice that these values mark the lowest accuracy of the numerical solution y _c (x, e) provided by the sdg code; in fact, for each value of e there are hundreds of values of M where the accuracy is much better than the indicated by ε_max.

Fig.3 shows the decimal logarithm of the maximum absolute error versus the eccentricity e. The initial seed – 0 iterations – provides an error which, in the worst cases, ranges from 10⁻⁴ to 10⁻¹². Since the calculations associated with Fig.3 have been carried out in double precision, the residual obtained after one iteration practically saturates the capacity of the machine providing the solution with 15 or 16 significant figures. It should be underlined that one additional iteration does not improve the accuracy of the numerical solution; figure shows that the lines corresponding to one and two iterations are, practically, the same.

Figure 3.

The maximum total absolute error ε_max versus the eccentricity e, for different number of iterations. Calculations in double precision and using 12 polynomials.

Fig.4 is the same as Fig.3 but with the calculations carried out in quadruple precision. It is clear that with one iteration we obtain more than 20 significant figures except for values of e close to unity, where the number of significant digits is about 17–20. With two iterations the sdg code saturates the precision capacities of the machine, providing 34 significant figures.

Figure 4.

The maximum total absolute error ε_max versus the eccentricity e, for different number of iterations. Calculations in quadruple precision and using 12 polynomials.

6.3 Improving the code

Note that the number of polynomials used for the discretization of the interval [0, π] of the mean anomaly M =x does not affect to the speed of the algorithm but it has an important impact on the accuracy. Instead of 12 we modified the code to use more polynomials in order to obtain a better initial seed. However, the number of polynomial cannot be increased arbitrarily because the intervals [x _i , x _i + 1] becomes too small and the polynomials becomes numerically unstable. Thus, when we use 24 polynomials the two first polynomials p ₁(x) and p ₂(x) associated with the two first intervals [x ₁, x ₂] and [x ₂, x ₃] become numerically unstable for some values of the eccentricity e. If these two polynomials are replaced by the first polynomial of the 12 family, the sdg code runs without problems and the accuracy of the algorithm is improved in a significant way.

Fig.5 is the same as Fig.2 but now the calculations have been carried out using 23 polynomials. The figure does not change if quadruple precision is used instead of double precision. The speed of the code is not affected for the increase of polynomials. Obviously, by using double precision the code is faster than using quadruple precision.

Figure 5.

Results with the sdg code (23 polynomials) using the modified Newton–Raphson (MNR) iteration scheme.

Fig.6 is the same as Fig.3 but now 23 polynomials have been used to provide the initial seed. The initial seed – 0 iterations – provides an error which, in the worst cases, ranges from 10⁻⁷ to 10⁻¹⁵. Since the calculations associated with Fig.6 have been carried out in double precision, the residual obtained after one iteration practically saturates the capacity of the machine providing the solution with 15 or 16 significant figures. It should be underlined that one additional iteration does not improve the accuracy of the numerical solution; figure shows that the lines corresponding to one and two iterations are, practically, the same. By comparing Figs3 and6 a clear increase in the accuracy of the initial seed is detected, when the 23 polynomials family is used. However, there are no differences in the accuracy of the results after one iteration.

Figure 6.

The maximum total absolute error ε_max versus the eccentricity e, for different number of iterations. Calculations in double precision and using 23 polynomials.

Fig.7 is the same as Fig.4 but now 23 polynomials have been used. Since the calculations associated with Fig.7 have been carried out in quadruple precision, the residual obtained after one iteration practically saturates the capacity of the machine providing the solution with 34 significant figures except in the region 0.75 <e < 1 where the number of significant digits ranges from 27 to 34. One additional iteration improves the accuracy of the numerical solution for all values of e.

Figure 7.

The maximum total absolute error ε_max versus the eccentricity e, for different number of iterations. Calculations in quadruple precision and using 23 polynomials.

7 COMPARISON OF RESULTS

In our calculations, instead of 12 approximating polynomials we use 23. In effect, after check the behaviour of the algorithm by using several number of polynomials we came to the conclusion that (1) the differences are no important when using 12, 16, 18 or 23 polynomials, and (2) by using 23 the algorithm behaviour is slightly smoother from a global point of view. All the calculations have been carried out in a workstation with Intel(R) Core(TM) i7 CPU 860 @ 2.80 GHz, 2.79 GHz microprocessor in a Windows 7 64 bits operative system and with the same Intel c/c++ compiler.

The Kepler equation for the elliptical case has been solved with the method proposed in this paper by using different Newton–Raphson algorithms. Moreover, the same analysis has been performed with typical solvers of the KE:

• the MNR;

• the Conway method;

• the CNR;

• the method described in Gooding & Odell (1988);

• the method described in Nijenhuis (1991);

• the method described in Fukushima (1996);

• the method described in Mortari & Elipe (2014).

In each case, the KE has been solved as described in Section6.1.

Fig.5 shows, in the plane (e, M), the number of iterations associated with the sdg code when using the MNR iteration scheme. Notice that in the 98.737 per cent of cases we reach the solution with only one iteration (white colour) ¹ and in the 1.263 per cent no iteration is necessary (see Table4). These results justify the option of work with quadruple precision which does not slow the calculations due to the very small number of iterations required. The CPU time invested in the ≈4 × 10⁶ of times that we solved the KE was 41.012 s and the average number of iterations was 0.987. In the figure the zero iteration cases are localized in lines that, approximately, mark the extremes of the intervals associated with the approximating polynomials.

Table4 shows as well the number of iterations associated with the sdg code when using the Conway and CNR methods. Notice that in small number of cases located close to the singular corner – the 0.106 and 0.138 per cent, respectively – two iterations are required. The CPU time invested was 40.779 s with Conway and 40.591 s with CNR, while the average number of iterations was 0.988 and 0.989, respectively. Moreover, Table4 shows the number of iterations associated with the approaches described in Gooding & Odell (1988), Nijenhuis (1991), Fukushima (1996) and Mortari & Elipe (2014) when using double precision. Notice that all of them require more than three iterations for some particular cases, reaching a maximum value of 18, 6, 9 and 4 iterations, respectively. The CPU time invested was 41.2 s with Gooding code, 42.136 s with Nijenhuis code, 59.872 s with Fukushima code and 44.173 s with Mortari code. The average number of iterations was 2.205, 4.5, 3.495 and 2.284, respectively.

Table 4.

Number of iteration (percentage) obtained with the different Newton–Raphson algorithms checked in double precision.

Scheme	i = 0	i = 1	i = 2	i = 3	i ≥ 4
MNR	1.263	98.737	0	0	0
Conway	1.263	98.631	0.106	0	0
CNR	1.263	98.599	0.138	0	0
Gooding	0.499	0	98.144	0	1.357
Nijenhuis	0	0	0	50.024	49.976
Fukushima	0.099	0	1	52.605	46.296
Mortari	0.099	0.050	71.162	28.686	0.003

Scheme	i = 0	i = 1	i = 2	i = 3	i ≥ 4
MNR	1.263	98.737	0	0	0
Conway	1.263	98.631	0.106	0	0
CNR	1.263	98.599	0.138	0	0
Gooding	0.499	0	98.144	0	1.357
Nijenhuis	0	0	0	50.024	49.976
Fukushima	0.099	0	1	52.605	46.296
Mortari	0.099	0.050	71.162	28.686	0.003

Table 4.

Number of iteration (percentage) obtained with the different Newton–Raphson algorithms checked in double precision.

Scheme	i = 0	i = 1	i = 2	i = 3	i ≥ 4
MNR	1.263	98.737	0	0	0
Conway	1.263	98.631	0.106	0	0
CNR	1.263	98.599	0.138	0	0
Gooding	0.499	0	98.144	0	1.357
Nijenhuis	0	0	0	50.024	49.976
Fukushima	0.099	0	1	52.605	46.296
Mortari	0.099	0.050	71.162	28.686	0.003

Scheme	i = 0	i = 1	i = 2	i = 3	i ≥ 4
MNR	1.263	98.737	0	0	0
Conway	1.263	98.631	0.106	0	0
CNR	1.263	98.599	0.138	0	0
Gooding	0.499	0	98.144	0	1.357
Nijenhuis	0	0	0	50.024	49.976
Fukushima	0.099	0	1	52.605	46.296
Mortari	0.099	0.050	71.162	28.686	0.003

Table5 is the same as Table4 but now in quadruple precision. Notice that in this case the sdg code keeps again the smallest number of iterations. The methods described in Gooding & Odell (1988), Nijenhuis (1991) and Fukushima (1996) do not experiment any change, keeping the same percentages as in double precision. The maximum number of iterations reached was the same as in double precision excepting Gooding code with 25 iterations and Mortari code with 17 iterations. The CPU time invested was 83.190 s (MNR), 87.840 s (Conway) and 81.840 s (CNR) with sdg code, 64.335 s with Gooding code, 63.071 s with Nijenhuis code, 260.926 s with Fukushima code and 129.824 s with Mortari code. When quadruple precision is applied, the sdg code does not become the fastest algorithm, but provides the most precise solution with more than 30 significant digits in practically the whole plane (e, M). Indeed, the Figs8 and9 show the number of significant digits obtained with the different codes when using the MNR method in double and quadruple precision. As it is shown, only the sdg code experiments and improvement of 20 significant figures for almost the complete plane when changing from double to quadruple precision. However, the results obtained with the remaining codes highlight the fact that they were developed to work in double precision but not in quadruple precision, especially in the case of Gooding, Nijenhuis and Fukushima codes. As a consequence, the results justify again the option to work with quadruple precision when using the sdg code.

Figure 8.

Number of significant figures, approximately, obtained with the code of Gooding & Odell (1988) (upper), Nijenhuis (1991) (medium) and Fukushima (1996) (lower). On the left, double precision; on the right, quadruple precision.

Figure 9.

Number of significant figures, approximately, obtained with the Mortari & Elipe (2014) (upper) and sdg code (lower). On the left, double precision; on the right, quadruple precision.

Table 5.

Number of iteration (percentage) obtained with the different Newton–Raphson algorithms checked in quadruple precision.

Scheme	i = 0	i = 1	i = 2	i = 3	i ≥ 4
MNR	0.161	99.839	0	0	0
Conway	0.161	34.865	64.974	0	0
CNR	0.161	35.588	66.251	0	0
Gooding	0.499	0	98.144	0	1.357
Nijenhuis	0	0	0	50.024	49.976
Fukushima	0.099	0	1	52.605	46.296
Mortari	0.099	0.049	0.00015	87.036	12.815

Scheme	i = 0	i = 1	i = 2	i = 3	i ≥ 4
MNR	0.161	99.839	0	0	0
Conway	0.161	34.865	64.974	0	0
CNR	0.161	35.588	66.251	0	0
Gooding	0.499	0	98.144	0	1.357
Nijenhuis	0	0	0	50.024	49.976
Fukushima	0.099	0	1	52.605	46.296
Mortari	0.099	0.049	0.00015	87.036	12.815

Table 5.

Number of iteration (percentage) obtained with the different Newton–Raphson algorithms checked in quadruple precision.

Scheme	i = 0	i = 1	i = 2	i = 3	i ≥ 4
MNR	0.161	99.839	0	0	0
Conway	0.161	34.865	64.974	0	0
CNR	0.161	35.588	66.251	0	0
Gooding	0.499	0	98.144	0	1.357
Nijenhuis	0	0	0	50.024	49.976
Fukushima	0.099	0	1	52.605	46.296
Mortari	0.099	0.049	0.00015	87.036	12.815

Scheme	i = 0	i = 1	i = 2	i = 3	i ≥ 4
MNR	0.161	99.839	0	0	0
Conway	0.161	34.865	64.974	0	0
CNR	0.161	35.588	66.251	0	0
Gooding	0.499	0	98.144	0	1.357
Nijenhuis	0	0	0	50.024	49.976
Fukushima	0.099	0	1	52.605	46.296
Mortari	0.099	0.049	0.00015	87.036	12.815

8 APPLICATIONS OF THE SDG CODE

In the introduction of this paper (Wisdom & Holman1991) we can read: long-term integrations are playing an increasingly important role in investigations in dynamical astronomy. The reason is twofold. First, numerical exploration is an essential tool in the study of complex dynamical systems which can exhibit chaotic behaviour, and there has been a growing realization of the importance of chaotic behaviour in dynamical astronomy.

In this sense, the code developed in this paper is very useful in many situations, especially in the numerical integration of perturbed problems which appears in dynamical astronomy. In particular, when the integration requires to solve the KE at every step, that means, a large number of times. Not only in Keplerian perturbed problems, but in more involved situations like the n-body problem.

We must underline that this method works regardless of whether or not the eccentricity maintains a constant value. In some cases, the eccentricity is constant, like, for example in the traditional Encke's method (Encke & Schäfli1854) where the reference orbit is a Keplerian one. The Encke's method is used in some long-term integration, like the LONGSTOP project (Milani et al.1986).

In other cases, the value of the eccentricity changes with time. This situation appears, for example, when the integration is performed with the Deprit's method, based on a Hansen ideal frame (Deprit1975,1976). This integration scheme applies for elliptical orbits and requires to solve the KE in each integration step.

The KE must be solved as well in symplectic integration codes as for example in the Wisdom–Holman method (Wisdom & Holman1991), where the orbital elements experiment small changes at every step.

To avoid the propagation of any kind of error involved in these numerical procedures it makes necessary to solve the KE with the highest possible accuracy. As we have shown, our algorithm, which can be used in double or quadruple precision, satisfies this condition assessing the convergence independent of the numerical method selected.

Finally, we want to highlight that in strong perturbed problems is quite easy to enter in the singular corner. Thanks to the asymptotic expansion, our code is very reliable and accurate and it avoids some chaotic behaviour of numerical nature that appears in that region with other solvers.

9 CONCLUSIONS

Several conclusions can be drawn from the results obtained of our numerical research.

• When the starting seed of the Newton–Raphson algorithm is very good, convergence is always assured as confirmed the thoroughly analysis carried out in our computers.

• First of all we should underline the stability and reliability of our scheme combined with the Newton–Raphson algorithm in its different versions. Notice that the CNR provides the solution in the 98.599 per cent (double precision) of cases with only one iteration. This is remarkable result that we do not found in other algorithms used in the Astrodynamics community. Improving the algorithm by using the Conway or the MNR permits to obtain the right solution with only one iteration in the 98.631 and 98.737 per cent of cases, respectively.

• The low number of iterations permits to use quadruple precision if you like, because the speed of the calculations is not jeopardized. Therefore, we obtain the benefit of a greater accuracy with a minimal cost.

• The global algorithm solves successfully the solution of the KE in the singular corner, M ≪ 1 and e ≈ 1. The asymptotic expansions used to generate the initial seed assure the reliability and convergence of the Newton–Raphson iterative scheme. This is another remarkable results since the convergence in this special region not always is assured by some of the algorithms usually considered in the literature.

• When quadruple precision is used in numerical simulations, the absolute error of the numerical solution provided by the sdg code is clearly under the 10⁻²⁵ after one iteration. In any case, after two iterations the error is always lower than 10⁻³⁰, including the region e > 0.75 (the singular corner).

This work was carried out in the framework of the research project entitled Dynamic Analysis, Advanced Orbital Propagation and Simulation of Complex Space Systems (ESP2013-41634-P) supported by the State Secretariat for Research, Development and Innovation of the Spanish Ministry of Economy and Competitiveness. Authors thank to Professor T. Fukushima for his support codifying the algorithm of his paper (Fukushima1996).

†

Present address: E.T.S.I. Aeronáutica y del Espacio, Pza. Cardenal Cisneros 3, E-28040 Madrid, Spain.

‡

Present address: E.T.S.I. Aeronáutica y del Espacio, Pza. Cardenal Cisneros 3, E-28040 Madrid, Spain.

1

The background white colour in the figures corresponds to cases in which only one iteration is needed to reach the solution with the prescribed tolerance.

REFERENCES

Battin

R.

,

1999

,

An Introduction to the Mathematics and Methods of Astrodynamics

, revised edn,

AIAA

,

Reston, VA

Colwell

P.

,

1993

,

Solving Kepler's Equation over Three Centuries

,

Willmann-Bell

,

Richmond, VA

Conway

B. A.

,

1986

,

Celest. Mech.

,

39

,

199

Danby

J.

,

1992

,

Fundamentals of Celestial Mechanics

,

Willman-Bell

,

Richmond, VA

Deprit

A.

,

1975

,

J. Res. Natl. Bureau Standards, Sec. B: Math. Sci.

,

79B

,

1

Deprit

A.

,

1976

,

Celest. Mech. Dyn. Astron.

,

13

,

253

Encke

J. F.

,

Schäfli

L.

,

1854

,

Über die Allgemeinen Störungen der Planeten. Berliner Astronomisches Jahrbuch fr 1857

,

Königlich Akademie der Wissenschaften

,

Berlin

Fukushima

T.

,

1996

,

Celest. Mech. Dyn. Astron.

,

66

,

309

Gooding

R.

,

Odell

A.

,

1988

,

Celest. Mech.

,

44

,

267

Markley

F. L.

,

1995

,

Celest. Mech. Dyn. Astron.

,

63

,

101

Milani

A.

,

Nobili

A.

,

Fox

K.

,

Carpino

M.

,

1986

,

Nature

,

319

,

386

Mortari

D.

,

Elipe

A.

,

2014

,

Celest. Mech. Dyn. Astron.

,

118

,

1

Nijenhuis

A.

,

1991

,

Celest. Mech. Dyn. Astron.

,

51

,

319

Pál

A.

,

2009

,

MNRAS

,

396

,

1737

Rambaut

A.

,

1906

,

MNRAS

,

66

,

520

Serafin

R.

,

1986

,

Celest. Mech.

,

38

,

111

Siewert

C.

,

Burniston

E.

,

1972

,

Celest. Mech.

,

6

,

294

Stumpf

L.

,

1999

,

Celest. Mech. Dyn. Astron.

,

74

,

95

Wisdom

J.

,

Holman

M.

,

1991

,

AJ

,

102

,

1528

APPENDIX A: The fifth degree polynomial |$\boldsymbol {p_1(x,e)}$|

The polynomial p ₁(x, e) is given by the expression

\begin{equation*} p_1(x,e)=a^{(1)}_{0}(e)+a^{(1)}_{1}(e)x+a^{(1)}_{2}(e)x^{2}+a^{(1)}_{3}(e)x^{3}+a^{(1)}_{4}(e)x^{4}+a^{(1)}_{5}(e)x^{5}, \end{equation*}

where the coefficients are defined as a function of the eccentricity by

\begin{equation*} a^{(1)}_{5}(e)=\frac{432.0\,e \left( 0.5393817\,{e}^{3}- 0.8743466\,{e}^{2 }+ 0.2650564\,e+ 0.0453931 \right) }{ d^{(1)}_{5}(e) }, \end{equation*}

\begin{equation*} a^{(1)}_{4}(e)=\frac{72.0\,e \left( 0.4657683\,{e}^{3}- 0.402960\,{e}^{2}- 0.125744\,e- 0.002438 \right) }{ d^{(1)}_{4}(e) }, \end{equation*}

\begin{equation*} a^{(1)}_{3}(e)=\frac{ 12.0\, \left( 2.1335977\,{e}^{3}- 7.174617\,{e}^{2}+ 8.438043\,e- 3.4460552 \right) e}{d^{(1)}_{3}(e)}, \end{equation*}

\begin{equation*} a^{(1)}_{2}(e)=0, \end{equation*}

\begin{equation*} a^{(1)}_{1}(e)=\frac{-1.0}{- 1.0+e}, \end{equation*}

\begin{equation*} a^{(1)}_{0}(e)=0, \end{equation*}

with

\begin{eqnarray*} d^{(1)}_{5}(e) &=&\left( 3739.002236\,{e}^{3}-2549.554211\,{e}^{4}+ 926.7532457\,{e}^{5}\right.\\ &&- 140.2961154\,{e}^{6} \left.- 3082.919872\,{e}^{2}+ 1355.064876\,e- 248.0502135 \right)\\ && \left( -1.0+e \right)\left( 9.0\,{e}^{2}- 18.84955592\,e+ 9.869604404 \right), \end{eqnarray*}

\begin{eqnarray*} d^{(1)}_{4}(e) &=&\left(3.141592654- 3.0\,e \right) \left( 231.3573259\,e- 323.4979098\,{e} ^{2}\right. \\ &&\left. +\, 200.9177486\,{e}^{3}- 46.76537182\,{e}^{4}- 62.01255338 \right)\\ &&\left( - 1.0+e \right) \left( 6.928203232\,e- 3.0\,{e}^{2}- 4.0 \right), \end{eqnarray*}

\begin{eqnarray*} d^{(1)}_{3}(e)&=&\left( - 27.0\,{e}^{3}+ 84.82300166\,{e}^{2}- 88.82643964\,e+ 31.00627669 \right)\\ &&\times \left( 5.196152424\,{e}^{3}- 18.0\,{e}^{2}+ 20.78460970\,e- 8.0\right) \left( - 1.0+e \right). \end{eqnarray*}

© 2017 The Authors Published by Oxford University Press on behalf of the Royal Astronomical Society

© 2017 The Authors Published by Oxford University Press on behalf of the Royal Astronomical Society

ashcroftanderser.blogspot.com

Source: https://academic.oup.com/mnras/article/467/2/1702/2929272

Share this post