Monte Carlo

import matplotlib.pyplot as plt
from numpy import random, pi
import time

Na figura, se a área do círculo é $$A_C = \pi r^2$$ a área do quadrado é $$A_Q = 4r^2$$ Então $\frac{A_C}{A_Q}=\frac{\pi}{4}$ logo $$\pi= 4\frac{A_C}{A_Q}$$

ts=time.time()
N = 1e7
circle_x = []
circle_y = []
square_x = []
square_y = []
pim =[4]
x = random.uniform(-1,1,int(N))
y = random.uniform(-1,1,int(N))
for i in range(1,int(N)):
    if (x[i]**2 + y[i]**2) <= 1:
        circle_x.append(x[i])
        circle_y.append(y[i])
    else:
        square_x.append(x[i])
        square_y.append(y[i])
    pim.append(4*len(circle_x)/i)
Pi = 4*len(circle_x)/N
print('Pi ~ ',Pi,'quando é ',pi)
tf=time.time()
print('Tempo %.2f s'%(tf-ts))
plt.plot(circle_x,circle_y,'g.')
plt.plot(square_x,square_y,'r.')
plt.grid()
plt.rcParams["figure.figsize"] = [4, 4]
plt.rcParams["figure.dpi"] = 100.

Pi ~  3.1415548 quando é  3.141592653589793
Tempo 19.72 s

plt.plot(range(int(N)),pim)
plt.semilogx([1,N],[pi,pi])
plt.ylim(3.1,3.2)
plt.xlim(100,)
plt.grid()

plt.plot(range(int(N)),pim)
plt.semilogx([1,N],[pi,pi])
plt.ylim(3.14,3.15)
plt.xlim(1e4,)
plt.grid()