Continued fraction approximation for $\pi$

Question

How to code this formula in Mathematica to approximate $\pi$? Thanks for your help!

$$\frac4{\pi}=1+\cfrac{1^2}{2+\cfrac{3^2}{2+\cfrac{5^2}{2+\cfrac{7^2}{2+\cdots}}}}$$

score 14 · Answer 1 · answered Jan 31 '17 at 03:05

14

1 + ContinuedFractionK[(2 n - 1)^2, 2, {n, 1, Infinity}]
(*  4/π  *)

Pick a termination point less than Infinity to get an approximation.

answered Jan 31 '17 at 03:05

Michael E2

235,386
17
334
747

score 7 · Answer 2 · edited Apr 13 '17 at 12:55

Michael's method of using ContinuedFractionK[] is the canonical way. If you want to look at the forward recursion method manually, you can use repeated matrix multiplication to implement it:

al = Range[1, 45, 2]^2; (* partial numerators *)
bl = ConstantArray[2, Length[al]]; (* partial denominators *)

1 + Divide @@@ FoldList[{{0, 1}, #2}.#1 &, IdentityMatrix[2],
                        Transpose[{al, bl}]][[All, 2]]
   {1, 3/2, 15/13, 105/76, 315/263, 3465/2578, 45045/36979, 45045/33976, 765765/622637,
    14549535/11064338, 14549535/11757173, 334639305/255865444, 1673196525/1346255081,
    5019589575/3852854518, 145568097675/116752370597, 4512611027925/3473755390832,
    4512611027925/3610501179557, 4512611027925/3481569435902,
    166966608033225/133330680156299, 166966608033225/129049485078524,
    6845630929362225/5457995496252709, 294362129962575675/227848175409504262,
    294362129962575675/234389556075339277, 13835020108241056725/10721947005578370344}

The other methods I mentioned in this answer can also be adapted to this case.

Every time they add another 1000 functions, I feel it just takes some of the fun out of programming. :) — Michael E2, Jan 31 '17 at 20:01
On the other hand, one usually hopes that the new built-ins outdo the stuff we can come up with, in either performance or generality. — J. M.'s missing motivation, Feb 01 '17 at 08:08

Karsten7 · Answer 3 · 2017-01-31T04:19:16.333

6

pi[n_] := 1 + Fold[(2*#2 - 1)^2/(2 + #1) &, 2*n - 1, Reverse[Range[n]]]

pi /@ Range[10] // N

{1.33333, 1.26316, 1.2766, 1.27176, 1.27401, 1.27279, 1.27353, 1.27305, 1.27338, 1.27314}

DiscretePlot[pi[n], {n, 10}, AxesOrigin -> {0, 4/Pi}, PlotRange -> All, FillingStyle -> Red]

1 + Fold[Inactivate[(2*#2 - 1)^2, Power]/(2 + #1) &, 2*4 - 1, Reverse[Range[4]]]

edited Jan 31 '17 at 04:19

answered Jan 31 '17 at 04:09

Karsten7

27,448
5
73
134

score 5 · Answer 4 · edited Jan 31 '17 at 19:37

A quick hack :) (someone good in Mathematica can make this more functional )

rest[n_] := If[n < 1024, (n + 2)^2/(2 + rest[n + 2]), n]
1 + 1/(2 + rest[1]) // N

Mathematica graphics

4/Pi // N

Mathematica graphics

Showing speed of convergence

rest[n_, max_] := If[n < max, (n + 2)^2/(2 + rest[n + 2, max]), n];
data = Table[{i, 1 + 1/(2 + rest[1, i])}, {i, 3, 100}];
Show[ListLinePlot[data], Plot[4/Pi, {x, 0, 100}, PlotStyle -> Red], 
 PlotRange -> All]

Mathematica graphics

score 4 · Answer 5 · answered Jan 31 '17 at 03:08

4

Lord Brouncker's Formula is the same as:

(Pi^2)/8 = 1/1^2 + 1/3^2 + 1/5^2 + ...

so for something simple in Mathematica

N[Sqrt[8 Sum[1/n^2, {n, 1, 99999, 2}]], 6]

3.14159

answered Jan 31 '17 at 03:08

Young

7,495
1
20
45

Continued fraction approximation for $\pi$

5 Answers5