数学建模————统计问题之仿真（四）

杨利霞

数学建模————统计问题之仿真（四）
1 w4 h, J( g5 v; C# I5 }5 P  仿真，顾名思义，就是利用计算机模拟研究对象，对于那些用数学公式或者规则描述的系统，计算机可以将其通过数值模拟出来，还能实现可视化。就好比我们看的小说一样，创造一个世界，需要有初始的人或物质，再加上法则（规则），那么这个世界就会逐步成型，仿真也是如此，我们需要给这个模拟世界一个初始的状态（包含应有的数据），然后告诉他运转的规则。
5 t2 Q, J# c& t) u8 x0 r  N) e


$ F, E* g  E# v* u+ @. o: a+ \       真实的系统往往存在着很多不确定因素 ，比如：要模拟某条道路的交通，我们就得知道路上的行车的情况，除了基本的交通规则之外，我们需要的是车辆的模拟。一般来说，我们都会给车流量一个分布，这样我们就相当于有了一个车辆生成器，然后通过行车规则，就可以完整的模拟出整条道路的交通。
0 m  `( }. O9 @- z3 g$ U       不过，大多数时候我们都只是设定几个交通状况指标，然后仿真不同时间的情况，就可以实现交通状况的数值模拟。当然，有时候为了论文（观赏）效果，还可以将整条道路分成很多个小块，当车经过时就让小块发亮，这样就可以看到整个交通的运行情况，这种方法我们叫做元胞自动机。
3 b9 w; m8 B4 U, |
  {5 s" ?9 s1 b6 c* Y
        既然是模拟系统，那么就需要一个系统的推进方式，我们依此可以将仿真分为时间步长法和事件步长法。时间步长法即将每经过一定时间步长就仿真一次活动，然后推进下去，而事件步长法即每发生一件事情就推进一次，当然这个步长也可以看做是每两个事件之间的时间。
) F0 Q) v6 e+ d; ]+ ^
       上面介绍的仿真方法都讲究推进，也就是说是动态的 ，除此之外还有静态仿真。静态仿真比较有名的是蒙特卡洛模拟，下面给大家展示一道百度校招笔试题：
       在平面上有一组间距为d的平行线，将一根长度为l(l<d)的针任意掷在这个平面上，求此针与平行线中任意一根相交的概率，用高等数学(微积分、概率的方法)求解，基于布丰投针的结论，任选一种编程语言(C/C++, matlab, Python, Java)，写出模拟投针实验(程序中允许把一个理想的π作为常量使用)，求解圆周率。
注：前面的高等数学部分可以求解，已证明这个概率=2l / πd，另外针中点到相邻平行线的距离x≤l/2sinφ，l是针的长度，φ是针与平行线的夹角。
4 B5 ^8 _0 l+ I  |& w
" [' C) i" N- P       现在我们知道了规则，那就是x≤l/2sinφ，为了模拟各种情况，我们现在需要做的就是对未知量x和未知夹角φ进行随机模拟，然后计算符合规则的概率，最后依次计算圆周率。
( J3 [; W. S6 l4 b" K% Q

$ {& X8 ]" U+ n5 p
& o. [7 B' _/ A
2 V) C) a: O0 y6 {/ I- _. z  Fclc;clear;close all;
4 }) W' U! v8 Md = 2;%设定平行线之间的距离
l = 1;%设定针的长度
; j; @; b- b2 z" j1 an = 1000;%设定投针个数
! X$ c9 J5 o; F3 lbeta = 0 : 0.002 : pi;
& {4 W9 d% a, nplot(beta, l/2*sin(beta), 'k-')%绘出l/2*sin(φ)曲线
axis([0, pi, 0, d/2])%横坐标范围设在0~pi,纵坐标范围设在0~d/2
title('蒲丰投针实验')
9 K; X: Q# `# S9 t6 \/ `hold on
beta = rand(1, n) * pi;%随机生成n个角度（0~180度）
x = (d/2) * rand(1, n);%在平行线中线以下生成n个针中心
5 b1 Z* A( h; M6 _  }m = 0;
for i = 1 : n
) B3 ]- r3 l# K1 ~  ~    if x(i) <= l/2 * sin(beta(i))
7 F$ T5 U; M/ w6 v        m = m + 1;%符合条件就增计数
! j1 H. q  C$ e. z: b1 s: s        plot(beta(i), x(i), '.r')%将符合条件的针以红点形式画在图中
$ j, J3 Q' w5 J        pause(0.00001)
    else
8 L( H! ~2 @; j2 X& ]        plot(beta(i),x(i),'.b')%将不符合条件的针以蓝点形式画在图中
       pause(0.00001)
    end
' k$ K0 z$ r. s& e# X. v3 Rend
8 {) t, F" k0 np = m/n;%计算概率
* B" P( m+ w0 m1 Tpai = 2*l / (d*p);%计算圆周率
! a# ?% j' m% A. q! ^: Ydisp(['圆周率为：',num2str(pai)])

8 s4 U$ [$ D' k% ^  A
! @  w+ `2 x& b
9 f$ P9 I2 j3 d结果如下：



) M( |$ f5 p; M' B" M
; z' h9 M2 u6 r/ a3 q! f3 m; Z
1 A$ R3 Z& }& S( G
# V5 |3 C- `* D4 @4 X, W
2 t) Q5 M8 p, a: z" R