Notes on codes, projects and everything

Estimating the value of pi with a line of Python

I came across a video on Youtube on Pi day. Coincidently it was about estimating the value of Pi produced by Matt Parker aka standupmaths. While I am not quite interested in knowing the best way to estimate Pi, I am quite interested in the algorithm he showed in the video however. Specifically, I am interested to find out how easy it is to implement in Python.

So the idea presented in the video is practically a series, where each term is shown as below:

1, {-\frac{1}{3}}, \frac{1}{5}, {-\frac{1}{7}}, \frac{1}{9}, \ldots

Besides the alternating positive and negative sign, each term’s denominator is an odd number as well,

2k + 1, \text{where } k = 0, 1, 2, 3, \ldots

By summing up the series (it is a divergent series) up to n elements, and multiply it by 4, we should be able to estimate a value close to Pi. So in short, the value we are looking for is

\pi \approx 4 \times \sum_{k=0}^{n} \frac{{-1}^k}{2k + 1}, \text{where } n \in  \mathbb Z_{\ge 0}

Surprisingly, the above statement can be represented in a line of python code (assuming n = 100) without importing any package.

4 * sum((1 if k % 2 == 0 else -1) / ((2 * k) + 1) for k in range(100))

It is also noted in the video, where the estimated value of Pi is always alternating below and above the real value of Pi. In order to have an idea of how it looks in real life, I have expanded the code as follows:

import matplotlib.pyplot as plt
from itertools import count
from math import pi
import pandas
from sympy import *

def f(k):
    return Rational(1 if k % 2 == 0 else -1, (2 * k) + 1)

def item():
    for k in count():
        yield k, f(k)

total = 0
history = {'k': [], 'pi': []}

for k, item in item():
    total += item

    print('k = {}: {}'.format(k, (4 * total if k < 6 else 4 * total.evalf())))
    history['k'].append(k)
    history['pi'].append(float(4 * total.evalf()))

    if input('Continue? [Y|n]: ').strip().lower() == 'n':
        break

pandas.Series(history['pi'], index=history['k']).plot()

plt.show()

So running it (only tested with Python3, which is my default code target from now on) by allowing it to loop until k = 20 yields the result below:

k = 0: 4
Continue? [Y|n]:
k = 1: 8/3
Continue? [Y|n]:
k = 2: 52/15
Continue? [Y|n]:
k = 3: 304/105
Continue? [Y|n]:
k = 4: 1052/315
Continue? [Y|n]:
k = 5: 10312/3465
Continue? [Y|n]:
k = 6: 3.28373848373848
Continue? [Y|n]:
k = 7: 3.01707181707182
Continue? [Y|n]:
k = 8: 3.25236593471888
Continue? [Y|n]:
k = 9: 3.04183961892940
Continue? [Y|n]:
k = 10: 3.23231580940559
Continue? [Y|n]:
k = 11: 3.05840276592733
Continue? [Y|n]:
k = 12: 3.21840276592733
Continue? [Y|n]:
k = 13: 3.07025461777918
Continue? [Y|n]:
k = 14: 3.20818565226194
Continue? [Y|n]:
k = 15: 3.07915339419743
Continue? [Y|n]:
k = 16: 3.20036551540955
Continue? [Y|n]:
k = 17: 3.08607980112383
Continue? [Y|n]:
k = 18: 3.19418790923194
Continue? [Y|n]:
k = 19: 3.09162380666784
Continue? [Y|n]:
k = 20: 3.18918478227759
Continue? [Y|n]: n

pi
Matplotlib output

The estimation does seem to be approaching the real value of Pi as n inceases. The reason I imported SymPy was to try seeing how the fraction grows, however it becomes increasingly difficult to make sense of after n > 5 hence I dropped the idea and show the result in decimal numbers instead.

So there it is, a quick and fun and short exercise in Python.

leave your comment

name is required

email is required

have a blog?

This blog uses scripts to assist and automate comment moderation, and the author of this blog post does not hold responsibility in the content of posted comments. Please note that activities such as flaming, ungrounded accusations as well as spamming will not be entertained.

Click to change color scheme