matlab模拟投针实验以及rand和unifrnd函数的效率

目录：

1.投针实验简述
2.Matlab模拟投针实验原理
3.四种方式代码实现
4.rand函数和unifrnd函数执行效率的比较总结即函数使用方法
不想看我啰嗦可以直接看总结奥嘿嘿。

简单说明一下投针实验：

设一组间距为a的等距平行线，向其投下长度为 l 的质地均匀针若干。当针的程度l小于平行线间距a的时候，针接触平行线的概率满足：
$$p=\frac{2l}{a\pi }$$

由上述公式，可以得到π的近似值：
$$\pi =\frac{2l}{ap}$$
这是18世纪法国数学家布封提出的圆周率π的计算方法。

matlab模拟投针实验原理

如上图所示，一根针的位置只需要中点与平行线间距 x 和针与平行线夹角 φ 就可以确定。
那么x和φ的取值范围分别是：
$$0\le x\le \frac{a}{2}$$
$$0\le \phi \le \pi$$
要满足针与平行线接触，只需满足：
$$x\le \frac{l}{2}\sin \phi$$
那么使用matlab模拟的基本思想就是：
随机生成针的位置x和夹角φ模拟针的落点情况，若针与平行线接触，就记一次数。循环n次模拟n次投针，计算这n次针接触平行线的概率。然后根据公式计算π

代码实现

这里使用四个产生随机数的方法，对比四者的效率。

unifrnd函数一次性生成:

%% unifrnd函数一次性生成。
tic
clc,clear;
n = 100000;
a = 1;         %平行线间距离a
l = 0.8;      %针的长度 小于平行线间距a即可
fi = unifrnd(0, pi, 1, n); %夹角φ
x = unifrnd(0, a/2, 1, n);  %中心到平行线间距
sum = 0;
for j = 1:100
    tol = 0;%计数君
    for i = 1:n
        if x(i) <= l / 2 * sin(fi(i))
            tol = tol + 1;
        end
    end
    p = tol / n;% 针和平行线相交出现的频率
    mypi = (2 * l) / (a * p);  %根据公式计算得到pi
    sum = sum + mypi;
end
avg = sum / 100;
disp(['100实验计算得到的pi=',num2str(avg)]);
toc

循环100次试验结果：
100实验计算得到的pi=3.1319
时间已过 0.509256 秒。

rand函数一次性生成

%%  rand函数一次性生成。
tic
clc,clear;
n = 100000;
a = 1;         %平行线间距离a
l = 0.8;      %针的长度 小于平行线间距a即可
x = rand(1, n) * a / 2; %生成n个x
fi = rand(1, n)*pi; %n个fi
sum = 0;
for j = 1:100
    tol = 0;%计数君
    for i = 1:n
        if x(i) <= l / 2 * sin(fi(i))
            tol = tol + 1;
        end
    end
    p = tol / n;% 针和平行线相交出现的频率
    mypi = (2 * l) / (a * p);  %根据公式计算得到pi
    sum = sum + mypi;
end
avg = sum / 100;
disp(['100实验计算得到的pi=',num2str(avg)]);
toc

循环100次试验的结果：
100实验计算得到的pi=3.1414
时间已过 0.322584 秒。

unifrnd函数循环生成

%投针实验的rand、unifrnd耗时对比：
% unifrnd函数循环生成
tic              %计时
clc,clear
n = 100000;      %重复次数
a = 1;         %平行线间距离a
l = 0.8;      %针的长度 小于平行线间距a即可
count = 0;       %计数君
for i = 1:n
    fi = unifrnd(0, pi);  %夹角φ
    x = unifrnd(0, a/2);  %中心到平行线间距
    if x <= l / 2 * sin(fi)
        count = count + 1;
    end
end
p = count / n;  %计算概率
pai = 2 * l / (a * p);      %公式求π
disp(['计算得到的π为：',num2str(pai)]);
toc

结果：（由于循环100次实在太慢，这里只做了单次的）
计算得到的π为：3.1599
时间已过 1.011857 秒。

rand函数循环生成

%投针实验的rand、unifrnd耗时对比：
% rand函数循环生成
tic              %计时
clc,clear
n = 100000;      %重复次数
a = 1;         %平行线间距离a
l = 0.8;      %针的长度 小于平行线间距a即可
count = 0;       %计数君
for i = 1:n
    fi = rand(1) * pi;  %夹角φ
    x = rand(1) * (a / 2);  %中心到平行线间距
    if x <= l / 2 * sin(fi)
        count = count + 1;
    end
end
p = count / n;  %计算概率
pai = 2 * l / (a * p);      %公式求π
disp(['计算得到的π为：',num2str(pai)]);
toc

单次结果：
计算得到的π为：3.1305
时间已过 0.141556 秒。
循环100次结果：
100次试验得到的π均值为：3.1415
时间已过 14.956320 秒。

rand函数和unifrnd函数执行效率总结：

1. rand函数快于unifrnd函数（因为unifrnd不是matlab内置函数）
2. 随机数一次性生成矩阵再循环调用矩阵快于每次循环的时候生成
   功能上rand完全可以替代unifrnd。所以尽可能使用rand吧。
   附上两种函数的使用方法：
   unifrnd(m,n)产生[m,n]范围内的一个随机浮点数
   rand(m,n)产生[0,1]之间的均匀分布m×n的矩阵
   rand(n)产生[0,1]之间的均匀分布n行方阵
   a + rand(m,n)(b-a)产生[a,b]之间的均匀分布m×n的矩阵
   % a + rand(m,n)(b-a)等价于unifrnd(a,b,m,n)

模型学习自清风老师的数模教程