Landen transformation for incomplete elliptic integral approximation

Question

While reading about Landen trandformation in the book of Murray R. Spiegel: Mathematical Handbook of Formulas and Tables

The formulas from the book

I tried to program that in Python as follows:

from numpy import pi, log, tan, sqrt
import numpy as np
from scipy.special import ellipk

def k(n, k0=1/sqrt(2)): 
     r = k0 
     if n == 0: 
         return k0 
     for i in range(1,n+1): 
         r = 2*np.sqrt(r)/(1.0+r) 
     return r 

k0 = 1/sqrt(2)
phi = pi/2

N = 10

product = np.prod([k(i) for i in range(1, N+1)])

approx_result = sqrt(product/k0)*log(tan(pi/4 + phi/2))

expected_result = ellipk(k0)

print("Calculated result", approx_result)
print("Expected result:", expected_result)

I got the following:

Calculated result 44.06430550154578
Expected result: 2.0859736981949366

As you can see, the calculated result using Landen transformation is wrong. Could you please tell me what I am missing here?

I appreciate any help.

Answer 1

There is an unfortunate error in your source. In (34.9), on the left hand side it used $\Phi$ to denote the upper bound of integration in the integral defining $F$ but on the right hand side it was used to denote the limit of the sequence $\phi_1, \phi_2, \ldots.$ The left hand side of (34.9) should be changed to $F(k, \phi)$ to make everything correct.

As it is currently, warning bells should ring because $F(k, \pi/2)$ is finite for $|k|\leq 1$ yet taking $\Phi = \pi/2$ on the right hand side results in a $\tan(\pi/2)$ term which is not defined. Because of how floating point arithmetic approximates real numbers, your program evaluates the $\tan(\pi/2)$ term as approximately $1.633 \times \ 10^{16}$ (see here) so the log term in your "approx_result" evaluates to about $37.3,$ which is why your approximate result ended up being about $44.06.$

If you correct the error I mentioned in the first paragraph I believe your code can be made to produce correct results. If you are interested not just about correctness but in making your code more numerically efficient, let me know and I can point some things out.

EDIT:

You may want to check out [1] from page 12 onwards.

[1] : www.mygeodesy.id.au/documents/Elliptic%20Integrals%20and%20Landen's%20Transformation.pdf