2

Suppose I have given the value of a complete elliptic integral of the second kind $E(\varphi | m)$, with a positive value $z$, and the value $m$. How can I solve for $\varphi$ in the following equation? \begin{equation} E(\varphi | m) = z \end{equation} I am aware that there are already answers to this question here or here. However as someone who has only a surface level understanding of math, these are very confuseing to me.

I would like to have an approximate result, dependent on $m$ and $z$, using a function or series expansion that I can then use in Matlab or Python. Thanks in advance!

haifisch123
  • 453

1 Answers1

4

The incomplete elliptic integral of the 2nd kind is $$E(\varphi | m) = \int_0^\varphi \sqrt{1-m\sin^2(\theta)} \, d\theta$$

The complete integral is defined as $E(\frac\pi2 | m)$.

There doesn't seem to be a simple algorithm to invert this integral. But I wouldn't be surprised if there is such an algorithm using the arithmetic-geometric mean (AGM), which can be used to rapidly evaluate all the elliptic integrals and the Jacobi elliptic functions. The AGM converges quadratically, so these integrals can be swiftly calculated to high precision. The elliptic integrals and functions are not usually supplied in standard math libraries, but they are found in advanced math libraries. And they are easy enough to code in only a few lines (using the AGM).

However, we can easily invert this elliptic integral using Newton's method. Of course, this assumes that we can actually evaluate the integral itself. ;)

Newton's method says that given an initial approximation $x_0$ for the zero of a function $f(x)$, we can often find a better approximation $x_1$:

$$x_1 = x_0 - \frac{f(x_0)}{f'(x_0)}$$ where $f'(x)$ is the derivative of $f(x)$, with respect to $x$.

Newton's method isn't perfect. We must supply an initial approximation, preferably one near a root of $f$. If the derivative is small near that root, eg, in the neighbourhood of a double root, then Bad Things happen. Please see the Wikipedia article for further details.

Fortunately, this elliptic integral is well-behaved for typical values of the parameters, eg when $\varphi$ is real and when $|m|\le 1$. We don't need to worry too much about finding a good initial approximation: even zero can work, although it does save time if we do have a better initial value.

Let $u = E(\varphi | m)$. Given an initial approximation $\varphi_0$,

$$\varphi_{i+1} = \varphi_i - \frac{u - E(\varphi_i | m)}{\sqrt{1-m\sin^2(\varphi_i)}}$$

Once Newton's method gets close to the root, it converges quadratically. That is, the number of correct bits / decimals doubles on each iteration.

Here's a brief demo in Sage, a free mathematics software system based on Python. This code uses arbitrary precision arithmetic, but the syntax in this demo is mostly identical to plain Python. Note that Sage permits the use of ^ for exponentiation.

Code

""" Invert an elliptic integral of the 2nd kind
    Uses Newton's method
    Written by PM 2Ring 2024.01.21
"""

bits = 80
RF = RealField(bits)
eps = RF.epsilon().sqrt()
print(f"{bits=}, {eps=:.3e}")


Invert an elliptic integral of the 2nd kind


phi is the initial approximation


def invert_elliptic_e(u, m, phi, eps=1e-12):
    for i in range(99):
        q = elliptic_e(phi, m)
        dphi = (u - q) / sqrt(1 - m * sin(phi)^2)
        phi += dphi
        print(i, phi, dphi)
        if abs(dphi) < eps:
            break
    return phi


Test


m = RF(0.7)
phi = RF(1.234)
print(f"   {phi=}")


u = elliptic_e(phi, m)
print(f" {u=}")


Initial approximation


Good for small phi


#approx1 = u + 1/6 * m * u^3 + 1/120 * m(13m - 4) * u^5


Interpolate from complete integral


alpha = RF(pi/2)
uc = elliptic_e(alpha, m)
print(f"{uc=}")
approx = alpha * u / uc
print(f"{approx=}")


newphi = invert_elliptic_e(u, m, approx, eps=eps)
print(f"{newphi=}")

Output

bits=80, eps=1.286e-12
   phi=1.2340000000000000000000
 u=1.0495259481188760945322
uc=1.2416705679458227508715
approx=1.3277205296960812349754
0 1.2316872646292563899235 -0.096033265066824845051867
1 1.2339984464246334166675 0.0023111817953770267439825
2 1.2339999999993000845721 1.5535746666679046024772e-6
3 1.2340000000000000000000 6.9991542788168335651453e-13
newphi=1.2340000000000000000000

In this demo, the initial approximation is a simple linear interpolation based on the complete integral. For small $\varphi$, $u$ itself is a good approximation, especially when $m$ is small. Another option is to use a few terms of the series given by zeraoulia rafik on MathOverflow. It performs very well for small $\varphi$, but it over-estimates for larger $\varphi$. I have a few terms of that series in my code, but I commented it out.

Here's a live version of my demo, running on the SageMathCell server.

If you don't want to use Sage, another option in standard Python for advanced mathematics with arbitrary precision is mpmath.

Incidentally, this elliptic integral arises in the computation of Geodesics on an ellipsoid, and Wikipedia has a lot of useful info and links on this topic. I found the articles by C.F.F. Karney particularly helpful (and Karney is a major contributor to those Wikipedia articles).

PM 2Ring
  • 5,544