On the integral solution of elliptic Kepler’s equation

In a recent paper, Philcox, Goodman and Slepian obtain an explicit solution of the elliptic Kepler’s equation (KE) as a quotient of two contour integrals along a Jordan curve C=C(M,e)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathcal {C} = \mathcal {C}(M,e)$$\end{document} that contains the unique real solution of KE but not includes other complex zeros of KE in its interior. The aim of this paper is to study the main issues that arise in the practical implementation of this integral solution. Thus, after a study of the complex zeros of KE, several families of Jordan contours C=C(M,e)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathcal {C} = \mathcal {C}(M,e)$$\end{document} that are suitable for this integral solution are proposed. Since contours with minimal length turn out to be the more accurate for numerical purposes, several families that minimize their length are constructed. Secondly, the approximation of the contour integrals by the composite trapezoidal rule is considered. Recall that this rule is employed in the fast Fourier transform and, in spite of its lower order, displays a spectral convergence as a function of the number of nodes, which implies a very fast convergence. Finally, the results of some numerical experiments are presented to show that such a combination of appropriate contours with the composite trapezoidal rule leads to a powerful numerical method to solve KE with any desired accuracy for all values of eccentricity.


Introduction
The elliptic Kepler's equation which relates the position in an elliptic orbit with eccentricity e ∈ [0, 1) of the two-body problem with the physical time, defines the eccentric anomaly E in terms of the mean anomaly M, and this equation has been one of the most studied transcendental equations in Celestial Mechanics along several centuries (Colwell 1993) and as remarked by R. Battin (Battin 1999, Chapter 4) is connected with the development of many mathematical topics such as Bessel functions, Fourier series, Lagrange expansion theorem and numerical approximation of functions. Thus, the solution of (1) can be written as a sine Fourier series of M in the form where J k (·) is the Bessel function of the first kind of order k. Also Lagrange's approach to solving Kepler's equation (1) by Newton and other high-order iterative methods. In these methods for given e and M, and starting from an approximation E 0 close to E(M; e) a sequence (E k ) k≥1 is defined which converges to the exact solution as fast as possible (see, e.g., Odell and Gooding (1986), Markley (1995), Fukushima (1997), Feinstein and McLaughlin (2006), Davis et al. (2010), Mortari and Elipe (2014), Calvo et al. (2013), Calvo et al. (2017), Raposo-Pulido and Peláez (2017), to mention just a few works).
In this paper we study a different approach to solve Kepler's equation proposed in a recent paper of Philcox et al. (2021) in which they obtain an explicit formula E = E(M; e), solution of Kepler's equation (1), following a technique used by Ullisch (2020) for solving the classical goat's problem. This technique is based on the following theorem of Jackson (1916): Theorem 1.1 Let U ⊆ C be an open simply connected subset and f : U −→ C a nonzero analytic function. For every simple zero z 0 ∈ U of f , there is a closed curve C in U such that . Note: C does not contain other zeros of f than z 0 .
Thus, Philcox et al. (2021) applying the above theorem shows that the solution E = E(M, e) of elliptic Kepler's equation (1) can be given explicitly as the quotient of the two contour integrals where f is the analytical complex valued function of z given by and C = C(M, e) is a Jordan curve that encloses the unique real solution of Eq. (1), z = E(M, e) for (e, M) ∈ D = (0, 1) × (0, π), and f (z) = 0 for all z in C and its interior, except z = E. Hereafter we will say that if a family of contours C = C(M; e), with (e, M) ∈ D, satisfies these conditions, it is admissible for the integral solution (2). In particular, these authors take as C 0 = C(M, e) the family of circles centered at the midpoint μ 0 of M and M + e and radius ρ 0 = e/2, i.e., showing that this family of Jordan curves is admissible for (2). With this choice and approximating the contour integrals of Eq. (2) by some Riemann sums, it is possible to provide explicit approximate solutions of elliptic Kepler's equation for all (e, M) ∈ D that are faster than standard iterative methods that provide approximate solutions to the transcendental Kepler equation (1). We already proved (Calvo et al. 2022) that for solving a similar (but simpler) formula emanating from the collapse's radial evolution in time , the length of the Jordan curve and the numerical quadrature formula used to compute the above line integrals have a big influence in the final result in both aspects of speed and accuracy. Thus, in this paper we study the two main issues that are relevant in the practical application of this integral solution (2), namely the choice of the Jordan curve and the numerical method to compute the line integrals in Eq. (2).
A first point is the choice of the Jordan contour C = C M,e that is admissible for the integrals of (2). To make clear this point in Sect. 2, we study the complex zeros of f (z; M, e) = 0 for M ∈ (0, π) and e ∈ [0, 1). It turns out that (2) is independent on the choice of the contour provided that it satisfies the above suitability condition but in the numerical approximation of these integrals their accuracy is much better for curves with small length. Then in Sect. 3 we focus our attention in the construction of suitable contours that minimize the length of their boundary. Several families are proposed with circular and elliptic contours.
In Sect. 4 we deal with the numerical quadrature method to compute the line integrals of Eq. (2). We study the application of the composite trapezoidal rule to approximate the integrals of (2) along suitable circular and elliptic contours. The choice of this quadrature rule is motivated because in spite of the lower order of the trapezoidal rule (second order), it is very reliable when it is used in the fast Fourier transform (FFT) and even more important, it exhibits spectral convergence with the number N of nodes, i.e., the error in the approximation of the integral of a periodic function (here the numerator or the denominator of (2)) by this quadrature rule with N nodes behaves as O(exp(−α N )) with some positive constant α, as we prove in Appendix A. In our case, the quadratures in (2) depend on the parameters (e, M) ∈ D and therefore α = α(M; e) also depends on these parameters. This spectral behavior of the numerator and denominator of (2) is inherited by the approximation E N to E. Observe that with a standard quadrature rule with order p the error behaves as O(N − p ) and therefore spectral convergence is much stronger than any potential convergence. We will show also that the factor α(M; e) → 0 when (M, e) → (1, 0) that is the singular point of the solution E = E(M; e) of Kepler's equation and consequently spectral convergence fails for values of (M, e) close to (1, 0).
In Sect. 5 we examine several points with the purpose to make the implementation of the quadrature rules as efficient as possible. Observe that the integrals of numerator and denominator of Eq. (2) could be obtained by a FFT by using a standard routine, but this is not the most efficient way because we only need two Fourier coefficients and the remaining are not necessary. Hence, we have made a direct approach to calculate these coefficients. Finally, we present the results of some numerical experiments showing that very high accuracy can be attained with a moderate computational cost.

The complex zeroes of f (z; M, e)
Here we consider Kepler's equation (KE): f (z; M, e) = 0 with complex z = x + iy, where x and y are reals and we study their solutions for the values of the parameters (e, M) ∈ D.
It is worth to remark that as shown in Wintner (page 216) the solution of KE can be considered as the determination of the inverse function z = z(M; e) of the meromorphic function M = z − e sin z, for (e, M) ∈ D, and taking into account that z(M; e) is a multivalued analytic function for each M ∈ (0, π) there is an infinite set of solutions. In this section, our aim is to describe this set of solutions in order to find suitable Jordan curves for the integral solution (2).
First of all, the complex KE can be written in the form therefore, z = x + i y is a zero of f if and only if their real and imaginary parts satisfy Clearly, for e = 0 the equations (6) have the unique solution y = 0, x = M; then we will consider e ∈ (0, 1). On the other hand if y = 0, equations (6) reduce to the real KE M = x − e sin x that has also a unique solution x ∈ [0, π] for all e ∈ (0, 1), M ∈ [0, π]. Moreover, if (x, y) is a solution of (6), (x, −y) is also a solution of these equations, i.e., in the complex plane the solutions are symmetric with respect to the real axis, and if z ∈ C is a solution then its conjugate z is also solution. Because of this, hereafter we will consider only the case y > 0. Under this assumption, the first equation of (6) can be written in the equivalent form cos x = y e sinh y , and introducing the notation g(y) ≡ y sinh y , the equation (7) can be written as Clearly, for all e = e 0 ∈ (0, 1), the function cos x will be defined only for the values of y such that |g(y)| ≤ e 0 , and taking into account that when y > 0, the function g(y) is strict monotonic decreasing with g(0) = 1 and g(+∞) = 0, for all e 0 ∈ (0, 1) there is a unique y 0 > 0 such that g(y 0 ) = e 0 and Eq. (9) is defined only for y ≥ y 0 . Moreover, inasmuch cos x is an even and 2π-periodic function, (9) defines x as a function of y by the multivalued function Next we will substitute the function (10) into the second equation of (6) to determine the values of y that satisfy this equation. The right-hand side of this equation along the functions (10) will be denoted by We study firstly the function ψ + k (y; e 0 ). It can be seen that it is a monotonic function of y for y ≥ y 0 . Moreover, for y = y 0 , g(y 0 ) = e 0 and ϕ + k (y 0 ; k) = 2kπ, which implies that ψ + k (y 0 ; e 0 ) = 2kπ, k ∈ Z, and lim Consequently, by the mean value theorem, for all M ∈ (0, π] the equation has no solution for k ≤ 0, whereas for each k > 0 there exists a unique solution of (12) depending on M that will be denoted by y + 0,k = y + 0,k (M; e 0 ), and these solutions satisfy is a monotonic increasing sequence. Hence, for each e 0 ∈ (0, 1), M ∈ (0, π], the complex equation (5) possesses the complex roots Secondly, we consider the equations Now ψ − k (y; e 0 ) are monotonic increasing functions of y for y ≥ y 0 with ψ − k (y 0 ; e 0 ) = 2 k π and ψ − k (y; e 0 ) → +∞ when y → +∞. Therefore, for k ≥ 1, Eq. (14) does not have solutions and for all k ≤ 0, there is a unique solution y = y − 0,k > y 0 depending on M ∈ (0, π] such that y 0 < y − 0,0 < y − 0,1 < . . . Consequently, the complex equation (5) has also the complex roots From this analysis we conclude • For all e ∈ [0, 1), M ∈ (0, π] there is a unique real equation. • If z ∈ C is a solution of Eq. (5), then z is also solution.
. . , consequently, there are no complex roots of Kepler's Equation in the complex band {z ∈ C ; 0 ≤ Re z ≤ π } and any Jordan curve contained in this band is suitable for the application of the integral formula Eq. (2). In Fig. 1 we display the locus of the complex zeros z − 0 of KE for eccentricity e = 0.5 and M ∈ (0, π). Observe that in general, for all e 0 ∈ (0, 1) and M ∈ (0, π), the two branches of z − 0 of non-real solutions remain in Re (z) < 0 (as shown in Philcox et al. (2021)), but when e 0 → 1 and M → 0, since y → 0, they become arbitrarily close to the origin.

Remark 2
The above study can be repeated for the limit Kepler's equation and for all k = 0, 1, 2, . . . the complex solutions z k are not contained in the complex strip Re z ∈ (0, π] for all M ∈ (0, π).

The Jordan contours C = C(M, e)
In theory, the exact solution E = E(M; e) of Eq. (2) is independent of the Jordan contour C = C(M, e) of the integrals provided that for each e ∈ [0, 1) and M ∈ (0, π) it contains the exact solution E(M; e) and f (z) = 0 for all z in the interior, except z = E(M; e). However, the practical approximation of these integrals depends on the curve C and, as we will see next, curves with smaller lengths lead to essential improvements in the accuracy when the integrals are approximated by quadrature rules. In this section we will consider several Jordan contours C that are designed with the purpose to minimize their length.
First of all, we consider admissible circular contours in which each C is a circle of radius ρ = ρ(M, e) and center μ = μ(M, e) so that for all (e, M) ∈ D, the curve C(μ, ρ) contains the exact solution of Kepler equation: E = E(M; e) and, moreover, are admissible for (2). For this Jordan curve (16) the integral solution (2) becomes with A crucial point for the choice of the contours (16) is that as follows from the exact solution of KE: E = E(M; e) is a monotonic increasing function of M and also is convex downwards. Then, by choosing appropriate upper and lower bounds of E = E(M; e) in M ∈ (0, π) we may derive admissible circular contours (16). Thus, taking as a lower bound that leads to the contour that is the contour used by Philcox et al. (2021). It contains the exact solution E(M; e) of (1) and by Proposition (2.1) is admissible for (2). Note that in (19) (20) that ρ 0 = e/2 is in fact independent of the mean anomaly M and, as we will see next, this fact has some computational advantages because it remains constant along the solution of the same elliptic orbit. Clearly, C 0 is a good choice for small eccentricities because the length of all contours is πe, but when e is close to 1 the length of the contours becomes close to π. Furthermore, when M → 0 and e → 1, the contour C 0 tends to the circle of center μ 0 = 1/2 and radius ρ 0 = 1/2 that is tangent to the imaginary axis at the origin, and recall that E = M = 0 with e = 1 is a singular point of Kepler's equation, and this fact introduces additional difficulties in the numerical calculations. Now with the purpose to have admissible circular contours with radius ρ(e) < e/2, we first note that we may split the interval of mean anomaly (0, π) into two (or more) subintervals and then to obtain admissible circular contours in each subinterval such that the corresponding radius in each subinterval is smaller than the above ρ 0 = e/2.
and let E L , E R be their corresponding eccentric anomalies. We take as lower bound E − of E = E(M; e) in the interval (M L , M R ) the chord between the two end points of it, and as upper bound E + , the tangent at some point (M * , E * ) of E = E(M; e) that is parallel to (21). Since that defines E * and M * by and then the upper bound is Now, we take as radius ρ, Note that α and ρ are independent of M. For the center μ, there results that is an affine function of M.
In Fig. 2 we present for this value of the eccentricity (e = 0.5) the upper (E + ) and lower bounds (E − ) of the eccentric anomaly in the above two subintervals. Besides, on the right part of this figure, we plot the respective radii (ρ 1 , ρ 2 ) of the Jordan circles and we can see that both are smaller than the radius ρ 0 = 0.5 of the circle C 0 ; hence, the length of each new circle is smaller than the length of C 0 .
Another possibility to construct contours for integrals (2) with small lengths is to replace the circular contours of center μ and radius ρ by elliptic contours with the same center, semimajor axis ρ, and semi-minor axis ε ρ with some 0 < ε < 1. It follows from Proposition (2.1) that if a circular contour is suitable for (2) the corresponding elliptic contour defined in this way will be also suitable. This contour can be written in the parametric form as and now the integral expression (2) becomes with

Fig. 2 Left) Lower bounds E − (in blue) and upper bounds E + (in orange) of the eccentric anomaly E = E(M; e) (in red)
when M is in the subintervals I 1 = (0, π/2 − e) and I 2 = (π/2 − e, π) for e = 0.5, as given in Eqs. (21) and (22). Right The radius ρ of the Jordan circle in the interval I 1 (in red), in the interval I 2 (in blue), and ρ 0 = e for the circle C 0 is given in Eq. (20)

Approximation of the line integrals along circular and elliptic contours
For circular contours, the solution E = E(M; e) is given by Eqs. (17), (18) and G(θ ) is a 2π-periodic function of θ that has a Fourier expansion with coefficients G(n) given by therefore, the exact solution (17) can be written also in the form A straightforward approach (used by Philcox et al. (2021)) to approximate (29) is to compute the FFT of G(θ ) by using a standard FFT solver in which, for a given number N ( 1) of points, we get approximations G N (k) to G(k) for |k| ≤ N /2. With these approximations we have the approximate solution E N (M; e) of (29) The main advantage of this approach is its simplicity and reliability because it is enough to use well-established FFT solvers and N sufficiently big to get accurate approximations. However, FFT solvers provide all coefficients G N (k) with |k| ≤ N /2, but for Eq. (30) we only need the corresponding to k = −2 and k = −1. Because of this, we focus on the approximations to G(−2) and G(−1) and we will use the same quadrature rule as the one used in the FFT, namely the composite trapezoidal rule that, in spite of its second order, is a reliable quadrature rule and, more important, the discrete approximations G N (k) have a spectral convergence to G(k). In fact, it will be seen in the Appendix 1 that for all |k| ≤ N /2 there exist α = α(e) so that with some constant C k . In other words, the discrete Fourier coefficients G N (k) converge exponentially to their corresponding exact values G(k). Note that in a standard quadrature rule with order p this convergence behaves asymptotically as 1/(N p ), (N → ∞).
and this implies that the logarithm of the error behaves linearly with N , whereas in the standard case, the logarithm of the error behaves as log N .
To simplify the calculation of G(k), k = {−2, −1}, we observe that for the particular function G(θ ) given by (18) there are some symmetries that allow us to reduce the computational cost. In fact, taking into account that for k = {−2, −1} (33)

Re exp(i kθ)G(θ ) = Re exp(−i kθ)G(−θ) ,
In the composite trapezoidal rule, given an (even) positive integer K the integral is approximated by with θ j = j π/K , (j = 0, 1, . . . , K ) and then (33) is approximated by Note that (36) is equivalent to the approximation provided by the FFT with N = 2K and the symmetry of G reduces the computational cost by a factor of 1/2.
A second remark is that in the computation of if ρ is independent of the mean anomaly and μ = β 0 + β 1 M is an affine function of M (here β 0 and β 1 may depend on e), the evaluations of ρ exp(i θ j ), j = 0, 1, . . . , K , can be stored and used for all evaluations of the mean anomaly along the orbit, and then, the calculation of G(θ ) −1 only requires the computation of sin μ and cos μ.

Numerical experiments
We will compute for circular and elliptic contours the absolute errors in the numerical solution denoted by AErr(K ) = |E(M; e) − E K (M; e)|, for (e, M) ∈ D. Moreover, in the case of elliptic contours we will denote by AErr(K , ε) the corresponding errors for ε ∈ (0, 1). We made many experiments to show the accuracy of our method, but here we present only a few of them, although the approach performs properly in all examples chosen. In all samples we take values of M ∈ (0, π) and a fixed value of the eccentricity e = 0.9, a higher value that may lead to convergence problems in some iterative methods for M close to zero. In the grid to apply the composite trapezoidal rule, we use three different values of K , namely, K = 8 (Fig. 3, top), K = 16 (Fig. 3, center), and K = 32 (Fig. 3, bottom).
For each value of K , we present the log 10 AErr(K ) corresponding to the seven contour Jordan used, from top to bottom on each plot, the circle C 0 , and ellipses with ε = 1/2, 1/4, 1/8, 1/16, 1/32 and a very eccentric ellipse with ε = 10 −3 ), close to a straight segment.
We can see in Fig. 3 that the smaller the length of the Jordan curve is the more accuracy we get. This pattern is the same for every number K of points in the grid. However, the accuracy dramatically increases with K ; thus, for K = 32, we reach from 20 significant digits for M ≈ 0 to 40 digits at the end of the interval (M ≈ π). Anyway, results for K = 8 are also very good; we reach from 10 to 20 significant digits.
We also note that there are several discontinuities in the plots. These are due to the fact that the numerical solution E K (M; e) oscillates around the exact solution E(M; e), and then, the discontinuities of log 10 AErr(K ) appear at the crossing points. Observe that the number of discontinuities increases with K , showing a similar behavior to the best approximation of continuous functions.
Another important issue is how fast is our method when compared with the one of Philcox et al. (2021). The answer can be found in Figure 4, where we represent the CPU time for three cases depending on the number of nodes taken in the grid. We consider three different methods, all of them sharing the same Jordan circle C used by Philcox et al. (2021), but we make the numerical computation of the line integrals (17) in three different ways: first (in gray), by using the FFT as in ; second (in red), by means of the composite trapezoidal rule; and finally (in red), with the optimized trapezoidal rule, that is, taking into consideration the symmetry and periodicity properties of the involved functions. We can see that the composite trapezoidal rule gives for even number of nodes a similar result in speed to the FFT method, the latter is slower for odd nodes. Clearly, optimized trapezoidal rule is the fastest, and the slope of its growing is more smooth than in the former cases.

Conclusions
The integral expression of the transcendental elliptic Kepler's Equation (KE) as quotient of two contours integrals around an appropriate Jordan curve is not only a mathematically Fig. 3 In log 10 scale, the absolute errors for eccentricity e = 0.9 and M ∈ (0, π) with K = 8, 16, 32 (from top to bottom), and for each one computed using seven elliptical contour C 0 , ranging from a circle (ε = 1) to an almost straight segment (ε = 10 −3 ) beautiful explicit expression of the solution but also can be used for practical purposes. The aim of this paper is to examine the main issues (accuracy and computational cost) that determines their practical application. Since the definition of the contour integrals depends on a Jordan curve that contains inside only the unique real zero of the KE, a complete study of the complex zeroes of this equation is carried out in Sect. 2.
Moreover, as Jordan contours with smaller lengths improve the accuracy of the numerical approximation of integrals along them, several circular and elliptic Jordan contours have been proposed. Such contours allow us to express the integrals as the Fourier coefficients of a 2π-periodic function for the frequencies k = −2 and k = −1. In this context, taking into Fig. 4 Comparison of the CPU time depending on the number of points in the grid, and for three used methods: in gray, the corresponding to the FFT, in red for the composite trapezoidal rule, and in green for the optimized trapezoidal rule, which is much faster than the other two account that the composite trapezoidal rule is used in the well-experimented FFT, we use this rule with some additional simplifications due to the symmetries of the integrals proposed to approximate these Fourier coefficients with great accuracy. The spectral accuracy with the number of nodes is obtained in Appendix A, and the results of a number of numerical experiments are presented to show the accuracy as a function of the number of nodes in the quadrature rule. Special attention has been paid to the accuracy for values near the singularity M = 0, e = 1 of KE.
Besides, it is worth to remark that as shown in Fig. 3, the use of ellipses as Jordan contours with parameter ε close to zero leads to a great improvement in the accuracy of the numerical approximation of the curvilinear integrals with the same computational cost as in the case of circular contours. Then, elliptical contours are highly recommended in computing the integrals of KE, particularly in the case of orbits with high eccentricities and for values of M close to zero.
Finally, the above study on the explicit integral solutions of elliptic Kepler's equation can be extended to other similar transcendental equations. In particular, for the limit Kepler's equation with e = 1, as noted in Remark 2, the study of its complex zeroes allows us to consider also elliptic Jordan contours for the efficient solution of the corresponding integrals. Moreover, the case of the hyperbolic Kepler's equation e sinh F − F = M (e > 1, M > 0) in the unknown F will be subject of a forthcoming paper.
Funding Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Appendix: A On the convergence of G h (k) to G(k)
An interesting study of the convergence of the composite trapezoidal rule was carried out by Johnson (2021), and a first pessimistic approximation leads to however, in our case G(θ ) is periodic and this is not an accurate estimate.
A more refined tool to analyze the convergence of G h (k) to G(k) is the Euler-Maclaurin formula (see, e.g., Hairer et al. (1992)) that with our notations becomes where B j are the Bernoulli numbers defined as the coefficients of the series expansion with their first values For the periodic function (θ) = exp(2i θ)G(θ ) with period 2π, it turns out that all derivatives satisfy ( j) (0) = ( j) (2π), and therefore for all N and, consequently, we have arbitrarily large order of convergence when N → ∞. Similar result holds for = exp(i θ)G(θ ). However, as remarked by Johnson (2021), Fourier analysis will allow us to give sharper bounds on the asymptotic convergence. To do that, we will use the statements (b) and (c) of Theorem 3, page 33 of Trefethen's book (2000) (also known as Paley Wiener theorem) that in our notation can be formulated as for every m ≥ 0. (c) If there exist positive constants a and c such that G(θ ) can be extended to an analytic function u(θ ) in the complex θ set with |u(θ )| ≤ c for all θ in Eq. (40), then for every > 0.
Under the conditions of the last statement for h = 2π/N , the discrete wavenumbers G h (k) approximate the corresponding continuous wavenumber G(k) with an error of asymptotic order and then we say that they exhibit spectral approximation with factor a/2. In the application of the statement (c) to G h (−2) and G h (−1) observe that G(θ ) = G(θ ; M, e) = f (z = μ + ρ exp(iθ)) −1 depends on the parameters (e, M) ∈ D through the center μ and the radius ρ of the circular contour C(μ, ρ). This implies that the constant a of (41) In the application of the statement (b) to the (2π)-periodic function S(θ ), we get that for every m > 0, or equivalently, taking into account that N = 2π/h and, in particular, the exact integrals in the numerator and denominator of (43) are approximated by their DFT counterparts with the asymptotic accuracy To apply the stronger result of Theorem A.1 (c), we must examine the analytic behavior of S(θ ) (42) for θ ∈ C in a set Re θ ∈ [0, 2π], and |Im θ | < a.
According to the study of Sect. 2, when M → 0 and e → 1, there is a branch of complex zeroes of G(θ ) −1 that tends to θ = 0, and similarly, when M → 2π and e → 1, that branch tends to θ = 2π. More precisely, for each e 0 ∈ (0, 1) there exist y * = y * (e 0 ) > 0 that is the unique solution of g(y) = y sinh y = e 0 so that G(θ ) −1 has no zeros for |Im (θ )| < y * and therefore G(θ ) is analytic in the set (47) with a = y * . Some particular values are e 0 0.3 0.5 0.9 y * 2.99734 2.17732 0.803436 Hence, by choosing a < y * (e 0 ), |G(θ ; M, e 0 )| will be uniformly bounded in (47) for all M ∈ [0, π] and the statement (c) implies the spectral accuracy of the discrete wavenumbers with error O (exp(−a N)) with N → ∞.
In conclusion, for a given eccentricity e 0 ∈ (0, 1) when the exact value of the numerator of (43), that in Fourier notation is S(−2), is approximated by the composite trapezoidal rule with (N + 1) nodes, i.e., S h (−2) in terms of the discrete Fourier transform, with h = 2π/N , the error satisfies an asymptotic error estimate with a coefficient a = y * (e 0 ) = g −1 (e 0 ) which depends only on e 0 and does not depend on M; therefore, with some constant C 2 . Similarly for the denominator of (43), we have  (49) and (50) show that the composite trapezoidal rule applied to the corresponding integrals exhibits an spectral accuracy of both numerator and denominator with the same constant factor a < y * (e 0 ).
Moreover, for the quotient, by putting = exp(−a/2 N ), by (49)  hence, for the quotient we have also spectral accuracy with the same coefficient. Note that such spectral accuracy could be improved if the equality S(−1)C 1 − S(−2)C 2 = 0 holds. In Fig. 5 we display the log 10 AErr(K ) for the mean anomaly M = 0.1 (upper plot) and M = 3.0 (lower plot), as a function of K (N = 2K ), with 2 ≤ K ≤ 32 and for a discrete range of values of the eccentricity between e = 0.1 and e = 0.99.
We see that log 10 AErr(K ) exhibits a linear behavior of K with a negative slope −α = −α(e) < 0 that increases with the eccentricity and then log 10 [AErr(K )] = log(AErr(K )) log 10 ∼ −α(e) K which implies AErr(K ) ∼ exp[−α(e)K log 10], which shows that AErr(K ) behaves as an exponential function of K with a negative coefficient that depends on the eccentricity so that it is smaller when the eccentricity gets closer to one. In Figure 5 (upper plot) we observe that for the mean anomaly M = 0.1 we have the same behavior as in the case M = 3.0 (lower plot), but as M is close to zero and (M = 0, e = 1) is a singularity of KE, the spectral behavior is affected by this fact and α(e) is closer to zero than in the above case.
This confirms numerically the spectral behavior of the error that depends on the eccentricity. Note that for the small eccentricity e = 0.1 the negative slope implies that arbitrarily small errors a attained even with moderate values of K . However, for the eccentricity e = 0.99 it is necessary to take large values of K to get moderate errors. Moreover, this is particularly relevant for values of M close to zero as shows the case M = 0.1.