matlab 蒙特卡罗计算pi值

蒙特卡罗法计算pi值是比较基础的入门应用之一，网上流传的大部分代码包括百度百科上的代码都是使用for循环完成迭代的，运算速度非常慢，这里我们提供一个向量化运算的方式，以加快运算速度。接触编程久了后，会发现运算速度是一个非常重要的问题，尤其是对于matlab这类高层语言，想提高matlab的运算速度，可以数学建模公会之前的一篇总结——优化matlab运行速度的方案的教程

方法一：蒲丰投针

l = 0.6; %针的长度
a = 1; %平行线的宽度
n = 100000; %做n次投针试验
y = unifrnd(0, a/2, 1, n); %在[0, a/2]内服从均匀分布随机产生n个数
theta = unifrnd(0, pi, 1, n); %在[0, pi]内服从均匀分布随机产生n个数
k=sum(y<(l/2)*sin(theta));%记录针与平行线相交的次数
pival = (2*l*n)/(a*k);
disp(pival)

第5行代码中，直接生成了theta的100000组数据，采用matlab自带的sin函数可以同时计算这100000组数据的sin值。实际上，matlab大多数自带的函数都是支持向量化运算的，大多数情况下可以避免使用for循环。

方法二：通过圆的面积计算

n = 100000; %做n次投点试验
x=2*rand(1,n)-1;
y=2*rand(1,n)-1;
pival=4*mean(x.^2+y.^2<1);
disp(pival)

在方法二中，我们在第四行代码中同样使用了向量化操作，不用每生成一个点判断一次是否在圆形当中。

可以重复以上两种方法100次，观察pi值计算的平均误差和误差的标准差。

方法一的误差：

numrepeat=100;
pivalset=zeros(1,numrepeat);
for i=1:numrepeat
    l = 0.6; %针的长度
    a = 1; %平行线的宽度
    n = 100000; %做n次投针试验
    y = unifrnd(0, a/2, 1, n); %在[0, a/2]内服从均匀分布随机产生n个数
    theta = unifrnd(0, pi, 1, n); %在[0, pi]内服从均匀分布随机产生n个数
    k=sum(y<(l/2)*sin(theta));%记录针与平行线相交的次数
    pival = (2*l*n)/(a*k);
    pivalset(i)=pival;
end
err1=abs(pi-pivalset);
disp(mean(err1))
disp(std(err1))

方法二的误差：

numrepeat=100;
pivalset=zeros(1,numrepeat);
for i=1:numrepeat
    n = 100000; %做n次投点试验
    x=2*rand(1,n)-1;
    y=2*rand(1,n)-1;
    pival=4*mean(x.^2+y.^2<1);
    pivalset(i)=pival;
end
err2=abs(pi-pivalset);
disp(mean(err2))
disp(std(err2))

了解更多matlab应用及教程，欢迎关注公众号“数学建模公会”——