数学建模清风微信公众号的习题答案(挑战篇-蒙特卡罗思想、枚举法和网格搜索法)

以下是个人结合在微信公众号是所学到的知识取做的，如有不正确或不足之处，欢迎指正！

蒙特卡罗模拟计算$\pi$

Q11.蒙特卡罗模拟是一种以概率和统计理论为基础的计算方法，它能通过随机数来解决很多计算问题。蒙特卡罗方法将所求解的问题同一定的概率模型相联系，用计算机实现统计模拟或抽样，以获得问题的近似解。本题将利用蒙特卡罗模拟来计算圆周率。

有同学可能会好奇为什么这个关系式会成立，实际上这里用到了概率论与数理统计中几何概率和大数定律的思想，我们这里不深入探究，大家有个印象即可。

现在请大家使用MATLAB生成N个随机点（N取10000），每个随机点的横纵坐标都分别在区间[0,1]上均匀分布，接下来统计出位于1/4圆内的随机点的数量n，然后套用上面的公式即可计算出一个近似的圆周率。将你的N分别变成之前的10倍、100倍和1000倍，你计算出来的圆周率的精度有何变化？

解:随着N取得越大，则越来越逼近$\pi$。但取得太大，数量太大导致计算机运算时间更长，需要更大的内存。

clear,clc
x=rand(1,10000);
y=rand(1,10000);
N=length(x);
scatter(x,y,'.')
hold on
t=0:pi/200:pi/2;
x1=cos(t);
y1=sin(t);
plot(x1,y1,'r','Linewidth',2)      %画1/4圆
for i=1:length(x)
    if sqrt(x(i)^2+y(i)^2)<=1
        X(i)=x(i);
        Y(i)=y(i);        %在1/4圆内的点
    else
        X1(i)=x(i);
        Y1(i)=y(i);       %不在1/4圆内的点

    end
end
scatter(X,Y,'b.')         %在1/4圆内的点的散点图
scatter(X1,Y1,'k.')
ind=find(X~=0);
n=length(ind);
format long
pai=4*n/N
e=abs(pai-pi)             %这是误差

pai =

   3.143200000000000


e =

   0.001607346410207

蒙特卡罗模拟求定积分

解：学过微积分的都知道，$f(x)$在$[a,b]$上的积分就是曲线$y=f(x)$在积分区间是与$x$轴围成的曲边梯形的面积，即图中的$S1$部分。用蒙特卡洛的思想就是：曲线下方与$x$轴围成的曲边梯形的面积($S1$)与以$a，b边为矩形的面积($$S1+S2$$)之比，即$$$P=\frac{S1}{S1+S2}=\frac{落在区域S1上的点数}{落在矩形上全部点数}$$
设落在区域$S1$上的点数为$n$,落在矩形上全部点数为$N$,则
$$\frac{S1}{S1+S2}=\frac{落在区域S1上的点数}{落在矩形上全部点数}=\frac{n}{N}$$即$${\int }_a^{b}f(x)\mathrm{d}x=S1=\frac{n(S1+S2)}{N}$$

clear,clc
x=2*rand(1,10000);
y=2*rand(1,10000);
S=2*2;
N=length(x);
y1=x./(exp(x)-1);
for i=1:length(x)
    if y(i)<=y1(i)
        Y(i)=y(i);
    end
end
ind=find(Y~=0);
n=length(ind);
disp('积分的估计值为：')
jf=S*n/N
积分的估计值为：

jf =

    1.2164

蒙特卡罗模拟计算概率

Q13.清风订了一份报纸，送报人可能在早上6：30至7：30之间把报纸送到他家门口，而清风出门的时间在早上7：00到8：00之间（假设送报人到达的时间和清风出发的时间在对应的时间区间内都是均匀分布的）。请使用蒙特卡罗模拟来估计清风在离开家之前能拿到报纸的概率。注意：这是概率论中几何概型的一道经典题目，准确的概率是7/8.

解：如下图所示，在红线框内即是清风在离开家之前能拿到报纸的概率。通过分析，我们知道它是$S1$,则有$$P=\frac{S1}{S1+S2}=\frac{落在区域S1上的点数}{落在小正方形上全部点}$$这里为了方便计算，我们记落在小三角形$S1$上的$n$,全部点为$N$,则有$$P=\frac{S1}{S1+S2}=\frac{N-n}{N}$$

我们把小正方形按同一比例大小方缩到0~1之间，则有

clear,clc
x=rand(1,10000);
y=rand(1,10000);
t=0.5:0.01:1;
y1=t-0.5;
plot(t,y1,'r','Linewidth',2)
hold on
scatter(x,y,'k.')
for i=1:length(x)
    if x(i)>=0.5&y(i)<=x(i)-0.5     %找在S1和S2的点
        X(i)=x(i);Y(i)=y(i);
    else
        X1(i)=x(i);Y1(i)=y(i);
    end
end
ind=find(X~=0);
ind1=find(X1~=0);
scatter(X,Y,'b.')
scatter(X1,Y1,'k.')
n=length(ind);
N=length(x);
format long
p=(N-n)/N
e=abs(p-7/8)   %误差

p =

   0.876100000000000


e =

   0.001100000000000

枚举法与网格搜索法

Q18.枚举法与网格搜索法：

解：枚举法和网格搜索法不知道的同学可以自行去查一下。
(1)这个很简单，只要给的$x$向量步长足够的小，那么得到的结果就会更精确。但不能太小，前面说了，会耗时间，占内存，只是要理解就可以。

clear,clc
 x=0:0.0001:1;
y1=2*x;
y2=exp(x);
Y=y1-y2;
plot(x,y1,'b',x,y2,'r')
legend('\ity=2*x','\ity=e(x)')
disp('最大值为：')
M=max(Y)
[ind,location]=find(Y==max(Y))
x(location)
最小值为：

M =

  -0.613705641106080


ind =

     1


location =

        6932


ans =

   0.693100000000000

大家可以看到曲线$y=2x$始终位于曲线$y=e^x$的下方，所以最大值固然是个负数，所以通过图形来看，图像离得越近，上面他们的差值就越大。

(2)这里可以去了解一下到数值分析中的共轭梯度法里面的最速下降法，取$y$的相反数$f=-y=-(2x-e^x)$即$f=e^x-2x$。求$y$的最大值即求$f$的最小值，感兴趣的同学可以去了解一下。最速下降法之后我会单独讲。
(3)这里说一下曼哈顿距离¹。
$$d(x,y)=\sum _{i=1}^n|x_i-y_i|$$

这个图中，红、蓝、黄这些轨迹都是曼哈顿距离，那条直线是欧氏距离。举个例子，如果点$A(x1,y1),B(x2,y2)$,则它们的曼哈顿距离为$$d(A,B)=|x1-x2|+|y1-y2|$$

clear,clc
x=[2 5 9 12 15];
y=[8 -4 -9 3 -1];
x1=0:0.01:20;
for i=1:length(x1)
s=0;
for j=1:5
 d=abs(x1(i)-x(j))+abs(0-y(j));
 s=s+d;
end
 s1(i)=s;
end
s1;
[m,n]=min(s1)
x0=x1(n)

m =

    45


n =

   901


x0 =

     9

通过分析可以看到，该桥修在$x(9,0)$处，使得它到5个村庄的总的曼哈尔顿距离最短。最短为45。

---

1. 内容来自: von Neumann ↩︎