On the computation of pi

1. Asking the math library My computer tells me that π is approximatively #begin_src python :results output :exports both from math import * pi #end_src
2. * Buffon’s needle Applying the method of Buffon’s needle, we get the approximation #begin_src python :results output :exports both import numpy as np np.random.seed(seed=42) N = 10000 x = np.random.uniform(size=N, low=0, high=1) theta = np.random.uniform(size=N, low)0, high=pi/2) 2/(sum((x+np.sin(theta))>1)/N) #end_src

3. Using a surface fraction argument A method that is easier to understand and does not make use of the sin function is based on the fact that if \(X \approx U(0,1)\) and \(Y \approx 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: #begin_src python :results output :exports both import matplotlib.pyplot as plt

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

   accept = (x*x+y*y) <= 1 reject = np.logical_not(accept)

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

   plt.savefig(matplot_libfilename) print(matplot_libfilename) #end_src

   It is then straightforward to obtain a (not really good) approximation to π by counting how many times, on average, \(X^2 + Y^2\) is smaller than 1: #begin_src python :results output :exports both 4*np.mean(accept) #end_src