基于蒙特卡洛法的圆周率 pi 的近似解
方法简析
简单来说,当我们无法求某个问题的精确解时,可以进行随机抽样,根据统计试验求该问题的近似解。
\(\pi\) 的近似求解
已知:圆的面积 \(S(r)=\pi r^2\)(别问我为什么),于是
\[ \pi = S(1)
\]
构造单位圆的外接正方形,如下图
先向正方形内随机撒点(均匀分布),当点的个数达到一定程度时,从统计分析的角度,一定有
\[ \frac{S_{圆}}{S_{正方形}} \approx \frac{圆内点的个数}{正方形内点的总数量}
\]
因此
\[ \pi = S(1) \approx 4\times\frac{圆内点的个数}{正方形内点的总数量}
\]
下面用Python实现
import random
import matplotlib.pyplot as plt
def IsInC(x,y):
'''
判断点(x,y)是否在圆内
'''
if (x-1)**2 + (y-1)**2 < 1:
return True
return False
if __name__ == "__main__":
## 参数定义
cx = 1 # 圆心横坐标
cy = 1 # 圆心纵坐标
## 随机模拟10000个点
n = 10000
## 存放圆内的点数
k = 0
## 随机模拟过程
plt.figure()
for i in range(n):
# 随机生成点坐标
x = random.uniform(0,2)
y = random.uniform(0,2)
# 若点在圆内,设置颜色为蓝色;在圆外则设为红色
if IsInC(x,y):
plt.plot(x,y,'b.')
k += 1
else:
plt.plot(x,y,'r.')
plt.gca().set_aspect(1) # 保证横纵坐标单位长度相同
plt.xlabel("x")
plt.ylabel("y")
plt.show()
Pi = (k/n)*4
## 打印pi的近似值
print('pi的近似值为:',Pi)
输出:
pi的近似值为: 3.1592