In the excellent $\pi$ approximation video, Katie Steckles asked for $\pi$ approximations. I teach a first year techniques module (mostly calculus and a little complex numbers and linear algebra). This year I have changed a few bits in my module; in particular I gave some of my more numerical topics to the numerical methods module and took in return some of the more analytic bits from that module. This gives both modules greater coherence, but it means I have lost one of my favourite examples, from the Taylor series topic, which uses a Maclaurin series to approximate $\pi$.

This uses a formula for $\pi$ due to John Machin (1680–1751) (for which a derivation can be found):

$$\pi =16{\tan }^{-1}\left(\frac{1}{5}\right)–4{\tan }^{-1}\left(\frac{1}{239}\right)\text{.}$$

First, we need a Maclaurin series for ${\tan }^{-1}$. That would be:

$$f(x)=f(0)+f^{\prime }(0)x+\frac{f”(0)}{2!}{x}^2+\frac{f^{(3)}(0)}{3!}{x}^3+\dots$$

To find this, we need to know the derivative of $f(x)={\tan }^{-1}(x)$, which I claim to be $\frac{1}{x^2+1}$.

(To see this, let $x=\tan (\theta )$ in $\int \frac{1}{x^2+1}\mathrm{d}x$, remembering $\frac{\mathrm{d}x}{\mathrm{d}\theta }={\sec }^2(\theta )$ and ${\tan }^2(\theta )+1={\sec }^2(\theta )$.)

So, back to our Maclaurin series, the relevant derivatives are: $f(x)={\tan }^{-1}(x)$, $f^{\prime }(x)=(x^2+1{)}^{-1}$, $f”(x)=-2x(x^2+1{)}^{-2}$, and so on (I’m waving my arms here because the quotient rule is involved at this point and it gets messy!).

Then the function values end up as: $f(0)=0$, $f^{\prime }(0)=1$, $f”(0)=0$, $f^{(3)}(0)=-2!$, $f^{(4)}(0)=0$, $f^{(5)}(0)=4!$, $f^{(6)}(0)=0$, $f^{(7)}(0)=-6!$, etc.

So

$$\begin{array}{rl}{\tan }^{-1}(x)& =0+x+\frac{0}{2!}{x}^2+\frac{-2!}{3!}{x}^3+\frac{0}{4!}{x}^4+\frac{4!}{5!}{x}^5+\frac{0}{6!}{x}^6+\frac{-6!}{7!}{x}^6+\dots \\ & =x–\frac{x^3}{3}+\frac{x^5}{5}–\frac{x^7}{7}+\dots \end{array}$$

I’m happy, for an approximation, to say ${\tan }^{-1}(x)\approx x–\frac{x^3}{3}+\frac{x^5}{5}–\frac{x^7}{7}$, so that

$${\tan }^{-1}\left(\frac{1}{5}\right)\approx \left(\frac{1}{5}\right)–\frac{{\left(\frac{1}{5}\right)}^3}{3}+\frac{{\left(\frac{1}{5}\right)}^5}{5}–\frac{{\left(\frac{1}{5}\right)}^7}{7}\approx 0.197395504761905$$

and

$${\tan }^{-1}\left(\frac{1}{239}\right)\approx \left(\frac{1}{239}\right)–\frac{{\left(\frac{1}{239}\right)}^3}{3}+\frac{{\left(\frac{1}{239}\right)}^5}{5}–\frac{{\left(\frac{1}{239}\right)}^7}{7}\approx 0.004184076002075\text{.}$$

Finally,

$$\pi \approx 16\times 0.197395504761905–4\times 0.004184076002075=3.141591772182177\text{.}$$

I think it is neat to get agreement with the first five decimal places from only four terms.

The first time I did this example in a lecture, I started by joking “this is a long and complicated example. When I get to the end, I fully expect a round of applause”. When I finished, somewhat embarrassingly, I received one — along with ironic whoops from the back row!

To take this a little further, I wrote this quick Python code.

import decimal
import math

for loop in range(1,12):
    pivalue=0
    firstterm=0
    secondterm=0
    for i in range(0, loop):
        firstterm = firstterm + decimal.Decimal((-1)**i * (1/5**(2*i+1))/(2*i+1))
        secondterm = secondterm + decimal.Decimal((-1)**i * (1/239**(2*i+1))/(2*i+1))
    pivalue = decimal.Decimal(16 * firstterm - 4 * secondterm)
    print("Using {} terms: {:.15f}".format(loop,pivalue))
    
print('math.pi:        {:.15f}'.format(math.pi))

This gives the following values, showing that this finds 15 digits of $\pi$ by the time eleven terms of the sequence are computed.

Using 1 terms: 3.183263598326360
Using 2 terms: 3.140597029326060
Using 3 terms: 3.141621029325035
Using 4 terms: 3.141591772182177
Using 5 terms: 3.141592682404400
Using 6 terms: 3.141592652615309
Using 7 terms: 3.141592653623555
Using 8 terms: 3.141592653588602
Using 9 terms: 3.141592653589836
Using 10 terms: 3.141592653589792
Using 11 terms: 3.141592653589793
math.pi:        3.141592653589793

Apparently Machin used his formula to compute 100 digits of $\pi$, but to do that I’d need to get my head around increasing Python’s decimal places. Or get a lot more free time and calculate it by hand!