蒙特卡洛方法近似计算圆周率

一、原理图

正方形区域内有1/4圆区域，向区域内随机大量掷点。根据概率论知识，落在每一点的概率相等，当n无限大时，落在1/4区域内的点占总点数量的比例即1/4圆占正方形面积的比例，圆周率PI=4（N1/N）[N1:落在1/4圆中点的数量，N:总点数]。

二、c++试验

#include <iostream>
#include <random>
#include <cmath>
#include <iomanip>  // 为了设置输出精度

using namespace std;
typedef long long ll;


int main()
{
    int N = 10; // 点阵的规模N*N
    int M = 100;// 掷点数量
    for (int N = 10; N <= 10000; N *= 10)
    {
        for (int M = N*N; M <= 100000000; M *= 10)
        {
            int cnt = 3;
            while (cnt--)
            {
                std::random_device rd;  // 真随机数生成器
                std::mt19937 gen(rd()); // 使用Mersenne Twister算法生成高质量随机数
                std::uniform_real_distribution<> dis(0.0, N); // 生成0到N之间的均匀分布的浮点数

                long double PI; // 近似圆周率
                ll N1 = 0; // 落在半圆内的点数量
                for (ll i = 0; i < M; i++)
                {
                    long double x = dis(gen);
                    long double y = dis(gen);
                    long double d = sqrt(x * x + y * y);
                    if (d <= N)  N1++;
                }
                PI = (N1 / (long double)M) * 4;
                cout<<"N="<<N<<",M="<<M << setprecision(15) << "时近似计算圆周率π = " << PI << endl; // 设置输出精度
            }
        }
    }
}

改变M、N的值，观察输出圆周率的精度（M最大取了100000000，再大循环太折磨了...）

结果如下：

N=10,M=100时近似计算圆周率π = 3.2
N=10,M=100时近似计算圆周率π = 3.04
N=10,M=100时近似计算圆周率π = 3.28
N=10,M=1000时近似计算圆周率π = 3.076
N=10,M=1000时近似计算圆周率π = 3.264
N=10,M=1000时近似计算圆周率π = 3.12
N=10,M=10000时近似计算圆周率π = 3.1464
N=10,M=10000时近似计算圆周率π = 3.158
N=10,M=10000时近似计算圆周率π = 3.1344
N=10,M=100000时近似计算圆周率π = 3.13672
N=10,M=100000时近似计算圆周率π = 3.13812
N=10,M=100000时近似计算圆周率π = 3.13512
N=10,M=1000000时近似计算圆周率π = 3.14016
N=10,M=1000000时近似计算圆周率π = 3.142876
N=10,M=1000000时近似计算圆周率π = 3.143608
N=10,M=10000000时近似计算圆周率π = 3.1416004
N=10,M=10000000时近似计算圆周率π = 3.1419076
N=10,M=10000000时近似计算圆周率π = 3.1425656
N=10,M=100000000时近似计算圆周率π = 3.14145444
N=10,M=100000000时近似计算圆周率π = 3.14154308
N=10,M=100000000时近似计算圆周率π = 3.14189244
N=100,M=10000时近似计算圆周率π = 3.1204
N=100,M=10000时近似计算圆周率π = 3.1444
N=100,M=10000时近似计算圆周率π = 3.1388
N=100,M=100000时近似计算圆周率π = 3.14372
N=100,M=100000时近似计算圆周率π = 3.14164
N=100,M=100000时近似计算圆周率π = 3.14344
N=100,M=1000000时近似计算圆周率π = 3.13904
N=100,M=1000000时近似计算圆周率π = 3.141212
N=100,M=1000000时近似计算圆周率π = 3.14156
N=100,M=10000000时近似计算圆周率π = 3.1418392
N=100,M=10000000时近似计算圆周率π = 3.1423576
N=100,M=10000000时近似计算圆周率π = 3.1415224
N=100,M=100000000时近似计算圆周率π = 3.14161424
N=100,M=100000000时近似计算圆周率π = 3.14168548
N=100,M=100000000时近似计算圆周率π = 3.14168172
N=1000,M=1000000时近似计算圆周率π = 3.142296
N=1000,M=1000000时近似计算圆周率π = 3.140968
N=1000,M=1000000时近似计算圆周率π = 3.14106
N=1000,M=10000000时近似计算圆周率π = 3.142186
N=1000,M=10000000时近似计算圆周率π = 3.1408512
N=1000,M=10000000时近似计算圆周率π = 3.1420132
N=1000,M=100000000时近似计算圆周率π = 3.14153392
N=1000,M=100000000时近似计算圆周率π = 3.14163044
N=1000,M=100000000时近似计算圆周率π = 3.14182124
N=10000,M=100000000时近似计算圆周率π = 3.14151084
N=10000,M=100000000时近似计算圆周率π = 3.14142316
N=10000,M=100000000时近似计算圆周率π = 3.14153156

从你的实验结果可以看出，随着 M 的增加，计算出的圆周率 π 的精度逐渐提高，结果也更加接近实际值 (  3.14159265 )。同时，增加 N 也可以在一定程度上提高精度。