数学建模之计算机仿真.pptx

《数学建模之计算机仿真.pptx》由会员分享，可在线阅读，更多相关《数学建模之计算机仿真.pptx（87页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、一个问题我们做一个实验:把一个硬币掷一万次,统计两个面出现的次数。这样做很简单但却需要大师时间，有没有一咱较快的办法把这个实验完成呢？第1页/共87页第2页/共87页掷硬币仿真流程图掷硬币仿真流程图第3页/共87页概述计算机科学技术的迅猛发展，给许多学科带来了巨大的影响计算计算机科学技术的迅猛发展，给许多学科带来了巨大的影响计算机不但使问题的求解变得更加方便、快捷和精确，而且使得解决实机不但使问题的求解变得更加方便、快捷和精确，而且使得解决实际问题的领域更加广泛计算机适合于解决那些规模大、难以解析际问题的领域更加广泛计算机适合于解决那些规模大、难以解析化以及不确定的数学模型例如对于一些带随机因

2、素的复杂系统，化以及不确定的数学模型例如对于一些带随机因素的复杂系统，用分析方法建模常常需要作许多简化假设，与面临的实际问题可能用分析方法建模常常需要作许多简化假设，与面临的实际问题可能相差甚远，以致解答根本无法应用，这时仿真几乎成为人们的唯一相差甚远，以致解答根本无法应用，这时仿真几乎成为人们的唯一选择在历届的美国和中国大学生的数学建模竞赛（选择在历届的美国和中国大学生的数学建模竞赛（MCM）中，学）中，学生们经常用到计算机仿真方法去求解、检验等计算机仿真生们经常用到计算机仿真方法去求解、检验等计算机仿真(computer simulation)是建模过程中较为重要的一类方法是建模过程中较为

3、重要的一类方法 第4页/共87页第5页/共87页第6页/共87页第7页/共87页第8页/共87页第9页/共87页计算机仿真可以解决以下5类问题:(1)(1)难以用数学公式表示的系统难以用数学公式表示的系统,或者没有建立和求解的有效方法或者没有建立和求解的有效方法.(2)(2)虽然可以用解析的方法解决问题虽然可以用解析的方法解决问题,但数学的分析与计算过于复杂但数学的分析与计算过于复杂,这时这时计算机仿真可能提供简单可行的求解方法计算机仿真可能提供简单可行的求解方法.(3)(3)希望能在较短的时间内观察到系统发展的全过程希望能在较短的时间内观察到系统发展的全过程,以估计某些参数对系以估计某些参数

4、对系统行为的影响统行为的影响.(4)(4)难以在实际环境中进行试验和观察时难以在实际环境中进行试验和观察时,计算机仿真是唯一可行的方法计算机仿真是唯一可行的方法,如太空飞行的研究如太空飞行的研究.(5)(5)需要对系统或过程进行长期运行比较需要对系统或过程进行长期运行比较,从大量方案中寻找最优方案从大量方案中寻找最优方案.第10页/共87页计算机仿真案例1模型建立：模型建立：由于本题要求使从搅拌中心到各个工地运输混凝土由于本题要求使从搅拌中心到各个工地运输混凝土的总的吨公里数最少，所以，该问题的目标函数是的总的吨公里数最少，所以，该问题的目标函数是 第11页/共87页求解方法：求解方法：1、高

5、数中的方法、高数中的方法2、数值计算方法、数值计算方法3、计算机仿真：、计算机仿真：离散化，遍历！离散化，遍历！第12页/共87页计算机仿真案例2例例2（赶火车过程仿真）（赶火车过程仿真）一列火车从一列火车从A站经过站经过B站站开往开往C站，某人每天赶往站，某人每天赶往B站乘这趟火车。已知站乘这趟火车。已知火车从火车从A站到站到B站的运行时间是均值为站的运行时间是均值为30min、标、标准差为准差为2min的正态随机变量。火车大约在下午的正态随机变量。火车大约在下午1点离开站。火车离开时刻的频率分布和这个人点离开站。火车离开时刻的频率分布和这个人到达站时刻的频率分布如下表所示。问他能赶到达站时

6、刻的频率分布如下表所示。问他能赶上火车的概率有多大？上火车的概率有多大？出发时刻出发时刻1:001:051:10到达时刻到达时刻1:281:301:321:34频率频率0.70.20.1频率频率0.30.40.20.1第13页/共87页u仿真过程：仿真过程：u1、生成火车的发车时间、运行时间，从而达得到、生成火车的发车时间、运行时间，从而达得到其到达其到达B站的时间。站的时间。u2、生成此人达到、生成此人达到B站的时间。站的时间。u3、如果此人到达站的时间早于火车到达时间，、如果此人到达站的时间早于火车到达时间，则算赶上火车一次。则算赶上火车一次。u4、将上述过程重复一万次，统计赶上火车的频率

7、、将上述过程重复一万次，统计赶上火车的频率作为所求概率。作为所求概率。分析：这个问题用概率论的方法求解十分困难，分析：这个问题用概率论的方法求解十分困难，它涉及此人到达时刻、火车离开它涉及此人到达时刻、火车离开A站的时刻、火站的时刻、火车运行时间几个随机变量。车运行时间几个随机变量。我们可以用计算机仿真的方法来解决。我们可以用计算机仿真的方法来解决。第14页/共87页第15页/共87页概述计算机仿真在计算机中运行实现计算机仿真在计算机中运行实现,不怕破坏不怕破坏,易修改易修改,可重用可重用,安全经安全经济济,不受外界条件和场地空间的限制不受外界条件和场地空间的限制.仿真仿真分为分为静态静态仿真

8、仿真(static simulation)和和动态动态仿真仿真(dynamic simulation).数值积分中的蒙特卡洛方法是典型的静态数值积分中的蒙特卡洛方法是典型的静态仿真仿真动态动态仿仿真真又分为又分为连续系统连续系统仿真仿真和和离散系统离散系统仿真仿真连续系统是指状态变量随连续系统是指状态变量随着时间连续变化的系统，例如传染病的检测与预报系统着时间连续变化的系统，例如传染病的检测与预报系统.离散系统是离散系统是指系统状态变量只在有限的时间点或可数的时间点上发生变化的系指系统状态变量只在有限的时间点或可数的时间点上发生变化的系统统,例如排队系统例如排队系统.第16页/共87页概述仿真

9、系统，必须设置一个仿真时钟仿真系统，必须设置一个仿真时钟(simulate clock)，它能将时间从一个时刻向另一个时刻进行推进，并且它能将时间从一个时刻向另一个时刻进行推进，并且能随时反映系统时间的当前值其中，模拟时间推进能随时反映系统时间的当前值其中，模拟时间推进方式有两种方式有两种:时间步长法时间步长法(均匀间隔时间推进法均匀间隔时间推进法,连续系连续系统常用统常用)和事件步长法和事件步长法(下次事件推进法下次事件推进法,离散系统常用离散系统常用)第17页/共87页主要内容一一:准备知识准备知识:随机数的产生随机数的产生二二:随机变量的模拟随机变量的模拟三三:连续系统的模拟连续系统的模

10、拟-时间步长法时间步长法四四:离散系统的模拟离散系统的模拟-事件步长法事件步长法五：蒙特卡洛方法五：蒙特卡洛方法第18页/共87页一:准备知识:随机数的产生由于仿真研究的实际系统要受到多种随机因素的作用和影响由于仿真研究的实际系统要受到多种随机因素的作用和影响,在在仿真过程中必须处理大量的随机因素仿真过程中必须处理大量的随机因素.要解决此问题的前提是要解决此问题的前提是确确定随机变量的类型定随机变量的类型和选择和选择合适的随机数产生的方法合适的随机数产生的方法.对随机现象进行模拟对随机现象进行模拟,实质是要给出随机变量的模拟实质是要给出随机变量的模拟,也就是说也就是说要利用计算机随机产生一系列

11、数值要利用计算机随机产生一系列数值,使它们使它们服从一定的概率分布服从一定的概率分布服从一定的概率分布服从一定的概率分布,称这些数值为称这些数值为随机数随机数.最基本最基本,最常用的是最常用的是(0,1)(0,1)区间内均匀分布的随机数区间内均匀分布的随机数.其他分布的其他分布的随机数均可利用它来产生随机数均可利用它来产生.第19页/共87页1:1:产生模拟随机数的计算机命令产生模拟随机数的计算机命令在在MATLAB中中,可以直接产生满足各种分布的随机数可以直接产生满足各种分布的随机数,命令如下命令如下:常见的分布函数常见的分布函数 MATLAB语句语句 均匀分布均匀分布U0,1 R=rand

12、(m,n)均匀分布均匀分布Ua,b R=unifrnd(a,b,m,n)指数分布指数分布 E()R=exprnd(,m,n)正态分布正态分布N(mu,sigma)R=normrnd(mu,sigma,m,n)标准正态分布标准正态分布N(0,1)R=randn(m,n)二项分布二项分布B(n,p)R=binornd(n,p,m,n1)泊松分布泊松分布 P()R=poissrnd(,m,n)以上语句均产生以上语句均产生m n 的矩阵的矩阵第20页/共87页2:案例分析例1:unifrnd(2,3)unifrnd(1,32,1,4)normrnd(1,2)normrnd(1,2,2,3)rand(2

13、,3)randn(2,3)lanslanslanslans=lans=lans=l 第21页/共87页2:案例分析例例2:2:敌空战部队对我方港口进行空袭敌空战部队对我方港口进行空袭,其到达规律服从其到达规律服从泊松分布泊松分布,平均每分钟到达平均每分钟到达4 4架飞机架飞机.(1)(1)模拟敌机在模拟敌机在3 3分钟内到达目标区域的数量分钟内到达目标区域的数量,以及在第以及在第1,2,31,2,3分钟内各到达几架飞机分钟内各到达几架飞机;(2)(2)模拟在模拟在3 3分钟内每架飞机的到达时刻分钟内每架飞机的到达时刻.分析分析:(1)n1=poissrnd(4),n2=poissrnd(4),

14、:(1)n1=poissrnd(4),n2=poissrnd(4),n3=poissrnd(4),n=n1+n2+n3n3=poissrnd(4),n=n1+n2+n3 (2)(2)由排队论知识由排队论知识,敌机到达规律服从泊松分布等敌机到达规律服从泊松分布等价于敌机到达港口的间隔时间服从参数为价于敌机到达港口的间隔时间服从参数为1/41/4的指数分的指数分布布,故可由指数分布模拟每架飞机的到达时刻故可由指数分布模拟每架飞机的到达时刻.注：如果单位时间发生的次数（如到达的人数）服从参数为r的泊松分布,则任连续发生的两次时间的间隔时间序列服从参数为r的指数分布!泊松分布的期望是，根据到泊松分布和

15、指数分布的关系，可以推出指数分布的期望是1/。第22页/共87页2:案例分析cleart=0;j=0;%到达的飞机数 while t3 j=j+1 t=t+exprnd(1/4)end第23页/共87页二:随机变量的模拟在很多实际问题中，我们需要模拟服从一定分布的随机变量在很多实际问题中，我们需要模拟服从一定分布的随机变量,来进行计算来进行计算和预测。和预测。利用均匀分布的随机数可以产生具有任意分布的随机变量的样本利用均匀分布的随机数可以产生具有任意分布的随机变量的样本,从而可从而可以对随机变量的取值情况进行模拟以对随机变量的取值情况进行模拟.1 1 连续型随机变量的模拟连续型随机变量的模拟具

16、有给定分布的连续型随机变量可以利用在区间具有给定分布的连续型随机变量可以利用在区间(0,1)(0,1)上均匀分布的随机上均匀分布的随机数来模拟数来模拟,最常用的方法是逆变换法最常用的方法是逆变换法.结论结论:若随机变量若随机变量Y Y有连续的分布函数有连续的分布函数F(y),F(y),则则Z Z与与Y Y有相同的分布有相同的分布.第24页/共87页1 连续型随机变量的模拟连续型随机变量的模拟 一般说来，具有给定分布的连续型随机变量可以利用在区间(0,1)上均匀分布的随机数来模拟最常用的方法是反函数法反函数法由概率论的理论可以证明，若随机变量Y有连续的分布函数F(y)，而X是区间(0,1)上均匀

17、分布的随机变量，令 ，则Z与Y有相同的分布由此，若已知Y的概率密度为 ，由 可得反函数法反函数法第25页/共87页如果给定区间(0,1)上均匀分布的随机数 ,则具有给定分布Y的随机数 可由方程 解出.例:模拟服从参数为 的指数分布时,由 可得 注：指数分布的密度函数第26页/共87页2 离散型随机变量的模拟设随机变量X的分布律为:有相同的发生的概率有相同的发生的概率.因此我们可以用随机变量因此我们可以用随机变量R落在落在小区间内的情况来模拟离散的随机变量小区间内的情况来模拟离散的随机变量X的取值情况的取值情况.第27页/共87页 因此我们可以用随机变量R落在小区间内的情况来模拟离散的随机变量X

18、的取值情况具体执行的过程是：产生一个(0,1)上均匀分布的随机数 r(简称随机数)，若p(n-1)r p(n)则理解为发生事件(X=xn)于是就可以模拟随机变量的取值情况第28页/共87页2 离散型随机变量的模拟例例 3:3:随机变量随机变量 表示每分钟到达银行柜台的表示每分钟到达银行柜台的顾客数顾客数.X.X的分布列见下表的分布列见下表,试模拟试模拟1010分钟内顾客到达柜台分钟内顾客到达柜台的情况的情况.表表1 101 10分钟内顾客到达柜台的情况分钟内顾客到达柜台的情况 Xk 0 1 2Xk 0 1 2 pk 0.4 0.3 pk 0.4 0.3 分析分析:因为每分钟到达柜台的人数是随机

19、的因为每分钟到达柜台的人数是随机的,所以可用计算所以可用计算机随机生成一组机随机生成一组(0,1)(0,1)的数据的数据,由由X X的概率分布情况的概率分布情况,可认可认为随机数在为随机数在(0,0.4)(0,0.4)范围内时没有顾客光顾范围内时没有顾客光顾,在在0.4,0.7)0.4,0.7)时时,有一个顾客光顾有一个顾客光顾,在在0.7,1)0.7,1)时时,有两个顾客光顾有两个顾客光顾.从而有从而有MATLABMATLAB程序程序:第29页/共87页2 离散型随机变量的模拟r=rand(1,10);for i=1:10;if r(i n(i)=0;elseif 0.40.1)for i=

20、1:2:7%循环求每个点的下一个时间点的坐标d=sqrt(x(i)-x(i+2)2+(x(i+1)-x(i+3)2);%与追击点的距离x(i)=x(i)+v*dt*(x(i+2)-x(i)/d;%下一个时间点的x坐标x(i+1)=x(i+1)+v*dt*(x(i+3)-x(i+1)/d;%下一个时间点的y坐标plot(x(i),x(i+1),.)endx(9)=x(1);x(10)=x(2);endhold3 MATLAB实现第44页/共87页3 MATLAB实现第45页/共87页例6 水池含盐量问题某水池有某水池有2000m2000m3 3水水,其中含盐其中含盐2kg,2kg,以以6m6m3

21、 3/min/min的速率向水池内注入含盐为的速率向水池内注入含盐为3 3的盐水的盐水,同时又以同时又以4m4m3 3/min/min的速率从水池流出搅拌均匀的盐水的速率从水池流出搅拌均匀的盐水.试用计算机试用计算机仿真该水池内盐水的变化过程仿真该水池内盐水的变化过程,并每隔并每隔10min10min计算水池中水的体积计算水池中水的体积,含盐量含盐量,含盐率含盐率.欲使池中盐水含盐率达到欲使池中盐水含盐率达到3 3,需经过多长时间需经过多长时间?分析分析:这是一个连续系统这是一个连续系统,首先要将系统离散化首先要将系统离散化,在一些离散点上进行考察在一些离散点上进行考察,这些离散点的间隔就是时

22、间步长这些离散点的间隔就是时间步长.可取步长为可取步长为1min,1min,即隔即隔1min1min考察一次系考察一次系统的状态统的状态,并相应地记录和分析并相应地记录和分析.在注入和流出的作用下在注入和流出的作用下,池中水的体积与池中水的体积与含盐量含盐量,含盐率均随时间变化含盐率均随时间变化,初始时刻含盐率为初始时刻含盐率为3 3,以后每分钟注入含盐率以后每分钟注入含盐率为为3 3的水的水6m6m3 3,流出混合均匀的盐水为流出混合均匀的盐水为4m4m3 3,当池中水的含盐率达到当池中水的含盐率达到3 3时时,仿真过仿真过程结束程结束.第46页/共87页例6 水池含盐量问题记记T T时刻的

23、体积为时刻的体积为w(w(m m3 3)，水的含盐量为，水的含盐量为s(kg)s(kg)，水的含盐率为，水的含盐率为r=s/w(kg/r=s/w(kg/m m3 3)，每隔，每隔1min1min池水的动态变化过程如下：每分钟水的体池水的动态变化过程如下：每分钟水的体积增加积增加6-4=2(6-4=2(m m3 3)；每分钟向池内注入盐；每分钟向池内注入盐60.5=3(kg)60.5=3(kg)；每分钟向池；每分钟向池外流出盐外流出盐4r(kg)4r(kg)；每分钟池内增加盐；每分钟池内增加盐3-4r(kg).3-4r(kg).本例还可以用微分方程建立数学模型，并求出它的解析解，这个解本例还可以

24、用微分方程建立数学模型，并求出它的解析解，这个解析解就是问题的精确解，有兴趣的读者可以按照这个思路求出该问析解就是问题的精确解，有兴趣的读者可以按照这个思路求出该问题的精确解，考察相应时刻精确解与仿真解的差异，还可以进一步题的精确解，考察相应时刻精确解与仿真解的差异，还可以进一步调整仿真过程的时间步长，通过与精确解的比较来研究时间步长的调整仿真过程的时间步长，通过与精确解的比较来研究时间步长的大小对仿真度的影响。大小对仿真度的影响。第47页/共87页MATLAB实现clearh=1;%时间步长为1s0=2;%初始含盐2kgw0=2000;%初始水池有水2000m3r0=s0/w0;%初始浓度s

25、(1)=s0+0.5*6*h-4*h*r0;%一分钟后的含盐量w(1)=w0+2*h;%一分钟后水池中的盐水体积r(1)=s(1)/w(1);%一分钟后的浓度 t(1)=h;y(1)=(2000000+3000000*h+3000*h2+h3)/(1000+h)2;for i=2:200 t(i)=i*h;s(i)=s(i-1)+0.5*6*h-4*h*r(i-1);%第i步后的含盐量 w(i)=w(i-1)+2*h;%第i步后的盐水体积 r(i)=s(i)/w(i);%第i步后的盐水浓度 y(i)=(2000000+3000000*t(i)+3000*t(i)2+t(i)3)/(1000+t

26、(i)2;m=floor(i/10);第48页/共87页MATLAB实现 if i/10-m0.2%若第i步后的盐水浓度大于 t02=i*h;r02=r(i);break endendt02,r0210*tm,sm,rm%表示逆subplot(1,2,1),plot(t,s,blue);hold onsubplot(1,2,2),plot(t,y,red);第49页/共87页应用举例-面积计算例例7:面积计算 求由曲线x2+y2=16,x2/36+y2=1,以及(x-2)2+(y+1)2=9所围成图形的面积。用Matlab语句作出图形得区域图 t=0:0.01:2*pi;x=sin(t);y=

27、cos(t);plot(4*x,4*y,6*x,y,2+3*x,3*y-1)axis(-6,6,-6,6)第50页/共87页应用举例-面积计算第51页/共87页应用举例-面积计算此区域形状复杂，理论分析困难，可以用计算机仿真实现。将可能的区域等分，考察每个小区域是否在此区域中，将在此区域中的小面积相加即可其仿真图如下：第52页/共87页应用举例-面积计算Matlab程序为x=-2:0.01:6;y=-2:0.01:2;s=0;h=0.01;for i=1:800 for j=1:400 xx=-2+i*h;yy=-2+j*h;if xx2+yy2=16 if xx2/36+yy2=1 if(x

28、x-2)2+(yy+1)2=9 s=s+h2;end end end endends运行后给出面积的值。第53页/共87页四:离散系统的模拟-事件步长法 离散系统离散系统(discrete system)(discrete system)是指系统状态只是指系统状态只在有限的时间点或可数的时间点上有随机事件在有限的时间点或可数的时间点上有随机事件驱动的系统例如排队系统（驱动的系统例如排队系统（queue systemqueue system），），显然状态量的变化只是在离散的随机时间点上显然状态量的变化只是在离散的随机时间点上发生假设离散系统状态的变化是在一个时间发生假设离散系统状态的变化是在一

29、个时间点上瞬间完成的点上瞬间完成的 常用的是事件步长法（下次事件推进法）其常用的是事件步长法（下次事件推进法）其过程是：置模拟时钟的初值为过程是：置模拟时钟的初值为0 0，跳到第一个事，跳到第一个事件发生的时刻，计算系统的状态，产生未来事件发生的时刻，计算系统的状态，产生未来事件并加入到队列中去，跳到下一事件，计算系件并加入到队列中去，跳到下一事件，计算系统状态，统状态，重复这一过程直到满足某个终，重复这一过程直到满足某个终止条件为止止条件为止第54页/共87页例7 收款台前的排队过程假设假设:(1)(1)顾客到达收款台是随机的顾客到达收款台是随机的,平均时间间隔为平均时间间隔为0.5min,

30、0.5min,即间隔时间服从即间隔时间服从lamda=2lamda=2的指数分布的指数分布.(2)(2)对不同的顾客收款和装袋的时间服从正态分布对不同的顾客收款和装袋的时间服从正态分布N(1,1/3).N(1,1/3).试模拟试模拟2020位顾客到收款台前的排队情况位顾客到收款台前的排队情况,我们关心的问我们关心的问题是每个顾客的平均等待时间题是每个顾客的平均等待时间,队长及服务员的工作效率队长及服务员的工作效率第55页/共87页例7 收款台前的排队过程分析分析:单服务台结构的排队系统有两类原发事件即单服务台结构的排队系统有两类原发事件即到来和离去到来和离去,顾客到来的顾客到来的后继事件是后继

31、事件是顾客接受服务顾客接受服务,顾客离去的后继事件是顾客离去的后继事件是服务台寻找服务服务台寻找服务,这四类这四类事件各自的子程序框图如图事件各自的子程序框图如图1 1所示所示.假设假设:t(i):t(i)为第为第i i位顾客位顾客到达时刻到达时刻;t;t2 2(i)(i)为第为第i i位顾客受到位顾客受到服务的时间服务的时间(随机随机变量变量);T(i);T(i)为第为第i i位顾客位顾客离去时刻离去时刻.将第将第i i位顾客到达作为事件发生位顾客到达作为事件发生:t(i+1)-t(i)=r(i)(:t(i+1)-t(i)=r(i)(随机变量随机变量););平衡关系平衡关系:当当t(i+1)

32、t(i+1)T(i)T(i)时时,T(i+1)=t(i+1)+t,T(i+1)=t(i+1)+t2 2(i+1);(i+1);否则否则,T(i+1)=,T(i+1)=T(i)+tT(i)+t2 2(i+1).(i+1).第56页/共87页图1 到来事件子程序 系统人数+1产生下一个顾客到来时刻调用接受服务事件的程序第57页/共87页图 2 离去事件子程序 系统人数-1置服务台”闲”已服务人数+1调用寻找服务事件子程序第58页/共87页图3:接受服务事件子程序 服务台空置服务台”忙”产生服务结束时刻登记到事件表排队人数+1否是第59页/共87页图4:寻找服务事件子程序 排队人数1置服务台”忙”产

33、生服务结束时刻排队人数-1排队人数+1否是第60页/共87页五：蒙特卡洛方法在用传统方法难以解决的问题中，有很大一部分可以用概率在用传统方法难以解决的问题中，有很大一部分可以用概率模型进行描述由于这类模型含有不确定的随机因素，分析模型进行描述由于这类模型含有不确定的随机因素，分析起来通常比确定性的模型困难有的模型难以作定量分析，起来通常比确定性的模型困难有的模型难以作定量分析，得不到解析的结果，或者是虽有解析结果，但计算代价太大得不到解析的结果，或者是虽有解析结果，但计算代价太大以至不能使用在这种情况下，可以考虑采用以至不能使用在这种情况下，可以考虑采用Monte CarloMonte Car

34、lo方方法。法。Monte CarloMonte Carlo方法是计算机模拟的基础，它的名字来源于方法是计算机模拟的基础，它的名字来源于世界著名的赌城世界著名的赌城摩纳哥的蒙特卡洛，其历史起源于摩纳哥的蒙特卡洛，其历史起源于17771777年法国科学家蒲丰提出的一种计算圆周年法国科学家蒲丰提出的一种计算圆周的方法的方法随机投随机投针法，即著名的针法，即著名的蒲丰投针问题。蒲丰投针问题。第61页/共87页18世纪，法国数学家蒲丰和勒可莱尔提出的“投针问题”，记载于蒲丰1777年出版的著作中：“在平面上画有一组间距为a的平行线，将一根长度为L（L=a/2）的针任意掷在这个平面上，求此针与平行线中任

35、一条相交的概率。”蒲丰本人证明了，这个概率是：p=2L/（a）为圆周率第62页/共87页利用这个公式可以用概率的方法得到圆周率的近似值。下面是一些资料试验者试验者时间时间投掷次数投掷次数相交次数相交次数圆周率估计圆周率估计值值Wolf1850年500025323.1596Smith1855年32041218.53.1554C.De Morgan1860年600382.53.137Fox1884年10304893.1595Lazzerini1901年340818083.1415929Reina1925年25208593.1795第63页/共87页五：蒙特卡洛方法Monte CarloMonte

36、Carlo方法的基本思想是首先建立一个概率模型，使所方法的基本思想是首先建立一个概率模型，使所求问题的解正好是该模型的参数或其他有关的特征量然后求问题的解正好是该模型的参数或其他有关的特征量然后通过模拟一统计试验，即多次随机抽样试验（确定通过模拟一统计试验，即多次随机抽样试验（确定m m和和n n），），统计出某事件发生的百分比只要试验次数很大，该百分比统计出某事件发生的百分比只要试验次数很大，该百分比便近似于事件发生的概率这实际上就是概率的统计定义便近似于事件发生的概率这实际上就是概率的统计定义利用建立的概率模型，求出要估计的参数蒙特卡洛方法属利用建立的概率模型，求出要估计的参数蒙特卡洛方法

37、属于试验数学的一个分支于试验数学的一个分支第64页/共87页五：蒙特卡洛方法蒙特卡洛方法适用范围很广泛，它既能求解确定性的问题，也能求解随机性的问题以及科学研究中的理论蒙特卡洛方法适用范围很广泛，它既能求解确定性的问题，也能求解随机性的问题以及科学研究中的理论问题例如利用蒙特卡洛方法可以近似地计算定积分，即产生数值积分问题问题例如利用蒙特卡洛方法可以近似地计算定积分，即产生数值积分问题第65页/共87页蒙特卡洛法求圆周率clearn=50000X=rand(n,1);Y=rand(n,1);k=0;for i=1:n;if X(i)2+Y(i)2a;fprintf(error：针长必须小于等于

38、%dn,a);fprintf(请重新调用函数pinpi(k,d,l)n);pi_value=0;else for i=1:ka*rand(1)pinpi(100000,4,3)投针法求得pi=3.143797e+000ans=第68页/共87页任意曲边梯形面积的近似计算一个古老的问题：用一堆石头测量一个水塘的面积应该怎样做呢？测量方法如下：假定水塘位于一块面一个古老的问题：用一堆石头测量一个水塘的面积应该怎样做呢？测量方法如下：假定水塘位于一块面积已知的矩形农田之中如图积已知的矩形农田之中如图8 82 2所示随机地向这块农田扔石头使得它们都落在农田内被扔到农田中所示随机地向这块农田扔石头使得它

39、们都落在农田内被扔到农田中的石头可能溅上了水，也可能没有溅上水，估计被的石头可能溅上了水，也可能没有溅上水，估计被“溅上水的溅上水的”石头量占总的石头量的百分比试想如何石头量占总的石头量的百分比试想如何利用这估计的百分比去近似计算该水塘面积？利用这估计的百分比去近似计算该水塘面积？第69页/共87页任意曲边梯形面积的近似计算第70页/共87页任意曲边梯形面积的近似计算结合图结合图8 82 2中的图形中的图形(1)(1)分析，只要已知各种参数及函数（分析，只要已知各种参数及函数（a a，b b，H H，f(x)f(x)），有以下两种方法可近似计算水塘面积），有以下两种方法可近似计算水塘面积 1

40、1随机投点法随机投点法 1 1）赋初值：试验次数）赋初值：试验次数n=0n=0，成功次数，成功次数m=0m=0；规定投点试验的总次数；规定投点试验的总次数N N；2 2）随机选择）随机选择m m个数对个数对xi,yi,1ixi,yi,1im m，其中，其中a a xixi b,b,00yiyi H H，置，置 n nn nl l；3 3）判断）判断nNnN，若是，转，若是，转4 4，否则停止计算；，否则停止计算；4 4）判断条件）判断条件 (表示一块溅水的石头表示一块溅水的石头)是否成立，若成立则置是否成立，若成立则置m=m+1m=m+1，转，转2 2，否则转，否则转2 2；5 5）计算水塘面

41、积的近似值）计算水塘面积的近似值S=HS=H(b-a)m/(b-a)m/N N 第71页/共87页任意曲边梯形面积的近似计算2 2平均值估计法平均值估计法1 1）产生）产生a,ba,b区间的均匀随机数区间的均匀随机数x x,i=1,2,i=1,2,N,N 2)2)计算计算f(xif(xi），），i=1,2,Ni=1,2,N3 3）计算）计算该方法的特点是估计函数该方法的特点是估计函数f(x)f(x)在在a,ba,b上的平均值，面积近似等于该平均值乘以上的平均值，面积近似等于该平均值乘以(b-a)(b-a)第72页/共87页例9 库存问题在物资的供应过程中在物资的供应过程中,由于到货和销售不可能

42、做到同步同量由于到货和销售不可能做到同步同量,故故总要保持一定的库存储备总要保持一定的库存储备.如果库存过多如果库存过多,就会造成积压浪费以就会造成积压浪费以及保管费的上升及保管费的上升;如果库存过少如果库存过少,会造成缺货会造成缺货.如何选择库存和订如何选择库存和订货策略货策略,就是一个需要研究的问题就是一个需要研究的问题.库存问题有多种类型库存问题有多种类型,一般都一般都比较复杂比较复杂,下面讨论一种简单的情形下面讨论一种简单的情形.某电动车行的仓库管理人员采取一种简单的订货策略某电动车行的仓库管理人员采取一种简单的订货策略,当库存当库存降低到降低到P P辆电动车时就向厂家辆电动车时就向厂

43、家第73页/共87页例9 库存问题订货订货,每次订货每次订货Q Q辆辆,如果某一天的需求量超过了库存量如果某一天的需求量超过了库存量,商店就商店就有销售损失和信誉损失有销售损失和信誉损失,但如果库存量过多但如果库存量过多,就会导致资金积压就会导致资金积压和保管费增加和保管费增加.若现在有如下表所示的两种库存策略若现在有如下表所示的两种库存策略,试比较选试比较选择一种策略以使总费用最少择一种策略以使总费用最少.表表 订货方案订货方案 方案方案 重新订货点重新订货点P/P/辆辆 重新订货量重新订货量Q/Q/辆辆 方案方案1 125 1501 125 150 方案方案2 150 2502 150 2

44、50第74页/共87页例9 库存问题这个问题的已知条件是这个问题的已知条件是:(1)(1)从发出订货到收到货物需隔从发出订货到收到货物需隔3 3天天.(2)(2)每辆电动车保管费为元每辆电动车保管费为元/天天,每辆电动车的缺货损失为元每辆电动车的缺货损失为元/天天,每次的订货费为每次的订货费为7575元元.(3)(3)每天电动车需求量是每天电动车需求量是0 0到到9999之间均匀分布的随机数之间均匀分布的随机数.(4)(4)原始库存为原始库存为110110辆辆,假设第一天没有发出订货假设第一天没有发出订货.第75页/共87页例9 库存问题分析分析:这一问题用解析法讨论比较麻烦这一问题用解析法讨

45、论比较麻烦,但用计算机按天仿真仓但用计算机按天仿真仓库货物的变动情况却很方便库货物的变动情况却很方便.我们以我们以3030天为例天为例,依次对这两种方依次对这两种方案进行仿真案进行仿真,最后比较各方案的总费用最后比较各方案的总费用,从而就可以作出决策从而就可以作出决策.计算机仿真时的工作流程是早上到货计算机仿真时的工作流程是早上到货,全天销售全天销售,晚上订货晚上订货,输入输入一些常数和初始数据后一些常数和初始数据后,以一天为时间步长进行仿真以一天为时间步长进行仿真.首先检查首先检查这一天是否为预订到货日期这一天是否为预订到货日期,如果是如果是,则原有库存量加则原有库存量加Q,Q,并把预并把预

46、定到货量清为定到货量清为0;0;如果不是如果不是,则库存量不变则库存量不变.接着用计算机中的接着用计算机中的随随机函数仿真随机需求量机函数仿真随机需求量,若库存量大于需求量若库存量大于需求量,则新的库存量减则新的库存量减去需求量去需求量;反之反之,则新库存量变为则新库存量变为0,0,并且要在总费用上加缺货损并且要在总费用上加缺货损失失,然后检查实际然后检查实际第76页/共87页例9 库存问题库存量加上预定到货量是否小于重新订货点库存量加上预定到货量是否小于重新订货点P,P,如果是如果是,则需要重则需要重新订货新订货,这时就加一次订货费这时就加一次订货费.如此重复运行如此重复运行3030天天,即

47、可得所需费即可得所需费用总值用总值.由此比较这两种方案的总费用由此比较这两种方案的总费用,即得最好方案即得最好方案.第77页/共87页cleardays=30;P=125,150;%重新订货点P，单位辆Q=150,250;%重新订货量Q，单位辆cost=0,0;%花费初值arrivalinterval=2;%到货时间间隔为2天storagefee=0.5;lossfee=1.6;bookfee=75;%各种费用storage0=110;booknumber=0;arrivedate=0;nr=rand(days,1);%随机产生30天的随机数for i=1:2%两种方案分别仿真%下面是第一天的

48、情况 storage(1)=storage0;n=round(99*nr(1);%30天中每天的需求量,0-99之间的均匀分布随机数 sale=n;remain=storage(1)-n;%计算库存 if remain=P(i);%第i种方案，库存少于重新订货量，就需要订货 booknumber=Q(i);%订货数量 arrivedate=4;%到货日期 orderfee=bookfee;%订单费用 else orderfee=0;%未出现缺货情况，订单费用为0 end storage(1)=remain;%更新货存 cost(i)=cost(i)+remain*storagefee+orde

49、rfee;%总费用累加，第一天计算结束 for j=2:days%第2天到第30天仿真 dh=j;if abs(dh-arrivedate)=n%库存足够 sale=n;%成交量 remain=storage(j)-n;%库存量 shortagenumber=0;%缺货数为0 else%否则，库存不够的处理 sale=storage(j);%只能销售已有货存的量，有缺货 remain=0;%库存卖完 shortagenumber=n-storage(j);%缺货量 end storage(j)=remain;%更新库存 if remain+booknumber=P(i);%库存量加上预定到货量

50、小于重新订货点P，重新订货 booknumber=Q(i);%订货数 arrivedate=dh+arrivalinterval;%到货日期 orderfee=bookfee;%订单费用 else orderfee=0;%无需重新订货，未发生费用 end cost(i)=cost(i)+remain*storagefee+shortagenumber*lossfee+orderfee;%总费用增加 end;mincost=min(cost);%求两个方案30天总费用最小值 end cost mincost第79页/共87页例10 赶火车过程仿真一列火车从一列火车从A A站经站经B B站开往站开