Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

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

Fast Switch and Spline Scheme for Accurate Inversion of Nonlinear Functions: The New First Choice... Numerically obtaining the inverse of a function is a common task for many scienti c 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 speci c 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 de ned 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  10 factor for the computation of W in the range W (y) 2 [0; 10], and by a factor 2  10 for Kepler's problem. In our tests, this 4 3 scheme is much faster than Newton-Raphson's method, by a factor in the range 10 to 10 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 rad, the algorithm outperforms Newton's method for all values of the number of points N  2. arXiv:1812.02273v2 [physics.comp-ph] 17 Dec 2018 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 [1, 2], the computation of Lambert W function [3, 4], and the solution of Kepler's equation for the orbital motion of a body in a gravitational eld [5, 6]. 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 speci c 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 [7] or some of its variants [8{11]. 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 rst 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 [12{14]. 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 nding the zeros of the functions F (x)  f (x) y. If a good initial guess x of the true y 0 value x = f (y) of the zero is available, the zeros of F can be computed by recursively y f (x ) F (x ) k applying the equation x = x ; i.e. x = x + . k+1 k 0 k+1 k F (x ) y k f (x ) 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 speci c spline routine that makes the FSSI more accurate and much 2 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 minimiza- tion methods) for computing values of inverse functions, especially if a good rst approxima- tion is dicult to obtain. Moreover, the FSSI is shown to be superior to Newton-Raphson when the values f (y) of the inverse function are required for many di erent 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 di erent instants, depending on the eccentricity e and the re- quested precision. For example, for e = 0:8 and accuracy  10 rad, the FSSI algorithm is already faster than Newton's even when the computation is done at N = 2 points, and 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 2 [x ; x ]. The function min max 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 (y) in the co-domain. The FSSI consists of a two step approach. First, when the input function is given an- alytically, the values of f on a given grid of points x , for j = 1; ; n, are computed to obtain the matrix (x ; y ), where y = f (x ). The matrix (y ; x ), obtained by switching the j j j j j j arrays, gives then the exact values g(y ) = x of the inverse function on the grid y . From j j j this modi ed matrix, the cubic spline interpolant S(y) of (y ; g(y )) is computed by using j j 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 speci ed on a grid by a set of tuples (x ; y ). In a high-level programming language such as Python, this tuple array is represented j j as (x; y), and the FSSI interpolant S(y) can be obtained by calling a cubic spline routine of 3 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 . Input Output analytic function Vector values and OR Domain of inverse: Error Speed: Evaluate Switch domains for spline Calculate coefficients for spline of Calculate breakpoints Evaluate on co-domain Build polynomial on grid Figure 1: Flow diagram of the FSSI method for obtaining the function inverse f . The diagram indicates the key steps of the method, as well as how it is interfaced to the input and output. Figure 1 shows the ow 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 () and its derivative f (), 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 is determined by the number of points n + 1 of the input grid, and which is used for sampling N points fY ; ; Y g in the function's co-domain. In subsequent sections, we 1 N show that the FSSI algorithm has an error bound proportional to jkxj , for k a constant (described in the text), and with a maximal time complexity of O(1)+ 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 . This is the case for Kepler's equation and for all the other examples that we 4 0 will consider hereafter. Therefore, we can pro t from the knowledge of f to design a speci c 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 di erence 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 coecients of the spline [7], where n is the number of grid intervals. An exception is Akima's cubic spline [15, 16], which is fast, diagonal and regular. The speci c spline that is designed here is based on a similar idea to Akima's, but it is signi cantly 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 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 di erent in the 0 0 problem of the inversion of a function f (x). In this case, the derivatives g (y ) = 1=f (x ) j j can be given on a grid, while the values g(y) are not known. Let us build the speci c spline S(y) piecewise in each interval, S(y) = S (y) for y < y < j j y , where j takes the values j = 0; ; n 1. If we de ne the arrays y  (y ; ; y ) j+1 0 0 n1 and y  (y ; ; y ), obtained by removing the last and the rst point of the y array, 1 1 n respectively, then the j-th interval can also be written as y < y < y . In this interval, the 0 1 j j cubic spline can be de ned as S (y) = c y y ; (1) j q 0 j j q=0 where for each value of q = 1; ; 4 the coecients c can also be thought as the n compo- nents of an array c . Since the values of the derivative f (x ) of the input function are known on the grid points, we can construct an array d whose n components are the derivatives of the inverse function g on the points y = f (x ), j j d  g (y ) = ; for j = 0; ; n: (2) j j f (x ) As we did for y, which was used to generate the arrays y and y by removing one end 1 2 point, it is convenient to de ne similar arrays of n components also from x and d, namely 5 x  (x ; ; x ), x  (x ; ; x ), d  (d ; ; d ), and d  (d ; ; d ). With 0 0 n1 1 1 n 0 0 n1 1 1 n this convention, we have to choose the spline coecients that lead to the best approximation of the inverse function g(y). The most natural choice is to force S to coincide with the known values of the inverse function, x and x , at the end points, and to ask the same for the 0 1 j j derivatives d and d . In other words, the conditions to be imposed in each interval are, 0 1 j j 0 0 S (y ) = x ; S (y ) = x ; S (y ) = d ; S (y ) = d ; (3) j 0 0 j 1 1 0 0 1 1 j j j j j j j j j j where S (y) is given by equation 1. For every xed value of j, these conditions give a system of four equations for the four unknown coecients c , c , c and c , which is decoupled from the similar systems of 0 1 2 3 j j j j equations corresponding to di erent 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 coecients [7], 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 xed j can be solved analytically in a straight- forward 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 vw+2z of components in such a way that an equality like, e.g., u = between vectors having v w +2z j j j the same number of elements represents the equations u = for every j. With this convention, equations (3) give the solution c = x ; 0 0 c = d ; 1 0 (2 d + d ) (y y ) 3 (x x ) 0 1 0 1 0 1 c = ; (y y ) 0 1 (d + d ) (y y ) 2 (x x ) 0 1 0 1 0 1 c = : (4) (y y ) 0 1 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 00 00 adding the additional conditions S (y ) = d and S (y ) = d to equation (3) we can also 0 0 1 1 j j j j j j 6 design a diagonal quintic spline, and if also the third derivatives are known the additional 000 000 conditions S (y ) = d and S (y ) = d allow for the construction of a septic spline. We 0 0 1 1 j j j j j j 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 speci c 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 fth 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 di erences 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. A. Derivation of an upper bound on the error of the FSSI by using the known analysis of cubic spline interpolation Following Ref. [7, 17], 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 jg(y) S(y)j  M ; (5) (4) 4 where M = max g (y) , and  = max (y y ) . j+1 j y yy 0jn1 max min 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, 00 3 (3) 00 (4) 15f (x) 10f (x)f (x) f (x) M = max + : (6) 0 7 0 6 0 5 x xx 0 n f (x) f (x) f (x) Equations (5) and (6), along with the de nition of , can be used to obtain an upper limit on the error. As shown in examples below, the actual errors are several orders of 7 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. B. Ab initio derivation of an improved estimation of the error for the FSSI Let us assume that the function g(y) is in nitely di erentiable. Thus, it can be expanded P (q) 1 g (y ) in a Taylor series g(y) = (yy ) around one of the points of the grid y = f (x ), j j j q=0 q! chosen to be the closest grid point to y, so that jyy j  jyy j and jyy j  jyy j. j j+1 j j1 We also assume that the cubic spline interpolation is made in such a way so that it is equivalent to the rst terms of this Taylor expansion, expanded around the same point and P (q) 3 g (y ) truncated beyond cubic order, S(y) = (y y ) . Therefore, the di erence is q=0 q! (q) (4) (4) g (y ) g (y  ) g (y ) j j j n 4 4 jg(y) S(y)j = (y y ) = (y y )  (y y ) ; (7) j j j n! 4! 4! q=4 where y  is an unknown intermediate point between y and y , and the last approximation is j j expected to hold for suciently small values of jyy j, which is the case when a suciently high number of grid points is chosen. On one hand, the exact equality in equation (7) can be translated in the following bound, (4) g (y ) y y 1 j+1 j jg(y) S(y)j  max max = M ; (8) y y y 0jn1 min max 4! 2 384 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 y is due to the j+1 j y y j+1 j fact that y was chosen as the closest grid point to y, so that jy y j  . j j On the other hand, equation (7) can be elaborated further and expressed in terms of the function f (x), 00 3 (3) 00 (4) 1 15f (x ) 10f (x )f (x ) f (x ) j j j j 4 jg(y) S(y)j  + [f (x )(x x )] ; j j 0 7 0 6 0 5 4! f (x ) f (x ) f (x ) j j j where x = g(y) and y y  f (x )(x x ). The last approximation is expected to be j j j accurate over the entire interval provided the following condition f (x) max jx x j max  1 (9) j+1 j 0jn1 x xx n 2f (x) 8 is satis ed. x x j+1 j Assuming that the point x is the one closest to x, so that jx x j  max , j j 0jn1 the following estimation is obtained jg(y) S(y)j / max jx x j  (10) j+1 j 0jn1 00 3 (3) 00 (4) 15f (x) 10f (x)f (x) f (x) max + : 0 3 0 2 0 x xx 0 n f (x) f (x) f (x) As described in the next section, this error approximation is usually much smaller, pos- sibly 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 satis es 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 grid, we can substitute 4 (x x ) x x n 0 n 0 max jx x j = in the condition (9), and max jx x j = in the j+1 j j+1 j n n 0jn1 0jn1 estimation of the error (10). Not only is this the simplest choice, if the grid is not given 4 (x x ) n 0 otherwise, but it is usually also the best option. In fact, when jx x j = for j+1 j 4 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), de ned over a domain x  x  x . For the numerical computation, we developed a python 0 n 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 rst example, in which the exact inverse function g(y) is known, we also provide an additional evaluation of the error by computing the di erence S(y)g(y). In both cases, we use a grid Y that contains 10 times as many points as the original grid y = f (x ). In fact, by construction S(y) is exactly equal j j to g(y){within machine errors{over the original grid y , so it is important to ensure that S(y) is compared with g(y) in between the grid points y . It is true that even if they were equal in number the points Y , chosen to be equally spaced, would not coincide in general 9 0 with the y , whose spacing is variable and roughly proportional to f (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 di erence 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 speci c 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 [18, 19]. 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 speci c spline. An exception is Akima routine [16], which is less accurate by three orders of magnitude in the bulk of the interval. This is an additional reason for preferring our speci c spline, besides the fact that it is the fastest one. A. Exponential The rst 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 , from the bounds expressions of (5) and (6), can be readily computed. The results are, 4x 4x M = max 6e = 6e ; (11) x xx 0 n and 4 4 x x x x n 0 n 0 4x 4x = max e = e ; (12) x xx n 0 n n so that the bound on the error as computed from equations (5) and (6) is 6 x x n 0 4(x x ) n 0 jg(y) S(y)j  e : (13) 384 n 10 On the other hand, the analytic estimation we derived in Equation (10) gives 6 x x n 0 jg(y) S(y)j / : (14) 384 n 4(x x ) n 0 Therefore, this analytic error approximation of FSSI is smaller by a factor e as compared to the limit that was derived in equation (13) by applying the standard cubic spline error bound. For example, if x x = 10, then our error estimation is a factor exp(40), n 0 i.e. 17 orders of magnitude, smaller than what could be expected from the literature. In order to bene t by this accuracy improvement, the grid must be chosen in such a way that the condition (9) is satis ed, i.e. x x f (x) x x n 0 n 0 =  1: (15) n 2f (x) 2n 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 2 [0; 10], the condition becomes n  5, so that values of n of the order of 50 or larger could be sucient. This is also what we have observed by performing numerical computations for di erent values of n. In general, for n & 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 2 [0; 10] using n = 10 grid points. In this case, our theoretical prediction of Equation (14) gives jg(y) S(y)j / 1:6  10 , which is in excellent agreement with the numerical computation over the entire interval. The results for this case also con rm 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. B. 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) [3, 4]. 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 11 Figure 2: Result of the FSSI applied to the function f (x) = exp(x) (top left) over the domain x 2 [0; 10]. The FSSI interpolant S(y) is shown for n = 10 grid points (top right), together with 1 1 two independent evaluations of the numerical errors: i) jS(y) f (y)j, where f (y) = ln(y) (bottom left); ii) jS(y) S(f (S(y)))j (bottom right). of M and  from equations (5) and (6) for the theoretical bound are, x x x x x x x x 15 (e x + 2e ) 10 (e x + 3e ) (e x + 2e ) e x + 4e M = max + ; (16) 7 6 5 x x x x x x x xx 0 n (e x + e ) (e x + e ) (e x + e ) which is a monotonically decreasing function for x > 1, and x x n 0 x x 4 = max je x + e j ; (17) x xx n 0 n which increases monotonically. Therefore, the bound (5) becomes jg(y) S(y)j  (x = x ) M (x = x ): (18) n 0 20 x x n 0 For example, over the domain x 2 [0; 10] this gives jg(y) S(y)j / 5:7 10 . On the other hand, our analytical estimation from Equation (10) becomes 1 x x n 0 jg(y) S(y)j /  (19) 384 n x x x x x x x x 15 (e x + 2e ) 10 (e x + 3e ) (e x + 2e ) e x + 4e max + : 3 2 x x x x x x x xx 0 n e x + e (e x + e ) (e x + e ) 12 The function to be maximized in equation (19) monotonically decreases for x > 1, so that its maximum is achieved for x = x . Thus, over the domain x 2 [0; 10] our estimate of x x n 22 the error (19) gives jg(y) S(y)j / 0:17 , which is 3 10 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 x x x + 2 0 n n  max = 10; (20) x xx 2 n x + 1 which is a surprisingly low value, for such a huge variation of f . Figure 3: Numerical result of the FSSI applied to the function f (x) = x exp(x) (top left) over the domain x 2 [0; 10]. The FSSI interpolant S(y) is shown for n = 10 grid points (top right), together 1 1 with two independent evaluations of the numerical errors: i) jS(y)f (y)j, where f (y) = W (y) as computed with scipy.special.lambertw routine (bottom left); ii) jS(y) S(f (S(y)))j (bottom right). Figure 3 shows the numerical result of the FSSI inversion of f (x) = x exp(x) in the domain x 2 [0; 10] using n = 10 grid points. Two independent evaluations of the nu- 1 1 merical errors are given: i) jS(y) f (y)j, where f (y) = W (y) as computed with scipy.special.lambertw routine; ii) jS(y) S(f (S(y)))j. The fact that they agree with each other provides con rmation concerning our treatments of the errors. Moreover, in this case 13 5 our theoretical prediction of Equation (19) gives jg(y) S(y)j / 1:7  10 , and the nu- merical error not only agrees with it, but it is even much smaller, by almost an order of magnitude, jg(y) S(y)j < 2:5 10 over the entire interval. In this case, FSSI is numerical more accurate by an astonishing factor 4  10 than what could be expected by naively applying the general results for cubic splines, as in equation (18). C. Kepler's equation Kepler's equation for an elliptical orbital motion of eccentricity e can be written as y = x e sin x; (21) where y and x represent the so-called mean and eccentric anomaly, respectively [5, 6]. The 2t former is the time elapsed since periapsis, as measured in radians, y = , where T is the period of the orbit. The eccentric anomaly x is related to the angle  between the position vectors at periapsis and at time t, with origin in the center of gravity, through the equation 1 + e x = 2 arctan tan : (22) 1 e 2 A fundamental problem in orbital dynamics [5, 6] is to obtain the time dependence of the angle  describing the position of the orbiting body at time t, which requires the inversion of the function y = f (x)  x e sin x. Taking into account that the orbit is periodic, and that for x 2 [; 2] we have f (x) = 2 f (2 x), it is sucient to consider only the interval x 2 [0; ] to obtain the behavior for all values of x. The corresponding co-domain is then y 2 [0; ] [5, 6]. The inverse function x = g(y) will yield the eccentric anomaly as a function of the mean anomaly, and thus the evolution (t) will be obtained. This is usually done in an ecient way using Newton's method with the rst guess x = y + e=2 [5, 6, 8]. Here, FSSI is considered as an alternative to Newton-based methods for solving Kepler's equation. In this case, the values of M and  in the theoretical bound of equations (5) and (6) are 3 2 15e sin x 10e sin x cos x e sin x M = max + + (23) 7 6 5 0x (1 e cos x) (1 e cos x) (1 e cos x) and = max j1 e cos xj : (24) 0x 14 As a concrete example, the case of e = 0:8 is considered. Thus, the maximum values 4 4 are M = 10275:1, which is obtained for x = 0:166, and  = j1 + ecj = 10:4976 , n n obtained for x = . Therefore the bound (5) becomes 2:7 10 jg(y) S(y)j . : (25) However, the expression from our analytic estimation from Equation (10) becomes 3 2 1  15e sin x 10e sin x cos x e sin x jg(y) S(y)j / max + + : (26) 3 2 0x 384 n (1 e cos x) (1 e cos x) 1 e cos x For e = 0:8, the expression in the j j bracket has a maximum value 21.586 obtained for x=0.214657, therefore we obtain 5:5 jg(y) S(y)j / : (27) Thus, our estimation for the theoretical error (27) in this case is 2  10 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 e sin x n  max ' 2: (28) 2 0x 1 e cos x 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 2 [0; ] using n = 10 and n = 10 grid points, respectively. In these 4 8 cases, our theoretical prediction of Equation (27) gives jg(y) S(y)j / 5:5 (10 or 10 ), 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 di erence of the FSSI interpolation with the values of g (y) obtained with Newton's method, and the other given by the di erence S(y) S(f (S(y))). The fact that these evaluations of the error also agree with each other is a further con rmation of the validity of our error analysis. By comparing gures 4 and 5, we also see that the accuracy scales with n , as equation (27) predicts, and that our estimation for the error is reliable even for just n = 10 grid points. 15 Figure 4: Numerical result of the FSSI applied to the function f (x) = x 0:8 sin x over the domain x 2 [0; ] (top left), corresponding to Kepler's equation for an elliptical orbit of eccentricity 0.8. The FSSI interpolant S(y) is shown for n = 10 grid points (top right), together with two independent evaluations of the numerical errors: i) jS(y) g (y)j, where g (y) is computed with Newton's N N method (bottom left); ii) jS(y) S(f (S(y)))j (bottom right). VI. NUMERICAL COMPARISONS WITH NEWTON-BASED METHODS Apart from the numerical calculations for error analysis, we carried out numerical com- parisons 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 imple- mented 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 di erence 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, 16 Figure 5: Numerical result of the FSSI applied to the function f (x) = x 0:8 sin x over the domain x 2 [0; ] (top left), corresponding to Kepler's equation for an elliptical orbit of eccentricity 0.8. The FSSI interpolant S(y) is shown for n = 10 grid points (top right), together with two independent evaluations of the numerical errors: i) jS(y) g (y)j, where g (y) is computed with Newton's N N method (bottom left); ii) jS(y) S(f (S(y)))j (bottom right). the di erence in the execution times being usually below  30% [11]. 17 Algorithm 1 Benchmark for FSSI 1: procedure bench FSSI(Y; x; f (x); f (x)) 2: y = f (x) 1 Algorithm 2 Benchmark for Newton 3: d = f (x) 1: procedure 4: c = x[: 1] bench Newton(Y; f (x); f (x); tol) 5: c = d[: 1] 2: for Y in Y do 6: d1 = d[1 :] 3: X = g (Y ) k 0 k 7: xx = c x[1 :] Y f (X ) k k 4:  = f (X ) 8: y0 = y[: 1] 5: while  > tol do 9: y1 = y[1 :] Y f (X ) k k 6: = f (X ) 10: yy = y0 y1 7: X = X + k k 11: yy2 = yy yy 8:  = j j 12: yd1 = yy d1 9: end while 13: yd0 = yy c1 10: end for 2yd0+yd13xx 14: c = yy2 11: return X . =f (Y) yd0+yd12xx 15: c = yy2yy 12: end procedure 16: call P: X = P((c ; c ; c ; c ); Y) 3 2 1 0 17: return X . =f (Y) 18: end procedure In algorithm 2, we have called g (Y ) the initial guess for Newton's method as a function 0 k 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 fi = 0; ng or over fi = 1; ng, for lowercase arrays, or over fi = 1; Ng, for uppercase arrays. An exception are the expressions v[1 :] and v[: 1], which mean the removal of the rst or the last element from v, respectively. In Python, the piecewise polynomial function P , corresponding to equation (1), can be obtained in terms of the breakpoints and the coecients c using the subroutine PPoly [20], so that P = PPoly. Another possibility is to write an explicit subroutine for com- puting the polynomial, in which the insertion points j are located by binary search using 18 scipy.searchsorted [21]. The two possibilities are shown below: Subroutine P for FSSI Search in Python 1: function P((c ; c ; c ; c ); Y) 3 2 1 0 Subroutine P for FSSI PPoly in Python 2: j = numpy.searchsorted(y1; Y) 1: function P((c ; c ; c ; c ); Y ) 3: P1 = Y y 3 2 1 0 0 2: P = scipy.PPoly 4: P2 = P1 P1 3: end function 5: S = c +c P1+c P2+c P2P1 o 1 2 3 j j j j 6: return S 7: end function A discrete analysis of the algorithms 1 and 2 shows that the FSSI executes in constant time O(1), because once the spline coecients 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, nite cache sizes and the search of the breakpoints overtake this behavior, so that the algorithm follows a linear time dependence O(N ) [22, 23]. In other words, the execution time can be written as t ' N +  and FSSI PPoly t ' N + , for the python implementations of FSSI with PPoly or Search- FSSI Search sorted, respectively. As we show below, t < t for large N , typically FSSI PPoly FSSI Search N & 10 , and t > t for lower values of N . We can then merge the two FSSI PPoly FSSI Search python routines for P, algorithms FSSI Search and FSSI PPoly, by choosing the fastest one with an if statement, e.g. if N > 10 do PPoly, else do the routine with searchsorted. The execution time for this combined routine is t ' + N , i.e. it behaves as O(1) + 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 O(mN ) for all values of N . To obtain an empirical execution time comparison between methods, we ran the bench- marks 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 19 (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 3 4 Lambert W problem, we used a tolerance 2 10 =n for scipy.lambertw and Newton; For the Kepler solution, we used a tolerance 6=n for Newton and Pynverse [24] quasi- Newton method. x +x 0 n For Lambert W, we chose the simplest rst guess, g (Y ) = = 5. Of course, 0 k better choices may be found, but we want to use this case to compare FSSI and Newton in the absence of a good rst 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 e = 0:8 and we use a very good rst guess, g (Y ) = Y + , as was mentioned in section V. 0 k k Figure 6 shows empirical execution time comparisons between di erent numerical algo- rithms 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 O(N ) behavior is negligible as compared with the O(1) behavior; however for very large N , when the O(N ) part dominates, the linear coecient is several orders of magnitude smaller than those of the other methods available. In all the cases, Pynverse [24] (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 gure 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 ts to the data. For example, for Kepler's problem with n = 50, corresponding to tolerance 10 rad (which can be a sucient accuracy 8 4 for orbit determination in many cases), we found t ' 2:1 10 N + 1:9 10 FSSI PPoly 8 5 and t ' 5:3  10 N + 7:2  10 . These values of the coecients have been FSSI Search obtained by separate ts to the low N data, for  and , and to the high N data, for  and 20 Newton FSSI_s Newton pynverse FSSI_s FSSI_p FSSI_p pynverse Newton Newton FSSI_s FSSI_p FSSI_s FSSI_p N evaluations N evaluations FSSI_s FSSI_s FSSI_s FSSI_s Newton Newton FSSI_p FSSI_p scipy.lambertw scipy.lambertw scipy.lambertw Newton Newton scipy.lambertw FSSI_p FSSI_p N evaluations N evaluations Figure 6: Numerical comparisons of FSSI and other methods for the solution of Kepler's equation (top) and computation of the Lambert W function (bottom). FSSI p and FSSI s stand for the algorithm using PPoly or Searchsorted subroutines, respectively. , in order to get the best estimates in these regimes, so that the approximation is slightly 3 4 worse for 10 . N . 10 . In any case, as shown in gure 6, FSSI PPoly is faster than FSSI Search for N & 10 and slower for N . 10 . By choosing the best of the two variants with an if statement, we 8 5 obtain a combined behavior t ' 2:1 10 N + 7:2 10 . This should be compared FSSI with the execution time for Newton-Raphson method, t = 4:210 N . We nd that Newton t > t for every N  2, and that FSSI is  2 10 faster than Newton-Raphson Newton FSSI for large N . 4 16 Similarly, for Kepler's problem with n = 10 , corresponding to tolerance 6 10 rad, we 8 3 5 obtain a combined behavior t ' 2:110 N +1:110 while t = 5:010 N . FSSI Newton We nd that t > t for every N & 20, and that FSSI is still  2 10 faster Newton FSSI than Newton-Raphson for large N . For Lambert W with n = 50, corresponding to tolerance 3 10 , we obtain a combined 8 5 4 behavior t ' 2:1 10 N + 6:2 10 while t = 1:3 10 N . We nd that FSSI Newton t > t for every N , and that FSSI is  6 10 faster than Newton-Raphson for Newton FSSI Execution Time (sec) Execution Time (sec) Execution Time (sec) Execution Time (sec) large N . This shows that, in the lack of a good rst guess, FSSI can be better than Newton- Raphson method for every value of N . Of course, for small N the speci c, semi-analytic routine scipy.lambertw [25], having t = 3:8 10 N outperforms the FSSI, scipy LambertW but surprisingly the opposite is true in the large N regime, in which FSSI is  20 times faster than scipy.lambertw. 4 13 Finally, for Lambert W with n = 10 , corresponding to tolerance 2  10 , we obtain 8 3 4 a combined behavior t ' 2:1  10 N + 1:1  10 while t = 1:3  10 N . FSSI Newton We nd that t > t for every N & 8, and that FSSI is  7  10 faster than Newton FSSI Newton-Raphson for large N . This shows that, in the lack of a good rst guess, FSSI is much better than Newton-Raphson method for every value of N . Again, for large N , t < t = 3:8 10 N by a factor  20. FSSI scipy LambertW 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 di erent 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 expres- sions 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, 22 4 by a  10 factor for W 2 [0; 10] or by a factor 2 10 for Kepler problem. This scheme is superior to, and much faster than, Newton-Raphson method when the latter is dicult to apply, when no good rst guess is available, or when the values of 22 the inverse function are required on an entire interval or in a large number of di erent points. When applied to Kepler's problem (e.g. with eccentricity e = 0:8), FSSI becomes faster than Newton's methods for N greater than a few points, and is  2 10 times faster for large N . If the requested accuracy is of the order of 10 rad, which is a low enough value for most applications, the speed of the FSSI algorithm is faster than Newton's for N  2. The N dependence of the scheme can be described as O(1) +O(N ). For a wide range of N values, the O(N ) behavior is negligible as compared with the O(1) behavior; however for very large N , when the O(N ) part dominates, the linear coecient  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 rst choice for solving Kepler's equation. VIII. ACKNOWLEDGEMENTS We thank A. Paredes, D. Gonz alez-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 on e Ordenaci on Universitaria (Xunta de Galicia). [1] Toshio Fukushima. Numerical computation of inverse complete elliptic integrals of rst and second kinds. Journal of Computational and Applied Mathematics, 249:37{50, 2013. [2] John P. Boyd. Four ways to compute the inverse of the complete elliptic integral of the rst kind. Computer Physics Communications, 196:13{18, 2015. [3] R. M. Corless, G. H. Gonnet, D. E. G. Hare, D. J. Je rey, and D. E. Knuth. On the lambert w function. Advances in Computational Mathematics, 5:329{359, 1996. [4] Darko Veberi c. Lambert w function for applications in physics. Computer Physics Communications, 183:2622{2628, 2012. 23 [5] J. E. Prussing and B. A. Conway. Orbital Mechanics. Oxford University Press, 2 edition, [6] Howard D. Curtis. Orbital Mechanics for Engineering Students. Elsevier, 3 edition, 2014. [7] John H. Mathews and Kurtis D. Fink. Numerical Methods Using MATLAB. Pearson, 4 edition, 2004. [8] J.M.A. Danby and T.M. Burkardt. The solution of kepler's equation, i. Celestial Mechanics, 31:95{107, 1983. [9] J.M.A. Danby and T.M. Burkardt. The solution of kepler's equation, iii. Celestial Mechanics, 40:303{312, 1987. [10] J. Gerlach. Accelerated convergence in newton's method. SIAM Rev., 36:272{276, 1994. [11] M. Palacios. Kepler equation and accelerated newton method. Journal of Computational and Applied Mathematics, 138:335{346, 2002. [12] B. A. Conway. An improved algorithm due to laguerre for the solution of kepler's equation. Celestial Mechanics, 39:199{211, 1986. [13] E. D. Charles and J. B. Tatum. The convergence of newton{raphson iteration with kepler's equation. Celestial Mechanics and Dynamical Astronomy, 69:357{372, 1998. [14] Laura Stumpf. Chaotic behaviour in the newton iterative function associated with kepler's equation. Celestial Mechanics and Dynamical Astronomy, 74:95{109, 1999. [15] Hiroshi Akima. A new method of interpolation and smooth curve tting based on local procedures. Journal of the Association for Computing Machinery, 17(4):589{602, 1970. [16] SciPy, 2016. https://docs.scipy.org/doc/scipy-0.18.1/reference/generated/scipy.interpolate.Akima1DInterpolator.html [17] P. Sonneveld. Errors in cubic spline interpolation. Journal of Engineering Mathematics, 3(2):107{111, 1969. [18] Carl de Boor. A Practical Guide to Splines. Springer-Verlag New York, 1978. [19] SciPy, 2016. https://docs.scipy.org/doc/scipy-0.18.1/reference/generated/scipy.interpolate.CubicSpline.html [20] SciPy, 2014. https://docs.scipy.org/doc/scipy-0.16.0/reference/generated/scipy.interpolate.PPoly.html [21] Numpy, 2014. https://docs.scipy.org/doc/numpy/reference/generated/numpy.searchsorted.html 24 [22] Oded Goldreich. Computational complexity: A conceptual perspective. Cambridge University Press, 1 edition, 2008. [23] Paul Zimmermann Richard Brent. Modern Computer Arithmetic. Cambridge University Press, 2010. [24] Alvaro Sanchez-Gonzalez, 2016. https://github.com/alvarosg/pynverse [25] SciPy, 2014. https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.special.lambertw.html http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Astrophysics arXiv (Cornell University)

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

Tommasini, Daniele; Olivieri, David N.
Astrophysics , Volume 2018 (1812) – Dec 5, 2018
Download PDF
25 pages

Loading next page...
 
/lp/arxiv-cornell-university/fast-switch-and-spline-scheme-for-accurate-inversion-of-nonlinear-xqGyVwf7bi

References

References for this paper are not available at this time. We will be adding them shortly, thank you for your patience.

ISSN
0096-3003
eISSN
ARCH-3330
DOI
10.1016/j.amc.2019.124677
Publisher site
See Article on Publisher Site

Abstract

Numerically obtaining the inverse of a function is a common task for many scienti c 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 speci c 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 de ned 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  10 factor for the computation of W in the range W (y) 2 [0; 10], and by a factor 2  10 for Kepler's problem. In our tests, this 4 3 scheme is much faster than Newton-Raphson's method, by a factor in the range 10 to 10 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 rad, the algorithm outperforms Newton's method for all values of the number of points N  2. arXiv:1812.02273v2 [physics.comp-ph] 17 Dec 2018 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 [1, 2], the computation of Lambert W function [3, 4], and the solution of Kepler's equation for the orbital motion of a body in a gravitational eld [5, 6]. 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 speci c 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 [7] or some of its variants [8{11]. 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 rst 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 [12{14]. 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 nding the zeros of the functions F (x)  f (x) y. If a good initial guess x of the true y 0 value x = f (y) of the zero is available, the zeros of F can be computed by recursively y f (x ) F (x ) k applying the equation x = x ; i.e. x = x + . k+1 k 0 k+1 k F (x ) y k f (x ) 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 speci c spline routine that makes the FSSI more accurate and much 2 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 minimiza- tion methods) for computing values of inverse functions, especially if a good rst approxima- tion is dicult to obtain. Moreover, the FSSI is shown to be superior to Newton-Raphson when the values f (y) of the inverse function are required for many di erent 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 di erent instants, depending on the eccentricity e and the re- quested precision. For example, for e = 0:8 and accuracy  10 rad, the FSSI algorithm is already faster than Newton's even when the computation is done at N = 2 points, and 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 2 [x ; x ]. The function min max 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 (y) in the co-domain. The FSSI consists of a two step approach. First, when the input function is given an- alytically, the values of f on a given grid of points x , for j = 1; ; n, are computed to obtain the matrix (x ; y ), where y = f (x ). The matrix (y ; x ), obtained by switching the j j j j j j arrays, gives then the exact values g(y ) = x of the inverse function on the grid y . From j j j this modi ed matrix, the cubic spline interpolant S(y) of (y ; g(y )) is computed by using j j 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 speci ed on a grid by a set of tuples (x ; y ). In a high-level programming language such as Python, this tuple array is represented j j as (x; y), and the FSSI interpolant S(y) can be obtained by calling a cubic spline routine of 3 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 . Input Output analytic function Vector values and OR Domain of inverse: Error Speed: Evaluate Switch domains for spline Calculate coefficients for spline of Calculate breakpoints Evaluate on co-domain Build polynomial on grid Figure 1: Flow diagram of the FSSI method for obtaining the function inverse f . The diagram indicates the key steps of the method, as well as how it is interfaced to the input and output. Figure 1 shows the ow 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 () and its derivative f (), 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 is determined by the number of points n + 1 of the input grid, and which is used for sampling N points fY ; ; Y g in the function's co-domain. In subsequent sections, we 1 N show that the FSSI algorithm has an error bound proportional to jkxj , for k a constant (described in the text), and with a maximal time complexity of O(1)+ 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 . This is the case for Kepler's equation and for all the other examples that we 4 0 will consider hereafter. Therefore, we can pro t from the knowledge of f to design a speci c 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 di erence 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 coecients of the spline [7], where n is the number of grid intervals. An exception is Akima's cubic spline [15, 16], which is fast, diagonal and regular. The speci c spline that is designed here is based on a similar idea to Akima's, but it is signi cantly 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 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 di erent in the 0 0 problem of the inversion of a function f (x). In this case, the derivatives g (y ) = 1=f (x ) j j can be given on a grid, while the values g(y) are not known. Let us build the speci c spline S(y) piecewise in each interval, S(y) = S (y) for y < y < j j y , where j takes the values j = 0; ; n 1. If we de ne the arrays y  (y ; ; y ) j+1 0 0 n1 and y  (y ; ; y ), obtained by removing the last and the rst point of the y array, 1 1 n respectively, then the j-th interval can also be written as y < y < y . In this interval, the 0 1 j j cubic spline can be de ned as S (y) = c y y ; (1) j q 0 j j q=0 where for each value of q = 1; ; 4 the coecients c can also be thought as the n compo- nents of an array c . Since the values of the derivative f (x ) of the input function are known on the grid points, we can construct an array d whose n components are the derivatives of the inverse function g on the points y = f (x ), j j d  g (y ) = ; for j = 0; ; n: (2) j j f (x ) As we did for y, which was used to generate the arrays y and y by removing one end 1 2 point, it is convenient to de ne similar arrays of n components also from x and d, namely 5 x  (x ; ; x ), x  (x ; ; x ), d  (d ; ; d ), and d  (d ; ; d ). With 0 0 n1 1 1 n 0 0 n1 1 1 n this convention, we have to choose the spline coecients that lead to the best approximation of the inverse function g(y). The most natural choice is to force S to coincide with the known values of the inverse function, x and x , at the end points, and to ask the same for the 0 1 j j derivatives d and d . In other words, the conditions to be imposed in each interval are, 0 1 j j 0 0 S (y ) = x ; S (y ) = x ; S (y ) = d ; S (y ) = d ; (3) j 0 0 j 1 1 0 0 1 1 j j j j j j j j j j where S (y) is given by equation 1. For every xed value of j, these conditions give a system of four equations for the four unknown coecients c , c , c and c , which is decoupled from the similar systems of 0 1 2 3 j j j j equations corresponding to di erent 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 coecients [7], 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 xed j can be solved analytically in a straight- forward 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 vw+2z of components in such a way that an equality like, e.g., u = between vectors having v w +2z j j j the same number of elements represents the equations u = for every j. With this convention, equations (3) give the solution c = x ; 0 0 c = d ; 1 0 (2 d + d ) (y y ) 3 (x x ) 0 1 0 1 0 1 c = ; (y y ) 0 1 (d + d ) (y y ) 2 (x x ) 0 1 0 1 0 1 c = : (4) (y y ) 0 1 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 00 00 adding the additional conditions S (y ) = d and S (y ) = d to equation (3) we can also 0 0 1 1 j j j j j j 6 design a diagonal quintic spline, and if also the third derivatives are known the additional 000 000 conditions S (y ) = d and S (y ) = d allow for the construction of a septic spline. We 0 0 1 1 j j j j j j 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 speci c 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 fth 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 di erences 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. A. Derivation of an upper bound on the error of the FSSI by using the known analysis of cubic spline interpolation Following Ref. [7, 17], 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 jg(y) S(y)j  M ; (5) (4) 4 where M = max g (y) , and  = max (y y ) . j+1 j y yy 0jn1 max min 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, 00 3 (3) 00 (4) 15f (x) 10f (x)f (x) f (x) M = max + : (6) 0 7 0 6 0 5 x xx 0 n f (x) f (x) f (x) Equations (5) and (6), along with the de nition of , can be used to obtain an upper limit on the error. As shown in examples below, the actual errors are several orders of 7 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. B. Ab initio derivation of an improved estimation of the error for the FSSI Let us assume that the function g(y) is in nitely di erentiable. Thus, it can be expanded P (q) 1 g (y ) in a Taylor series g(y) = (yy ) around one of the points of the grid y = f (x ), j j j q=0 q! chosen to be the closest grid point to y, so that jyy j  jyy j and jyy j  jyy j. j j+1 j j1 We also assume that the cubic spline interpolation is made in such a way so that it is equivalent to the rst terms of this Taylor expansion, expanded around the same point and P (q) 3 g (y ) truncated beyond cubic order, S(y) = (y y ) . Therefore, the di erence is q=0 q! (q) (4) (4) g (y ) g (y  ) g (y ) j j j n 4 4 jg(y) S(y)j = (y y ) = (y y )  (y y ) ; (7) j j j n! 4! 4! q=4 where y  is an unknown intermediate point between y and y , and the last approximation is j j expected to hold for suciently small values of jyy j, which is the case when a suciently high number of grid points is chosen. On one hand, the exact equality in equation (7) can be translated in the following bound, (4) g (y ) y y 1 j+1 j jg(y) S(y)j  max max = M ; (8) y y y 0jn1 min max 4! 2 384 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 y is due to the j+1 j y y j+1 j fact that y was chosen as the closest grid point to y, so that jy y j  . j j On the other hand, equation (7) can be elaborated further and expressed in terms of the function f (x), 00 3 (3) 00 (4) 1 15f (x ) 10f (x )f (x ) f (x ) j j j j 4 jg(y) S(y)j  + [f (x )(x x )] ; j j 0 7 0 6 0 5 4! f (x ) f (x ) f (x ) j j j where x = g(y) and y y  f (x )(x x ). The last approximation is expected to be j j j accurate over the entire interval provided the following condition f (x) max jx x j max  1 (9) j+1 j 0jn1 x xx n 2f (x) 8 is satis ed. x x j+1 j Assuming that the point x is the one closest to x, so that jx x j  max , j j 0jn1 the following estimation is obtained jg(y) S(y)j / max jx x j  (10) j+1 j 0jn1 00 3 (3) 00 (4) 15f (x) 10f (x)f (x) f (x) max + : 0 3 0 2 0 x xx 0 n f (x) f (x) f (x) As described in the next section, this error approximation is usually much smaller, pos- sibly 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 satis es 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 grid, we can substitute 4 (x x ) x x n 0 n 0 max jx x j = in the condition (9), and max jx x j = in the j+1 j j+1 j n n 0jn1 0jn1 estimation of the error (10). Not only is this the simplest choice, if the grid is not given 4 (x x ) n 0 otherwise, but it is usually also the best option. In fact, when jx x j = for j+1 j 4 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), de ned over a domain x  x  x . For the numerical computation, we developed a python 0 n 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 rst example, in which the exact inverse function g(y) is known, we also provide an additional evaluation of the error by computing the di erence S(y)g(y). In both cases, we use a grid Y that contains 10 times as many points as the original grid y = f (x ). In fact, by construction S(y) is exactly equal j j to g(y){within machine errors{over the original grid y , so it is important to ensure that S(y) is compared with g(y) in between the grid points y . It is true that even if they were equal in number the points Y , chosen to be equally spaced, would not coincide in general 9 0 with the y , whose spacing is variable and roughly proportional to f (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 di erence 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 speci c 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 [18, 19]. 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 speci c spline. An exception is Akima routine [16], which is less accurate by three orders of magnitude in the bulk of the interval. This is an additional reason for preferring our speci c spline, besides the fact that it is the fastest one. A. Exponential The rst 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 , from the bounds expressions of (5) and (6), can be readily computed. The results are, 4x 4x M = max 6e = 6e ; (11) x xx 0 n and 4 4 x x x x n 0 n 0 4x 4x = max e = e ; (12) x xx n 0 n n so that the bound on the error as computed from equations (5) and (6) is 6 x x n 0 4(x x ) n 0 jg(y) S(y)j  e : (13) 384 n 10 On the other hand, the analytic estimation we derived in Equation (10) gives 6 x x n 0 jg(y) S(y)j / : (14) 384 n 4(x x ) n 0 Therefore, this analytic error approximation of FSSI is smaller by a factor e as compared to the limit that was derived in equation (13) by applying the standard cubic spline error bound. For example, if x x = 10, then our error estimation is a factor exp(40), n 0 i.e. 17 orders of magnitude, smaller than what could be expected from the literature. In order to bene t by this accuracy improvement, the grid must be chosen in such a way that the condition (9) is satis ed, i.e. x x f (x) x x n 0 n 0 =  1: (15) n 2f (x) 2n 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 2 [0; 10], the condition becomes n  5, so that values of n of the order of 50 or larger could be sucient. This is also what we have observed by performing numerical computations for di erent values of n. In general, for n & 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 2 [0; 10] using n = 10 grid points. In this case, our theoretical prediction of Equation (14) gives jg(y) S(y)j / 1:6  10 , which is in excellent agreement with the numerical computation over the entire interval. The results for this case also con rm 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. B. 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) [3, 4]. 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 11 Figure 2: Result of the FSSI applied to the function f (x) = exp(x) (top left) over the domain x 2 [0; 10]. The FSSI interpolant S(y) is shown for n = 10 grid points (top right), together with 1 1 two independent evaluations of the numerical errors: i) jS(y) f (y)j, where f (y) = ln(y) (bottom left); ii) jS(y) S(f (S(y)))j (bottom right). of M and  from equations (5) and (6) for the theoretical bound are, x x x x x x x x 15 (e x + 2e ) 10 (e x + 3e ) (e x + 2e ) e x + 4e M = max + ; (16) 7 6 5 x x x x x x x xx 0 n (e x + e ) (e x + e ) (e x + e ) which is a monotonically decreasing function for x > 1, and x x n 0 x x 4 = max je x + e j ; (17) x xx n 0 n which increases monotonically. Therefore, the bound (5) becomes jg(y) S(y)j  (x = x ) M (x = x ): (18) n 0 20 x x n 0 For example, over the domain x 2 [0; 10] this gives jg(y) S(y)j / 5:7 10 . On the other hand, our analytical estimation from Equation (10) becomes 1 x x n 0 jg(y) S(y)j /  (19) 384 n x x x x x x x x 15 (e x + 2e ) 10 (e x + 3e ) (e x + 2e ) e x + 4e max + : 3 2 x x x x x x x xx 0 n e x + e (e x + e ) (e x + e ) 12 The function to be maximized in equation (19) monotonically decreases for x > 1, so that its maximum is achieved for x = x . Thus, over the domain x 2 [0; 10] our estimate of x x n 22 the error (19) gives jg(y) S(y)j / 0:17 , which is 3 10 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 x x x + 2 0 n n  max = 10; (20) x xx 2 n x + 1 which is a surprisingly low value, for such a huge variation of f . Figure 3: Numerical result of the FSSI applied to the function f (x) = x exp(x) (top left) over the domain x 2 [0; 10]. The FSSI interpolant S(y) is shown for n = 10 grid points (top right), together 1 1 with two independent evaluations of the numerical errors: i) jS(y)f (y)j, where f (y) = W (y) as computed with scipy.special.lambertw routine (bottom left); ii) jS(y) S(f (S(y)))j (bottom right). Figure 3 shows the numerical result of the FSSI inversion of f (x) = x exp(x) in the domain x 2 [0; 10] using n = 10 grid points. Two independent evaluations of the nu- 1 1 merical errors are given: i) jS(y) f (y)j, where f (y) = W (y) as computed with scipy.special.lambertw routine; ii) jS(y) S(f (S(y)))j. The fact that they agree with each other provides con rmation concerning our treatments of the errors. Moreover, in this case 13 5 our theoretical prediction of Equation (19) gives jg(y) S(y)j / 1:7  10 , and the nu- merical error not only agrees with it, but it is even much smaller, by almost an order of magnitude, jg(y) S(y)j < 2:5 10 over the entire interval. In this case, FSSI is numerical more accurate by an astonishing factor 4  10 than what could be expected by naively applying the general results for cubic splines, as in equation (18). C. Kepler's equation Kepler's equation for an elliptical orbital motion of eccentricity e can be written as y = x e sin x; (21) where y and x represent the so-called mean and eccentric anomaly, respectively [5, 6]. The 2t former is the time elapsed since periapsis, as measured in radians, y = , where T is the period of the orbit. The eccentric anomaly x is related to the angle  between the position vectors at periapsis and at time t, with origin in the center of gravity, through the equation 1 + e x = 2 arctan tan : (22) 1 e 2 A fundamental problem in orbital dynamics [5, 6] is to obtain the time dependence of the angle  describing the position of the orbiting body at time t, which requires the inversion of the function y = f (x)  x e sin x. Taking into account that the orbit is periodic, and that for x 2 [; 2] we have f (x) = 2 f (2 x), it is sucient to consider only the interval x 2 [0; ] to obtain the behavior for all values of x. The corresponding co-domain is then y 2 [0; ] [5, 6]. The inverse function x = g(y) will yield the eccentric anomaly as a function of the mean anomaly, and thus the evolution (t) will be obtained. This is usually done in an ecient way using Newton's method with the rst guess x = y + e=2 [5, 6, 8]. Here, FSSI is considered as an alternative to Newton-based methods for solving Kepler's equation. In this case, the values of M and  in the theoretical bound of equations (5) and (6) are 3 2 15e sin x 10e sin x cos x e sin x M = max + + (23) 7 6 5 0x (1 e cos x) (1 e cos x) (1 e cos x) and = max j1 e cos xj : (24) 0x 14 As a concrete example, the case of e = 0:8 is considered. Thus, the maximum values 4 4 are M = 10275:1, which is obtained for x = 0:166, and  = j1 + ecj = 10:4976 , n n obtained for x = . Therefore the bound (5) becomes 2:7 10 jg(y) S(y)j . : (25) However, the expression from our analytic estimation from Equation (10) becomes 3 2 1  15e sin x 10e sin x cos x e sin x jg(y) S(y)j / max + + : (26) 3 2 0x 384 n (1 e cos x) (1 e cos x) 1 e cos x For e = 0:8, the expression in the j j bracket has a maximum value 21.586 obtained for x=0.214657, therefore we obtain 5:5 jg(y) S(y)j / : (27) Thus, our estimation for the theoretical error (27) in this case is 2  10 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 e sin x n  max ' 2: (28) 2 0x 1 e cos x 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 2 [0; ] using n = 10 and n = 10 grid points, respectively. In these 4 8 cases, our theoretical prediction of Equation (27) gives jg(y) S(y)j / 5:5 (10 or 10 ), 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 di erence of the FSSI interpolation with the values of g (y) obtained with Newton's method, and the other given by the di erence S(y) S(f (S(y))). The fact that these evaluations of the error also agree with each other is a further con rmation of the validity of our error analysis. By comparing gures 4 and 5, we also see that the accuracy scales with n , as equation (27) predicts, and that our estimation for the error is reliable even for just n = 10 grid points. 15 Figure 4: Numerical result of the FSSI applied to the function f (x) = x 0:8 sin x over the domain x 2 [0; ] (top left), corresponding to Kepler's equation for an elliptical orbit of eccentricity 0.8. The FSSI interpolant S(y) is shown for n = 10 grid points (top right), together with two independent evaluations of the numerical errors: i) jS(y) g (y)j, where g (y) is computed with Newton's N N method (bottom left); ii) jS(y) S(f (S(y)))j (bottom right). VI. NUMERICAL COMPARISONS WITH NEWTON-BASED METHODS Apart from the numerical calculations for error analysis, we carried out numerical com- parisons 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 imple- mented 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 di erence 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, 16 Figure 5: Numerical result of the FSSI applied to the function f (x) = x 0:8 sin x over the domain x 2 [0; ] (top left), corresponding to Kepler's equation for an elliptical orbit of eccentricity 0.8. The FSSI interpolant S(y) is shown for n = 10 grid points (top right), together with two independent evaluations of the numerical errors: i) jS(y) g (y)j, where g (y) is computed with Newton's N N method (bottom left); ii) jS(y) S(f (S(y)))j (bottom right). the di erence in the execution times being usually below  30% [11]. 17 Algorithm 1 Benchmark for FSSI 1: procedure bench FSSI(Y; x; f (x); f (x)) 2: y = f (x) 1 Algorithm 2 Benchmark for Newton 3: d = f (x) 1: procedure 4: c = x[: 1] bench Newton(Y; f (x); f (x); tol) 5: c = d[: 1] 2: for Y in Y do 6: d1 = d[1 :] 3: X = g (Y ) k 0 k 7: xx = c x[1 :] Y f (X ) k k 4:  = f (X ) 8: y0 = y[: 1] 5: while  > tol do 9: y1 = y[1 :] Y f (X ) k k 6: = f (X ) 10: yy = y0 y1 7: X = X + k k 11: yy2 = yy yy 8:  = j j 12: yd1 = yy d1 9: end while 13: yd0 = yy c1 10: end for 2yd0+yd13xx 14: c = yy2 11: return X . =f (Y) yd0+yd12xx 15: c = yy2yy 12: end procedure 16: call P: X = P((c ; c ; c ; c ); Y) 3 2 1 0 17: return X . =f (Y) 18: end procedure In algorithm 2, we have called g (Y ) the initial guess for Newton's method as a function 0 k 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 fi = 0; ng or over fi = 1; ng, for lowercase arrays, or over fi = 1; Ng, for uppercase arrays. An exception are the expressions v[1 :] and v[: 1], which mean the removal of the rst or the last element from v, respectively. In Python, the piecewise polynomial function P , corresponding to equation (1), can be obtained in terms of the breakpoints and the coecients c using the subroutine PPoly [20], so that P = PPoly. Another possibility is to write an explicit subroutine for com- puting the polynomial, in which the insertion points j are located by binary search using 18 scipy.searchsorted [21]. The two possibilities are shown below: Subroutine P for FSSI Search in Python 1: function P((c ; c ; c ; c ); Y) 3 2 1 0 Subroutine P for FSSI PPoly in Python 2: j = numpy.searchsorted(y1; Y) 1: function P((c ; c ; c ; c ); Y ) 3: P1 = Y y 3 2 1 0 0 2: P = scipy.PPoly 4: P2 = P1 P1 3: end function 5: S = c +c P1+c P2+c P2P1 o 1 2 3 j j j j 6: return S 7: end function A discrete analysis of the algorithms 1 and 2 shows that the FSSI executes in constant time O(1), because once the spline coecients 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, nite cache sizes and the search of the breakpoints overtake this behavior, so that the algorithm follows a linear time dependence O(N ) [22, 23]. In other words, the execution time can be written as t ' N +  and FSSI PPoly t ' N + , for the python implementations of FSSI with PPoly or Search- FSSI Search sorted, respectively. As we show below, t < t for large N , typically FSSI PPoly FSSI Search N & 10 , and t > t for lower values of N . We can then merge the two FSSI PPoly FSSI Search python routines for P, algorithms FSSI Search and FSSI PPoly, by choosing the fastest one with an if statement, e.g. if N > 10 do PPoly, else do the routine with searchsorted. The execution time for this combined routine is t ' + N , i.e. it behaves as O(1) + 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 O(mN ) for all values of N . To obtain an empirical execution time comparison between methods, we ran the bench- marks 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 19 (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 3 4 Lambert W problem, we used a tolerance 2 10 =n for scipy.lambertw and Newton; For the Kepler solution, we used a tolerance 6=n for Newton and Pynverse [24] quasi- Newton method. x +x 0 n For Lambert W, we chose the simplest rst guess, g (Y ) = = 5. Of course, 0 k better choices may be found, but we want to use this case to compare FSSI and Newton in the absence of a good rst 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 e = 0:8 and we use a very good rst guess, g (Y ) = Y + , as was mentioned in section V. 0 k k Figure 6 shows empirical execution time comparisons between di erent numerical algo- rithms 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 O(N ) behavior is negligible as compared with the O(1) behavior; however for very large N , when the O(N ) part dominates, the linear coecient is several orders of magnitude smaller than those of the other methods available. In all the cases, Pynverse [24] (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 gure 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 ts to the data. For example, for Kepler's problem with n = 50, corresponding to tolerance 10 rad (which can be a sucient accuracy 8 4 for orbit determination in many cases), we found t ' 2:1 10 N + 1:9 10 FSSI PPoly 8 5 and t ' 5:3  10 N + 7:2  10 . These values of the coecients have been FSSI Search obtained by separate ts to the low N data, for  and , and to the high N data, for  and 20 Newton FSSI_s Newton pynverse FSSI_s FSSI_p FSSI_p pynverse Newton Newton FSSI_s FSSI_p FSSI_s FSSI_p N evaluations N evaluations FSSI_s FSSI_s FSSI_s FSSI_s Newton Newton FSSI_p FSSI_p scipy.lambertw scipy.lambertw scipy.lambertw Newton Newton scipy.lambertw FSSI_p FSSI_p N evaluations N evaluations Figure 6: Numerical comparisons of FSSI and other methods for the solution of Kepler's equation (top) and computation of the Lambert W function (bottom). FSSI p and FSSI s stand for the algorithm using PPoly or Searchsorted subroutines, respectively. , in order to get the best estimates in these regimes, so that the approximation is slightly 3 4 worse for 10 . N . 10 . In any case, as shown in gure 6, FSSI PPoly is faster than FSSI Search for N & 10 and slower for N . 10 . By choosing the best of the two variants with an if statement, we 8 5 obtain a combined behavior t ' 2:1 10 N + 7:2 10 . This should be compared FSSI with the execution time for Newton-Raphson method, t = 4:210 N . We nd that Newton t > t for every N  2, and that FSSI is  2 10 faster than Newton-Raphson Newton FSSI for large N . 4 16 Similarly, for Kepler's problem with n = 10 , corresponding to tolerance 6 10 rad, we 8 3 5 obtain a combined behavior t ' 2:110 N +1:110 while t = 5:010 N . FSSI Newton We nd that t > t for every N & 20, and that FSSI is still  2 10 faster Newton FSSI than Newton-Raphson for large N . For Lambert W with n = 50, corresponding to tolerance 3 10 , we obtain a combined 8 5 4 behavior t ' 2:1 10 N + 6:2 10 while t = 1:3 10 N . We nd that FSSI Newton t > t for every N , and that FSSI is  6 10 faster than Newton-Raphson for Newton FSSI Execution Time (sec) Execution Time (sec) Execution Time (sec) Execution Time (sec) large N . This shows that, in the lack of a good rst guess, FSSI can be better than Newton- Raphson method for every value of N . Of course, for small N the speci c, semi-analytic routine scipy.lambertw [25], having t = 3:8 10 N outperforms the FSSI, scipy LambertW but surprisingly the opposite is true in the large N regime, in which FSSI is  20 times faster than scipy.lambertw. 4 13 Finally, for Lambert W with n = 10 , corresponding to tolerance 2  10 , we obtain 8 3 4 a combined behavior t ' 2:1  10 N + 1:1  10 while t = 1:3  10 N . FSSI Newton We nd that t > t for every N & 8, and that FSSI is  7  10 faster than Newton FSSI Newton-Raphson for large N . This shows that, in the lack of a good rst guess, FSSI is much better than Newton-Raphson method for every value of N . Again, for large N , t < t = 3:8 10 N by a factor  20. FSSI scipy LambertW 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 di erent 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 expres- sions 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, 22 4 by a  10 factor for W 2 [0; 10] or by a factor 2 10 for Kepler problem. This scheme is superior to, and much faster than, Newton-Raphson method when the latter is dicult to apply, when no good rst guess is available, or when the values of 22 the inverse function are required on an entire interval or in a large number of di erent points. When applied to Kepler's problem (e.g. with eccentricity e = 0:8), FSSI becomes faster than Newton's methods for N greater than a few points, and is  2 10 times faster for large N . If the requested accuracy is of the order of 10 rad, which is a low enough value for most applications, the speed of the FSSI algorithm is faster than Newton's for N  2. The N dependence of the scheme can be described as O(1) +O(N ). For a wide range of N values, the O(N ) behavior is negligible as compared with the O(1) behavior; however for very large N , when the O(N ) part dominates, the linear coecient  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 rst choice for solving Kepler's equation. VIII. ACKNOWLEDGEMENTS We thank A. Paredes, D. Gonz alez-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 on e Ordenaci on Universitaria (Xunta de Galicia). [1] Toshio Fukushima. Numerical computation of inverse complete elliptic integrals of rst and second kinds. Journal of Computational and Applied Mathematics, 249:37{50, 2013. [2] John P. Boyd. Four ways to compute the inverse of the complete elliptic integral of the rst kind. Computer Physics Communications, 196:13{18, 2015. [3] R. M. Corless, G. H. Gonnet, D. E. G. Hare, D. J. Je rey, and D. E. Knuth. On the lambert w function. Advances in Computational Mathematics, 5:329{359, 1996. [4] Darko Veberi c. Lambert w function for applications in physics. Computer Physics Communications, 183:2622{2628, 2012. 23 [5] J. E. Prussing and B. A. Conway. Orbital Mechanics. Oxford University Press, 2 edition, [6] Howard D. Curtis. Orbital Mechanics for Engineering Students. Elsevier, 3 edition, 2014. [7] John H. Mathews and Kurtis D. Fink. Numerical Methods Using MATLAB. Pearson, 4 edition, 2004. [8] J.M.A. Danby and T.M. Burkardt. The solution of kepler's equation, i. Celestial Mechanics, 31:95{107, 1983. [9] J.M.A. Danby and T.M. Burkardt. The solution of kepler's equation, iii. Celestial Mechanics, 40:303{312, 1987. [10] J. Gerlach. Accelerated convergence in newton's method. SIAM Rev., 36:272{276, 1994. [11] M. Palacios. Kepler equation and accelerated newton method. Journal of Computational and Applied Mathematics, 138:335{346, 2002. [12] B. A. Conway. An improved algorithm due to laguerre for the solution of kepler's equation. Celestial Mechanics, 39:199{211, 1986. [13] E. D. Charles and J. B. Tatum. The convergence of newton{raphson iteration with kepler's equation. Celestial Mechanics and Dynamical Astronomy, 69:357{372, 1998. [14] Laura Stumpf. Chaotic behaviour in the newton iterative function associated with kepler's equation. Celestial Mechanics and Dynamical Astronomy, 74:95{109, 1999. [15] Hiroshi Akima. A new method of interpolation and smooth curve tting based on local procedures. Journal of the Association for Computing Machinery, 17(4):589{602, 1970. [16] SciPy, 2016. https://docs.scipy.org/doc/scipy-0.18.1/reference/generated/scipy.interpolate.Akima1DInterpolator.html [17] P. Sonneveld. Errors in cubic spline interpolation. Journal of Engineering Mathematics, 3(2):107{111, 1969. [18] Carl de Boor. A Practical Guide to Splines. Springer-Verlag New York, 1978. [19] SciPy, 2016. https://docs.scipy.org/doc/scipy-0.18.1/reference/generated/scipy.interpolate.CubicSpline.html [20] SciPy, 2014. https://docs.scipy.org/doc/scipy-0.16.0/reference/generated/scipy.interpolate.PPoly.html [21] Numpy, 2014. https://docs.scipy.org/doc/numpy/reference/generated/numpy.searchsorted.html 24 [22] Oded Goldreich. Computational complexity: A conceptual perspective. Cambridge University Press, 1 edition, 2008. [23] Paul Zimmermann Richard Brent. Modern Computer Arithmetic. Cambridge University Press, 2010. [24] Alvaro Sanchez-Gonzalez, 2016. https://github.com/alvarosg/pynverse [25] SciPy, 2014. https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.special.lambertw.html

Journal

Astrophysics – arXiv (Cornell University)

Published: Dec 5, 2018

References

  • Numerical computation of inverse complete elliptic integrals of first and second kinds
    Fukushima
  • Four ways to compute the inverse of the complete elliptic integral of the first kind
    Boyd
  • On the lambert w function
    Corless
  • Precise and fast computation of lambert w-functions without transcendental function evaluations
    Fukushima
  • Lambert w function for applications in physics
    Veberič