Mohamed M. Hafez
Department of Mechanical and Aerospace Engineering, University of California, Davis, CA, USA.
DOI: 10.4236/ijaa.2026.164021   PDF    HTML   XML   2 Downloads   27 Views   Citations

Abstract

This study compares numerical methods for the elliptic, parabolic, and hyperbolic forms of Kepler’s equation. Newton-type iterations are combined with initial approximations constructed by bisection and local inverse interpolation using linear, quadratic, ordinary cubic, and piecewise cubic Hermite formulas. The central contribution is a systematic construction of stretched grids guided by the singular behaviour of the inverse mapping as $e\to 1^{-}$ and $e\to 1^{+}$ . Regular, geometrically stretched, and analytically optimised grids are compared for elliptic and hyperbolic cases. Continuation, Wynn’s epsilon algorithm, and homotopy are also examined as alternative means of obtaining initial approximations. Numerical results show that concentrating nodes near the parabolic limit substantially improves local Hermite interpolation and usually leaves only a few final Newton corrections.

Keywords

Kepler’s Equation, Newton-Type Iterative Methods, Inverse Interpolation, Stretched Grids, Singular Asymptotics, Bisection Methods, Numerical Results

Share and Cite:

1. Introduction

Kepler’s equation arises in many applications of science and engineering, particularly in orbital mechanics and celestial dynamics. Although it has a simple scalar form, its numerical solution has attracted sustained attention for several centuries. Colwell [1] reviewed its historical development and cited more than four hundred related publications up to 1993. The continuing appearance of new solution procedures demonstrates that Kepler’s equation remains an instructive and practically important nonlinear problem. General treatments of the orbital and celestial mechanics setting, including discussions of Kepler’s equation and its applications, are available in the books by Meeus [2], Danby [3], Fitzpatrick [4], Hintz [5], and Giorgilli et al. [6].

The present work was motivated in part by the recent study of Bekir [7], in which useful linear and quadratic empirical formulae were used to construct initial approximations for Newton’s method in the elliptic case. The resulting approximations allowed the solution to be obtained in no more than a few Newton iterations over the parameter ranges considered. On the other hand, the purpose of the present work is to develop a systematic approach based on standard numerical procedures. This viewpoint leads to a broader investigation of interpolation-based initial approximations and their behaviour near the difficult limiting regimes.

Kepler’s equation is a scalar nonlinear equation whose solution depends on the mean anomaly and the eccentricity. Its apparent simplicity makes it an effective example through which several standard ideas in numerical analysis can be introduced and compared. These include bracketing by bisection; Newton, Halley, and Chebyshev iterations; inverse interpolation; convergence order; and the treatment of nearly singular behaviour. The relationships among higher-order iterations for a single nonlinear equation are discussed in the classical work of Householder [8].

In this paper, linear, quadratic, ordinary cubic, and piecewise cubic Hermite inverse approximations are examined as systematic procedures for constructing initial guesses. Uniform grids are first considered, followed by geometrically stretched and analytically optimised grids. Inverse interpolation and spline-based inversion of Kepler’s equation have previously been studied, notably by Tommasini and Olivieri [9]. The distinguishing feature of the present approach is the systematic construction of non-uniform stretched grids guided by the asymptotic behaviour of the solution near the parabolic limit. As $e\to 1^{-}$ in the elliptic case and $e\to 1^{+}$ in the hyperbolic case, the inverse mapping develops increasingly rapid variation near the singular region. A uniform grid does not distribute its nodes efficiently in this limit. We therefore concentrate the interpolation nodes in the region where they are most needed, using stretching strategies derived from the corresponding singular behaviour. This provides a unified and physically motivated treatment of the elliptic and hyperbolic regimes while maintaining the parabolic case as their common limiting configuration.

This asymptotically guided grid construction is not merely an implementation detail; it is a central element of the present strategy. The inverse interpolation framework supplies the local approximation, whereas the stretched grid determines how the available information is distributed so that accuracy is maintained as the parabolic limit is approached. Both elliptic and hyperbolic forms of Kepler’s equation are studied, and the accuracy and efficiency of the resulting algorithms are evaluated numerically.

Several alternative approaches are also considered. Continuation methods are introduced in the spirit of solution-path and fixed-point techniques [10]. Wynn’s epsilon algorithm is used as a nonlinear sequence transformation, following the general framework discussed by Weniger [11]. Finally, algebraic and derivative homotopy formulations are examined in connection with the homotopy solution of the elliptic and hyperbolic Kepler equations proposed by Fitz-Coy and Jang [12]. These alternatives illustrate different ways of transforming, approaching, or accelerating the solution of a nonlinear equation.

Here, Kepler’s equation is used as a unified example to study how analysis, approximation, iteration, and numerical experimentation can be combined in the solution of a nonlinear equation. More advanced approaches requiring complex-variable analysis are beyond the scope of the present work.

2. Kepler’s Equation and Its Three Regimes

Consider the two-body problem of orbital mechanics, in which the relative motion of two bodies—such as a planet and a star, or an artificial satellite and Earth—is described by a conic section. According to the eccentricity $e$ , the orbit is elliptic for $0\le e<1$ , parabolic for $e=1$ , and hyperbolic for $e>1$ , as illustrated in Figure 1. Kepler’s equation relates time, represented by a mean anomaly, to an angular or auxiliary anomaly that locates the moving body on its orbit.

Figure 1. The elliptic, parabolic, and hyperbolic orbital regimes, classified by eccentricity.

Geometrically, the parabolic orbit is the continuous boundary between the two other regimes: with the periapsis distance held fixed, the elliptic orbit approaches the parabola as $e\to 1^{-}$ , while the hyperbolic orbit approaches the same parabola as $e\to 1^{+}$ . This two-sided limiting behaviour is shown in Figure 2. The corresponding time-of-flight equations require a suitable limiting change of variables; direct substitution of $e=1$ into the elliptic or hyperbolic Kepler equation does not by itself produce Barker’s equation.

2.1. Elliptic Case

For an elliptic orbit, Kepler’s equation is

$M=E-e\sin E,0\le e<1,$ (1)

where $M$ is the mean anomaly and $E$ is the eccentric anomaly. The geometric definition of $E$ is shown in Figure 3. The ellipse is stretched vertically into its circumscribing auxiliary circle. The circular sector has area $a^{2}E/2$ , while the focus triangle has area $a^{2}e\sin E/2$ . Their difference is $a^2\left(E-e\sin E\right)/2$ , which leads directly to Equation (1) [13].

Figure 2. Two-sided convergence to the parabolic orbit for a fixed periapsis distance $q$ . The curves $e=0.99$ and $e=1.01$ nearly coincide with the parabola near periapsis.

Figure 3. Geometric construction of the elliptic eccentric anomaly $E$ and true anomaly $\nu$ . Original Figure based on the geometric definitions summarised by [13].

For fixed $e$ , the function $M\left(E\right)$ is monotone because $\text{d}M/\text{d}E=1-e\cos E>0$ . Consequently, a prescribed $M$ determines a unique $E$ .

2.2. Parabolic Case

At $e=1$ , the elliptic eccentric anomaly becomes degenerate, and the time-of-flight relation is expressed conveniently by Barker’s equation [14]-[16]. The precise limiting relationship between the elliptic and hyperbolic Kepler equations and Barker’s equation is derived in Appendix. Let $\nu$ denote the true anomaly, measured at the focus from the direction of periapsis, and define the parabolic anomaly

$D=\text{tan}\frac{\nu }{2}.$ (2)

Barker’s equation may then be written as

$t-\tau =\frac{1}{2}\sqrt{\frac{p^3}{\mu }}\left(D+\frac{D^3}{3}\right),$ (3)

where $t$ is the current time, $\tau$ is the time of periapsis passage, $p$ is the semilatus rectum, and $\mu =G\left(m_1+m_2\right)$ is the gravitational parameter. Introducing $M_p=2\sqrt{\mu /p^3}\left(t-\tau \right)$ gives

$M_p=D+\frac{D^3}{3}.$ (4)

Figure 4 shows the corresponding half-angle construction [17]. Since $\text{d}{M}_p/\text{d}D=1+D^2>0$ , every prescribed $M_p$ determines a unique $D$ , and $\nu =2\arctan D$ .

Figure 4. Geometric definition of the parabolic anomaly $D=\text{tan}\left(\nu /2\right)$ . Original Figure based on the parabolic formulation summarised by [17].

2.3. Hyperbolic Case

For a hyperbolic orbit, the hyperbolic anomaly is denoted by $H$ . The hyperbolic form of Kepler’s equation is

$M=e\sinh H-H,e>1.$ (5)

The anomaly $H$ is a hyperbolic angle, represented by twice the area of an auxiliary hyperbolic sector rather than an ordinary Euclidean angle. In Figure 5, subtracting the auxiliary sector $a^{2}H/2$ from the triangle $a^{2}e\sinh H/2$ gives $a^2\left(e\sinh H-H\right)/2$ and hence Equation (5) [18].

Figure 5. Geometric construction of the hyperbolic anomaly $H$ and true anomaly $\nu$ . Original Figure based on the hyperbolic definitions summarised by [18].

For $H\ge 0$ , $\text{d}M/\text{d}H=e\cosh H-1>0$ , so a prescribed mean anomaly determines a unique hyperbolic anomaly. These constructions provide a unified introduction to the three orbital regimes and reveal the monotonicity and uniqueness of the inverse relation before numerical iteration is introduced.

3. Iterative Methods

Let

$f\left(x\right)=0$ (6)

represent a nonlinear algebraic equation whose solution is sought by an iterative method. We denote an exact root by $x^{\text{*}}$ and the current approximation by $x_n$ .

3.1. Fixed-Parameter Iteration and an Artificial-Time Interpretation

The simplest fixed-parameter iteration may be written as

$x_{n+1}=x_n+\alpha f\left(x_n\right),$ (7)

where $\alpha \ne 0$ is fixed. Equation (7) may also be viewed as a forward Euler discretisation, in an artificial time $\tau$ , of

$\frac{\text{d}x}{\text{d}\tau }=\alpha f\left(x\right),$ (8)

with the artificial time step absorbed into $\alpha$ . Starting from an initial guess $x_0$ , and with a suitable choice of $\alpha$ , the iteration seeks a steady state $x^{\text{*}}$ satisfying $f\left(x^{\text{*}}\right)=0$ .

To see the local convergence condition, let $e_n=x_n-x^{\text{*}}$ . Linearisation of (7) about a simple root gives

$e_{n+1}=\left[1+\alpha f^{\prime }\left(x^{\text{*}}\right)\right]e_n+O\left(e_n^2\right).$ (9)

Thus, local convergence requires

$\left|1+\alpha f^{\prime }\left(x^{\text{*}}\right)\right|<1.$ (10)

When it converges, the fixed-parameter iteration is generally only linearly convergent. Faster convergence can be obtained by allowing the effective parameter to depend on the current iterate.

3.2. Taylor-Series Basis

The following methods can be motivated by expanding $f$ about the current approximation $x_n$ , assumed to be sufficiently close to the simple root $x^{\text{*}}$ . With

$\text{Δ}x=x^{\text{*}}-x_n,$ (11)

we have

$\begin{array}{l}0=f\left(x^{*}\right)=f\left(x_n\right)+f^{\prime }\left(x_n\right)\Delta x+\frac{1}{2}{f}^{″}\left(x_n\right){\left(\Delta x\right)}^2\\ +\frac{1}{6}{f}^{‴}\left(x_n\right){\left(\Delta x\right)}^3+\frac{1}{24}{f}^{\left(4\right)}\left(x_n\right){\left(\Delta x\right)}^4+\cdots .\end{array}$ (12)

In the formulas below, $f$ and all its derivatives are evaluated at $x_n$ unless an argument is displayed explicitly.

3.3. Newton’s Method

Newton’s method retains only the first two terms of Equation (12) and neglects the higher-order terms, assuming that $\text{Δ}x$ is small. Hence

$0\approx f+f^{\prime }\Delta x,$ (13)

which gives the Newton correction

$\Delta x_N=-\frac{f}{f^{\prime }}.$ (14)

The resulting iteration is

$x_{n+1}=x_n-\frac{f\left(x_n\right)}{f^{\prime }\left(x_n\right)}.$ (15)

For a sufficiently smooth function, a simple root, and an initial guess close enough to that root, Newton’s method converges quadratically.

3.4. Chebyshev’s Method

To obtain faster convergence, Chebyshev retains the first three terms of Equation (12). Direct retention would produce a quadratic equation for $\text{Δ}x$ . Chebyshev avoids solving it by approximating ${\left(\Delta x\right)}^2$ with the square of the Newton correction:

${\left(\Delta x\right)}^2\approx {\left(\Delta x_N\right)}^2.$ (16)

Substitution into the truncated Taylor expansion gives

$\Delta x_C=-\frac{f}{f^{\prime }}-\frac{f^{″}{f}^2}{2{\left(f^{\prime }\right)}^3},$ (17)

and therefore

$x_{n+1}=x_n+\Delta x_C.$ (18)

Chebyshev’s method has cubic local convergence.

3.5. Halley’s Method

Halley uses the same three-term Taylor approximation but replaces ${\left(\Delta x\right)}^2$ by the product of the Newton correction and the unknown Halley correction:

${\left(\Delta x\right)}^2\approx \Delta x_N\Delta x_H.$ (19)

The resulting equation is linear in $\text{Δ}{x}_H$ and yields

$\text{Δ}{x}_H=-\frac{2f{f}^{\prime }}{2{\left(f^{\prime }\right)}^2-f{f}^{″}}=-\frac{f/f^{\prime }}{1-f{f}^{″}/\left[2{\left(f^{\prime }\right)}^2\right]}.$ (20)

Thus,

$x_{n+1}=x_n+\text{Δ}{x}_H.$ (21)

Halley’s method also has cubic local convergence. Both Chebyshev’s and Halley’s methods require evaluation of $f^{″}$ , whereas Newton’s method requires only $f^{\prime }$ . Thus, to avoid solving a quadratic equation for $\text{Δ}x$ , Chebyshev and Halley introduced appropriate simplifications without sacrificing cubic convergence. Similarly, Danby avoided solving the cubic polynomial equation for $\text{Δ}x$ by using a sequence of nested corrections.

3.6. Danby’s Fourth-Order Method

Danby and Burkardt [19] retained additional information from the Taylor expansion and constructed nested corrections of progressively higher order. In a sign convention consistent with Equation (14), define

$\Delta x_1=-\frac{f}{f^{\prime }},$ (22)

$\Delta x_2=-\frac{f}{f^{\prime }+\frac{1}{2}\Delta x_1{f}^{″}},$ (23)

$\Delta x_3=-\frac{f}{f^{\prime }+\frac{1}{2}\Delta x_2{f}^{″}+\frac{1}{6}{\left(\Delta x_2\right)}^2{f}^{‴}}.$ (24)

The new approximation is

$x_{n+1}=x_n+\Delta x_3.$ (25)

The corrections $\text{Δ}{x}_1$ , $\text{Δ}{x}_2$ , and $\text{Δ}{x}_3$ generate, respectively, second-, third-, and fourth-order local convergence. The last formula is commonly known as Danby’s method. It requires derivatives through $f^{‴}$ . Danby [3] provides the broader celestial mechanics setting in which this iteration is used.

3.7. Householder’s General Family

Householder developed a systematic family of higher-order iterations using a Padé construction; see his monograph [8] and related paper [20]. The result can be expressed compactly in terms of derivatives of the reciprocal function $1/f$ :

$x_{n+1}=x_n+d \frac{{\left(1/f\right)}^{\left(d-1\right)}\left(x_n\right)}{{\left(1/f\right)}^{\left(d\right)}\left(x_n\right)},d=1,2,3,\cdots ,$ (26)

where the superscript in parentheses denotes differentiation. For a sufficiently smooth function and a simple root, Equation (26) has local convergence of order $d+1$ . In particular,

• $d=1$ gives Newton’s second-order method;

• $d=2$ gives Halley’s third-order method;

• $d=3$ gives a fourth-order method requiring derivatives through $f^{‴}$ ;

• $d=4$ gives a fifth-order method requiring derivatives through $f^{\left(4\right)}$ .

The compact notation in Equation (26) should not obscure its computational cost: increasing the formal convergence order requires additional derivatives and more algebra per iteration. Consequently, a higher-order method is not automatically more efficient. Its advantage depends on the cost and reliability of the derivatives, the quality of the initial approximation, and the number of function evaluations saved.

3.8. Summary of Iterative Methods

This progression gives a useful systematic interpretation of several familiar root-finding procedures. The fixed-parameter iteration is generally linear; Newton’s method uses a variable correction and is quadratic; Chebyshev’s and Halley’s methods use the second derivative and are cubic; Danby’s nested correction uses the third derivative and is fourth order; and Householder’s formula extends this construction to arbitrary order. These orders describe local behaviour near a simple root. They do not guarantee convergence from a distant initial guess or near a multiple root, and the denominators appearing in the formulas must remain safely separated from zero. Several other, more sophisticated methods for solving nonlinear equations are available; however, they will not be discussed here.

4. Initial Guesses

Newton-type methods generally require a sufficiently accurate initial guess to ensure rapid and reliable convergence. Several studies have examined this issue in detail, particularly for Kepler’s equation. Esmaelzadeh and Ghadiri [21] compared several starting approximations for the Newton-Raphson method and evaluated their performance in terms of stability, number of iterations, and computational time. Charles and Tatum [22] examined the convergence of Newton’s iteration for Kepler’s equation, while Mikkola [23] developed a widely used cubic starting approximation. Mathar [24] proposed additional improved first estimates in an arXiv preprint revised in 2025.

The role of the starting approximation becomes especially important near the parabolic limit, where the variation of the inverse mapping becomes strongly non-uniform. A practical strategy is first to identify a small interval that contains the solution and then to use a local approximation from that interval as the initial value for a higher-order iteration.

5. Bisection and Local Approximation Methods

The classical bisection method provides a reliable procedure for locating a solution. Let the root $x^{\text{*}}$ be bracketed by $\left[a_0,b_0\right]$ such that

$f\left(a_0\right)f\left(b_0\right)<0.$ (27)

At iteration $n$ , define

$c_n=\frac{a_n+b_n}{2}.$ (28)

If $f\left(a_n\right)f\left(c_n\right)<0$ , set

$a_{n+1}=a_n,b_{n+1}=c_n;$ (29)

otherwise set

$a_{n+1}=c_n,b_{n+1}=b_n.$ (30)

After $n$ bisections,

$b_n-a_n=\frac{b_0-a_0}{2^n},\left|x^{\text{*}}-c_n\right|\le \frac{b_0-a_0}{2^{n+1}}.$ (31)

Bisection is globally reliable once a valid bracket has been established, although its convergence is only linear. It can therefore locate a small region containing the root before a higher-order iteration or interpolation formula is applied.

Approximation and interpolation methods provide another approach. For Kepler’s equation, the function and all required derivatives are analytically available at every grid point. It is therefore unnecessary to construct a highly accurate approximation over the entire domain. Once the interval containing the desired solution is identified, a local approximation can be constructed from neighbouring function and derivative values.

6. Inverse Interpolation

Inverse interpolation is particularly attractive for Kepler’s equation. Instead of approximating $M=M\left(E\right)$ and then solving the approximation for $E$ , we interchange the variables and directly approximate

$E=E\left(M\right).$ (32)

This switching of variables has been used in spline-based inversion methods; a particularly relevant treatment is given by Tommasini and Olivieri [9]. Closely related interpolation-based solvers were developed by Raposo-Pulido and Peláez for the elliptic [25] and hyperbolic [26] Kepler equations. Their work emphasises accurate local starting approximations, high-order endpoint interpolation, and special treatment of the nearly singular region.

For the elliptic equation,

$M\left(E\right)=E-e\sin E,0\le e<1,$ (33)

the mapping is monotone because

$\frac{\text{d}M}{\text{d}E}=1-e\cos E>0.$ (34)

For neighbouring points $E_i$ and $E_{i+1}$ , define

$M_i=E_i-e\sin E_i,M_{i+1}=E_{i+1}-e\sin E_{i+1}.$ (35)

Suppose that the interval containing the required solution has been located. The interpolation stencil is then formed from nearby tabulated pairs $\left(M_j,E_j\right)$ . At an endpoint of the tabulated range, the stencil is shifted so that only available nodes are used. The following formulas provide a natural progression in approximation order.

6.1. Linear Inverse Interpolation

The simplest inverse approximation uses two neighbouring nodes:

$E^{\left[1\right]}\left(M\right)=E_i+\frac{E_{i+1}-E_i}{M_{i+1}-M_i}\left(M-M_i\right).$ (36)

For a sufficiently smooth inverse function, its error is

$E\left(M\right)-E^{\left[1\right]}\left(M\right)=\frac{E^{″}\left(\xi \right)}{2!}\left(M-M_i\right)\left(M-M_{i+1}\right),$ (37)

where $\xi$ lies in the smallest interval containing the interpolation nodes and the target value $M$ .

6.2. Quadratic Inverse Interpolation

Using three local nodes $M_{i-1}$ , $M_i$ , and $M_{i+1}$ gives the quadratic Lagrange approximation

$\begin{array}{c}{E}^{\left[2\right]}\left(M\right)=E_{i-1}\frac{\left(M-M_i\right)\left(M-M_{i+1}\right)}{\left(M_{i-1}-M_i\right)\left(M_{i-1}-M_{i+1}\right)}\\ +E_i\frac{\left(M-M_{i-1}\right)\left(M-M_{i+1}\right)}{\left(M_i-M_{i-1}\right)\left(M_i-M_{i+1}\right)}\\ +E_{i+1}\frac{\left(M-M_{i-1}\right)\left(M-M_i\right)}{\left(M_{i+1}-M_{i-1}\right)\left(M_{i+1}-M_i\right)}.\end{array}$ (38)

The corresponding error has the form

$E\left(M\right)-E^{\left[2\right]}\left(M\right)=\frac{E^{\left(3\right)}\left(\xi \right)}{3!}\left(M-M_{i-1}\right)\left(M-M_i\right)\left(M-M_{i+1}\right).$ (39)

6.3. Cubic Inverse Interpolation

Ordinary cubic interpolation uses four function values. With the local stencil $I_3=\left\{i-1,i,i+1,i+2\right\}$ , define

$L_j\left(M\right)={\displaystyle \prod }_{\begin{array}{c}k\in I_3\\ k\ne j\end{array}}\frac{M-M_k}{M_j-M_k}.$ (40)

The cubic Lagrange approximation is then

$E^{\left[3\right]}\left(M\right)={\displaystyle \sum }_{j\in I_3}{E}_j{L}_j\left(M\right),$ (41)

with error

$E\left(M\right)-E^{\left[3\right]}\left(M\right)=\frac{E^{\left(4\right)}\left(\xi \right)}{4!}{\displaystyle \prod }_{j\in I_3}\left(M-M_j\right).$ (42)

This formula remains completely local and requires no derivatives, but it uses data from four distinct grid points.

6.4. Piecewise Cubic Hermite Inverse Interpolation

An alternative cubic approximation uses the two endpoint values and the two endpoint derivatives on a single interval. Define

$s=\frac{M-M_i}{M_{i+1}-M_i},\Delta M_i=M_{i+1}-M_i,0\le s\le 1.$ (43)

The inverse derivatives are known analytically:

$E{\text{'}}_i=\frac{1}{1-e\cos E_i},{E^{\prime }}_{i+1}=\frac{1}{1-e\cos E_{i+1}}.$ (44)

The cubic Hermite approximation on the interval is

$\begin{array}{c}{E}^{\left[H\right]}\left(M\right)=\left(2s^3-3s^2+1\right)E_i+\left(s^3-2s^2+s\right)\Delta M_i{E^{\prime }}_i\\ +\left(-2s^3+3s^2\right)E_{i+1}+\left(s^3-s^2\right)\Delta M_i{E^{\prime }}_{i+1}.\end{array}$ (45)

Its error is

$E\left(M\right)-E^{\left[H\right]}\left(M\right)=\frac{E^{\left(4\right)}\left(\xi \right)}{4!}{\left(M-M_i\right)}^2{\left(M-M_{i+1}\right)}^2.$ (46)

Thus ordinary cubic interpolation uses four inverse-function values, whereas piecewise cubic Hermite interpolation uses two values and two analytically available inverse derivatives. Adjacent intervals share the same nodal value and first derivative, so the resulting piecewise interpolant is $C^1$ . Its second derivative, however, is not generally continuous across an interval boundary. Consequently, this construction is not a classical cubic spline, which normally has $C^2$ continuity. The cubic Hermite formula directly produces an approximation to the eccentric anomaly, is local to one interval, and requires neither the solution of an interpolating polynomial nor the global linear system required by a cubic spline. This makes it especially natural for Kepler’s equation. If greater interpolation accuracy is required, the next natural choice is a piecewise quintic Hermite interpolant. In addition to the function values and first derivatives used by the cubic Hermite formula, it matches the second derivatives at the two endpoints. Adjacent intervals then share values, first derivatives, and second derivatives, giving a $C^2$ piecewise interpolant. It is still distinct from a conventional quintic spline, which generally imposes higher global smoothness.

The same construction applies to the hyperbolic Kepler equation,

$M\left(H\right)=e\sinh H-H,e>1,$ (47)

for which

$\frac{\text{d}H}{\text{d}M}=\frac{1}{e\cosh H-1},$ (48)

and to Barker’s parabolic equation,

$M_p=D+\frac{D^3}{3},\frac{\text{d}D}{\text{d}{M}_p}=\frac{1}{1+D^2}.$ (49)

Thus, in all three regimes, the forward function and the derivatives required for inverse interpolation are available analytically at every grid point.

The linear, quadratic, ordinary cubic, and cubic Hermite formulas apply in the same way after replacing $\left(E,M\right)$ by $\left(H,M\right)$ or $\left(D,M_p\right)$ and using the appropriate tabulated values. It is not necessary to construct a high-accuracy global approximation over the complete domain. After the interval containing the solution has been identified, an accurate local inverse approximation can be formed using only neighbouring grid values. The efficiency of this procedure depends strongly on grid placement. Near $e=1$ , the proposed stretched grids concentrate nodes according to the singular behaviour of the inverse solution. This choice also follows directly from the error formulas above: reducing the local spacing reduces the products of the interpolation factors precisely where the higher inverse derivatives become large. The singular point itself is excluded from the interpolation stencil. The stretched grid thereby complements the inverse interpolation formula and is a central distinguishing feature of the present strategy.

7. Asymptotic Behaviour near the Parabolic Limit

The difficulty near $e=1$ can be understood directly from the local asymptotic forms of the elliptic and hyperbolic equations. The analysis also shows that the two regimes possess the same leading singular behaviour. Only real variables are required.

7.1. Elliptic Limit: $e\to 1^{-}$

Write

$e=1-{\delta }_{-},{\delta }_{-}>0,$ (50)

and consider the region of small eccentric anomaly $E$ . Expanding the sine in the elliptic Kepler equation gives

$\begin{array}{c}M=E-e\sin E\\ ={\delta }_{-}E+\frac{1-{\delta }_{-}}{6}{E}^3-\frac{1-{\delta }_{-}}{120}{E}^5+O\left(E^7\right).\end{array}$ (51)

Consequently, the common leading terms are

$M={\delta }_{-}E+\frac{E^3}{6}+O\left({\delta }_{-}{E}^3+E^5\right).$ (52)

At the origin, the inverse slope is

${\frac{\text{d}E}{\text{d}M}|}_{M=0}=\frac{1}{1-e}=\frac{1}{{\delta }_{-}},$ (53)

which becomes unbounded as $e\to 1^{-}$ . Thus the inverse solution changes rapidly near $M=0$ , even though the forward equation remains smooth.

7.2. Hyperbolic Limit: $e\to 1^{+}$

Similarly, write

$e=1+{\delta }_{+},{\delta }_{+}>0.$ (54)

For small hyperbolic anomaly $H$ ,

$\begin{array}{c}M=e\sinh H-H\\ ={\delta }_{+}H+\frac{1+{\delta }_{+}}{6}{H}^3+\frac{1+{\delta }_{+}}{120}{H}^5+O\left(H^7\right),\end{array}$ (55)

and hence

$M={\delta }_{+}H+\frac{H^3}{6}+O\left({\delta }_{+}{H}^3+H^5\right).$ (56)

The corresponding inverse slope at the origin is

${\frac{\text{d}H}{\text{d}M}|}_{M=0}=\frac{1}{e-1}=\frac{1}{{\delta }_{+}},$ (57)

and it becomes unbounded as $e\to 1^{+}$ . The signs of the fifth-order terms in the elliptic and hyperbolic expansions are different, but their linear and cubic leading terms have the same form.

7.3. Common Transition Scaling

Equations (52) and (56) can be written in the common notation

$M=\delta A+\frac{A^3}{6}+O\left(\delta A^3+A^5\right),$ (58)

where $\left(A,\delta \right)=\left(E,{\delta }_{-}\right)$ on the elliptic side and $\left(A,\delta \right)=\left(H,{\delta }_{+}\right)$ on the hyperbolic side. The linear and cubic terms balance when

$A=O\left({\delta }^{1/2}\right),M=O\left({\delta }^{3/2}\right).$ (59)

More precisely, introduce

$A=\sqrt{6\delta }y,\mu =\frac{M}{\sqrt{6}{\delta }^{3/2}}.$ (60)

Then the transition equation becomes

$\mu =y+y^3+O\left(\delta \right).$ (61)

This scaled relation is regular and describes the narrow region in which the linear and cubic contributions are comparable.

E − sin(E) = M(62)

Therefore,

sinh(H) − H = M(63)

For fixed positive M, the limiting elliptic and hyperbolic equations as e→1 are given by (62) and (63), respectively. Their solutions are generally different. The common cube-root approximation A ~ (6M)1/3, where A = E or H, applies in the joint limit M→0 and δ/M2/3→0. Barker’s equation is obtained through the physical-time rescaling described in the Appendix.

7.4. Inverse Derivatives and Interpolation

Let $c=6^{1/3}$ . Differentiating the local cube-root approximation $A\left(M\right)~c{M}^{1/3}$ gives the following leading-order asymptotic derivatives:

$\begin{array}{l}{A}^{\prime }\left(M\right)~\frac{c}{3}{M}^{-2/3},A^{″}\left(M\right)~-\frac{2c}{9}{M}^{-5/3},\\ A^{\left(3\right)}\left(M\right)~\frac{10c}{27}{M}^{-8/3},A^{\left(4\right)}\left(M\right)~-\frac{80c}{81}{M}^{-11/3}.\end{array}$ (64)

These expressions explain why a uniform grid in $M$ becomes inefficient near the parabolic limit. The derivatives entering the interpolation-error formulas grow rapidly as $M\to 0^{+}$ .

For piecewise cubic Hermite interpolation, the leading local error is proportional to $\left|A^{\left(4\right)}\left(M\right)\right|{\left(\text{Δ}M\right)}^4$ . An asymptotic error monitor is therefore

$\rho \left(M\right)={\left|A^{\left(4\right)}\left(M\right)\right|}^{1/4}~C{M}^{-11/12},$ (65)

where $C$ is an unimportant positive constant. Equidistributing this monitor on a truncated interval $M_{\text{min}}\le M\le M_{\text{max}}$ gives

$M_i={\left[M_{\text{min}}^{1/12}+\frac{i}{N}\left(M_{\text{max}}^{1/12}-M_{\text{min}}^{1/12}\right)\right]}^{12},i=0,\cdots ,N.$ (66)

Here $M_{\text{min}}>0$ , so the singular point itself is excluded. In the formal limit $M_{\text{min}}\to 0^{+}$ , Equation (66) reduces to $M_i/M_{\text{max}}~{\left(i/N\right)}^{12}$ .

The same reasoning gives exponents 6 and 9 for linear and quadratic inverse interpolation, respectively, while ordinary cubic and cubic Hermite interpolation both give exponent 12. These exponents are not empirical fitting parameters; they follow from the cube-root singularity and the corresponding interpolation-error formulas. In practical calculations, the asymptotic distribution may be replaced by a simpler truncated geometric grid, but its stretching direction and required concentration are determined by the analysis above.

8. Construction of Stretched Grids

The preceding asymptotic formulas can now be converted into practical node distributions. Let an inverse interpolant have a local error of the form

${ℰ}_r\left(M\right)\simeq C_r{A}^{\left(r\right)}\left(M\right){\left(\Delta M\right)}^r,$ (67)

where $r=2$ for linear interpolation, $r=3$ for quadratic interpolation, and $r=4$ for ordinary cubic or cubic Hermite interpolation. A suitable grid monitor is

${\rho }_r\left(M\right)={\left|A^{\left(r\right)}\left(M\right)\right|}^{1/r}.$ (68)

The nodes are selected by equidistributing the cumulative monitor:

${\displaystyle {\int }_{M_a}^{M_i}{\rho }_r\left(M\right)\text{d}M}=\frac{i}{N}{\displaystyle {\int }_{M_a}^{M_b}{\rho }_r\left(M\right)\text{d}M},i=0,\cdots ,N.$ (69)

This makes ${\rho }_r\left(M_i\right)\text{Δ}{M}_i$ approximately constant and therefore balances the leading interpolation error among the intervals.

8.1. Closed-Form Asymptotic Grids

At the parabolic limit, $A\left(M\right)~c{M}^{1/3}$ , and consequently

$\left|A^{\left(r\right)}\left(M\right)\right|~C_r{M}^{1/3-r},{\rho }_r\left(M\right)~{\tilde{C}}_r{M}^{-1+1/\left(3r\right)}.$(70)

Substitution in Equation (69) gives the explicit truncated grid

$M_i={\left[M_a^{1/\left(3r\right)}+\frac{i}{N}\left(M_b^{1/\left(3r\right)}-M_a^{1/\left(3r\right)}\right)\right]}^{3r}.$ (71)

The stretching exponents can also be derived separately and more directly from the error formula for each interpolant. This makes clear why a different grid is associated with each approximation order.

8.1.1. Linear Interpolation

The local error is proportional to

$\left|A^{″}\left(M\right)\right|{\left(\text{Δ}M\right)}^2.$ (72)

Since $\left|A^{″}\left(M\right)\right|~C{M}^{-5/3}$ , approximately equal error in every interval requires

$M^{-5/3}{\left(\Delta M\right)}^2\simeq \text{constant},\Delta M\propto M^{5/6}.$ (73)

Treating the node index $i$ as a continuous variable gives

$\frac{\text{d}M}{\text{d}i}=K{M}^{5/6}.$ (74)

Integration yields $M^{1/6}\propto i$ , and hence

$M_i\propto {\left(\frac{i}{N}\right)}^6.$ (75)

8.1.2. Quadratic Interpolation

Here the leading error is proportional to

$\left|A^{\left(3\right)}\left(M\right)\right|{\left(\text{Δ}M\right)}^3.$ (76)

Using $\left|A^{\left(3\right)}\left(M\right)\right|~C{M}^{-8/3}$ gives

$M^{-8/3}{\left(\Delta M\right)}^3\simeq \text{constant},\Delta M\propto M^{8/9}.$ (77)

Therefore

$\frac{\text{d}M}{\text{d}i}=K{M}^{8/9},M^{1/9}\propto i,$ (78)

and the quadratic grid is

$M_i\propto {\left(\frac{i}{N}\right)}^9.$ (79)

8.1.3. Cubic Interpolation

For ordinary cubic interpolation,

${ℰ}_{\text{C}}\propto \left|A^{\left(4\right)}\left(M\right)\right|{\left(\Delta M\right)}^4.$ (80)

Because $\left|A^{\left(4\right)}\left(M\right)\right|~C{M}^{-11/3}$ ,

$M^{-11/3}{\left(\Delta M\right)}^4\simeq \text{constant},\Delta M\propto M^{11/12}.$ (81)

It follows that

$\frac{\text{d}M}{\text{d}i}=K{M}^{11/12},M^{1/12}\propto i,$ (82)

and consequently

$M_i\propto {\left(\frac{i}{N}\right)}^{12}.$ (83)

The cubic Hermite error contains the same fourth inverse derivative and four powers of the local interval length. It therefore produces the same exponent 12, although its data consist of two function values and two endpoint derivatives rather than four function values.

Thus, the grids appropriate to the four interpolation formulas are

$\begin{array}{l}{M}_i^{\text{L}}={\left[M_a^{1/6}+s_i\left(M_b^{1/6}-M_a^{1/6}\right)\right]}^6,\\ M_i^{\text{Q}}={\left[M_a^{1/9}+s_i\left(M_b^{1/9}-M_a^{1/9}\right)\right]}^9,\\ M_i^{\text{C}}=M_i^{\text{H}}={\left[M_a^{1/12}+s_i\left(M_b^{1/12}-M_a^{1/12}\right)\right]}^{12},\end{array}$ (84)

where $s_i=i/N$ . The superscripts L, Q, C, and H denote linear, quadratic, ordinary cubic, and cubic Hermite interpolation, respectively. These distributions cluster nodes increasingly strongly near $M=0$ as the order of interpolation increases.

For $e\ne 1$ , the inverse derivative at $M=0$ is large but finite, so the endpoint $M_0=0$ may be included in the elliptic and hyperbolic tables. At the exact parabolic limit, the cube-root derivative is infinite and the Hermite table must instead start from a small positive value $M_a=M_{\text{min}}>0$ , unless the first interval is treated directly by the analytic cube-root expansion.

8.2. Exact Monitor for Finite Eccentricity

The limiting formulas provide a simple grid and the correct concentration as $e\to 1^{\pm }$ . For a specified finite eccentricity, an optimised grid can be constructed from the exact inverse derivatives. Let

$p=\frac{\text{d}M}{\text{d}A},q=\frac{{\text{d}}^{2}M}{\text{d}{A}^2},r_f=\frac{{\text{d}}^{3}M}{\text{d}{A}^3},s_f=\frac{{\text{d}}^{4}M}{\text{d}{A}^4}.$ (85)

The required inverse derivatives are

$A^{″}\left(M\right)=-\frac{q}{p^3},$

$A^{\left(3\right)}\left(M\right)=\frac{3q^2}{p^5}-\frac{r_f}{p^4},$

$A^{\left(4\right)}\left(M\right)=-\frac{s_f}{p^5}+\frac{10q{r}_f}{p^6}-\frac{15q^3}{p^7}.$ (86)

For the elliptic equation, $A=E$ and

$p=1-e\cos E,q=e\sin E,r_f=e\cos E,s_f=-e\sin E.$ (87)

For the hyperbolic equation, $A=H$ and

$p=e\cosh H-1,q=e\sinh H,r_f=e\cosh H,s_f=e\sinh H.$ (88)

These formulas provide ${\rho }_2$ , ${\rho }_3$ , or ${\rho }_4$ without numerical differentiation.

Moreover, the grid may be generated in anomaly space, where the forward mapping is explicitly known. Since $\text{d}M=p\text{d}A$ , define

$R_r\left(A\right)={\displaystyle {\int }_{A_a}^A{\rho }_r\left(M\left(\alpha \right)\right)p\left(\alpha \right)\text{d}\alpha }.$ (89)

Choose $A_i$ from

$R_r\left(A_i\right)=R_r\left(A_a\right)+\frac{i}{N}\left[R_r\left(A_b\right)-R_r\left(A_a\right)\right],$ (90)

and finally set

$M_i=E_i-e\sin E_i$ (91)

in the elliptic case, or

$M_i=e\sinh H_i-H_i$ (92)

in the hyperbolic case. Therefore, construction of the optimised table does not require repeated inversion of Kepler’s equation at arbitrary mean anomalies.

8.3. Practical Comparison Grids

For the numerical comparisons, the same interval and number of stored values will be used for all grids. A uniform grid is

$M_i=M_a+\frac{i}{N}\left(M_b-M_a\right).$ (93)

A simple one-parameter geometric grid is defined through

$\Delta M_i=\Delta M_1{g}^{i-1},\Delta M_1=\left(M_b-M_a\right)\frac{g-1}{g^N-1},g>1,$ (94)

with $M_i=M_a+{\displaystyle {\sum }_{j=1}^i\Delta M_j}$ . The optimised grid is obtained from the exact monitor above, while Equation (84) supplies its closed-form parabolic-limit approximation.

A classical $C^2$ cubic spline provides an alternative to the local piecewise cubic Hermite approximation, but its increased smoothness does not by itself guarantee a smaller pointwise error. In representative elliptic and hyperbolic cases on the exact-monitor grid, the local Hermite construction gave approximately one half of the maximum error of a clamped cubic spline. It also remains entirely local and avoids the global tridiagonal solution required to construct the spline. On an unsuitable uniform grid, however, a classical spline can substantially reduce the severe endpoint-slope overshoot of the Hermite approximation. The comparison therefore reinforces the main conclusion of this section: the decisive improvement comes from combining the local Hermite formula with a grid consistent with the singular asymptotic behaviour, rather than from increasing continuity alone. After Newton correction, the two predictors require nearly identical numbers of iterations.

9. Numerical Results

This section provides the numerical validation of the proposed approach. We apply the grids to inverse cubic Hermite interpolation and retain enough detail to make the calculations reproducible.

9.1. Numerical Protocol

The interval $0\le M\le \pi$ is divided into $N=64$ subintervals, giving 65 stored values. At every node the anomaly and its exact inverse derivative are stored. Reference values at the test points are obtained by 90 steps of bracketed bisection applied independently to the elliptic or hyperbolic Kepler equation. All calculations use radians and IEEE 754 double-precision arithmetic.

For the elliptic problem, each tabulated anomaly $E_i$ satisfies

$E_i-e\sin E_i=M_i,{E^{\prime }}_i=\frac{1}{1-e\cos E_i}.$ (95)

For the hyperbolic problem, the corresponding data are

$e\sinh H_i-H_i=M_i,{H^{\prime }}_i=\frac{1}{e\cosh H_i-1}.$ (96)

For all Newton-type corrections, convergence requires both the anomaly correction and the equation residual to be no greater than 10⁻¹⁴. This two-part test is important near the parabolic limit because the derivative of Kepler’s equation can be small. The tolerance is substantially smaller than the interpolation errors reported below.

Within each interval $\left[M_i,M_{i+1}\right]$ , the inverse Hermite polynomial is evaluated from the two stored anomalies and the two exact inverse slopes. If $A$ denotes either $E$ or $H$ , the pointwise predictor error is

${\epsilon }_j=A_H\left(M_j\right)-A_{\text{ref}}\left(M_j\right).$ (97)

For $n_t$ test values, the reported measures are

${\epsilon }_{\text{max}}=\underset{1\le j\le n_t}{\text{max}}\left|{\epsilon }_j\right|,{\epsilon }_{\text{RMS}}={\left(\frac{1}{n_t}{\displaystyle \sum }_{j=1}^{n_t}{\epsilon }_j^2\right)}^{1/2}.$ (98)

For the iteration and correction tables, the reported mean is the arithmetic average of the counts over all mean anomalies—and, where applicable, all eccentricities—in the stated test set. The maximum is the largest count observed in that set and is included to quantify worst-case performance, which may not be apparent from the mean alone.

Figure 6. Normalised node positions for $N=64$ . The asymptotic Hermite grid follows $M_i/\pi ={\left(i/N\right)}^{12}$ .

For the interpolation-error tables, the test set is the union of 300,001 values uniformly spaced on $\left[0,\pi \right]$ and 300,001 values geometrically spaced from 10⁻¹⁸ to $\pi$ . The second set resolves the first grid intervals, where the near-parabolic error is largest.

Four node distributions are examined under identical conditions: uniform, geometric, the closed-form asymptotic distribution, and the exact-monitor distribution. The geometric value $g=1.19$ is used in the representative comparison. The asymptotic Hermite grid uses $M_i=\pi {\left(i/N\right)}^{12}$ , while the exact-monitor grid equidistributes ${\left|A^{\left(4\right)}\left(M\right)\right|}^{1/4}$ .

Figure 6 compares the uniform grid, a geometric grid with $g=1.19$ , and the closed-form Hermite grid derived from the parabolic asymptotics. The logarithmic scale shows that the stretched grids place many more nodes near $M=0$ , where the inverse solution changes most rapidly.

9.2. Elliptic Application

The elliptic calculations cover $e=0.9$ through $e=0.99999$ . For the most difficult case, the inverse slope at the origin is

${\frac{\text{d}E}{\text{d}M}|}_{M=0}=\frac{1}{1-e}={10}^5.$ (99)

On a uniform grid, this large slope is used across the comparatively wide first interval and the Hermite polynomial develops a severe overshoot. The maximum interpolation error is $7.269\times {10}^2$ . With geometric stretching the maximum error falls to $3.906\times {10}^{-5}$ , and with the exact monitor ${\rho }_E={\left|E^{\left(4\right)}\right|}^{1/4}$ it falls further to $4.588\times {10}^{-6}$ . The full eccentricity sweep is reported in Table 1.

Table 1. Elliptic inverse-Hermite maximum errors on $0\le M\le \pi$ with $N=64$ . The geometric factor is selected from $1.00\le g\le 1.30$ in increments of 0.01.

$e$	Uniform	Best $g$	Geometric	Exact monitor
0.90000	$3.815\times {10}^{-3}$	1.08	$2.360\times {10}^{-7}$	$1.822\times {10}^{-7}$
0.99000	$4.604\times {10}^{-1}$	1.13	$1.246\times {10}^{-6}$	$8.293\times {10}^{-7}$
0.99900	$6.970\times {10}^0$	1.19	$5.032\times {10}^{-6}$	$1.865\times {10}^{-6}$
0.99990	$7.242\times {10}^1$	1.25	$1.352\times {10}^{-5}$	$3.313\times {10}^{-6}$
0.99999	$7.269\times {10}^2$	1.30	$3.906\times {10}^{-5}$	$4.588\times {10}^{-6}$

The exact-monitor error remains below $4.6\times {10}^{-6}$ throughout the five elliptic cases. Its advantage over the best constant-factor geometric grid increases from a factor of 1.30 at $e=0.9$ to 8.51 at $e=0.99999$ .

9.3. Hyperbolic Application

The calculation is repeated entirely in real variables for $e=1.00001$ . In this case

${\frac{\text{d}H}{\text{d}M}|}_{M=0}=\frac{1}{e-1}={10}^5.$ (100)

The uniform-grid maximum error is again $7.269\times {10}^2$ . The geometric and exact-monitor grids reduce it to $3.906\times {10}^{-5}$ and $4.692\times {10}^{-6}$ , respectively. The close agreement with the elliptic results confirms the common asymptotic structure on the two sides of $e=1$ . Table 2 gives the complete hyperbolic sweep.

Table 2. Hyperbolic inverse-Hermite maximum errors on $0\le M\le \pi$ with $N=64$ , using the same geometric-factor search as in the elliptic case.

$e$	Uniform	Best $g$	Geometric	Exact monitor
1.10000	$4.336\times {10}^{-3}$	1.08	$2.586\times {10}^{-7}$	$1.935\times {10}^{-7}$
1.01000	$4.627\times {10}^{-1}$	1.13	$1.483\times {10}^{-6}$	$8.169\times {10}^{-7}$
1.00100	$6.972\times {10}^0$	1.19	$5.982\times {10}^{-6}$	$1.981\times {10}^{-6}$
1.00010	$7.242\times {10}^1$	1.25	$1.606\times {10}^{-5}$	$3.417\times {10}^{-6}$
1.00001	$7.269\times {10}^2$	1.30	$3.906\times {10}^{-5}$	$4.692\times {10}^{-6}$

The optimised hyperbolic grid is more accurate than the best geometric grid for every tested eccentricity. Its advantage grows from 1.34 at $e=1.1$ to 8.33 at $e=1.00001$ , while its maximum error remains below $4.7\times {10}^{-6}$ .

Figure 7 displays the same comparison. Relative to the uniform grid, the exact-monitor grid improves the representative maximum error by more than eight orders of magnitude. It also improves the already stretched geometric grid by a factor of approximately eight. This calculation is the direct numerical application of the singular asymptotics: the grid is concentrated where the inverse derivatives, and hence the interpolation error, are largest.

Figure 7. Maximum Hermite interpolation errors for the representative near-parabolic elliptic and hyperbolic cases.

These examples use only the Hermite predictor; no Newton correction has been applied. A subsequent Newton-type correction may reduce the residual to the desired tolerance, but the stretched grid has already supplied a uniformly useful initial approximation in the difficult near-parabolic region.

9.4. Comparison of the Iterative Corrections

The actual convergence of the iterative methods discussed earlier was also examined. Because the greatest difficulty occurs for small mean anomaly and large eccentricity, the test interval was

${0.001}^{\circ }\le M\le 2^{\circ },$ (101)

sampled at 2000 equally spaced values for each of $e=0.9$ , 0.99, 0.999, 0.9999, and 0.99999. The angles were converted to radians before calculation.

To compare the convergence orders without favouring any method through an empirical initial formula, every iteration was started from the midpoint of the guaranteed elliptic bracket

$M\le E\le \text{min}\left(M+e,\pi \right),E_0=\frac{M+\text{min}\left(M+e,\pi \right)}{2}.$ (102)

One iteration means one complete update of the corresponding method. The calculation was stopped only when both

$\left|E_{n+1}-E_n\right|\le {10}^{-14},\left|E_{n+1}-e\sin E_{n+1}-M\right|\le {10}^{-14}$ (103)

were satisfied. For the fixed-parameter procedure we used $\alpha =-1$ in Equation (7), giving the familiar fixed-point form

$E_{n+1}=M+e\sin E_n.$ (104)

Table 3 compares the mean and maximum iteration counts for the elliptic problem.

Table 3. Mean iteration counts over 2000 equally spaced mean anomalies for each eccentricity, with maximum counts in parentheses, for the elliptic Kepler equation. All methods use the common systematic starter (102) and stopping test (103). The fifth-order entry is Householder’s formula (26) with $d=4$ .

Method	$e=0.9$	$e=0.99$	$e=0.999$	$e=0.9999$	$e=0.99999$
Fixed point	238.05 (277)	418.16 (2597)	424.20 (19788)	408.65 (21466)	405.88 (19970)
Newton	6.00 (6)	5.72 (8)	5.59 (11)	5.58 (11)	5.57 (11)
Chebyshev	4.30 (5)	4.29 (6)	4.23 (8)	4.22 (8)	4.22 (8)
Halley	4.00 (4)	4.06 (6)	4.05 (7)	4.05 (7)	4.05 (7)
Danby, order 4	3.63 (4)	3.54 (5)	3.45 (6)	3.44 (6)	3.44 (6)
Householder, order 5	3.50 (4)	3.18 (4)	3.16 (5)	3.16 (5)	3.16 (5)

The fixed-point results display the weakness of linear convergence. Its local contraction factor is approximately $e\cos E^{*}$ and approaches unity in the difficult near-parabolic region; consequently, individual cases require roughly twenty thousand updates. The higher-order methods are far less sensitive. Over the entire calculation, Newton requires no more than 11 updates, Chebyshev no more than 8, Halley no more than 7, Danby no more than 6, and fifth-order Householder no more than 5.

The table also illustrates why formal order alone should not be used as the only efficiency measure. Chebyshev and Halley both have cubic convergence, but Halley’s rational correction gives slightly smaller maximum counts here. Danby’s fourth-order nested correction and Householder’s fifth-order formula save additional iterations, although each update evaluates more derivatives and performs more algebra. For Kepler’s equation these derivatives are available analytically, so the additional work is modest and the comparison clearly reveals the benefit of higher local order.

Finally, these counts deliberately use only the bracket midpoint. In the proposed algorithm a local inverse interpolant replaces this midpoint and is already close to the root. The following calculations examine this combined predictor-corrector procedure in both orbital regimes.

9.5. Bisection and Inverse-Interpolation Predictors

A symmetric calculation was performed on each side of the parabolic limit. For the elliptic cases, 100 eccentricities were sampled uniformly over $0.9\le e\le 0.999$ ; for the hyperbolic cases, 100 values were sampled over $1.001\le e\le 1.1$ . For every eccentricity, 2000 values of $M$ were used on ${0.001}^{\circ }\le M\le 2^{\circ }$ , giving 200,000 equations in each regime. The elliptic bracket was

$M\le E\le \text{min}\left(M+e,\pi \right),$ (105)

while $0\le H\le 1$ brackets all the hyperbolic roots in this particular test region. After three or four bisections, linear or quadratic inverse interpolation supplied the initial value for Newton’s method. The Newton iteration was stopped at an equation residual of 10⁻¹⁰.

Table 4 reports the correction counts following bisection and local inverse interpolation in both regimes.

Table 4. Newton corrections following bisection and local inverse interpolation. The percentage is the fraction requiring no more than three Newton updates.

Regime	Predictor	Mean	Maximum	At most 3
Elliptic	3 bisections + linear	2.304	5	99.978%
Elliptic	3 bisections + quadratic	1.748	4	99.997%
Elliptic	4 bisections + quadratic	1.316	2	100.000%
Hyperbolic	3 bisections + linear	2.412	5	99.977%
Hyperbolic	3 bisections + quadratic	1.816	4	99.996%
Hyperbolic	4 bisections + quadratic	1.386	2	100.000%

The results on the two sides are remarkably close. Four bisections followed by quadratic inverse interpolation never require more than two Newton corrections in either set of 200,000 equations. This provides a completely systematic realisation of the sequence

$\text{guaranteed bracket}\to \text{local approximation}\to \text{rapid correction}.$ (106)

9.6. Refinement of the Hermite Predictor

The exact-monitor Hermite grid was next coupled to Newton correction for the representative cases $e=0.999$ and $e=1.001$ . The same 2000 mean anomalies were used, and corrections were continued until both the step and residual were below 10⁻¹⁴. Table 5 shows the effect of increasing the number of grid intervals.

Table 5. Newton corrections after exact-monitor cubic Hermite prediction.

Regime	$N$	Mean	Maximum	At most 1
Elliptic, $e=0.999$	32	1.975	2	2.55%
	64	1.881	2	11.95%
	128	1.371	2	62.90%
	256	0.999	1	100.00%
Hyperbolic, $e=1.001$	32	1.975	2	2.50%
	64	1.881	2	11.90%
	128	1.390	2	61.00%
	256	1.000	1	100.00%

With 256 intervals, every elliptic and hyperbolic test requires at most one Newton update. With 128 intervals, 62.9% of the elliptic cases and 61.0% of the hyperbolic cases require at most one update, and none requires more than two. The closely matched counts provide another direct numerical confirmation that the same near-parabolic grid strategy applies on both sides of $e=1$ .

9.7. Additional Grid and Timing Checks

To test the grid beyond the two near-parabolic eccentricities in Table 5, Table 6 reports a separate refinement study. For each row we equidistributed the exact ${\left|A^{\left(4\right)}\right|}^{1/4}$ monitor over $0\le M\le \pi$ by trapezoidal quadrature on 250,001 uniformly spaced anomaly values, using 2000 uniform test anomalies in ${10}^{-8}\le M\le \pi$ . Reference anomalies were obtained from 85 bisection steps. The entries are the maximum absolute error of the Hermite predictor before any correction. The quadrature and sample interval are specified because these values are not directly interchangeable with the more densely sampled errors in Table 1 and Table 2.

For a limited runtime check on the same 2000 anomaly values, we compared vectorised NumPy implementations of Newton correction starting from a simple initial value and from the $N=128$ Hermite predictor. The elliptic simple start was the midpoint of the bracket in Equation (102); the hyperbolic simple start was $H=1$ . On this execution environment, median time per 2000 solves over seven batches of 100 runs was 0.798 ms versus 0.392 ms for $e=0.999$ , and 0.315 ms versus 0.192 ms for $e=1.001$ , respectively. These timings include predictor evaluation and Newton updates but exclude grid construction, which took approximately 22 ms and 16 ms, respectively. They depend on the language, hardware, batching, and number of repeated solves: grid setup must be amortised over many queries. They do not measure Mikkola’s, Raposo-Pulido and Peláez’s, or Tommasini and Olivieri’s implementations. Consequently, the present results establish grid accuracy and illustrate within-implementation costs, but do not establish a speed advantage over those published solvers. A direct comparative implementation under common accuracy and timing criteria remains useful future work.

Table 6. Additional exact-monitor grid sensitivity: maximum uncorrected Hermite anomaly error in radians over 2000 test values.

Case/$N$	32	64	128	256
Elliptic $e=0.90$	$2.898\times {10}^{-6}$	$1.815\times {10}^{-7}$	$1.135\times {10}^{-8}$	$7.096\times {10}^{-10}$
Elliptic $e=0.99$	$1.250\times {10}^{-5}$	$7.846\times {10}^{-7}$	$4.910\times {10}^{-8}$	$3.070\times {10}^{-9}$
Elliptic $e=0.999$	$2.808\times {10}^{-5}$	$1.764\times {10}^{-6}$	$1.105\times {10}^{-7}$	$6.910\times {10}^{-9}$
Elliptic $e=0.9999$	$4.691\times {10}^{-5}$	$2.951\times {10}^{-6}$	$1.848\times {10}^{-7}$	$1.156\times {10}^{-8}$
Hyperbolic $e=1.0001$	$3.956\times {10}^{-5}$	$2.470\times {10}^{-6}$	$1.544\times {10}^{-7}$	$9.647\times {10}^{-9}$
Hyperbolic $e=1.001$	$2.371\times {10}^{-5}$	$1.481\times {10}^{-6}$	$9.251\times {10}^{-8}$	$5.782\times {10}^{-9}$
Hyperbolic $e=1.01$	$1.058\times {10}^{-5}$	$6.609\times {10}^{-7}$	$4.130\times {10}^{-8}$	$2.581\times {10}^{-9}$
Hyperbolic $e=1.10$	$2.891\times {10}^{-6}$	$1.686\times {10}^{-7}$	$1.034\times {10}^{-8}$	$6.146\times {10}^{-10}$

9.8. Alternative Methods

Three additional ideas were tested: continuation in eccentricity, Wynn’s nonlinear sequence transformation, and homotopy. They are useful for showing how Kepler’s equation can be embedded in a sequence of simpler problems, but the calculations below also clarify their limitations relative to direct inverse interpolation.

9.8.1. Continuation in Eccentricity

Following the general solution-path viewpoint of Garcia and Zangwill [10], implicit differentiation gives

$\frac{\text{d}E}{\text{d}e}=\frac{\sin E}{1-e\cos E},\frac{\text{d}H}{\text{d}e}=-\frac{\sinh H}{e\cosh H-1}.$ (107)

For each regime, the reported continuation statistics are computed over the same 2000 equally spaced mean anomalies on ${0.001}^{\circ }\le M\le 2^{\circ }$ . The elliptic continuation starts from the known solution $E=M$ at $e=0$ and ends at $e=0.999$ . The hyperbolic calculation starts from a converged solution at $e=1.1$ and ends at $e=1.001$ . At each of 16 parameter levels, a tangent prediction is followed by Newton correction. We compare uniform parameter levels with the stretched distribution

$s_j=1-{\left(1-\frac{j}{16}\right)}^3,e_j=e_0+s_j\left(e_{\text{target}}-e_0\right),$ (108)

which places smaller changes near the difficult target.

Table 7 compares uniform and stretched continuation steps for the elliptic and hyperbolic problems.

Table 7. Continuation with 16 eccentricity steps for each of 2000 mean anomalies. The starting solution is not included in the correction counts.

Regime	Parameter steps	Mean total	Maximum total	Maximum per step
Elliptic, $e=0.999$	Uniform	39.25	40	9
	Stretched	41.10	42	4
Hyperbolic, $e=1.001$	Uniform	45.18	48	6
	Stretched	38.78	42	4

Continuation is reliable, and stretching substantially reduces the largest number of corrections at an individual parameter level. In this test it slightly increases the elliptic total but reduces the hyperbolic total. Its total work is nevertheless larger than a direct solve when only one eccentricity is needed. Its natural advantage occurs when solutions at all intermediate eccentricities are themselves required.

9.8.2. Wynn’s Nonlinear Sequence Transformation

Let $A_k$ denote either $E\left(e_k\right)$ or $H\left(e_k\right)$ and form five eccentricities approaching the target according to

$e_k=e_{\text{*}}+\left(e_0-e_{\text{*}}\right)2^{-k},k=0,\cdots ,4.$ (109)

Wynn’s epsilon algorithm [11] is

${ϵ}_{-1}^{\left(k\right)}=0,{ϵ}_0^{\left(k\right)}=A_k,{ϵ}_{r+1}^{\left(k\right)}={ϵ}_{r-1}^{\left(k+1\right)}+\frac{1}{{ϵ}_r^{\left(k+1\right)}-{ϵ}_r^{\left(k\right)}}.$ (110)

The even columns provide transformed approximations. Five values give the estimate ${ϵ}_4^{\left(0\right)}$ . Table 8 compares this value with the last untransformed member $A_4$ .

Table 8. Errors before and after the five-term Wynn epsilon transformation.

Regime	Raw maximum	Wynn maximum	Raw RMS	Wynn RMS
Elliptic, $e_{\text{*}}=0.999$	$6.691\times {10}^{-2}$	$2.104\times {10}^{-1}$	$3.110\times {10}^{-2}$	$6.573\times {10}^{-2}$
Hyperbolic, $e_{\text{*}}=1.001$	$6.684\times {10}^{-2}$	$2.094\times {10}^{-1}$	$3.181\times {10}^{-2}$	$6.543\times {10}^{-2}$

For this direct five-value construction the transformation does not accelerate convergence; it increases both the maximum and RMS errors. This negative result is informative. The anomaly has a singular dependence on eccentricity near $e=1$ , and a short geometrically spaced sequence does not match that structure. Wynn’s transformation should therefore not be used here without an additional asymptotic change of variable.

9.8.3. Homotopy

The homotopy statistics are likewise computed over the same 2000 mean anomalies in each regime. Motivated by the homotopy treatment of Kepler’s equations by Fitz-Coy and Jang [12], for the elliptic equation we use

$F_E\left(E,t\right)=E-te\sin E-M=0,$ (111)

and for the hyperbolic equation

$F_H\left(H,t\right)=H+t\left(e\sinh H-2H\right)-M=0.$ (112)

At $t=0$ both problems have the known solution $A=M$ , while $t=1$ gives the required Kepler equation. Differentiating $F\left(A,t\right)=0$ produces a first-order ordinary differential equation for the anomaly. It is integrated by the explicit midpoint rule, using either uniform values of $t$ or

$t_j=1-{\left(1-\frac{j}{K}\right)}^3.$ (113)

The endpoint is then used as the starting value for Newton correction.

Table 9 reports the homotopy predictor errors and subsequent Newton correction counts.

Table 9. Midpoint homotopy prediction and subsequent Newton correction for each of 2000 mean anomalies.

Regime	(K, grid)	Max. predictor error	RMS error	Mean corrections	Max. corrections
Elliptic	(16, uniform)	$1.189\times {10}^{-1}$	$4.019\times {10}^{-2}$	3.617	8
	(16, stretched)	$1.893\times {10}^{-2}$	$1.209\times {10}^{-2}$	3.764	5
	(32, uniform)	$7.920\times {10}^{-2}$	$1.640\times {10}^{-2}$	3.236	7
	(32, stretched)	$6.126\times {10}^{-3}$	$3.440\times {10}^{-3}$	3.049	4
Hyperbolic	(16, uniform)	$1.185\times {10}^{-1}$	$3.983\times {10}^{-2}$	3.628	8
	(16, stretched)	$1.890\times {10}^{-2}$	$1.181\times {10}^{-2}$	3.779	5
	(32, uniform)	$7.905\times {10}^{-2}$	$1.630\times {10}^{-2}$	3.235	7
	(32, stretched)	$6.118\times {10}^{-3}$	$3.360\times {10}^{-3}$	3.048	4

The stretched homotopy grid gives a consistently more accurate predictor. With $K=32$ , it reduces the maximum predictor error by approximately a factor of 13 and reduces the largest Newton count from 7 to 4 in both regimes. Nevertheless, the direct exact-monitor Hermite predictor remains more accurate and requires fewer final corrections. The alternative calculations therefore reinforce the central conclusion: parameter stretching is useful, but stretching based directly on the singular inverse solution is the most effective strategy among those examined here.

9.9. Principal Numerical Findings

The calculations establish several points. First, a uniform grid becomes unusable close to $e=1$ because it applies the very large endpoint slope over an interval that is too wide. Second, a constant-factor geometric grid restores high accuracy, but its factor must increase as the transition layer narrows. Third, equidistribution of the analytically derived monitor improves the result for every eccentricity without introducing an empirical starter formula. Finally, the nearly identical elliptic and hyperbolic trends confirm that the common cube-root asymptotics correctly describe the approach to the parabolic limit from both sides. These results supply the direct numerical validation of the proposed interpolation-and-stretching strategy. The iteration counts verify the expected progression from linear fixed-point convergence to the quadratic, cubic, fourth-order, and fifth-order methods; when the Hermite predictor is used, only a small number of final corrections is normally required. Bisection followed by local inverse interpolation is also reliable in both regimes. Finally, continuation and homotopy benefit from stretched parameter steps, whereas the direct five-term Wynn transformation does not accommodate the near-parabolic singularity. Taken together, the results favour the proposed combination of analytical asymptotics, stretched inverse interpolation, and a short final correction.

10. Conclusion

Kepler’s equation provides a useful unified example for studying several ideas in numerical analysis. The comparisons presented here illustrate the importance of choosing an accurate initial approximation, particularly near the nearly singular regimes. The elliptic and hyperbolic calculations show that asymptotically guided stretched grids can maintain accurate local inverse interpolation as the parabolic limit is approached from either side. With a good Hermite predictor, Newton’s method normally requires only a small number of final corrections, so the quality and placement of the interpolation data can be more important in practice than increasing the formal order of the corrector. Continuation, nonlinear sequence transformation, and homotopy provide useful alternative viewpoints, while the numerical comparisons clarify their advantages and limitations. From an educational viewpoint, the problem connects nonlinear iteration, inverse interpolation, asymptotic analysis, grid adaptation, continuation, sequence transformation, and homotopy within a single familiar equation.

Acknowledgements

ChatGPT (OpenAI) was used to assist with manuscript preparation, language editing, checking mathematical derivations, and production of numerical tables. The author remains solely responsible for the ideas, interpretations, conclusions, accuracy, and research presented.

The author would like to thank Dr. Esmat Bekir (see Ref. [7]) for several discussions about Kepler’s equation.

Data and Code Availability

The numerical procedures and sampling specifications are described in the text. No public code repository accompanies this submission.

Appendix

A. The Parabolic Limit and Barker’s Exact Solution

The elliptic and hyperbolic Kepler equations are defined for $e<1$ and $e>1$ , respectively. Although the eccentricity may be set formally equal to unity in either equation, this substitution alone does not produce the physical parabolic Kepler equation. The parabolic limit must instead be taken while the periapsis distance and dimensional time are retained. This appendix shows that the correctly scaled limits from both sides of $e=1$ agree with Barker’s equation and its exact real solution.

A.1. Formal Substitution of $e=1$

For an elliptic orbit, Kepler’s equation is

$M=E-e\sin E,0\le e<1.$ (114)

Formal substitution of $e=1$ gives $M=E-\text{sin}E$ . For a hyperbolic orbit, using the positive magnitude of the semimajor axis, the corresponding equation is

$M_h=e\sinh H-H,e>1,$ (115)

which formally gives $M_h=\sinh H-H$ at $e=1$ . These endpoint equations are not identical, and neither is Barker’s equation. Their mean-motion scalings degenerate as $e$ approaches unity.

A.2. Common Local Behaviour near Periapsis

The local cube-root expansions and their implications for grid design were derived in Section 7. Here we retain only the dimensional scaling needed to connect those inverse mappings to the physical parabolic time law.

A.3. Limit from the Elliptic Side

Let $q$ denote the periapsis distance, $\mu$ the gravitational parameter, $\nu$ the true anomaly, and $\tau$ the time of periapsis passage. Introduce Barker’s variable $D=\text{tan}\left(\nu /2\right)$ . For an ellipse,

$D={\left(\frac{1+e}{1-e}\right)}^{1/2}\text{tan}\frac{E}{2},a=\frac{q}{1-e}.$ (116)

If $q$ and $D$ are held fixed while $e\to 1^{-}$ , then $E~{\left[2\left(1-e\right)\right]}^{1/2}D$ . Writing $\delta =1-e$ and expanding equation (114) gives

$M=E-\left(1-\delta \right)\sin E~\sqrt{2}{\delta }^{3/2}\left(D+\frac{D^3}{3}\right).$ (117)

Since $t-\tau ={\left(a^3/\mu \right)}^{1/2}M$ and $a=q/\delta$ , the factors in $\delta$ cancel, leaving

$t-\tau =\sqrt{\frac{2q^3}{\mu }}\left(D+\frac{D^3}{3}\right).$ (118)

A.4. Limit from the Hyperbolic Side

For a hyperbola, let $a_h>0$ denote the magnitude of the semimajor axis. Then

$D={\left(\frac{e+1}{e-1}\right)}^{1/2}\tanh \frac{H}{2},a_h=\frac{q}{e-1}.$ (119)

Holding $q$ and $D$ fixed as $e\to 1^{+}$ gives $H~{\left[2\left(e-1\right)\right]}^{1/2}D$ . Writing $\epsilon =e-1$ and expanding Equation (115),

$M_h=\left(1+\epsilon \right)\sinh H-H~\sqrt{2}{\epsilon }^{3/2}\left(D+\frac{D^3}{3}\right).$ (120)

Since $t-\tau ={\left(a_h^3/\mu \right)}^{1/2}{M}_h$ and $a_h=q/\epsilon$ , the factors in $\epsilon$ again cancel, giving precisely Equation (118).

A.5. Barker’s Equation and Its Exact Solution

Defining the dimensionless parabolic time

$B=\left(t-\tau \right)\sqrt{\frac{\mu }{2q^3}},$ (121)

Equation (118) becomes

$B=D+\frac{D^3}{3},D^3+3D-3B=0.$ (122)

The cubic is strictly increasing and therefore has one real solution for every real $B$ . A convenient exact real form is

$D=2\sinh \left[\frac{1}{3}\text{asinh}\left(\frac{3B}{2}\right)\right].$ (123)

An equivalent Cardano form is

$D=\sqrt[\begin{array}{l}3\\ \end{array}]{\frac{3B}{2}+\sqrt{1+\frac{9B^2}{4}}}+\sqrt[\begin{array}{l}3\\ \end{array}]{\frac{3B}{2}-\sqrt{1+\frac{9B^2}{4}}}.$ (124)

Once $D$ is known, the true anomaly follows exactly from $\nu =2\arctan D$ . Thus the correctly scaled elliptic and hyperbolic limits agree with one another and with Barker’s exact parabolic solution.