MATLAB蒙特卡洛方法求圆周率π（代码详解）

一 计算原理

        一个正方形内部内切一个圆，设圆的半径为r，则正方体的边长为2r。在该正方体内部，随机投掷n个均匀分布的点，统计投掷点在圆内的与总点数的比例，得出圆的近似面积，就是π的值。理论上，n越大，实验次数越多，计算得出的π的值就越精准。

二 所使用的计算公式

圆形参数方程：x = a=r*cosθ;y = b+r*sinθ；

两点之间距离公式：

三 完整代码

%蒙特卡洛算法求解pi的值

%初始化数据：p为投掷点的个数，(a,b)为圆心，t为在在圆内的投掷点个数
p = 10000;
r = 1;
a = 1;
b = 1;
t = 0;

%绘制圆
theta = 0:pi/20:2*pi;
x = a+r*cos(theta);
y = b+r*sin(theta);
plot(x,y);
title('蒙特卡洛算法求pi的值');
hold on%使图像持续存在

%绘制点  t为在圆内的随机点数
for i = 1:p%循环p次，生成p个点
    p_x = rand*2;%随机生成点p在x轴的坐标
    p_y = rand*2;%随机生成点p在y轴的坐标
    %如果该点在圆内，则为绿色且t的计量加一，不在则为红色
    if (p_x-1)^2 + (p_y-1)^2 < 1
        plot(p_x,p_y,'.','Colo',[0,0,1]);
        t = t+1;
    else
        plot(p_x,p_y,'.','Colo',[1,0,0]);
    end
end

%计算圆的近似面积,s为圆的近似面积，因为r = 1，则pi = s
s = (t/p)*4;
pi = s;

四 实验结果

1.随机投掷点数：10000

运行结果：

        π = 3.141600 

        误差：0.0000074

2.随机投掷点数：200000

运行结果：

        π = 3.142804

        误差：0.001247

3.随机投掷点数：500000

运行结果：

        π = 3.142476

        误差：0.0008834