1、数学建模之计算机仿真2概述l计算机科学技术的迅猛发展，给许多学科带来了巨大计算机科学技术的迅猛发展，给许多学科带来了巨大的影响计算机不但使问题的求解变得更加方便、快的影响计算机不但使问题的求解变得更加方便、快捷和精确，而且使得解决实际问题的领域更加广捷和精确，而且使得解决实际问题的领域更加广泛计算机适合于解决那些规模大、难以解析化以及泛计算机适合于解决那些规模大、难以解析化以及不确定的数学模型例如对于一些带随机因素的复杂不确定的数学模型例如对于一些带随机因素的复杂系统，用分析方法建模常常需要作许多简化假设，与系统，用分析方法建模常常需要作许多简化假设，与面临的实际问题可能相差甚远，以致解答根本
2、无法应面临的实际问题可能相差甚远，以致解答根本无法应用，这时仿真几乎成为人们的唯一选择在历届的美用，这时仿真几乎成为人们的唯一选择在历届的美国和中国大学生的数学建模竞赛（国和中国大学生的数学建模竞赛（MCM）中，学生）中，学生们经常用到计算机仿真方法去求解、检验等计算机们经常用到计算机仿真方法去求解、检验等计算机仿真仿真(computer simulation)是建模过程中较为重要的是建模过程中较为重要的一类方法一类方法 3计算机仿真的基本概念计算机仿真的基本概念 计算机仿真，是根据已知的信息和计算机仿真，是根据已知的信息和知识（如数学、物理规律），利用计算知识（如数学、物理规律），利用计算机
3、模拟现实情况和系统的演变过程，发机模拟现实情况和系统的演变过程，发现新的知识或规律，从而解决问题的一现新的知识或规律，从而解决问题的一种方法种方法4计算机仿真的基本概念计算机仿真的基本概念 独立于理论研究与实验研究的认识独立于理论研究与实验研究的认识世界的第三中方法世界的第三中方法5计算机仿真的基本概念计算机仿真的基本概念计算机仿真的特点计算机仿真的特点 代价小、时间短、可重复、参数代价小、时间短、可重复、参数设置灵活设置灵活6计算机仿真可以解决以下计算机仿真可以解决以下5类问题类问题:l(1)(1)难以用数学公式表示的系统难以用数学公式表示的系统,或者没有建立和求解的或者没有建立和求解的有效
4、方法有效方法.l(2)(2)虽然可以用解析的方法解决问题虽然可以用解析的方法解决问题,但数学的分析与计算但数学的分析与计算过于复杂过于复杂,这时计算机仿真可能提供简单可行的求解方法这时计算机仿真可能提供简单可行的求解方法.l(3)(3)希望能在较短的时间内观察到系统发展的全过程希望能在较短的时间内观察到系统发展的全过程,以估以估计某些参数对系统行为的影响计某些参数对系统行为的影响.l(4)(4)难以在实际环境中进行试验和观察时难以在实际环境中进行试验和观察时,计算机仿真是唯计算机仿真是唯一可行的方法一可行的方法,如太空飞行的研究如太空飞行的研究.l(5)(5)需要对系统或过程进行长期运行比较需
5、要对系统或过程进行长期运行比较,从大量方案中寻从大量方案中寻找最优方案找最优方案.7 计算机仿真的分类计算机仿真的分类l计算机仿真在计算机中运行实现计算机仿真在计算机中运行实现,不怕破坏不怕破坏,易修改易修改,可重用可重用,安全经济安全经济,不受外界条件和场地空间的限制不受外界条件和场地空间的限制.l仿真仿真分为静态分为静态仿真仿真(static simulation)和动态和动态仿真仿真(dynamic simulation).数值积分中的蒙特卡洛方法数值积分中的蒙特卡洛方法（统计模拟方法）是典型的静态（统计模拟方法）是典型的静态仿真仿真动态动态仿真仿真又分又分为连续系统为连续系统仿真仿真和
6、离散系统和离散系统仿真仿真连续系统是指状态连续系统是指状态变量随着时间连续变化的系统，例如传染病的检测与变量随着时间连续变化的系统，例如传染病的检测与预报系统预报系统.离散系统是指系统状态变量只在有限的时离散系统是指系统状态变量只在有限的时间点或可数的时间点上发生变化的系统间点或可数的时间点上发生变化的系统,例如排队系例如排队系统统.8l 仿真系统，必须设置一个仿真时钟仿真系统，必须设置一个仿真时钟(simulate clock)，它能将时间从一个时刻，它能将时间从一个时刻向另一个时刻进行推进，并且能随时反映向另一个时刻进行推进，并且能随时反映系统时间的当前值其中，模拟时间推进系统时间的当前值
7、其中，模拟时间推进方式有两种方式有两种:时间步长法时间步长法(均匀间隔时间推均匀间隔时间推进法进法,连续系统常用连续系统常用)和事件步长法和事件步长法(下次事下次事件推进法件推进法,离散系统常用离散系统常用)910主要内容主要内容l一一:准备知识准备知识:随机数的产生随机数的产生l二二:随机变量的模拟随机变量的模拟l三三:连续系统的模拟连续系统的模拟-时间步长法时间步长法l四四:离散系统的模拟离散系统的模拟-事件步长法事件步长法五：蒙特卡洛方法五：蒙特卡洛方法11一一:准备知识准备知识:随机数的产生随机数的产生l由于仿真研究的实际系统要受到多种随机因素的由于仿真研究的实际系统要受到多种随机因素
8、的作用和影响作用和影响,在仿真过程中必须处理大量的随机因在仿真过程中必须处理大量的随机因素素.要解决此问题的前提是确定随机变量的类型和要解决此问题的前提是确定随机变量的类型和选择合适的随机数产生的方法选择合适的随机数产生的方法.l对随机现象进行模拟对随机现象进行模拟,实质是要给出随机变量的模实质是要给出随机变量的模拟拟,也就是说要利用计算机随机产生一系列数值也就是说要利用计算机随机产生一系列数值,使它们服从一定的概率分布使它们服从一定的概率分布,称这些数值为随机数称这些数值为随机数.l最基本最基本,最常用的是最常用的是(0,1)(0,1)区间内均匀分布的随机区间内均匀分布的随机数数.其他分布的
9、随机数均可利用它来产生其他分布的随机数均可利用它来产生.121:1:产生模拟随机数的计算机命令产生模拟随机数的计算机命令l在在MATLAB中中,可以直接产生满足各种分布的随机可以直接产生满足各种分布的随机数数,命令如下命令如下:l常见的分布函数常见的分布函数 MATLAB语句语句 l均匀分布均匀分布U0,1 R=rand(m,n)l均匀分布均匀分布Ua,b R=unifrnd(a,b,m,n)l指数分布指数分布E()R=exprnd(,m,n)l正态分布正态分布N(mu,sigma)R=normrnd(mu,sigma,m,n)l标准正态分布标准正态分布N(0,1)R=randn(m,n)l二
10、项分布二项分布B(n,p)R=binornd(n,p,m,n)l泊松分布泊松分布 P()R=poissrnd(,m,n)l以上语句均产生以上语句均产生m n 的矩阵的矩阵132:案例分析案例分析l例1:unifrnd(2,3)lunifrnd(1,32,1,4)lnormrnd(1,2)l normrnd(1,2,2,3)lrand(2,3)lrandn(2,3)142:案例分析案例分析lans=2.8132lans=1.3057 5.3056 7.2857 7.1604lans=0.2527lans=l 2.7429 0.0219 2.7759l 2.2756 0.0992 -0.9560l
11、ans=l 0.6038 0.1988 0.7468l 0.2722 0.0153 0.4451lans=l -0.0945 -1.3089 -0.2440l -0.2141 0.8248 -0.1778l 152:案例分析案例分析l例例2:2:敌空战部队对我方港口进行空袭敌空战部队对我方港口进行空袭,其到达规律服从泊其到达规律服从泊松分布松分布,平均每分钟到达平均每分钟到达4 4架飞机架飞机.l(1)(1)模拟敌机在模拟敌机在3 3分钟内到达目标区域的数量分钟内到达目标区域的数量,以及在第以及在第1,2,31,2,3分钟内各到达几架飞机分钟内各到达几架飞机;l(2)(2)模拟在模拟在3 3分
12、钟内每架飞机的到达时刻分钟内每架飞机的到达时刻.分析分析:(1)n1=poissrnd(4),n2=poissrnd(4),:(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的指数分布的指数分布,故可由故可由指数分布模拟每架飞机的到达时刻指数分布模拟每架飞机的到达时刻.162:案例分析案例分析lclear
13、lt=0;lj=0;%到达的飞机数 lwhile t3l j=j+1l t=t+exprnd(1/4)lend17二二:随机变量的模拟随机变量的模拟l利用均匀分布的随机数可以产生具有任意分布的随机变量利用均匀分布的随机数可以产生具有任意分布的随机变量的样本的样本,从而可以对随机变量的取值情况进行模拟从而可以对随机变量的取值情况进行模拟.l1 1 连续型随机变量的模拟连续型随机变量的模拟l具有给定分布的连续型随机变量可以利用在区间具有给定分布的连续型随机变量可以利用在区间(0,1)(0,1)上上均匀分布的随机数来模拟均匀分布的随机数来模拟,最常用的方法是逆变换法最常用的方法是逆变换法.l结论结论
14、:若随机变量若随机变量Y Y有连续的分布函数有连续的分布函数F(y),F(y),l 则则Z Z与与Y Y有相同的分布有相同的分布.1(0,1),()XUZFX令181 1 连续型随机变量的模拟连续型随机变量的模拟l若已知Y的概率密度为l如果给定区间(0,1)上均匀分布的随机数 ,则具有给定分布Y的随机数 可由方程 l 解出.l例:模拟服从参数为 的指数分布时,由 可得 (),()(),yf yF Yf y dy有()iyirfyd y iy()1iiyyirfy dye11ln(1),ln.iiiiyryr 也可简化为ir192 离散型随机变量的模拟离散型随机变量的模拟l设随机变量X的分布律为
15、:(1,2,3,),iiP Xxp i 令(0)()()1(1)()0,1,2,3,(0,1)(,).(0,1),nnniinnppp npppRU,将作为分点 将区间分为一系列小区间对于均匀的随机变量则有(1)()()(1)(1)(),1,2,nnnnnnnPpRpppnpRpXx由 此 可 知 事 件和 事 件有相同的发生的概率有相同的发生的概率.因此我们可以用随机变量因此我们可以用随机变量R落在落在小区间内的情况来模拟离散的随机变量小区间内的情况来模拟离散的随机变量X的取值情况的取值情况.202 离散型随机变量的模拟离散型随机变量的模拟l例例 3:3:随机变量随机变量 表示每分钟到达银行
16、柜台的顾客表示每分钟到达银行柜台的顾客数数.X.X的分布列见下表的分布列见下表,试模拟试模拟1010分钟内顾客到达柜台的情况分钟内顾客到达柜台的情况.l 表表1 101 10分钟内顾客到达柜台的情况分钟内顾客到达柜台的情况l Xk 0 1 2Xk 0 1 2l pk 0.4 0.3 0.3 pk 0.4 0.3 0.3l分析分析:因为每分钟到达柜台的人数是随机的因为每分钟到达柜台的人数是随机的,所以可用计算机所以可用计算机随机生成一组随机生成一组(0,1)(0,1)的数据的数据,由由X X的概率分布情况的概率分布情况,可认为随机可认为随机数在数在(0,0.4)(0,0.4)范围内时没有顾客光顾
17、范围内时没有顾客光顾,在在0.4,0.7)0.4,0.7)时时,有一个有一个顾客光顾顾客光顾,在在0.7,1)0.7,1)时时,有两个顾客光顾有两个顾客光顾.l 从而有从而有MATLABMATLAB程序程序:0,1,2X 212 离散型随机变量的模拟离散型随机变量的模拟lr=rand(1,10);lfor i=1:10;l if r(i)0.4l n(i)=0;l elseif 0.4=r(i)&r(i)0.1)lfor i=1:2:7ld=sqrt(x(i)-x(i+1)2+(x(i+1)-x(i+3)2);lx(i)=x(i)+v*dt*(x(i+2)-x(i)/d;lx(i+1)=x(i
18、+1)+v*dt*(x(i+3)-x(i+1)/d;lplot(x(i),x(i+1),.)lendlx(9)=x(1);x(10)=x(2);lendlhold3 MATLAB实现3002468100123456789103 MATLAB实现31例例9.6 水池含盐量问题水池含盐量问题l某水池有某水池有2000m2000m3 3水水,其中含盐其中含盐2kg,2kg,以以6m6m3 3/min/min的速率向水池的速率向水池内注入含盐为内注入含盐为0.5kg/m0.5kg/m3 3的盐水的盐水,同时又以同时又以4m4m3 3/min/min的速率从的速率从水池流出搅拌均匀的盐水水池流出搅拌均匀
19、的盐水.试用计算机仿真该水池内盐水试用计算机仿真该水池内盐水的变化过程的变化过程,并每隔并每隔10min10min计算水池中水的体积计算水池中水的体积,含盐量含盐量,含含盐率盐率.欲使池中盐水含盐率达到欲使池中盐水含盐率达到0.2kg/m0.2kg/m3 3,需经过多长时间需经过多长时间?l分析分析:这是一个连续系统这是一个连续系统,首先要将系统离散化首先要将系统离散化,在一些离在一些离散点上进行考察散点上进行考察,这些离散点的间隔就是时间步长这些离散点的间隔就是时间步长.可取步可取步长为长为1min,1min,即隔即隔1min1min考察一次系统的状态考察一次系统的状态,并相应地记录和并相应
20、地记录和分析分析.在注入和流出的作用下在注入和流出的作用下,池中水的体积与含盐量池中水的体积与含盐量,含含盐率均随时间变化盐率均随时间变化,初始时刻含盐率为初始时刻含盐率为0.001kg/m0.001kg/m3 3,以后每以后每分钟注入含盐率为分钟注入含盐率为0.5kg/m0.5kg/m3 3的水的水6m6m3 3,流出混合均匀的盐水流出混合均匀的盐水为为4m4m3 3,当池中水的含盐率达到当池中水的含盐率达到0.2kg/m0.2kg/m3 3时时,仿真过程结束仿真过程结束.32例例9.6 水池含盐量问题水池含盐量问题l记记T T时刻的体积为时刻的体积为w(w(m m3 3)，水的含盐量为，水
21、的含盐量为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)；每分钟向；每分钟向池内注入盐池内注入盐6 60.5=3(kg)0.5=3(kg)；每分钟向池外流出盐；每分钟向池外流出盐4 4r(kg)r(kg)；每分钟池内增加盐；每分钟池内增加盐3-43-4r(kg).r(kg).l本例还可以用微分方程建立数学模型，并求出它的解本例还可以用微分方程建立数学模型，并求出它的解析解，这个解析解就是问题的精
22、确解，有兴趣的读者析解，这个解析解就是问题的精确解，有兴趣的读者可以按照这个思路求出该问题的精确解，考察相应时可以按照这个思路求出该问题的精确解，考察相应时刻精确解与仿真解的差异，还可以进一步调整仿真过刻精确解与仿真解的差异，还可以进一步调整仿真过程的时间步长，通过与精确解的比较来研究时间步长程的时间步长，通过与精确解的比较来研究时间步长的大小对仿真度的影响。的大小对仿真度的影响。33MATLAB实现实现lclearlh=1;%时间步长为1ls0=2;%初始含盐2kglw0=2000;%初始水池有水2000m3lr0=s0/w0;%初始浓度ls(1)=s0+0.5*6*h-4*h*r0;%一分
23、钟后的含盐量lw(1)=w0+2*h;%一分钟后水池中的盐水体积lr(1)=s(1)/w(1);%一分钟后的浓度 lt(1)=h;ly(1)=(2000000+3000000*h+3000*h2+h3)/(1000+h)2;lfor i=2:200l t(i)=i*h;l s(i)=s(i-1)+0.5*6*h-4*h*r(i-1);%第i步后的含盐量l w(i)=w(i-1)+2*h;%第i步后的盐水体积l r(i)=s(i)/w(i);%第i步后的盐水浓度l y(i)=(2000000+3000000*t(i)+3000*t(i)2+t(i)3)/(1000+t(i)2;l m=floor
24、(i/10);34MATLAB实现实现l if i/10-m0.2%若第i步后的盐水浓度大于0.2l t02=i*h;l r02=r(i);l break l endlendlt02,r02l10*tm,sm,rm%表示逆lsubplot(1,2,1),plot(t,s,blue);lhold onlsubplot(1,2,2),plot(t,y,red);35四四:离散系统的模拟离散系统的模拟-事件步长法事件步长法 l离散系统离散系统(discrete system)(discrete system)是指系统状态只在有是指系统状态只在有限的时间点或可数的时间点上有随机事件驱动的系限的时间点或
25、可数的时间点上有随机事件驱动的系统例如排队系统（统例如排队系统（queue systemqueue system），显然状态量），显然状态量的变化只是在离散的随机时间点上发生假设离散的变化只是在离散的随机时间点上发生假设离散系统状态的变化是在一个时间点上瞬间完成的系统状态的变化是在一个时间点上瞬间完成的 l常用的是事件步长法（下次事件推进法）其过程常用的是事件步长法（下次事件推进法）其过程是：置模拟时钟的初值为是：置模拟时钟的初值为0 0，跳到第一个事件发生，跳到第一个事件发生的时刻，计算系统的状态，产生未来事件并加入到的时刻，计算系统的状态，产生未来事件并加入到队列中去，跳到下一事件，计算系
26、统状态，队列中去，跳到下一事件，计算系统状态，重复这一过程直到满足某个终止条件为止重复这一过程直到满足某个终止条件为止36例例9.7 收款台前的排队过程收款台前的排队过程l假设假设:(1)(1)顾客到达收款台是随机的顾客到达收款台是随机的,平均时间间平均时间间隔为隔为0.5min,0.5min,即间隔时间服从即间隔时间服从lamda=2lamda=2的指的指数分布数分布.l(2)(2)对不同的顾客收款和装袋的时间服从对不同的顾客收款和装袋的时间服从正态分布正态分布N(1,1/3).N(1,1/3).l试模拟试模拟2020位顾客到收款台前的排队情况位顾客到收款台前的排队情况,我我们关心的问题是每
27、个顾客的平均等待时间们关心的问题是每个顾客的平均等待时间,队长及服务员的工作效率队长及服务员的工作效率37例例9.7 收款台前的排队过程收款台前的排队过程l分析分析:单服务台结构的排队系统有两类原发事件即到来和单服务台结构的排队系统有两类原发事件即到来和离去离去,顾客到来的后继事件是顾客接受服务顾客到来的后继事件是顾客接受服务,顾客离去的后顾客离去的后继事件是服务台寻找服务继事件是服务台寻找服务,这四类事件各自的子程序框图这四类事件各自的子程序框图如图如图1 1所示所示.l假设假设:t(i):t(i)为第为第i i位顾客到达时刻位顾客到达时刻;t;t2 2(i)(i)为第为第i i位顾客受到位
28、顾客受到服务的时间服务的时间(随机变量随机变量);T(i);T(i)为第为第i i位顾客离去时刻位顾客离去时刻.l将第将第i i位顾客到达作为件事发生位顾客到达作为件事发生:t(i+1)-t(i)=r(i)(:t(i+1)-t(i)=r(i)(随机随机变量变量););平衡关系平衡关系:当当t(i+1)t(i+1)T(i)T(i)时时,T(i+1)=t(i+1)+,T(i+1)=t(i+1)+t t2 2(i+1);(i+1);否则否则,T(i+1)=T(i)+t,T(i+1)=T(i)+t2 2(i+1).(i+1).38图图1 到来事件子程序到来事件子程序l 系统人数+1产生下一个顾客到来时
29、刻调用接受服务事件的程序39图图 2 离去事件子程序离去事件子程序 系统人数-1置服务台”闲”已服务人数+1调用寻找服务事件子程序40图图3:接受服务事件子程序接受服务事件子程序 服务台空置服务台”忙”产生服务结束时刻登记到事件表排队人数+1否是41图图4:寻找服务事件子程序寻找服务事件子程序 排队人数1置服务台”忙”产生服务结束时刻排队人数-1排队人数+1否是42五：蒙特卡洛方法五：蒙特卡洛方法l在用传统方法难以解决的问题中，有很大一部在用传统方法难以解决的问题中，有很大一部分可以用概率模型进行描述由于这类模型含分可以用概率模型进行描述由于这类模型含有不确定的随机因素，分析起来通常比确定性有
30、不确定的随机因素，分析起来通常比确定性的模型困难有的模型难以作定量分析，得不的模型困难有的模型难以作定量分析，得不到解析的结果，或者是虽有解析结果，但计算到解析的结果，或者是虽有解析结果，但计算代价太大以至不能使用在这种情况下，可以代价太大以至不能使用在这种情况下，可以考虑采用考虑采用Monte CarloMonte Carlo方法。方法。Monte CarloMonte Carlo方法方法是计算机模拟的基础，它的名字来源于世界著是计算机模拟的基础，它的名字来源于世界著名的赌城名的赌城摩纳哥的蒙特卡洛，其历史起源摩纳哥的蒙特卡洛，其历史起源于于17771777年法国科学家蒲丰提出的一种计算圆周
31、年法国科学家蒲丰提出的一种计算圆周的方法的方法随机投针法，即著名的随机投针法，即著名的蒲丰投针蒲丰投针问题。问题。43五：蒙特卡洛方法五：蒙特卡洛方法lMonte CarloMonte Carlo方法的基本思想是首先建立一个概方法的基本思想是首先建立一个概率模型，使所求问题的解正好是该模型的参数率模型，使所求问题的解正好是该模型的参数或其他有关的特征量然后通过模拟一统计试或其他有关的特征量然后通过模拟一统计试验，即多次随机抽样试验（确定验，即多次随机抽样试验（确定m m和和n n），统计），统计出某事件发生的百分比只要试验次数很大，出某事件发生的百分比只要试验次数很大，该百分比便近似于事件发生
32、的概率这实际上该百分比便近似于事件发生的概率这实际上就是概率的统计定义利用建立的概率模型，就是概率的统计定义利用建立的概率模型，求出要估计的参数蒙特卡洛方法属于试验数求出要估计的参数蒙特卡洛方法属于试验数学的一个分支学的一个分支44五：蒙特卡洛方法五：蒙特卡洛方法l蒙特卡洛方法适用范围很广泛，它既能求解确蒙特卡洛方法适用范围很广泛，它既能求解确定性的问题，也能求解随机性的问题以及科学定性的问题，也能求解随机性的问题以及科学研究中的理论问题例如利用蒙特卡洛方法可研究中的理论问题例如利用蒙特卡洛方法可以近似地计算定积分，即产生数值积分问题以近似地计算定积分，即产生数值积分问题45蒙特卡洛法求圆周率
33、蒙特卡洛法求圆周率lclearln=50000lX=rand(n,1);lY=rand(n,1);lk=0;lfor i=1:n;l if X(i)2+Y(i)2a;l fprintf(error：针长必须小于等于%dn,a);l fprintf(请重新调用函数pinpi(k,d,l)n);l pi_value=0;lelse lfor i=1:kl if a*rand(1)pinpi(100000,4,3)l投针法求得pi=3.143797e+000lans=l 3.1437972879508748任意曲边梯形面积的近似计算任意曲边梯形面积的近似计算l一个古老的问题：用一堆石头测量一个水塘的
34、一个古老的问题：用一堆石头测量一个水塘的面积应该怎样做呢？测量方法如下：假定水面积应该怎样做呢？测量方法如下：假定水塘位于一块面积已知的矩形农田之中如图塘位于一块面积已知的矩形农田之中如图8 82 2所示随机地向这块农田扔石头使得它们所示随机地向这块农田扔石头使得它们都落在农田内被扔到农田中的石头可能溅上都落在农田内被扔到农田中的石头可能溅上了水，也可能没有溅上水，估计被了水，也可能没有溅上水，估计被“溅上水的溅上水的”石头量占总的石头量的百分比试想如何利用石头量占总的石头量的百分比试想如何利用这估计的百分比去近似计算该水塘面积？这估计的百分比去近似计算该水塘面积？49任意曲边梯形面积的近似计
35、算任意曲边梯形面积的近似计算50任意曲边梯形面积的近似计算任意曲边梯形面积的近似计算l结合图结合图8 82 2中的图形中的图形(1)(1)分析，只要已知各种参数及函数分析，只要已知各种参数及函数（a a，b b，H H，f(x)f(x)），有以下两种方法可近似计算水塘面），有以下两种方法可近似计算水塘面积积 l1 1随机投点法随机投点法 1 1）赋初值：试验次数）赋初值：试验次数n=0n=0，成功次数，成功次数m=0m=0；规定投点试验；规定投点试验的总次数的总次数N N；l2 2）随机选择）随机选择m m个数对个数对xi,yi,1ixi,yi,1im m，其中，其中yaya xixi b,b
36、,00yiyi H H，置，置 n nn nl l；l3 3）判断）判断nNnN，若是，转，若是，转4 4，否则停止计算；，否则停止计算；l4 4）判断条件）判断条件 (表示一块溅水的石头表示一块溅水的石头)是否成立，是否成立，若成立则置若成立则置m=m+1m=m+1，转，转2 2，否则转，否则转2 2；l5 5）计算水塘面积的近似值）计算水塘面积的近似值S=HS=H(b-a)(b-a)m/m/N N ()iiyf x51任意曲边梯形面积的近似计算任意曲边梯形面积的近似计算l2 2平均值估计法平均值估计法l1 1）产生）产生a,ba,b区间的均匀随机数区间的均匀随机数x x,i=1,2,i=1
37、,2,N,Nl2)2)计算计算f(xif(xi），），i=1,2i=1,2,N,Nl3 3）计算）计算l该方法的特点是估计函数该方法的特点是估计函数f(x)f(x)在在a,ba,b上的平上的平均值，面积近似等于该平均值乘以均值，面积近似等于该平均值乘以(b-a)(b-a)1()NiibaSf xN52例例9 库存问题库存问题l在物资的供应过程中在物资的供应过程中,由于到货和销售不可能做到由于到货和销售不可能做到同步同量同步同量,故总要保持一定的库存储备故总要保持一定的库存储备.如果库存如果库存过多过多,就会造成积压浪费以及保管费的上升就会造成积压浪费以及保管费的上升;如果如果库存过少库存过少,
38、会造成缺货会造成缺货.如何选择库存和订货策略如何选择库存和订货策略,就是一个需要研究的问题就是一个需要研究的问题.库存问题有多种类型库存问题有多种类型,一般都比较复杂一般都比较复杂,下面讨论一种简单的情形下面讨论一种简单的情形.l 某电动车行的仓库管理人员采取一种简单的订货某电动车行的仓库管理人员采取一种简单的订货策略策略,当库存降低到当库存降低到P P辆电动车时就向厂家辆电动车时就向厂家53例例9 库存问题库存问题l订货订货,每次订货每次订货Q Q辆辆,如果某一天的需求量超过了库如果某一天的需求量超过了库存量存量,商店就有销售损失和信誉损失商店就有销售损失和信誉损失,但如果库存但如果库存量过
39、多量过多,就会导致资金积压和保管费增加就会导致资金积压和保管费增加.若现在若现在有如下表所示的两种库存策略有如下表所示的两种库存策略,试比较选择一种策试比较选择一种策略以使总费用最少略以使总费用最少.l 表表 订货方案订货方案 方案方案 重新订货点重新订货点P/P/辆辆 重新订货量重新订货量Q/Q/辆辆 方案方案1 125 1501 125 150 方案方案2 150 2502 150 25054例例9 库存问题库存问题l这个问题的已知条件是这个问题的已知条件是:l(1)(1)从发出订货到收到货物需隔从发出订货到收到货物需隔3 3天天.l(2)(2)每辆电动车保管费为每辆电动车保管费为0.50
40、0.50元元/天天,每辆电动每辆电动车的缺货损失为车的缺货损失为1.601.60元元/天天,每次的订货费为每次的订货费为7575元元.l(3)(3)每天电动车需求量是每天电动车需求量是0 0到到9999之间均匀分布的之间均匀分布的随机数随机数.l(4)(4)原始库存为原始库存为110110辆辆,假设第一天没有发出订假设第一天没有发出订货货.55例例9 库存问题库存问题l分析分析:这一问题用解析法讨论比较麻烦这一问题用解析法讨论比较麻烦,但用计算但用计算机按天仿真仓库货物的变动情况却很方便机按天仿真仓库货物的变动情况却很方便.我们以我们以3030天为例天为例,依次对这两种方案进行仿真依次对这两种
41、方案进行仿真,最后比较最后比较各方案的总费用各方案的总费用,从而就可以作出决策从而就可以作出决策.l计算机仿真时的工作流程是早上到货计算机仿真时的工作流程是早上到货,全天销售全天销售,晚上订货晚上订货,输入一些常数和初始数据后输入一些常数和初始数据后,以一天为以一天为时间步长进行仿真时间步长进行仿真.首先检查这一天是否为预订到首先检查这一天是否为预订到货日期货日期,如果是如果是,则原有库存量加则原有库存量加Q,Q,并把预定到货并把预定到货量清为量清为0;0;如果不是如果不是,则库存量不变则库存量不变.接着用计算机接着用计算机中的随机函数仿真随机需求量中的随机函数仿真随机需求量,若库存量大于需求
42、若库存量大于需求量量,则新的库存量减去需求量则新的库存量减去需求量;反之反之,则新库存量变则新库存量变为为0,0,并且要在总费用上加缺货损失并且要在总费用上加缺货损失,然后检查实际然后检查实际56例例9 库存问题库存问题l库存量加上预定到货量是否小于重新订货点库存量加上预定到货量是否小于重新订货点P,P,如如果是果是,则需要重新订货则需要重新订货,这时就加一次订货费这时就加一次订货费.如此如此重复运行重复运行3030天天,即可得所需费用总值即可得所需费用总值.由此比较这由此比较这两种方案的总费用两种方案的总费用,即得最好方案即得最好方案.57lclearldays=30;lP=125,150;
43、lQ=150,250;lcost=0,0;larrivalinterval=2;lstoragefee=0.5;lossfee=1.6;bookfee=75;lstorage0=110;booknumber=0;arrivedate=0;lnr=rand(days,1);lfor i=1:2l storage(1)=storage0;l n=round(99*nr(1);l sale=n;l remain=storage(1)-n;l if remain=P(i);l booknumber=Q(i);l arrivedate=4;l orderfee=bookfee;l elsel order
44、fee=0;l endl storage(1)=remain;l cost(i)=cost(i)+remain*storagefee+orderfee;l for j=2:daysl dh=j;l if abs(dh-arrivedate)=nl sale=n;l remain=storage(j)-n;l shortagenumber=0;l else l sale=storage(j);l remain=0;l shortagenumber=n-storage(j);l endl storage(j)=remain;l if remain+booknumber=P(i);l booknum
45、ber=Q(i);l arrivedate=dh+arrivalinterval;l orderfee=bookfee;l elsel orderfee=0;l endl cost(i)=cost(i)+remain*storagefee+shortagenumber*lossfee+orderfee;l end;l mincost=min(cost);l endl costl mincost59例例9.10 赶火车过程仿真赶火车过程仿真l一列火车从一列火车从A A站经站经B B站开往站开往C C站站,某人每天赶往某人每天赶往B B站乘站乘这趟火车这趟火车.已知火车从已知火车从A A站到站到B
46、 B站运行时间为均值站运行时间为均值30min,30min,标准差为标准差为2min2min的正态随机变量的正态随机变量.火车大约在火车大约在下午下午1 1点离开点离开A A站站,离开时刻的频率分布见表离开时刻的频率分布见表1.1.这个这个人到达人到达B B站时的频率分布表见表站时的频率分布表见表2.2.用计算机仿真火用计算机仿真火车开出车开出,火车到达火车到达B B站站,这个人到达这个人到达B B站的情况站的情况,并给并给出能赶上火车的仿真结果出能赶上火车的仿真结果.表表1 1 离开时刻的频率分布离开时刻的频率分布出发时刻出发时刻T 1:00 1:05 1:10T 1:00 1:05 1:1
47、0频率频率 0.7 0.2 0.10.7 0.2 0.160例例9.10 赶火车过程仿真赶火车过程仿真l表表2 2 人到达人到达B B站时的频率分布站时的频率分布l 到达时刻到达时刻T 1:28 1:30 1:32 1:34T 1:28 1:30 1:32 1:34l 频率频率:0.3 0.4 0.2 0.1:0.3 0.4 0.2 0.1l分析分析:引入以下变量引入以下变量:T:T1 1为火车从为火车从A A站开出的时刻站开出的时刻;T;T2 2为火车为火车从从A A站运行到站运行到B B站所需要的时间站所需要的时间;T;T3 3为此人到达为此人到达B B站的时刻站的时刻.则有表则有表3.3
48、.l 表表3 T3 T1,1,T T2 2的频率分布的频率分布lT T1,1,/min 0 5 10 T/min 0 5 10 T2 2/min 28 30 32 34/min 28 30 32 34l p 0.7 0.2 0.1 0.3 0.4 0.2 0.1 p 0.7 0.2 0.1 0.3 0.4 0.2 0.161例例9.10 赶火车过程仿真赶火车过程仿真l显然显然,这位旅客要赶上火车的条件是这位旅客要赶上火车的条件是T T3 3 T T1 1+T+T2 2,可以通过可以通过计算机模拟出这三个时间计算机模拟出这三个时间,再检验是否满足再检验是否满足T T3 3 T T1 1+T+T2
49、 2.如如果满足果满足,即能够赶上火车即能够赶上火车,若不然若不然,则不能赶上火车则不能赶上火车.火车运行时间的仿真程序火车运行时间的仿真程序x=randn(10000,1);for i=1:10000 y(i)=30+2*x(i);end62开车时间的仿真程序开车时间的仿真程序ls1=0;s2=0;s3=0;l x=rand(10000,1);l for i=1:10000l if x(i)0.9l s3=s3+1;l endl endls1/10000,1-s1/10000-s3/10000,s3/1000063人到达时刻的仿真程序人到达时刻的仿真程序ls1=0;s2=0;s3=0;s4=
50、0;l x=rand(10000,1);l for i=1:10000l if x(i)0.3l s1=s1+1;l elseif x(i)0.7l s2=s2+1;l elsel if x(i)0.9l s3=s3+1;l else l s4=s4+1;l endl endlendls1/10000,s2/10000,s3/10000,s4/1000064赶上火车的仿真程序赶上火车的仿真程序ls=0;lx1=rand(10000,1);lx2=rand(10000,1);lx3=randn(10000,1);lfor i=1:10000lif x1(i)0.7lT1=0;lelseif x1
