Consider the 2d iteration \begin{aligned} x_{n+1} &=\frac{x_n+y_n}2,\\ y_{n+1} &=\sqrt{x_ny_n}. \end{aligned} We assume that $x_0>y_0>0$. Hence, $x_n>x_{n+1}>y_{n+1}>y_n$ due to the AM-GM inequality. Since the sequence $x_n$ is bounded from below and decreases monotonically, there is a limit $x_\infty$. Also, $y_n$ grows monotonically and is bounded from above. Hence, $\lim_{n\rightarrow\infty}y_n=y_\infty$. We also know that $$ x_{n+1}-y_{n+1}= \frac{(\sqrt{x_n}-\sqrt{y_n})^2}2=\frac12\frac{\sqrt{x_n}-\sqrt{y_n}}{\sqrt{x_n}+\sqrt{y_n}}(x_n-y_n)<\frac12(x_n-y_n)<2^{-n-1}(x_0-y_0). $$ Therefore, we find that both sequences have the same limit $L(x_0, y_0)=x_\infty=y_\infty$.
I cannot find the fixed-point $(x^*, y^*)$ of the iteration above, since both equations yield only $x^*=y^*$. I wonder if there is an analytic expression for $L(x_0, y_0)$ or at least an aproximate expression.
In the paper by Borwein and Borwein "The arithmetic-geometric mean and fast computation of elementary functions", SIAM V26, 351 (1984), they claim that elliptic integral of the second kind $$ I(a,b)= \int_0^{\pi/2}\frac{d\theta}{\sqrt{a^2\cos^2\theta+b^2\sin^2\theta}}, $$ where $I(a,b)=I(L(a,b))=\frac{\pi}{2L(a,b)}$. By numerical integration, I got something else though. Edit: I fixed this numerical problem, it was a bug in my code. I added the code below. Unfortunately, it seems the code formatting tool of math stackexchange is not working properly.
import numpy as np
import scipy.integratefrom scipy.special import ellipk
a = 2.0
b = 1.0
x = a
y = b
###########################################################
def fun(x, a, b):
f = 1.0/np.sqrt((anp.cos(x))2+(bnp.sin(x))**2)
return f
###########################################################
def fun2(x, m):
f = 1.0/np.sqrt(1-m(np.sin(x))*2)
return f
##########################################################
integral, err = scipy.integrate.quad(fun,0., np.pi/2,args=(a,b,))
print('Original Gaussian quadrature integral', integral)
for i in range(100):
plt.plot(x, y, 'ro', ms=2.0)
x_n = (x+y)/2
y_n = np.sqrt(x*y)
x = x_n
L = x
print('iteration result', np.pi/(2*L))
m = 1-(b/a)**2
Im = ellipk(m)/a
print('ellipk', Im)
integral, err = scipy.integrate.quad(fun2, 0., np.pi/2,args=(m,))
print('Gaussian quadrature integral', integral/a)
- How to show that $I(a, b)=I\left(\frac{a+b}2, \sqrt{ab}\right)$?
According to the authors, this problem has a long history, dating back to Gauss and Legendre.
