Fast Switch and Spline Scheme for Accurate Inversion of Nonlinear Functions: The New First Choice Solution to Kepler's Equation

TL;DR

This paper introduces a fast, spline-based method for inverting monotonic functions, demonstrating superior speed and accuracy over traditional Newton-Raphson methods in solving equations like Kepler's and Lambert's W.

Contribution

The authors develop a switching and spline interpolation scheme with an ultra-fast routine, providing highly accurate inverse function computations with significantly improved speed.

Findings

The method achieves errors smaller than traditional spline analysis by up to 10^{-22}.

It is 10^{-4} to 10^{-3} times faster than Newton-Raphson in tested cases.

Outperforms Newton's method for Kepler's equation at high precision and large point sets.

Abstract

Numerically obtaining the inverse of a function is a common task for many scientific problems, often solved using a Newton iteration method. Here we describe an alternative scheme, based on switching variables followed by spline interpolation, which can be applied to monotonic functions under very general conditions. To optimize the algorithm, we designed a specific ultra-fast spline routine. We also derive analytically the theoretical errors of the method and test it on examples that are of interest in physics. In particular, we compute the real branch of Lambert's $W(y)$ function, which is defined as the inverse of $x \exp(x)$, and we solve Kepler's equation. In all cases, our predictions for the theoretical errors are in excellent agreement with our numerical results, and are smaller than what could be expected from the general error analysis of spline interpolation by many orders of magnitude, namely by an astonishing $3\times 10^{-22}$ factor for the computation of $W$ in the range $W(y)\in [0,10]$, and by a factor $2\times 10^{-4}$ for Kepler's problem. In our tests, this scheme is much faster than Newton-Raphson method, by a factor in the range $10^{-4}$ to $10^{-3}$ for the execution time in the examples, when the values of the inverse function over an entire interval or for a large number of points are requested. For Kepler's equation and tolerance $10^{-6}$ rad, the algorithm outperforms Newton's method for all values of the number of points $N\ge 2$.

Figures

Equations

$$S_j(y)=\sum _{q=0}^3{c}_{q_j}{\left(y-y_{0_j}\right)}^q,$$

$$d_j\equiv g^{\prime }(y_j)=\frac{1}{f^{\prime }(x_j)},\text{for }j=0,\cdots ,n.$$

$$S_j(y_{0_j})=x_{0_j},S_j(y_{1_j})=x_{1_j},S_j^{\prime }(y_{0_j})=d_{0_j},S_j^{\prime }(y_{1_j})=d_{1_j},$$

$${\mathbf{c}}_0={\mathbf{x}}_0,$$

$${\mathbf{c}}_1={\mathbf{d}}_0,$$

$${\mathbf{c}}_2=\frac{(2\ast {\mathbf{d}}_0+{\mathbf{d}}_1)\ast ({\mathbf{y}}_0-{\mathbf{y}}_1)-3\ast ({\mathbf{x}}_0-{\mathbf{x}}_1)}{({\mathbf{y}}_0-{\mathbf{y}}_1{)}^2},$$

$${\mathbf{c}}_3=\frac{({\mathbf{d}}_0+{\mathbf{d}}_1)\ast ({\mathbf{y}}_0-{\mathbf{y}}_1)-2\ast ({\mathbf{x}}_0-{\mathbf{x}}_1)}{({\mathbf{y}}_0-{\mathbf{y}}_1{)}^3}.$$

$$|g(y)-S(y)|\le \frac{1}{384}M\mu ,$$

$$M=\underset{x_0\le x\le x_n}{\max }\left|-\frac{15f^{\prime \prime }(x{)}^3}{f^{\prime }(x{)}^7}+\frac{10f^{(3)}(x)f^{\prime \prime }(x)}{f^{\prime }(x{)}^6}-\frac{f^{(4)}(x)}{f^{\prime }(x{)}^5}\right|.$$

$$|g(y)-S(y)|=\left|\sum _{q=4}^{\infty }\frac{g^{(q)}(y_j)}{n!}(y-y_j{)}^n\right|=\left|\frac{g^{(4)}({\bar{y}}_j)}{4!}(y-y_j{)}^4\right|\approx \left|\frac{g^{(4)}(y_j)}{4!}(y-y_j{)}^4\right|,$$

$$|g(y)-S(y)|\le \underset{y_{\min }\le \bar{y}\le y_{\max }}{\max }\left|\frac{g^{(4)}(\bar{y})}{4!}\right|\underset{0\le j\le n-1}{\max }{\left|\frac{y_{j+1}-y_j}{2}\right|}^4=\frac{1}{384}M\mu ,$$

$$|g(y)-S(y)|\approx \frac{1}{4!}\left|\left[-\frac{15f^{\prime \prime }(x_j{)}^3}{f^{\prime }(x_j{)}^7}+\frac{10f^{(3)}(x_j)f^{\prime \prime }(x_j)}{f^{\prime }(x_j{)}^6}-\frac{f^{(4)}(x_j)}{f^{\prime }(x_j{)}^5}\right]{\left[f^{\prime }(x_j)(x-x_j)\right]}^4\right|,$$

$$\underset{0\le j\le n-1}{\max }|x_{j+1}-x_j|\underset{x_0\le x\le x_n}{\max }\left|\frac{f^{\prime \prime }(x)}{2f^{\prime }(x)}\right|\ll 1$$

$$|g(y)-S(y)|⪅$$

$$\underset{x_0\le x\le x_n}{\max }\left|-\frac{15f^{\prime \prime }(x{)}^3}{f^{\prime }(x{)}^3}+\frac{10f^{(3)}(x)f^{\prime \prime }(x)}{f^{\prime }(x{)}^2}-\frac{f^{(4)}(x)}{f^{\prime }(x)}\right|.$$

$$M=\underset{x_0\le x\le x_n}{\max }6e^{-4x}=6e^{-4x_0},$$

$$\mu ={\left(\frac{x_n-x_0}{n}\right)}^4\underset{x_0\le x\le x_n}{\max }{e}^{4x}={\left(\frac{x_n-x_0}{n}\right)}^4{e}^{4x_n},$$

$$|g(y)-S(y)|\le \frac{6}{384}{\left(\frac{x_n-x_0}{n}\right)}^4{e}^{4(x_n-x_0)}.$$

$$|g(y)-S(y)|⪅\frac{6}{384}{\left(\frac{x_n-x_0}{n}\right)}^4.$$

$$\frac{x_n-x_0}{n}\left|\frac{f^{\prime \prime }(x)}{2f^{\prime }(x)}\right|=\frac{x_n-x_0}{2n}\ll 1.$$

$$M=\underset{x_0\le x\le x_n}{\max }\left|-\frac{15{\left(e^{x}x+2e^x\right)}^3}{{\left(e^{x}x+e^x\right)}^7}+\frac{10\left(e^{x}x+3e^x\right)\left(e^{x}x+2e^x\right)}{{\left(e^{x}x+e^x\right)}^6}-\frac{e^{x}x+4e^x}{{\left(e^{x}x+e^x\right)}^5}\right|,$$

$$\mu ={\left(\frac{x_n-x_0}{n}\right)}^4\underset{x_0\le x\le x_n}{\max }{\left|e^{x}x+e^x\right|}^4,$$

$$|g(y)-S(y)|\le \frac{1}{384}\mu (x=x_n)M(x=x_0).$$

$$|g(y)-S(y)|⪅\frac{1}{384}{\left(\frac{x_n-x_0}{n}\right)}^4\times$$

$$\underset{x_0\le x\le x_n}{\max }\left|-\frac{15{\left(e^{x}x+2e^x\right)}^3}{{\left(e^{x}x+e^x\right)}^3}+\frac{10\left(e^{x}x+3e^x\right)\left(e^{x}x+2e^x\right)}{{\left(e^{x}x+e^x\right)}^2}-\frac{e^{x}x+4e^x}{e^{x}x+e^x}\right|.$$

$$n\gg \frac{x_0-x_n}{2}\underset{x_0\le x\le x_n}{\max }\left|\frac{x+2}{x+1}\right|=10,$$

$$y=x-\text{e}\sin x,$$

$$\theta =2\arctan \left(\sqrt{\frac{1+\text{e}}{1-\text{e}}}\tan \frac{x}{2}\right).$$

$$M=\underset{0\le x\le \pi }{\max }\left|-\frac{15{\text{e}}^3{\sin }^{3}x}{(1-\text{e}\cos x{)}^7}+\frac{10{\text{e}}^2\sin x\cos x}{(1-\text{e}\cos x{)}^6}+\frac{\text{e}\sin x}{(1-\text{e}\cos x{)}^5}\right|$$

$$\mu ={\left(\frac{\pi }{n}\right)}^4\underset{0\le x\le \pi }{\max }{\left|1-\text{e}\cos x\right|}^4.$$

$$|g(y)-S(y)|\lesssim \frac{2.7\times 10^4}{n^4}.$$

$$|g(y)-S(y)|⪅\frac{1}{384}{\left(\frac{\pi }{n}\right)}^4\underset{0\le x\le \pi }{\max }\left|-\frac{15{\text{e}}^3{\sin }^{3}x}{(1-\text{e}\cos x{)}^3}+\frac{10{\text{e}}^2\sin x\cos x}{(1-\text{e}\cos x{)}^2}+\frac{\text{e}\sin x}{1-\text{e}\cos x}\right|.$$

$$|g(y)-S(y)|⪅\frac{5.5}{n^4}.$$

Full text

Algorithm for Inverse Function; Kepler Equation for Orbital Motion; Astrodynamics; Cubic Spline Interpolation; Newton-Raphson Iteration Method; Celestial Mechanics

I Introduction

Many problems in science and technology require the inversion of a known nonlinear function $f(x)$. Widely studied examples include the inversion of elliptic integrals Fukushima2013 ; Boyd2015 , the computation of Lambert W function Corless1996 ; Veberic2012 , and the solution of Kepler’s equation for the orbital motion of a body in a gravitational field Prussing2012 ; Curtis2014 .

In many cases, the inverse function cannot be found analytically, and numerical methods must be used. Besides possible special procedures that may be found for specific forms of $f(x)$, the most popular numerical inversion schemes are those based on the Newton-Raphson method for computing the zeros of a function Mathews1999 or some of its variants Danby1983 ; Danby1987 ; Gerlach1994 ; Palacios2002 . These schemes are largely universal, i.e. they can be applied to a wide class of functions and converge very rapidly, especially when the value of the inverse function at one given point is required, rather than on an entire interval. However, they require a reasonably good first guess in order to avoid problems of convergence, which may be a nontrivial issue in some cases, such as in Kepler’s problem for values of the eccentricity close to one Conway1986 ; Charles1998 ; Stumpf1999 .

The rationale behind using the Newton-Raphson method is based on the fact that solving the equation $y=f(x)$ for $x$ when the value of $y$ is given is equivalent to the problem of finding the zeros of the functions $F_{y}(x)\equiv f(x)-y$. If a good initial guess $x_{0}$ of the true value $x=f^{-1}(y)$ of the zero is available, the zeros of $F_{y}$ can be computed by recursively applying the equation $x_{k+1}=x_{k}-\frac{F_{y}(x_{k})}{F_{y}^{\prime}(x_{k})},$ i.e. ${\displaystyle x_{k+1}=x_{k}+{\frac{y-f(x_{k})}{f^{\prime}(x_{k})}}}$.

Here, a Fast Switch and Spline Inversion (FSSI) scheme is described that does not require an initial guess and can be applied under very general conditions provided the function $f$ is one-to-one. The basic idea underlying this method is remarkably simple, yet it can be turned into a very powerful and accurate tool, as shall be demonstrated. Surprisingly, to our knowledge, this scheme has not been explored in the published literature. Perhaps, this may be due to an underappreciation of its rate of convergence, given the known bounds on the precision of spline interpolation, and to the existence of standard alternatives such as Newton’s method.

After describing the FSSI method, we derive theoretically a set of analytical expressions of its error estimates, and show that they are much smaller than the limit that could be derived by merely applying the existing spline analysis to this case. To optimize the algorithm, we also designed a specific spline routine that makes the FSSI more accurate and much faster than using the known spline routines. We then test the scheme on several nonlinear functions, and demonstrate that in all cases our theoretical predictions for the errors are in excellent agreement with the numerical computations.

Based upon this error analysis and on the numerical computations, the FSSI is shown to be a valid alternative to the Newton-Raphson method (and similar quasi-Newton minimization methods) for computing values of inverse functions, especially if a good first approximation is difficult to obtain. Moreover, the FSSI is shown to be superior to Newton-Raphson when the values $f^{-1}(y)$ of the inverse function are required for many different $y$ points, or over an entire interval. In the case of Kepler’s equation for orbital motion, FSSI is faster than Newton and quasi-Newton methods when the position of the orbiting body must be known at more than a few different instants, depending on the eccentricity e and the requested precision. For example, for $\text{e}=0.8$ and accuracy $\sim 10^{-6}$ rad, the FSSI algorithm is already faster than Newton’s even when the computation is done at $N=2$ points, and $\sim 2000$ times faster for large $N$.

II The Fast Switch and Spline Inversion (FSSI) Scheme

In what follows the FSSI method is described. Let $f(x)$ be the input function, which is presumed to be single valued (monotonic) in a given domain $x\in\left[x_{\mathrm{min}},x_{\mathrm{max}}\right]$. The function $f(x)$ is assumed to be given analytically, but the case when it is known at discrete points shall also be considered. The goal of the method is to obtain a numerical approximation for the inverse function $g(y)=f^{-1}(y)$ in the co-domain.

The FSSI consists of a two step approach. First, when the input function is given analytically, the values of $f$ on a given grid of points $x_{j}$, for $j=1,\cdots,n$, are computed to obtain the matrix $(x_{j},y_{j})$, where $y_{j}=f(x_{j})$. The matrix $(y_{j},x_{j})$, obtained by switching the arrays, gives then the exact values $g(y_{j})=x_{j}$ of the inverse function on the grid $y_{j}$. From this modified matrix, the cubic spline interpolant $S(y)$ of $(y_{j},g(y_{j}))$ is computed by using a special routine that is designed in the next section. The resultant function $S(y)$ is the approximation of the inverse function.

FSSI can also be used when the function $f(x)$ is specified on a grid by a set of tuples $(x_{j},y_{j})$. In a high-level programming language such as Python, this tuple array is represented as $(x,y)$, and the FSSI interpolant $S(y)$ can be obtained by calling a cubic spline routine of the switched tuple array, $S=CubicSpline(y,x)$. In this way, the object S in a high-level computer language would act as a generator for points in the co-domain of $f$, giving the inverse $f^{-1}$.

Figure 1 shows the flow of the FSSI algorithm, that could be implemented in any high-level computer language. The central dotted box shows the two-step procedure of FSSI, while the outer boxes show possible interfacing between input and output. In particular, the input interface could accept two types of data: (1) a pointer to the analytic function $f(\cdot)$ and its derivative $f^{\prime}(\cdot)$, together with a grid of $n+1$ points $x$, or (2) the discrete tuple arrays $(x,y)$, for the case when the function or its derivative are not known explicitly. At the output, the procedure returns a generating function, whose precision as an approximation of $f^{-1}$ is determined by the number of points $n+1$ of the input grid, and which is used for sampling $N$ points $\{Y_{1},\cdots,Y_{N}\}$ in the function’s co-domain. In subsequent sections, we show that the FSSI algorithm has an error bound proportional to $|k\Delta x|^{4}$, for k a constant (described in the text), and with a maximal time complexity of $O(1)+\epsilon\,O(N)$, where, beyond some value of $N$, the second term dominates.

III Design of a specific ultra-fast spline for the FSSI scheme

In many problems, the function $f$ to be inverted is known analytically, along with its derivative $f^{\prime}$. This is the case for Kepler’s equation and for all the other examples that we will consider hereafter. Therefore, we can profit from the knowledge of $f^{\prime}$ to design a specific cubic spline $S(y)$ interpolant for the FSSI algorithm.

The resulting spline makes the FSSI algorithm more accurate and much faster than calling the spline routines that are currently available, which do not make use of the derivatives of the input function. This huge difference in speed is due to the fact that most spline routines require the numerical solution of a system of $4n$ coupled equations to compute the $4n$ coefficients of the spline Mathews1999 , where $n$ is the number of grid intervals. An exception is Akima’s cubic spline Akima ; ScipyAkima , which is fast, diagonal and regular.

The specific spline that is designed here is based on a similar idea to Akima’s, but it is significantly more accurate than the latter, usually by three orders of magnitude in the examples that we shall consider, and it is also faster. Its superior performance is due to the fact that it uses the derivative $f^{\prime}$ as an input, unlike Akima’s. Of course, the usual applications of splines are not meant for cases in which the function to be interpolated and its derivative are given analytically. However, the situation is completely different in the problem of the inversion of a function $f(x)$. In this case, the derivatives $g^{\prime}(y_{j})=1/f^{\prime}(x_{j})$ can be given on a grid, while the values $g(y)$ are not known.

Let us build the specific spline $S(y)$ piecewise in each interval, $S(y)=S_{j}(y)$ for $y_{j}<y<y_{j+1}$, where $j$ takes the values $j=0,\cdots,n-1$. If we define the arrays $\mathbf{y_{0}}\equiv(y_{0},\cdots,y_{n-1})$ and $\mathbf{y_{1}}\equiv(y_{1},\cdots,y_{n})$, obtained by removing the last and the first point of the $\mathbf{y}$ array, respectively, then the $j$-th interval can also be written as $y_{0_{j}}<y<y_{1_{j}}$. In this interval, the cubic spline can be defined as

[TABLE]

where for each value of $q=1,\cdots,4$ the coefficients $c_{q_{j}}$ can also be thought as the $n$ components of an array $\mathbf{c_{q}}$.

Since the values of the derivative $f^{\prime}(x_{j})$ of the input function are known on the grid points, we can construct an array $\mathbf{d}$ whose $n$ components are the derivatives of the inverse function $g$ on the points $y_{j}=f(x_{j})$,

[TABLE]

As we did for $\mathbf{y}$, which was used to generate the arrays $\mathbf{y_{1}}$ and $\mathbf{y_{2}}$ by removing one end point, it is convenient to define similar arrays of $n$ components also from $\mathbf{x}$ and $\mathbf{d}$, namely $\mathbf{x_{0}}\equiv(x_{0},\cdots,x_{n-1})$, $\mathbf{x_{1}}\equiv(x_{1},\cdots,x_{n})$, $\mathbf{d_{0}}\equiv(d_{0},\cdots,d_{n-1})$, and $\mathbf{d_{1}}\equiv(d_{1},\cdots,d_{n})$. With this convention, we have to choose the spline coefficients that lead to the best approximation of the inverse function $g(y)$. The most natural choice is to force $S_{j}$ to coincide with the known values of the inverse function, $x_{0_{j}}$ and $x_{1_{j}}$, at the end points, and to ask the same for the derivatives $d_{0_{j}}$ and $d_{1_{j}}$. In other words, the conditions to be imposed in each interval are,

[TABLE]

where $S_{j}(y)$ is given by equation 1.

For every fixed value of $j$, these conditions give a system of four equations for the four unknown coefficients $c_{0_{j}}$, $c_{1_{j}}$, $c_{2_{j}}$ and $c_{3_{j}}$, which is decoupled from the similar systems of equations corresponding to different values of $j$. As we have mentioned above, this is an important advantage, as compared with most of the other cubic spline interpolation methods, which must solve systems of $4n$ coupled equations to compute the coefficients Mathews1999 , with the exception of Akima’s. We can expect that this will make the FSSI algorithm using this spline much faster than using the alternative ones, and this is also what we have found numerically.

In fact, the system of equations (3) for a fixed $j$ can be solved analytically in a straightforward way, and then implemented numerically in a completely diagonal form. In order to write the solution in a compact form, we will use a convention for vector arrays that is common in computer languages like python: an equation for arrays is interpreted in terms of components in such a way that an equality like, e.g., $\mathbf{u}=\frac{\mathbf{v}*\mathbf{w}+2*\mathbf{z}}{\mathbf{s}}$ between vectors having the same number of elements represents the equations ${u}_{j}=\frac{{v}_{j}{w}_{j}+2{z}_{j}}{s_{j}}$ for every $j$. With this convention, equations (3) give the solution

[TABLE]

This result can be used when the function $f(x)$ and his derivative are known on the whole interval in which the inversion is required. If the second derivative is also known, by adding the additional conditions $S_{j}^{\prime\prime}(y_{0_{j}})=d_{0_{j}}$ and $S_{j}^{\prime\prime}(y_{1_{j}})=d_{1_{j}}$ to equation (3) we can also design a diagonal quintic spline, and if also the third derivatives are known the additional conditions $S_{j}^{\prime\prime\prime}(y_{0_{j}})=d_{0_{j}}$ and $S_{j}^{\prime\prime\prime}(y_{1_{j}})=d_{1_{j}}$ allow for the construction of a septic spline. We have done this in both cases for the FSSI algorithm, and checked in the examples that, for a given accuracy, the resulting versions of the method perform slightly worse than with the cubic spline as designed above. Therefore, the latter will be taken as the optimal specific spline for FSSI.

IV Computation of the theoretical error

In this section, the predicted theoretical error of the FSSI is developed for the case of an input function $f$ that is continuous and having continuous derivatives up to at least the fifth degree.

Note that this error analysis not only works for the cubic spline that we have designed in the previous section, but also holds when the FSSI method is implemented with most known cubic spline routines. The main differences between the use of a cubic spline routine or the other are the speed and the accuracy very close to the end points of the $y$ domain. In both these aspects, the FSSI performs better with the spline of section III than with the others.

IV.1 Derivation of an upper bound on the error of the FSSI by using the known analysis of cubic spline interpolation

Following Ref. Sonneveld1969 ; Mathews1999 , we can compute an upper bound for the error of the cubic spline $S(y)$, used to interpolate the function $g(y)$, from the formula

[TABLE]

where $M=\max\limits_{y_{\mathrm{min}}\leq y\leq y_{\mathrm{max}}}\left|g^{(4)}(y)\right|$, and $\mu=\max\limits_{0\leq j\leq n-1}(y_{j+1}-y_{j})^{4}$.

In our case, $g(y)$ is the inverse of the input function $f(x)$, therefore it is convenient to express this error in terms of $f(x)$ and its derivatives. The $M$ term becomes,

[TABLE]

Equations (5) and (6), along with the definition of $\mu$, can be used to obtain an upper limit on the error. As shown in examples below, the actual errors are several orders of magnitude smaller than this upper bound. In other words, the method converges much more rapidly than expected. Therefore, it is of great interest to obtain a more accurate, albeit approximate, analytical estimate of the error of FSSI , and check its consistency in examples. This is done in the next subsection.

IV.2 Ab initio derivation of an improved estimation of the error for the FSSI

Let us assume that the function $g(y)$ is infinitely differentiable. Thus, it can be expanded in a Taylor series $g(y)=\sum_{q=0}^{\infty}\frac{g^{(q)}(y_{j})}{q!}(y-y_{j})^{n}$ around one of the points of the grid $y_{j}=f(x_{j})$, chosen to be the closest grid point to $y$, so that $|y-y_{j}|\leq|y-y_{j+1}|$ and $|y-y_{j}|\leq|y-y_{j-1}|$.

We also assume that the cubic spline interpolation is made in such a way so that it is equivalent to the first terms of this Taylor expansion, expanded around the same point and truncated beyond cubic order, $S(y)=\sum_{q=0}^{3}\frac{g^{(q)}(y_{j})}{q!}(y-y_{j})^{q}$ . Therefore, the difference is

[TABLE]

where $\bar{y}_{j}$ is an unknown intermediate point between $y$ and $y_{j}$, and the last approximation is expected to hold for sufficiently small values of $|y-y_{j}|$, which is the case when a sufficiently high number of grid points is chosen.

On one hand, the exact equality in equation (7) can be translated in the following bound,

[TABLE]

which does not use the information that $g$ is the inverse function of $f$, and coincides with the limit of equations (5) and (6). The factor $2$ dividing the interval $y_{j+1}-y_{j}$ is due to the fact that $y_{j}$ was chosen as the closest grid point to $y$, so that $\left|{y-y_{j}}\right|\leq\left|\frac{y_{j+1}-y_{j}}{2}\right|$.

On the other hand, equation (7) can be elaborated further and expressed in terms of the function $f(x)$,

[TABLE]

where $x=g(y)$ and $y-y_{j}\approx f^{\prime}(x_{j})(x-x_{j})$. The last approximation is expected to be accurate over the entire interval provided the following condition

[TABLE]

is satisfied.

Assuming that the point $x_{j}$ is the one closest to $x$, so that $\left|{x-x_{j}}\right|\leq\max\limits_{0\leq j\leq n-1}\left|\frac{x_{j+1}-x_{j}}{2}\right|$, the following estimation is obtained

[TABLE]

As described in the next section, this error approximation is usually much smaller, possibly by many orders of magnitude, than the limit of equations (5) and (6) that we derived using existing literature on cubic spline interpolation. In fact, equation (10) can be expected to be a good approximation if the number of grid points $n$ is large enough and it satisfies the condition (9). The examples in the next section show that this is indeed the case.

Finally, we note that, for an equally spaced $x_{j}$ grid, we can substitute $\max\limits_{0\leq j\leq n-1}\left|x_{j+1}-x_{j}\right|=\frac{x_{n}-x_{0}}{n}$ in the condition (9), and $\max\limits_{0\leq j\leq n-1}\left|x_{j+1}-x_{j}\right|^{4}=\frac{(x_{n}-x_{0})^{4}}{n^{4}}$ in the estimation of the error (10). Not only is this the simplest choice, if the grid is not given otherwise, but it is usually also the best option. In fact, when $\left|x_{j+1}-x_{j}\right|^{4}=\frac{(x_{n}-x_{0})^{4}}{n^{4}}$ for every $j$, the grid values of the inverse function, where $g$ is known exactly, are also equally spaced. As a result, the error is distributed uniformly across the entire interval (as seen in the examples below).

V Examples

Here, the FSSI is applied to examples of interest involving nonlinear functions $f(x)$, defined over a domain $x_{0}\leq x\leq x_{n}$. For the numerical computation, we developed a python code that implements the FSSI as well as other Newton-based function inverse solvers.

Once the numerical interpolation $S(y)$ of the inverse function is obtained, the numerical errors are computed by evaluating $S(y)-S(f(S(y)))$. In the first example, in which the exact inverse function $g(y)$ is known, we also provide an additional evaluation of the error by computing the difference $S(y)-g(y)$. In both cases, we use a grid $Y_{k}$ that contains 10 times as many points as the original grid $y_{j}=f(x_{j})$. In fact, by construction $S(y)$ is exactly equal to $g(y)$–within machine errors–over the original grid $y_{j}$, so it is important to ensure that $S(y)$ is compared with $g(y)$ in between the grid points $y_{j}$. It is true that even if they were equal in number the points $Y_{k}$, chosen to be equally spaced, would not coincide in general with the $y_{j}$, whose spacing is variable and roughly proportional to $f^{\prime}(g(y))$; however, the election of ten times more points is more conservative. In this way, if the number of original grid points is large enough, a reasonable evaluation of the error is guaranteed. In fact, the examples show that when the exact analytic inverse function $g(y)$ is known, the difference $S(y)-g(y)$ has the same behaviour and magnitude of oscillations as $S(y)-S(f(S(y)))$. Moreover, the estimates are also in excellent agreement with the theoretical predictions for the error from equation (10).

In all the examples, we use an equally spaced input grid in $x$, which is expected to be the best choice in most cases, as we have discussed in the previous section. The spline routine used to implement the FSSI scheme is the specific one we have designed in section III. However, we have also checked that similar results are obtained by calling other splines routines that do not take the derivatives of $f$ as an input, such as Scipy cubic spline routines deBoor1978 ; SciPyCubicSpline . The errors in the bulk of the $y$ interval with most of those routines are very similar to each other, except very close to the boundary points, where they can be larger by an order of magnitude than those obtained using our specific spline. An exception is Akima routine ScipyAkima , which is less accurate by three orders of magnitude in the bulk of the interval. This is an additional reason for preferring our specific spline, besides the fact that it is the fastest one.

V.1 Exponential

The first function considered is $f(x)=\exp(x)$. Of course, in this case the exact inverse function is known analytically, $g(y)=\ln(y)$, thereby serving as a validation check of our scheme.

For this case, the quantities $M$ and $\mu$, from the bounds expressions of (5) and (6), can be readily computed. The results are,

[TABLE]

and

[TABLE]

so that the bound on the error as computed from equations (5) and (6) is

[TABLE]

On the other hand, the analytic estimation we derived in Equation (10) gives

[TABLE]

Therefore, this analytic error approximation of FSSI is smaller by a factor $e^{-4(x_{n}-x_{0})}$ as compared to the limit that was derived in equation (13) by applying the standard cubic spline error bound. For example, if $x_{n}-x_{0}=10$, then our error estimation is a factor $\exp(-40)$, i.e. 17 orders of magnitude, smaller than what could be expected from the literature. In order to benefit by this accuracy improvement, the grid must be chosen in such a way that the condition (9) is satisfied, i.e.

[TABLE]

If this condition on the number of grid points $n$ is met, our improved estimation of the error (14) can be expected to be a good approximation. For instance, if $x\in[0,10]$, the condition becomes $n\gg 5$, so that values of $n$ of the order of 50 or larger could be sufficient. This is also what we have observed by performing numerical computations for different values of $n$. In general, for $n\gtrsim 50$, equation (14) gives a correct estimate for the error over the entire interval.

Figure 2 shows the result of the FSSI for the inversion of $f(x)=\exp(x)$ over the domain $x\in[0,10]$ using $n=10^{2}$ grid points. In this case, our theoretical prediction of Equation (14) gives $|g(y)-S(y)|\lessapprox 1.6\times 10^{-6}$, which is in excellent agreement with the numerical computation over the entire interval. The results for this case also confirm the theoretical prediction that FSSI is 17 orders of magnitude more accurate than what could be expected by naively applying the general results for cubic splines, as in equation (13). An important feature of Figure 2 is that the error is distributed uniformly across the interval. As discussed previously, this is a consequence of choosing an equally spaced grid for $x$.

V.2 Lambert W function

Let $f(x)=x\exp(x)$, whose inverse function $g(y)$ in the real domain is the principal branch of Lambert’s W function, $W(y)$ Corless1996 ; Veberic2012 . In this case, the FSSI interpolation $S(y)$ can be compared with the values of $W(y)$ that are computed with other methods. The values of $M$ and $\mu$ from equations (5) and (6) for the theoretical bound are,

[TABLE]

which is a monotonically decreasing function for $x>-1$, and

[TABLE]

which increases monotonically. Therefore, the bound (5) becomes

[TABLE]

For example, over the domain $x\in[0,10]$ this gives $|g(y)-S(y)|\lessapprox 5.7\times 10^{20}\left(\frac{x_{n}-x_{0}}{n}\right)^{4}$.

On the other hand, our analytical estimation from Equation (10) becomes

[TABLE]

The function to be maximized in equation (19) monotonically decreases for $x>-1$, so that its maximum is achieved for $x=x_{0}$. Thus, over the domain $x\in[0,10]$ our estimate of the error (19) gives $|g(y)-S(y)|\lessapprox 0.17\left(\frac{x_{n}-x_{0}}{n}\right)^{4}$, which is $3\times 10^{-22}$ times smaller than the bound (18) that is obtained by applying the standard spline error analysis, as in equation (5). In this case, the condition (9) for the applicability of our approximation (19) becomes

[TABLE]

which is a surprisingly low value, for such a huge variation of $f$.

Figure 3 shows the numerical result of the FSSI inversion of $f(x)=x\exp(x)$ in the domain $x\in[0,10]$ using $n=10^{2}$ grid points. Two independent evaluations of the numerical errors are given: i) $|S(y)-f^{-1}(y)|$, where $f^{-1}(y)=W_{S}(y)$ as computed with scipy.special.lambertw routine; ii) $|S(y)-S(f(S(y)))|$. The fact that they agree with each other provides confirmation concerning our treatments of the errors. Moreover, in this case our theoretical prediction of Equation (19) gives $|g(y)-S(y)|\lessapprox 1.7\times 10^{-5}$, and the numerical error not only agrees with it, but it is even much smaller, by almost an order of magnitude, $|g(y)-S(y)|_{\text{numerical}}<2.5\times 10^{-6}$ over the entire interval. In this case, FSSI is more accurate by an astonishing factor $4\times 10^{-23}$ than what could be expected by naively applying the general results for cubic splines, as in equation (18).

V.3 Kepler’s equation

Kepler’s equation for an elliptical orbital motion of eccentricity e can be written as

[TABLE]

where $y$ and $x$ represent the so-called mean and eccentric anomaly, respectively Prussing2012 ; Curtis2014 . The former is the time elapsed since periapsis, as measured in radians, $y=\frac{2\pi t}{T}$, where $T$ is the period of the orbit. The eccentric anomaly $x$ is related to the angle $\theta$ between the position vectors at periapsis and at time $t$, with origin in the center of gravity, through the equation

[TABLE]

A fundamental problem in orbital dynamics Prussing2012 ; Curtis2014 is to obtain the time dependence of the angle $\theta$ describing the position of the orbiting body at time $t$, which requires the inversion of the function $y=f(x)\equiv x-\text{e}\sin x$. Taking into account that the orbit is periodic, and that for $x\in[\pi,2\pi]$ we have $f(x)=2\pi-f(2\pi-x)$, it is sufficient to consider only the interval $x\in[0,\pi]$ to obtain the behavior for all values of $x$. The corresponding co-domain is then $y\in[0,\pi]$ Prussing2012 ; Curtis2014 .

The inverse function $x=g(y)$ will yield the eccentric anomaly as a function of the mean anomaly, and thus the evolution $\theta(t)$ will be obtained. This is usually done in an efficient way using Newton’s method with the first guess $x_{0}=y+\text{e}/2$ Prussing2012 ; Curtis2014 ; Danby1983 .

Here, FSSI is considered as an alternative to Newton-based methods for solving Kepler’s equation. In this case, the values of $M$ and $\mu$ in the theoretical bound of equations (5) and (6) are

[TABLE]

and

[TABLE]

As a concrete example, the case of $\text{e}=0.8$ is considered. Thus, the maximum values are $M=10275.1$, which is obtained for $x=0.166$, and $\mu=\left|1+\text{ec}\right|^{4}\left(\frac{\pi}{n}\right)^{4}=10.4976\left(\frac{\pi}{n}\right)^{4}$, obtained for $x=\pi$. Therefore the bound (5) becomes

[TABLE]

However, the expression from our analytic estimation from Equation (10) becomes

[TABLE]

For $\text{e}=0.8$, the expression in the $|\;\;|$ bracket has a maximum value 21.586 obtained for x=0.214657, therefore we obtain

[TABLE]

Thus, our estimation for the theoretical error (27) in this case is $2\times 10^{-4}$ smaller than what could be expected by naively applying the known bounds on cubic spline interpolation. Here, the condition (9) for the applicability of our approximation (27) becomes

[TABLE]

As a result, for Kepler problem, the FSSI method and the estimation (27) start to be reliable for $n$ as small as the order of ten.

Figures 4 and 5 show the result of the FSSI for the inversion of $f(x)=x-0.8\sin{}x\;{}$ over the domain $x\in[0,\pi]$ using $n=10$ and $n=10^{2}$ grid points, respectively. In these cases, our theoretical prediction of Equation (27) gives $|g(y)-S(y)|\lessapprox 5.5\times(10^{-4}\;\text{or}\;10^{-8})$, respectively, in excellent agreement with our numerical computation over the entire interval.

Again, we provide two independent numerical computations of the error, one obtained by plotting the difference of the FSSI interpolation with the values of $g_{\rm N}(y)$ obtained with Newton’s method, and the other given by the difference $S(y)-S(f(S(y)))$. The fact that these evaluations of the error also agree with each other is a further confirmation of the validity of our error analysis.

By comparing figures 4 and 5, we also see that the accuracy scales with $n^{-4}$, as equation (27) predicts, and that our estimation for the error is reliable even for just $n=10$ grid points.

VI Numerical Comparisons with Newton-based methods

Apart from the numerical calculations for error analysis, we carried out numerical comparisons between FSSI and Newton-based methods (as well as the scipy.lambertw, for the case of Lambert $W$ calculation) for calculating the inverse of single-valued functions. As in the examples of the previous section, the FSSI and Newton-based methods were implemented in the Python programming language, respecting standard practice of minimizing loops and relying upon library function calls (that depend upon compiled code). When possible, we also tested accelerating all methods with Numba JIT compilation, however we found that no considerable difference in empirical execution times could be appreciated.

The algorithms 1 and 2 provide the steps of the FSSI method and the generalized Newton-Raphson method, respectively, used in the benchmark comparisons. This simple version of Newton-Raphson method has been shown to be almost as fast as more elaborate versions, the difference in the execution times being usually below $\sim 30\%$ Palacios2002 .

In algorithm 2, we have called $g_{0}(Y_{k})$ the initial guess for Newton’s method as a function of $Y$, which is to be chosen depending on the problem considered.

In algorithm 1, we have followed the conventions of section III for the arrays, which are indicated in boldface. Accordingly, the operations involving them are to be understood to be valid for the components, i.e. they run over $\{i=0,n\}$ or over $\{i=1,n\}$, for lowercase arrays, or over $\{i=1,N\}$, for uppercase arrays. An exception are the expressions $\mathbf{v}[1:]$ and $\mathbf{v}[:-1]$, which mean the removal of the first or the last element from $\mathbf{v}$, respectively.

In Python, the piecewise polynomial function $P$, corresponding to equation (1), can be obtained in terms of the breakpoints and the coefficients $\mathbf{c_{q}}$ using the subroutine PPoly PPoly , so that P = PPoly. Another possibility is to write an explicit subroutine for computing the polynomial, in which the insertion points $j$ are located by binary search using scipy.searchsorted searchsorted . The two possibilities are shown below:

Subroutine P for $\text{FSSI}\_\textsc{PPoly}$ in Python

1:function P($(\mathbf{c_{3}},\mathbf{c_{2}},\mathbf{c_{1}},\mathbf{c_{0}}),Y$)

2:     P = scipy.PPoly

3:end function

Subroutine P for $\text{FSSI}\_\textsc{Search}$ in Python

1:function P($(\mathbf{c_{3}},\mathbf{c_{2}},\mathbf{c_{1}},\mathbf{c_{0}}),\mathbf{Y}$)

2:     $j=\text{numpy.searchsorted}(\mathbf{y1},\mathbf{Y})$

3:     $\mathbf{P1}=\mathbf{Y}-\mathbf{y_{0_{j}}}$

4:     $\mathbf{P2}=\mathbf{P1}*\mathbf{P1}$

5:     $\mathbf{S}=\mathbf{c_{o_{j}}}+\mathbf{c_{1_{j}}}*\mathbf{P1}+\mathbf{c_{2_{j}}}*\mathbf{P2}+\mathbf{c_{3_{j}}}*\mathbf{P2}*\mathbf{P1}$

6:     return $\mathbf{S}$

7:end function

A discrete analysis of the algorithms 1 and 2 shows that the FSSI executes in constant time $\mathcal{O}(1)$, because once the spline coefficients are obtained with a grid given by $n$ points, all subsequent $N$ function evaluations are equivalent array access through the generating function. However, when $N$ is large, finite cache sizes and the search of the breakpoints overtake this behavior, so that the algorithm follows a linear time dependence $\mathcal{O}(N)$ Oded2008 ; Brent2010 .

In other words, the execution time can be written as $\Delta t_{\textsc{FSSI\_PPoly}}\simeq\epsilon N+\eta$ and $\Delta t_{\textsc{FSSI\_Search}}\simeq\beta N+\alpha$, for the python implementations of FSSI with PPoly or Searchsorted, respectively. As we show below, $\Delta t_{\textsc{FSSI\_PPoly}}<\Delta t_{\textsc{FSSI\_Search}}$ for large $N$, typically $N\gtrsim 10^{4}$, and $\Delta t_{\textsc{FSSI\_PPoly}}>\Delta t_{\textsc{FSSI\_Search}}$ for lower values of $N$. We can then merge the two python routines for P, algorithms FSSI_Search and FSSI_PPoly, by choosing the fastest one with an if statement, e.g. if $N>10^{4}$ do PPoly, else do the routine with searchsorted. The execution time for this combined routine is $\Delta t\simeq\alpha+\epsilon N$, i.e. it behaves as $\mathcal{O}(1)+\epsilon\,\mathcal{O}(N)$.

On the other hand, for the Newton minimization based methods, all evaluations of the function inverse occur with an average number of iterations, $m$ (as seen in the while loop of lines 6-10), therefore, these algorithms have a lower bound linear time behavior $\mathcal{O}(mN)$ for all values of $N$.

To obtain an empirical execution time comparison between methods, we ran the benchmarks for two cases: the calculation of the Lambert W function, and the solution of Kepler’s problem.

The details of the numerical comparison are as follows:

• •

Hardware: The numerical comparisons were carried out on a modest desktop computer (a 64 bit Intel i5-2400 CPU 3.10GHz, with 32GB memory, and with the Ubuntu/Linux operating system with 4.13.16 kernel).

• •

Tolerance: For each case the same level as the error of the FSSI in this case: For the Lambert W problem, we used a tolerance $2\times 10^{3}/n^{4}$ for scipy.lambertw and Newton; For the Kepler solution, we used a tolerance $6/n^{4}$ for Newton and Pynverse Pynverse quasi-Newton method.

• •

For Lambert W, we chose the simplest first guess, $g_{0}(Y_{k})=\frac{x_{0}+x_{n}}{2}=5$. Of course, better choices may be found, but we want to use this case to compare FSSI and Newton in the absence of a good first guess. On the other hand, we also penalize the FSSI method by taking the tolerance for Newton-based methods equal to the theoretical error of FSSI, which overestimates the numerical error by an order of magnitude as shown in section V.

• •

In the case of Kepler’s equation, we take $\text{e}=0.8$ and we use a very good first guess, $g_{0}(Y_{k})=Y_{k}+\frac{\text{e}}{2}$, as was mentioned in section V.

Figure 6 shows empirical execution time comparisons between different numerical algorithms and FSSI for calculating Lambert W and for solving Kepler’s equation. The results support the theoretical expectations described above. For the FSSI method, there is a wide range of $N$ values for which the $\mathcal{O}(N)$ behavior is negligible as compared with the $\mathcal{O}(1)$ behavior; however for very large $N$, when the $\mathcal{O}(N)$ part dominates, the linear coefficient $\epsilon$ is several orders of magnitude smaller than those of the other methods available.

In all the cases, Pynverse Pynverse (based on a quasi-newton optimization) is much slower than the other methods considered, which is not a surprise since it is meant to be universal, rather than fast. Therefore, we will limit our discussion to the comparison between FSSI and Newton-Raphson methods. As we see from figure 6, FSSI is not only universal, but it is also fast, and for large $N$ it is the fastest method.

These results have been used to obtain linear fits to the data. For example, for Kepler’s problem with $n=50$, corresponding to tolerance $10^{-6}$ rad (which can be a sufficient accuracy for orbit determination in many cases), we found $\Delta t_{\textsc{FSSI\_PPoly}}\simeq 2.1\times 10^{-8}N+1.9\times 10^{-4}$ and $\Delta t_{\textsc{FSSI\_Search}}\simeq 5.3\times 10^{-8}N+7.2\times 10^{-5}$. These values of the coefficients have been obtained by separate fits to the low $N$ data, for $\eta$ and $\alpha$, and to the high $N$ data, for $\epsilon$ and $\beta$, in order to get the best estimates in these regimes, so that the approximation is slightly worse for $10^{3}\lesssim N\lesssim 10^{4}$.

In any case, as shown in figure 6, FSSI_PPoly is faster than FSSI_Search for $N\gtrsim 10^{4}$ and slower for $N\lesssim 10^{4}$. By choosing the best of the two variants with an if statement, we obtain a combined behavior $\Delta t_{\textsc{FSSI}}\simeq 2.1\times 10^{-8}N+7.2\times 10^{-5}$. This should be compared with the execution time for Newton-Raphson method, $\Delta t_{\textsc{Newton}}=4.2\times 10^{-5}N$. We find that $\Delta t_{\textsc{Newton}}>\Delta t_{\textsc{FSSI}}$ for every $N\geq 2$, and that FSSI is $\sim 2\times 10^{3}$ faster than Newton-Raphson for large $N$.

Similarly, for Kepler’s problem with $n=10^{4}$, corresponding to tolerance $6\times 10^{-16}$ rad, we obtain a combined behavior $\Delta t_{\textsc{FSSI}}\simeq 2.1\times 10^{-8}N+1.1\times 10^{-3}$ while $\Delta t_{\textsc{Newton}}=5.0\times 10^{-5}N$. We find that $\Delta t_{\textsc{Newton}}>\Delta t_{\textsc{FSSI}}$ for every $N\gtrsim 20$, and that FSSI is still $\sim 2\times 10^{3}$ faster than Newton-Raphson for large $N$.

For Lambert W with $n=50$, corresponding to tolerance $3\times 10^{-4}$, we obtain a combined behavior $\Delta t_{\textsc{FSSI}}\simeq 2.1\times 10^{-8}N+6.2\times 10^{-5}$ while $\Delta t_{\textsc{Newton}}=1.3\times 10^{-4}N$. We find that $\Delta t_{\textsc{Newton}}>\Delta t_{\textsc{FSSI}}$ for every $N$, and that FSSI is $\sim 6\times 10^{3}$ faster than Newton-Raphson for large $N$. This shows that, in the lack of a good first guess, FSSI can be better than Newton-Raphson method for every value of $N$. Of course, for small $N$ the specific, semi-analytic routine scipy.lambertw SciPyLambertW , having $\Delta t_{\textsc{scipy\_LambertW}}=3.8\times 10^{-7}N$ outperforms the FSSI, but surprisingly the opposite is true in the large $N$ regime, in which FSSI is $\sim 20$ times faster than scipy.lambertw.

Finally, for Lambert W with $n=10^{4}$, corresponding to tolerance $2\times 10^{-13}$, we obtain a combined behavior $\Delta t_{\textsc{FSSI}}\simeq 2.1\times 10^{-8}N+1.1\times 10^{-3}$ while $\Delta t_{\textsc{Newton}}=1.3\times 10^{-4}N$. We find that $\Delta t_{\textsc{Newton}}>\Delta t_{\textsc{FSSI}}$ for every $N\gtrsim 8$, and that FSSI is $\sim 7\times 10^{3}$ faster than Newton-Raphson for large $N$. This shows that, in the lack of a good first guess, FSSI is much better than Newton-Raphson method for every value of $N$. Again, for large $N$, $\Delta t_{\textsc{FSSI}}<\Delta t_{\textsc{scipy\_LambertW}}=3.8\times 10^{-7}N$ by a factor $\sim 20$.

Note that the values of the FSSI execution times are almost equal for Kepler and Lambert problems with the same values of $n$ and $N$. The fact that the method performs at the same speed when applied to functions that are very different from each other is a further proof of its universality.

VII Conclusions

In this study, we described a scheme, called FSSI, based on switch and spline to invert monotonic functions under very general conditions. Moreover, we derived analytical expressions for the associated theoretical errors of this method, and tested it on examples that are of interest in physics, including the computation of Lambert W function and the solution of Kepler’s equation. As a summary, the FSSI method has several advantages over other more standard techniques for inverting functions:

• •

It is simple and universal and, unlike Newton methods, it does not require any initial guess.

• •

The error is much smaller than what could be expected from general spline analysis, by a $\sim 10^{-22}$ factor for $W\in[0,10]$ or by a factor $2\times 10^{-4}$ for Kepler problem.

• •

This scheme is superior to, and much faster than, Newton-Raphson method when the latter is difficult to apply, when no good first guess is available, or when the values of the inverse function are required on an entire interval or in a large number of different points.

• •

When applied to Kepler’s problem (e.g. with eccentricity $\text{e}=0.8$), FSSI becomes faster than Newton’s methods for $N$ greater than a few points, and is $\sim 2\times 10^{3}$ times faster for large $N$. If the requested accuracy is of the order of $10^{-6}$ rad, which is a low enough value for most applications, the speed of the FSSI algorithm is faster than Newton’s for $N\geq 2$.

• •

The $N$ dependence of the scheme can be described as $\mathcal{O}(1)+\epsilon\,\mathcal{O}(N)$. For a wide range of $N$ values, the $\mathcal{O}(N)$ behavior is negligible as compared with the $\mathcal{O}(1)$ behavior; however for very large $N$, when the $\mathcal{O}(N)$ part dominates, the linear coefficient $\epsilon$ is several orders of magnitude smaller than those of the other methods available.

For all these reasons, we believe that this method could become a competitive choice for inverting functions in a wide range of applications, and the first choice for solving Kepler’s equation.

VIII Acknowledgements

We thank A. Paredes, D. González-Salgado and H. Michinel for discussions. This work has been supported by grants FIS2017-83762-P from Ministerio de Economía y Competitividad (Spain), and grant GPC2015/019 from Consellería de Cultura, Educación e Ordenación Universitaria (Xunta de Galicia).