# Archimedes' Scheme to find $\pi$

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?

• That means RSolve can't solve it. – happy fish May 1 '16 at 6:50
• – user9660 May 1 '16 at 7:53

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)

• Oh I see. I think this is what I need. Thanks :) – Chen M Ling May 1 '16 at 11:49
• To make this much faster you can use a numeric value for any of the coefficients. So e.g. 2.0 instead of 2 and you won't need the N. – RunnyKine May 1 '16 at 12:21
• @RunnyKine Added! Thanks. My first version takes 10 seconds for n=12, yours is immediate! – mattiav27 May 1 '16 at 12:29
• Thanks @RunnyKine for making it faster :) – Chen M Ling May 1 '16 at 13:26
• i don't know that a numerical evaluation would be considered "proof".. – george2079 May 1 '16 at 17:17

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],
"\!$$\*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}],
{Range[0, 20],
{"\!$$\*SubscriptBox[\(a$$, $$n$$]\)",
"\!$$\*SubscriptBox[\(b$$, $$n$$]\)"}}
]


• Thanks so much for visualizing it. But I couldn't understand your code. Where did you use the recursion? :/ – Chen M Ling May 1 '16 at 11:48
• @ChenMLing the recursion is coded in f using Nest. I merely substituted $a_{n+1}$ on the RHS of $b_{n+1}$ equation. I suggest you try out Nest after you have looked at documentation. There are many ways to implement recursion. I just liked Nest. fn just "numericizes" (if you put exact input the expressions rapidly become more complex. :) – ubpdqn May 1 '16 at 11:54
• Oh I see. So we can use Nest as well. OK. I will recognize it as a new knowledge, as I am new in Mathematica. Thanks :) – Chen M Ling May 1 '16 at 13:25

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
*)

• Oh Great!! I will try this! Thanks @MichaelE2 – Chen M Ling May 1 '16 at 13:27
• @MichaelE2 enjoyed both the tongue in cheek and the lesson +1 :) 🖖 – ubpdqn May 1 '16 at 13:29