Numerical inversion of a general incomplete elliptic integral

摘要

We present a numerical method to invert a general incomplete elliptic integral with respect to its argument and/or amplitude. The method obtains a solution by bisection accelerated by half argument formulas and addition theorems to evaluate the incomplete elliptic integrals and Jacobian elliptic functions required in the process. If faster execution is desirable at the cost of complexity of the algorithm, the sequence of bisection is switched to allow an improvement by using Newton’s method, Halley’s method, or higher-order Schröder methods. In the improvement process, the elliptic integrals and functions are computed by using Maclaurin series expansion and addition theorems based on the values obtained at the end of the bisection. Also, the derivatives of the elliptic integrals and functions are recursively evaluated from their values. By adopting 0.2 as the critical value of the length of the solution interval to shift to the improvement process, we suppress the expected number of bisections to be as low as four on average. The typical number of applications of update formulas in the double precision environment is three for Newton’s method, and two for Halley’s method or higher-order Schröder methods. Whether the improvement process is added or not, our method requires none of the procedures to compute the incomplete elliptic integrals and Jacobian elliptic functions but only those to evaluate the complete elliptic integrals once at the beginning. As a result, it runs fairly quickly in general. For example, when using the improvement process, it is around 2–5 times faster than Newton’s method using Boyd’s starter (Boyd (2012) [25]) in inverting E(φ|m), Legendre’s incomplete elliptic integral of the second kind.