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}]

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. For this, we first need a list of ordered pairs:

pairs=Transpose[{Range[nMax],acc10}]
model=LinearModelFit[pairs,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 π.

Here is one way to do it:

computeMachin[numTerms_]:=Module[{sum=0},
Do[
sum = sum+(((-1))^(i+1))((16/((2i-1)(5^(2i-1))))-(4/((2i-1)(239^(2i-1))))),
{i,1,numTerms}
];

(* return the sum *)
sum
]

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

How to avoid numerical errors

If you get “precision limit reached” error, set $MaxExtraPrecision to a larger value, like this:

$MaxExtraPrecision=1000

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

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

The 5 above causes Mathematica to use 5 digits of precision for the accuracy calculation. You could use a different integer other than 5, but it doesn’t need to be a big integer.