Archimedes' Scheme to find $\pi$

Question

Given an Archimedes' scheme

$$a_{n+1}=\frac{2a_nb_n}{a_n+b_n}$$ $$b_{n+1}=\sqrt{a_{n+1}b_n}$$

with initial values $a_0=2\sqrt{3},b_0=3.$

How do I prove using Mathematica that this converges to $\pi$?

My attemp is to use RSolve or RSolveValue. But my when I input

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

It gives an output the same as input

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

What is wrong here?

Answer 1

You should use RecurrenceTable:

N@RecurrenceTable[{a[n + 1] == 2 a[n] b[n]/(a[n] + b[n]), 
   b[n + 1] == Sqrt[a[n + 1] b[n]], a[0] == 2 Sqrt[3], 
   b[0] == 3}, {a[n], b[n]}, {n, 1, 10}]
N@Pi

with result:

{{3.21539, 3.10583}, {3.15966, 3.13263}, {3.14609, 3.13935}, {3.14271, 3.14103}, {3.14187, 3.14145}, {3.14166, 3.14156}, {3.14161, 3.14158}, {3.1416, 3.14159}, {3.14159, 3.14159}, {3.14159, 3.14159}}

3.14159

You can also change the number of digit, with N[expr,n] (n is the desired number of digits). See documentation.

EDIT Following @RunnyKine suggestion:

RecurrenceTable[{a[n + 1] == 2. a[n] b[n]/(a[n] + b[n]), 
  b[n + 1] == Sqrt[a[n + 1] b[n]], a[0] == 2. Sqrt[3.], 
  b[0] == 3.}, {a[n], b[n]}, {n, 1, 20}]

is much faster (on my laptop, AbsoluteTiming gives 0.5 sec for the first version and 0.01 for the second)

Answer 2

This is a "visual proof" of the Archimedean limiting regular polygons. You could implement the recursion and it would progressively approach $\pi$. The proof lies in the "squeezing" argument. $\pi$ is transcendental (no solution to this recurrence in a closed algebraic expression). Whatever implementation of recursion with approach $\pi$ that is reassuring but not proof.

This is a slight modification. The $\pi$ approximation from half the circumference of the inscribed and circumscribed polygons. f is one way to implement recursion. The "geometric" method just calculates length of polygon side (euclidean distance of points):

f[n_, {a_, b_}] := 
 Nest[
   {2 #1 #2/(#1 + #2), #2 Sqrt[2 #1/(#1 + #2)]} & @@ # &, 
   {a, b}, n]
g[n_] := Module[{pi = Table[{Cos[j], Sin[j]}, {j, 0, 2 Pi, 2 Pi/n}],
   pe = Table[Sec[Pi/n] {Cos[j], Sin[j]}, {j, 0, 2 Pi, 2 Pi/n}]},
  {Graphics[{Circle[{0, 0}, 1], FaceForm[White], EdgeForm[Red], 
     Polygon[pi], FaceForm[None], EdgeForm[Blue], Polygon[pe],
     Line[{{0, 0}, #}] & /@ pe[[{1, 2}]],
     Text["\!\(\*SubscriptBox[\(β\), \(n\)]\)", 
      0.9 Mean[pi[[{1, 2}]]]],
     Text["\!\(\*SubscriptBox[\(α\), \(n\)]\)", 
      1.1 Mean[pe[[{1, 2}]]]]
     }], (Norm[#1 - #2]/2) & @@@ {pe[[{1, 2}]], pi[[{1, 2}]]}}]
dem[n_, r_] := Module[{s = g /@ PowerRange[n, n 2^(r - 1), 2], i},
  i = s[[1, 2]];
  MapIndexed[
   Flatten[{#2[[1]] - 1, #1[[1]], N[#1[[2]]], 
      N[n 2^(#2[[1]] - 1) #1[[2]]], N[f[#2[[1]] - 1, 2 n i]/2]}] &, 
   s]]

Visualizing:

TableForm[dem[6, 4], 
 TableHeadings -> {None, {"n", "Polygons", 
    "\!\(\*SubscriptBox[\(α\), \(n\)]\)", 
    "\!\(\*SubscriptBox[\(β\), \(n\)]\)", 
    "\!\(\*SubscriptBox[\(Pi\), \(geometric\)]\)\n \
\!\(\*SubscriptBox[\(a\), \(n\)]\)", 
    "\n \!\(\*SubscriptBox[\(b\), \(n\)]\)", 
    "\!\(\*SubscriptBox[\(Pi\), \(recursion\)]\)\n \
\!\(\*SubscriptBox[\(a\), \(n\)]\)", 
    "\n \!\(\*SubscriptBox[\(b\), \(n\)]\)"}}] 

Just for illustrative purposes (would have to very careful for arbitrary precision):

st[n_] := 
 Module[{a = RealDigits[n, 10, 20][[1]], 
   pi = RealDigits[N[Pi, 20]][[1]], lg, pos, a1},
  pos = FirstPosition[a - pi, _?(# != 0 &)];
  a1 = StringJoin@(ToString /@ a);
  Row[{Style[StringTake[a1, 1], Red, Bold], ".", 
    Style[StringTake[a1, 2 ;; pos[[1]]], Red, Bold], 
    StringTake[a1, pos[[1]] + 1 ;;]}]
  ]
fn[n_, {a_, b_}] := 
  Nest[
    N[{2 #1 #2/(#1 + #2), #2 Sqrt[2 #1/(#1 + #2)]}, 40] & @@ # &, 
    {a, b}, n]
res = fn[#, {4 Sqrt[3], 6}]/2 & /@ Range[0, 20];

So,

TableForm[
  Map[st, res, {2}], 
  TableHeadings -> 
   {Range[0, 20], 
    {"\!\(\*SubscriptBox[\(a\), \(n\)]\)", 
     "\!\(\*SubscriptBox[\(b\), \(n\)]\)"}}
]

Answer 3

You can use Mathematica to prove the induction step in a proof that $a_n, b_n \rightarrow L$ for some limit $L$ by showing that $a_n$ is decreasing and $b_n$ increasing toward each other and that the difference goes to zero. To show $L = \pi$ you would need a definition of $\pi$ you can relate to this sequence. For instance, if you define $\pi$ to be the limit $L$, then you're done. ;-) This is essentially what Archimedes did in deriving the sequence, so perhaps it is not such a flippant suggestion after all.

next = RecurrenceTable[
   {a[n + 1] == 2 a[n] b[n]/(a[n] + b[n]), b[n + 1] == Sqrt[a[n + 1] b[n]],
    a[0] == α, b[0] == β},       (* arbitrary starting points *)
   {a, b}, {n, 1, 1}] // First
diff = First@next - Last@next;
ratio = diff/(α - β);
"a[n] decreasing " ->  Simplify[First@next < α, α > β >= 3]
"b[n] increasing " ->  Simplify[Last@next > β, α > β >= 3]
"a[n]-b[n] --> 0 " ->  Simplify[0 < ratio < 1/4, α > β >= 3]
(*
  {(2 α β)/(α + β),  Sqrt[2] Sqrt[(α β^2)/(α + β)]}
  "a[n] decreasing " -> True
  "b[n] increasing " -> True
  "a[n]-b[n] --> 0 " -> True
*)