基于蒙特卡洛法的圆周率 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

参考

https://www.bilibili.com/video/BV1WS4y1x77p?spm_id_from=333.1007.top_right_bar_window_history.content.click

posted @ 2022-08-28 16:33  只会加减乘除  阅读(171)  评论(0编辑  收藏  举报