On the computation of π

1 Asking the math library

My computer tells me that π is approximatively

from math import *pi

3.141592653589793

2 Buffon's needle

By the method of Buffon's needle we get the approximation

import numpy as np

N = 10000
np.random.seed(seed=42)

x = np.random.uniform(size=N, low=0, high=1)
theta = np.random.uniform(size=N, low=0, high=pi/2)

P = sum(x + np.sin(theta) > 1) / N
2/P

3.128911138923655

3 Using Monte Carlo experiments

A method that is easier to understand and does not make use of the \(\sin\) function is based on the fact that if \(X \sim U(0,1)\) and \(Y \sim U(0,1)\), then \(P[X^2+Y^2 \leq 1] = \pi/4\) (see "Monte Carlo method" on Wikipedia). The following code uses this approach:

import matplotlib.pyplot as plt

np.random.seed(seed=42)
x = np.random.uniform(size=N, low=0, high=1)
y = np.random.uniform(size=N, low=0, high=1)

inside = (x**2 + y**2) <= 1
outside = np.logical_not(inside)

fig, ax = plt.subplots(1)
ax.scatter(x[inside], y[inside], c='b', alpha=0.2, edgecolor=None)
ax.scatter(x[outside], y[outside], c='r', alpha=0.2, edgecolor=None)
ax.set_aspect('equal')

plt.savefig(matplot_lib_filename)
print(matplot_lib_filename)

It is then straightforward to obtain an approximation to π by counting how many times, on average, \(X^2 + Y^2 \le 1\):

4 * np.mean(inside)

3.1556

4 Through Wallis' product

John Wallis' product for π states that the infinite product below converges to \(\pi / 2\). Check out the geometrical proof.

prod = 1
for i in range(2, N):
    prod *= i / (i-1) if i % 2 == 0 else (i-1) / i

prod * 2

3.141435562175509

5 As the root of a function

As we know, \(sin(\pi) = 0\), thus value for π can be found by approximating the root (one of them, that is) of \(sin(x)\). This can be done using Newton's method, for instace.

from scipy.optimize import fsolve
results = fsolve(sin, x0=3)
results[0]

3.141592653589793