12
$\begingroup$

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?

$\endgroup$
11
$\begingroup$

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)

$\endgroup$
  • $\begingroup$ Oh I see. I think this is what I need. Thanks :) $\endgroup$ – Chen M Ling May 1 '16 at 11:49
  • 1
    $\begingroup$ 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. $\endgroup$ – RunnyKine May 1 '16 at 12:21
  • $\begingroup$ @RunnyKine Added! Thanks. My first version takes 10 seconds for n=12, yours is immediate! $\endgroup$ – mattiav27 May 1 '16 at 12:29
  • $\begingroup$ Thanks @RunnyKine for making it faster :) $\endgroup$ – Chen M Ling May 1 '16 at 13:26
  • 3
    $\begingroup$ i don't know that a numerical evaluation would be considered "proof".. $\endgroup$ – george2079 May 1 '16 at 17:17
20
$\begingroup$

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\)]\)"}}] 

enter image description here

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\)]\)"}}
]

enter image description here

$\endgroup$
  • $\begingroup$ Thanks so much for visualizing it. But I couldn't understand your code. Where did you use the recursion? :/ $\endgroup$ – Chen M Ling May 1 '16 at 11:48
  • 1
    $\begingroup$ @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. :) $\endgroup$ – ubpdqn May 1 '16 at 11:54
  • $\begingroup$ 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 :) $\endgroup$ – Chen M Ling May 1 '16 at 13:25
10
$\begingroup$

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
*)
$\endgroup$
  • $\begingroup$ Oh Great!! I will try this! Thanks @MichaelE2 $\endgroup$ – Chen M Ling May 1 '16 at 13:27
  • 1
    $\begingroup$ @MichaelE2 enjoyed both the tongue in cheek and the lesson +1 :) 🖖 $\endgroup$ – ubpdqn May 1 '16 at 13:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.