Solución:
No ofrezco ninguna perspectiva analítica, pero espero que algunos experimentos numéricos sean útiles.
Una cosa sobre la que debo advertir: esta secuencia converge muy lentamente, de hecho, tan lentamente, que sin un software de precisión ilimitado (como Mathematica) es muy difícil encontrar incluso $ 3 $ dígitos del límite. Podría probar esto en Mathematica la próxima semana, pero hasta ahora he usado R para experimentos.
(Para el OP: R es gratis, puede descargarlo e instalarlo y copiar el siguiente script allí si desea hacer algunos experimentos usted mismo).
Para que la recurrencia sea más fácil de programar, la he reformulado de la siguiente manera:
$$ begin {cases} a_0 = a \ b_0 = b \ displaystyle s_0 = frac {a + b} {2} \ p_0 = (ab) ^ {1/3} end {cases} $ PS
$$ begin {cases} a_n = s_ {n-1} \ displaystyle b_n = p_ {n-1} a ^ { frac {1} {2n + 1}} \ displaystyle s_n = frac { n} {n + 1} s_ {n-1} + frac {a_n + b_n} {2n + 2} \ Displaystyle p_n = p_ {n-1} ^ { frac {2n + 1} {2n + 3}} left (a_n b_n right) ^ { frac {1} {2n + 3}} end {cases} $$
Puede que haya mejores formas de hacerlo, pero esta me vino a la mente. Todos los números parecen converger al mismo límite (aunque no tengo una prueba formal de ello).
Para $ n = 100 $ y para los datos iniciales que se utilizó el OP ($ a = 2, b = 4 $) los resultados son:
$$ a_n = 2.9337 dots \ b_n = 2.9320 dots $$
Podemos suponer que el límite es de aproximadamente $ 2.933 dots $, pero es difícil decir más sin usar una precisión ilimitada.
La convergencia es fácil de ver numéricamente, aunque estoy de acuerdo con el OP en que el hecho de que $ lim_ {n to infty} a_n = lim_ {n to infty} b_n $ queda por probar.
Aquí está el gráfico para $ a_n $ (rojo), $ b_n $ (azul), $ s_n $ (marrón) y $ p_n $ (verde) en el caso anterior, generado por R:
Aquí está el programa en R:
#Kind of AGM, but different
options(digits=22)
a <-2;
b <- 4;
a0 <- a;
b0 <- b;
s0 <- (a0+b0)/2;
p0 <- (a0*b0)^(1/3);
n <- 0;
N <- 100;
x <- 1:N;
y <- rep.int(a,N);
z <- rep.int(b,N);
u <- rep.int(s0,N);
w <- rep.int(p0,N);
while (n < N-1){
n <- n+1;
a1 <- s0;
b1 <- p0*a1^(1/(2*n+1));
s1 <- n*s0/(n+1)+(a1+b1)/(2*n+2);
p1 <- p0^((2*n+1)/(2*n+3))*(a1*b1)^(1/(2*n+3));
y[n+1] <- a1;
z[n+1] <- b1;
u[n+1] <- s1;
w[n+1] <- p1;
a0 <- a1;
b0 <- b1;
s0 <- s1;
p0 <- p1;
}
y2 <- (a+b)/2;
y1 <- (a*b*(a+b)/2)^(1/3);
plot(x,y,col="red",ylim=c(y1,y2),xlab="n",ylab="an,bn");
points(x,z,col="blue")
points(x,u,col="brown")
points(x,w,col="green")
a1
b1
Respecto a la prueba de convergencia. Tomando el límite de las relaciones de recurrencia que usé, tenemos:
Por $ n gg 1 $:
$$ begin {cases} a_n = s_ {n-1} \ b_n = p_ {n-1} \ s_n = s_ {n-1} \ p_n = p_ {n-1} end {cases} $$
Porque:
$$ lim_ {n to infty} x ^ {1 / n} = 1 $$
$$ lim_ {n to infty} frac {n} {n + 1} = 1 $$
$$ lim_ {n to infty} frac {x} {n} = 0 $$
Puede que esto no sea completamente riguroso, pero es lo suficientemente claro.