Solution to Simultaneous Arithmetic/Geometric Mean Recursion Relations

Question

I'm trying to solve the simultaneous convergent sequences of geometric/arithmetic means where $a_{n+1}=\frac{1}{2}(a_n+b_n)$ and $b_{n+1}=\sqrt{a_nb_n}$ and initial values are $a_0=1+x$ and $b_0=1-x$. I've tried an implementation using RSolve[] as follows but Mathematica just spits my input back at me without raising an error. Can someone point me in the right direction?

RSolve[{a[n + 1] == 1/2 (a[n] + b[n]), b[n + 1] == Sqrt[a[n] b[n]],
           a[0] == 1 + x, b[0] == 1 - x}, {a, b}, n]

EDIT/UPDATE: I have taken Sasha's advice and implemented ArithmeticGeometricMean[]. The solution $M(x)$ matches the iterative solution obtained from the last element of the following recurrence table.

RecurrenceTable[{a[n + 1] == 1/2 (a[n] + b[n]),b[n + 1] == Sqrt[a[n] b[n]], a[0] == 1 + x, b[0] == 1 - x}, {a,b}, {n, 1, 10}]

That is, the coefficients of the series expansion for the approximate iterative solution match those for the series expansion of the ArithmeticGeometricMean[] function up to order 10 at least. Further, and for completeness, the series expansion of $\frac{1}{M(x)}$ also matches that of the following integral, which is actually what I was ultimately after. $$\frac{1}{\pi}\int_0^\pi\frac{d\phi}{\sqrt{1-x^2\cos^2(\phi)}}$$ Thanks everyone for the help! I'm closing the question now.

Answer 1

If you are looking for the $a(x) = \lim_{n \to \infty} a_n(x) = \lim_{n \to \infty} b_n\left(x\right)$, then Mathematica provides this as ArithmeticGeometricMean function.

Moreover, it can be expressed in terms of the complete elliptic integral of the first kind:

In[12]:= FunctionExpand[ArithmeticGeometricMean[1 + x, 1 - x], 
 0 < x < 1]

Out[12]= \[Pi]/(2 EllipticK[x^2])

Answer 2

Clear[f1, f2, f3]

The recurrence relationship can be implemented using FixedPoint

f1[x_?NumericQ] :=
  FixedPoint[
     {(#[[1]] + #[[2]])/2, Sqrt[#[[1]]*#[[2]]]} &,
     {1 + x, 1 - x}][[1]] /; -1 < x < 1;

Or with NestWhile

f2[x_?NumericQ] :=
  NestWhile[
     {(#[[1]] + #[[2]])/2, Sqrt[#[[1]]*#[[2]]]} &,
     {1 + x, 1 - x}, Unequal @@ # &][[1]] /; -1 < x < 1;

Or, as suggested by Sasha, with ArithmeticGeometricMean

f3[x_] = FunctionExpand[
  ArithmeticGeometricMean[1 + x, 1 - x],
  -1 < x < 1]

(*  Pi/(2 EllipticK[x^2])  *)

Plot[{f1[x], f2[x], f3[x]}, {x, -1, 1},
 PlotStyle -> {
   Directive[Lighter[Red], AbsoluteThickness[4]],
   Directive[Blue, AbsoluteDashing[{10, 10}]],
   Directive[Green, AbsoluteDashing[{20, 15}]]},
 PlotPoints -> 25,
 PlotRange -> {0, 1.02},
 PlotLegends -> {FixedPoint, NestWhile, AGM}]