蒙特卡洛模拟求圆周率

红太狼 2024-03-16 13:25 58阅读 0赞

## 蒙特卡洛模拟求圆周率 ##

### 算法思路 ###

代码的基本思想

是利用蒙特卡洛方法（Monte Carlo method）来估计圆周率 π \\pi π。蒙特卡洛方法是一种以概率统计为基础的数值计算方法，通过随机采样得到结果或近似值。在这个程序中，我们生成了一个以原点为中心、半径为 r r r的圆。然后我们模拟投掷 n n n个随机点，如果该点在圆内，则代表当前的点可以绕圆周作一次完整的旋转；而若该点在圆外，则表示该时刻没有绕圆周运动。最后，估计出的圆内点数 p p p与总点数 n n n之比再乘以 4 4 4，即可得到圆周率的一个估计值。

代码主要包括以下几个部分：

1.  monte\_carlo\_pi函数

该函数在半径为r的单位圆内随机投掷n个点，并返回投掷到圆内的点的数量。投掷的过程实际就是根据均匀分布生成二元组的点(x,y), 然后判断该点是否在圆内，即是否满足 x 2 + y 2 ⩽ r 2 x^2+y^2\\leqslant r^2 x2\+y2⩽r2.

1.  主函数main

该函数使用monte\_carlo\_pi函数来模拟采样，各采样用不同的随机种子。然后，将得到的圆内点数之和除以总采样数量，得到圆面积与正方形面积之比  p n \\frac\{p\}\{n\} np 的近似值  p ^ \\hat\{p\} p^. 因为圆面积为 π r 2 4 \\frac\{\\pi r^2\}\{4\} 4πr2，正方形面积为 r 2 r^2 r2，所以  p n ≈ π 4 \\frac\{p\}\{n\} \\approx \\frac\{\\pi\}\{4\} np≈4π.

最后再乘以4即可得到近似的 π \\pi π的值，并根据样本标准差和置信区间计算估计偏差。

这个代码的用处是用蒙特卡罗方法来估计圆周率。该方法可以在很短的时间内得到较为精确的结果，在数值计算中经常被使用。

### 完整代码 ###

#include <stdio.h>
    #include <stdlib.h>
    #include <math.h>
    #include <time.h>
    
    /* 该函数模拟了在一个半径为r的圆内投掷n个点的情况，并返回样本中落在圆内的点的数量 */
    int monte_carlo_pi(int n, float r)
    {
        
        int i;
        float x, y, distance;
        int p = 0; /* 落在圆内的点的数量 */
    
        for (i = 0; i < n; i++) {
        
            /* 对每个点进行均匀随机分布，即随机生成一个二维点 (x, y)，且每个点被生成的概率相等 */
            x = (float) rand() / RAND_MAX * 2 * r - r;  /* 在[-r, r]范围内随机生成x坐标 */
            y = (float) rand() / RAND_MAX * 2 * r - r;  /* 在[-r, r]范围内随机生成y坐标 */
    
            /* 判断该点是否在圆内 */
            distance = sqrt(x*x + y*y);  /* 计算点到原点的距离 */
            if (distance <= r) {
           /* 如果距离小于等于r，则表示该点在圆内*/
                p++;  /* 投掷的点数加1 */
            }
        }
    
        return p;
    }
    
    /* 主函数，对输入的参数进行处理，调用monte_carlo_pi函数模拟多次投掷情况，并输出结果 */
    int main(void)
    {
        
        int i;  /* 循环计数器 */
        int n = 1000000;  /* 每次模拟中投掷的点的数量 */
        float r = 1.0;  /* 圆的半径 */
        double mu = 0, s = 0, p_mean, p_stddev, err;  /* p的均值、p的方差、样本均值、样本标准差、误差程度 */
    
        srand(time(NULL));  /* 设置随机数种子 */
    
        /* 利用不同的随机种子进行10次模拟 */
        for (i = 0; i < 10; i++) {
        
            /* 进行一次投掷，并累计落在圆内的点的数量 */
            int p = monte_carlo_pi(n, r);  /* 落在圆内的点的数量 */
            mu += p;  /* 各次投掷所得到的p的和，后面再除以10即为均值 */
            s += p * p;  /* 各次投掷所得到的p的平方的和，后面再除以10即为方差 */
        }
    
        /* 计算得到的p的平均值和标准差 */
        mu /= (n * 10); /* 根据大数定理，采样数量越多，样本均值越接近总体均值*/
    
        // 计算样本标准差
        s = s / (n * 10 - 1); // 样本方差：S^2=(∑(Xi-X)²)/(N-1)
        p_stddev = sqrt(s); // 样本标准差： S = √S^2
    
        err = 1.96 * p_stddev / sqrt(n * 10); /* 根据95%置信区间计算误差程度，1.96为正态分布的值 */
    
        /* 输出结果 */
        printf("样本数量：%d\n", n * 10);  /* 总共采样了多少个样本 */
        printf("圆周率的估计值：%lf\n", 4 * mu);  /* 各次投掷所得到的p的均值的均值 */
        printf("误差程度：+-%lf\n", 4 * err); /* 置信区间为平均值加减误差程度 */
    
        return 0;  /* 程序正常结束*/
    }

运行结果：

样本数量：10000000
    圆周率的估计值：3.141101
    误差程度：+-0.045014