蒙特卡洛法计算圆周率.md

蒙特卡洛法计算圆周率

#counting PI
import random
import matplotlib.pyplot as plt
import numpy as np

DARTS = 1024*1024 # 总共撒点的个数
counts = 0 # 落在圆内的点数
PI = 0 # PI
xs = [0,0]
ys = [0,0]

# 开始画左边的图：撒点估计曲线下方的面积
plt.subplot(121)
x = np.arange(0,1,0.001)
plt.ylim(0,1) # y轴坐标范围
plt.xlabel('x') # x轴标签
plt.ylabel('y') # y轴标签
plt.plot(x,(1-x**2)**0.5) # 绘制1/4圆曲线
plt.legend(loc=2) # 在左上角增加图例
plt.legend(['y = sqrt(1-x**2)']) # 图例的内容
plt.plot([0,0,1,1],[0,1,1,0],'r',linewidth=0.2) # 绘制撒点范围框

for i in range(DARTS):
    x = random.uniform(0,1)
    y = random.uniform(0,1)
    if x**2 + y**2 <= 1: # 点落在1/4圆内
        counts += 1
        plt.subplot(121)
        plt.plot(x,y,'g.')
    else: # 点落在1/4圆外
        plt.subplot(121)
        plt.plot(x,y,'r.')
    if i>0 :
        if counts>0:
            PI = 4*counts/i

        # 开始画右边的图：PI的计算值随投掷次数的关系
        plt.subplot(122)
        xs[0] = xs[1] # 上一个PI值与下一个PI值，通过xs与ys列表中的两个元素进行两点连线
        xs[1] = i
        ys[0] = ys[1]
        ys[1] = PI
        plt.ylim(0,4.5) # y轴坐标范围
        plt.xlabel('Number of try') # x轴标签
        plt.ylabel('Estimation of PI') # y轴标签
        plt.yticks(np.arange(0,4.5,0.5)) # y轴刻度线
        plt.title('PI:{:.10f}\ncount:{}'.format(PI,i)) # 图的标题动态更新
        plt.axhline(3.1415926,linewidth=0.05,color='r') # 绘制2.71828参考线
        plt.plot(xs,ys,'b--',linewidth=0.3) # 绘制e的计算值随撒点次数变化的曲线
        plt.ion() # 保持图像处于交互更新状态
        plt.pause(0.001) # 控制撒点速度