On the computation of pi
--
-
- Asking the math library
+++
Table of Contents
+ ++-1 1 Asking the math library
+++My computer tells me that π is approximatively -#beginsrc python :results output :exports both -from math import * +
+ +++from math import * pi -#endsrc +
+ - * Buffon’s needle
+++
2 2 * Buffon’s needle
+++Applying the method of Buffon’s needle, we get the approximation -#beginsrc 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) -#endsrc -
-Using a surface fraction argument - A method that is easier to understand and does not make use of the
+sinfunction 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: - #beginsrc python :results output :exports both -import matplotlib.pyplot as plt++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) +
++3 3. Using a surface fraction argument
++-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) +A method that is easier to understand and does not make use of the
+sinfunction 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:++ +import matplotlib.pyplot as plt -
+-accept = (x*x+y*y) <= 1 -reject = np.logicalnot(accept) -
+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) --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.setaspect(’equal’) -
+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_lib_filename') +print('matplot_lib_filename') ++matplot_lib_filename +
--plt.savefig(matplotlibfilename) -print(matplotlibfilename) -#endsrc -
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: -#beginsrc python :results output :exports both -4*np.mean(accept) -#endsrc -
- + +++ +4*np.mean(accept) +++3.112 +
+