Limit of a sequenced defined by arithmetic mean and geometric mean.

I don't offer any analytical insight, but I hope some numerical experiments would be useful.

One thing I have to warn about: this sequence converges very slowly, in fact, so slowly, that without unlimited precision software (like Mathematica) it's very hard to find even $3$ digits of the limit. I might try this in Mathematica next week, but so far I have used R for experiments.

(To the OP: R is free, you can download and install it, and copy the following script there if you want to do some experiments yourself).

To make the recurrence easier to program, I have reformulated it the following way:

$$\begin{cases} a_0=a \\ b_0=b \\ \displaystyle s_0=\frac{a+b}{2} \\ p_0 = (ab)^{1/3} \end{cases}$$

$$\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}$$

There might be better ways to do it, but this one came to mind. All the numbers seem to converge to the same limit (I have no formal proof for it though).

For $n=100$ and for the initial data the OP used ($a=2,b=4$) the results are:

$$a_n=2.9337 \dots \\ b_n=2.9320 \dots$$

We can assume that the limit is about $2.933 \dots$, but it's hard to say more without using unlimited precision.

The convergence is easy to see numerically, though I agree with the OP that the fact that $\lim_{n \to \infty} a_n = \lim_{n \to \infty} b_n$ remains to be proved.

Here's the plot for $a_n$ (red), $b_n$ (blue), $s_n$ (brown) and $p_n$ (green) in the above case, generated by R:

Here's the program in 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

---

Concerning the convergence proof. Taking the limit for the recurrence relations I used we have:

For $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}$$

Because:

$$\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$$

This may not be completely rigorous, but it's clear enough.