Machin-Like Formulas for π

Math 242 Modern Computational Mathematics

Madhava Series for π

The following module returns the sum of the first n terms of the Madhava series:

π=4(1-(1/3)+(1/5)-(1/7)+...)

madhava[n_]:=Module[
(* initialize accumulator variable *)
{sum=0},

(* main loop *)
Do[
sum+=((((-1))^(i+1))/(2i-1)),
{i,1,n}
];

(* return four times the sum *) 
4*sum
]

Note that this module returns an exact value. Also note that this module returns a single value, not a list of values:

madhava[10]

We can use the Table function to easily compute a list of values:

nMax=40;
piValues=Table[madhava[n],{n,1,nMax}]

N[piValues]

Let' s use the (Subscript[log, 10]) method to find the number of correct digits of π:

accuracy=N[-Log[10,Abs[Pi-piValues]]]

Now plot the accuracy values:

ListPlot[accuracy,AxesLabel->{"n","correct digits"}]

What is the shape of this plot? If it is shaped like a square-root function, then squaring it will make it linear. Let’s try that.

accSquared=accuracy^2
ListPlot[accSquared,AxesLabel->{"n","(accuracy^2)"}]

If the accuracy plot is a (Subscript[log, 10]) function, then raising 10 to the accuracy values will make a linear plot. Let’s try that.

acc10=10^accuracy
ListPlot[acc10,AxesLabel->{"n","(10^accuracy)"}]

Wow, that looks like the function y=x! Let’s compare our plot to the plot of y=x.

Show[
ListPlot[acc10,AxesLabel->{"n","(10^accuracy)"}],
Plot[x,{x,0,nMax},PlotStyle->Orange]
]

We could also use a linear regression to find the equation of the line that fits our data points. A linear regression isn’t necessary right now, but it could be useful in the future.

Transpose[{Range[nMax],acc10}]

model=LinearModelFit[%,x,x]

model["BestFit"]

This linear model is nearly the same as y=x.

Since the function y=x is a good fit for our transformed accuracy data, the function y=(Subscript[log, 10])(x) should be a good fit for our original accuracy data. Let’s find out.

Show[
ListPlot[accuracy,AxesLabel->{"n","correct digits"}],
Plot[Log[10,x],{x,1,nMax},PlotStyle->Orange]
]

About how many terms of the Madhava series would be required to compute 8 correct digits of π?

Machin's Formula

Here is Machin’s formula, also on page 28 in the text:

(π/4)=4 arctan((1/5))-arctan((1/239))

Convert Machin's formula to an infinite series. Then write a module that returns the sum of the first n terms of the series as an approximation of π.

computeMachin[n_]:=Module[
(* initialize accumulator variable *)
{sum=0},

(* main loop *)
Do[
sum+=4((((-1))^(i+1))/(2i-1))((((1/5)))^(2i-1))-((((-1))^(i+1))/(2i-1))((((1/239)))^(2i-1)),
{i,1,n}
];

(* return four times the sum *) 
4*sum
]

N[computeMachin[10],30]

How many correct digits of π are produced by each term of the series for Machin’s formula?

piValuesMachin=Table[computeMachin[n],{n,1,20}];
accuracyMachin=N[-Log[10,Abs[Pi-piValuesMachin]],5]

ListPlot[accuracyMachin,AxesLabel->{"n","number of correct digits"}]

It looks like the number of correct digits increases linearly with the number of terms. On average, it appears that each additional term results in about 1.4 more correct digits.

slope=(accuracyMachin[[20]]-accuracyMachin[[2]])/(20-2)

Note: how to avoid certain numerical errors

If you get “precision limit reached” error, try:

$MaxExtraPrecision=1000

If you get “Indeterminate” accuracy values, use a fixed precision inside of the accuracy computation, like this:

accuracy=N[-Log[10,Abs[Pi-piValues]],5]