我们知道以r为半径的圆的面积是
,以r为边长的正方形的面积是
,两者面积之比即为
的值。
具体算法是模拟一个边长为1的正方形,随机在其中生成n个点,当n趋向于无穷大时,整个正方形就被这n个点所填满(下图中蓝色区域)。
在这n个点之中,统计出落入以1为半径的扇形区域中的点的个数(下图中橙色区域),记为c。
于是扇形面积与正方形的面积之比可以用c/n表示,可得c/n=
,推得
=4*c/n。
import numpy as np
import matplotlib.pyplot as mp
#预设总样本个数,值越大结果越精确,由于后面要绘制散点图,这里先设一个小一点的值做测试
total_dots=100000
#随机生成total_dots个点
x=np.random.random(total_dots)
y=np.random.random(total_dots)
#计算在四分之一圆内的点的个数,用掩码使满足要求的值变为1
counts_in_quarter=np.zeros_like(x)
mask=np.sqrt(x**2+y**2)<=1
counts_in_quarter[mask]=1
#扇形与正方形的比值再乘以4即为π的值
res=counts_in_quarter.sum()/total_dots
pi=4*res
print(pi)
#用matplotlib绘制出样本分布的散点图,蓝色代表正方形,橙色代表扇形
mp.scatter(x,y,label=r'$r^2$',color='dodgerblue')
mp.scatter(x[mask],y[mask],label=r'$\frac{1}{4}\pi r^2$',alpha=0.5,color='orange')
mp.xlim(0,1)
mp.ylim(0,1)
mp.xticks(np.linspace(0,1,5))
mp.yticks(np.linspace(0,1,5))
mp.legend()
mp.tight_layout()
mp.show()