1、数学建模培训,计算机模拟,模拟的概念,模拟就是利用物理的、数学的模型来类比、模仿现实系统及其演变过程，以寻求过程规律的一种方法,模拟的基本思想是建立一个试验的模型，这个模型包含所研究系统的主要特点通过对这个实验模型的运行，获得所要研究系统的必要信息.,对实际系统及其过程用功能相似的实物系统去模仿例如，军事演习、船艇实验、沙盘作业等,物理模拟通常花费较大、周期较长，且在物理模型上改变系统结构和系数都较困难而且，许多系统无法进行物理模拟，如社会经济系统、生态系统等,模拟的方法,物理模拟,在一定的假设条件下，运用数学运算模拟系统的运行，称为数学模拟现代的数学模拟都是在计算机上进行的，称为计算机模拟,
2、数学模拟,计算机模拟可以反复进行，改变系统的结构和系数都比较容易,在实际问题中，面对一些带随机因素的复杂系统，用分析方法建模常常需要作许多简化假设，与面临的实际问题可能相差甚远，以致解答根本无法应用这时，计算机模拟几乎成为唯一的选择计算机模拟与数学建模的关系如图所示：,案例一 圆周率的近似计算 蒙特卡罗法,在正方形中随机撒一把豆子，通过考察落在其内切圆内黄豆的数目，用随机模拟的方法计算圆周率P的近似值（如图）,具体方案如下： 利用计算机产生两组0，1上的均匀随机数，x=RAND，y=RAND； 统计试验总次数N和落在内切圆内的点数N1（满足条件x2+y21的点（x，y）的个数）； 计算频率N1
3、 /N ，即为点落在圆内的概率的近似值； 设圆的面积为S，由几何概率公式得点落在内切圆部分的概率为P=S/ 4 S/ 4 = N1 /N S4N1/N ，即为圆的面积的近似值又S圆=r2=，=S4N1 /N ，即为圆周率的近似值,模拟过程流程图,cs=0 n=500 for i=1:na=rand(1,2);if a(1)2+a(2)2=1cs=cs+1end end 4*cs/n,Matlab程序：,从计算结果看：这种数据模拟算法，在实验次数较少时，算的近似值距离真实值误差较大，但这种数据模拟方法简单易行。在精度要求不高的情况下，这种取随机数进行数据模拟的方法还是有一定的实用价值，通过简单编
4、程我们可以模拟出许多数学现象。,依次取n=500,1000,100000取得的近似值分别为：3.0960；3.0680；3.1405,案例二 随机模拟抛硬币问题,硬币有正反两面，抛掷的均匀的硬币，正面和反面朝上的机会都是二分之一，我们利用matlab函数rand产生在（0,1）区间内均与分布的伪随机数的功能，设计算法，编写程序，模拟抛掷的均与的硬币这种随机行为。,算法：设x是（0，1）区间分布的随机数，定义x=0.5表示正面朝上，x=0.5，k=k+1;否则k不变；第5步 计算频率p=k/n,抛硬币的matlab程序,k=0 n=10 for i=1:nx=rand;if x=0.5k=k+1
5、;end end k p=k/n,案例三 随机模拟掷骰子问题,骰子呈正方体，有6面，点数分别是1、2、3、4、5、6，如果骰子是质地均匀的，那么投掷骰子掷出每一个点数的机会都是六分之一。我们要设计算法，编写程序，模拟投掷骰子的这种随机行为。,算法：设x是（0，1）区间分布的随机数，定义0x1/6表示掷出1点， 1/6x1/3表示掷出2点, 5/6x1表示掷出6点.输入 投骰子的总次数n输出 投出每一个点数的频率p（有六个分量）第1步 初始化计数器（有六个分量）；第2步 对i=1,2,.,n,循环进行第3,4步；第3步 x=rand;第4步 如果0x1/6 ，k(1)=k(1)+1;如果1/6x
6、1/3，k(2)=k(2)+1;.如果5/6x1，k(6)=k(6)+1;第5步 计算频率p=k/n,n=100; k=zeros(1,6); for i=1:nx=rand;if x1/6k(1)=k(1)+1;elseif x2/6k(2)=k(2)+1;elseif x3/6k(3)=k(3)+1;elseif x4/6k(4)=k(4)+1;elseif x5/6k(5)=k(5)+1;elsek(6)=k(6)+1;end end k p=k/n,投掷骰子的matlab程序,产生模拟随机数的计算机命令,在MATLAB软件中，可以直接产生满足各种分布的随机数，命令如下：,2产生mn阶，
7、均匀分布的随机数矩阵：rand (m, n)产生一个，均匀分布的随机数：rand,1产生mn阶a，b上均匀分布U（a，b）的随机数矩阵：unifrnd (a,b,m, n)产生一个a，b均匀分布的随机数：unifrnd (a,b),当只知道一个随机变量取值在（a，b）内，但不知道（也没理由假设）它在何处取值的概率大，在何处取值的概率小，就只好用U（a，b）来模拟它,例 1的计算机模拟,To MATLAB（rnd）,当研究对象视为大量相互独立的随机变量之和，且其中每一种变量对总和的影响都很小时，可以认为该对象服从正态分布,机械加工得到的零件尺寸的偏差、射击命中点与目标的偏差、各种测量误差、人的身
8、高、体重等，都可近似看成服从正态分布,若连续型随机变量X的概率密度函数为其中 0为常数，则称X服从参数为 的指数分布,指数分布的期望值为,指数分布可以用来表示独立随机事件发生的时间间隔，比如旅客进机场的时间间隔，排队服务系统中顾客到达率为常数时的到达间隔、故障率为常数时零件的寿命都服从指数分布,指数分布在排队论、可靠性分析中有广泛应用,注意：MATLAB中，产生参数为 的指数分布的命令为exprnd( ),例 顾客到达某商店的间隔时间服从参数为0.1的指数分布,指数分布的均值为1/0.1=10指两个顾客到达商店的平均间隔时间是10个单位时间.即平均10个单位时间到达1个顾客. 顾客到达的间隔时
9、间可用exprnd(10)模拟,设离散型随机变量X的所有可能取值为0,1,2,且取各个值的概率为 其中 0为常数，则称X服从参数为 的泊松分布,泊松分布在排队系统、产品检验、天文、物理等领域有广泛应用,泊松分布的期望值为,案例四 敌坦克分队对我方阵地实施突袭，其到达规律服从泊松分布，平均每分钟到达辆（1）模拟敌坦克在分钟内到达目标区的数量，以及在第、分钟内各到达几辆坦克（2）模拟在3分钟内每辆敌坦克的到达时刻,（1）用poissrnd(4)进行模拟,To MATLAB（poiss）,（2）坦克到达的间隔时间应服从参数为4的负指数分布，用exprnd（1/4）模拟,To MATLAB（time）
10、,案例四 企业库存系统的计算机模拟方法,在销售部门、工厂等领域中都存在库存问题, 库存太多造成浪费以及资金积压,库存少了不能满足需求也造成损失. 部门的工作人员需决定何时进货,进多少,使得所花费的平均费用最少,而收益最大,这就是库存问题.某企业生产易变质的产品. 当天生产的产品必须售出,否则就会变质. 该产品单位成本为2.5 元,单位产品售价为5 元. 企业为避免存货过多而造成损失,拟从以下2 种库存方案中选出一个较优的方案:方案甲:按前1 天的销售量作为当天的库存量;方案乙:按前2 天的平均销售量作为当天的库存量.,【问题的提出】,【模型建立】,假定市场对该产品的每天需求量是一个随机变量,但
11、从以往的统计分析得知它服从正态分布: N(135 , 22. 4) .,第一：获得市场对该产品的需求量的数据;,第二：计算出按照2 种不同方案经T 天后企业的利润值；,第三：比较大小,从中选出一个更优的方案.,计算机模拟的基本思路:,D :每天需求量; Q1 :方案甲当天的库存量; Q2 :方案甲当天的库存量; S1 :方案甲前1 天的销售量; S21 :方案乙前1 天的销售量; S22 :方案乙前2 天的销售量; S3 :方案甲当天实际销售量; S4 :方案乙当天实际销售量; L1 :方案甲当天的利润; L2 :方案乙当天的利润; TL1 :方案甲累计总利润; TL2 :方案甲累计总利润;
12、T :预定模拟天数.,引入下列符号:,模拟过程的流程图,【模型的求解】,利用Matlab 编程来实现这一过程, 建立如下M文件:,function TL1 , TL2 = kucun( T , S1 , S21 , S22) TL1 = 0; TL2 = 0 ; k = 1; while k TQ1 = S1 ; Q2 = ( S21 + S22) / 2 ;D = normrnd(135 ,22. 4) ;if D Q1S3 = Q1 ;elseS3 = D;endif D Q2S4 = Q2 ;elseS4 = D;endL1 = 5 *S3 - 2. 5 * Q1 ; L2 = 5 *S
13、4 - 2. 5 *Q2 ;TL1 = TL1 + L1 ; TL2 = TL2 + L2 ;k = k + 1 ; end S1 = S3 ; S22 = S21 ; S21 = S4 ;,案例五连续系统模拟实例: 追逐问题,状态随时间连续变化的系统称为连续系统对连续系统的计算机模拟只能是近似的，只要这种近似达到一定的精度，也就可以满足要求,追逐问题: 如图,正方形ABCD的4个顶点各有1人.在某一时刻,4人同时出发以匀速v=1m/s按顺时针方向追逐下一人,如果他们始终保持对准目标,则最终按螺旋状曲线于中心点O.试求出这种情况下每个人的行进轨迹.,1. 建立平面直角坐标系: A(x1,y1)
14、, B(x2,y2), C(x3,y3), D(x4,y4).,2. 取时间间隔为t,计算每一点在各个时刻的坐标.,4. 对每一个点，连接它在各时刻的位置,即得所求运动轨迹.,求解过程:,To MATLAB(chase),v=1; dt=0.05; x=0 0 10 10; y=0 10 10 0;for i=1:4plot(x(i),y(i),.),hold on endd=20; while(d0.1)x(5)=x(1);y(5)=y(1);for i=1:4d=sqrt(x(i+1)-x(i)2+(y(i+1)-y(i)2);x(i)=x(i)+v*dt*(x(i+1)-x(i)/d;y
15、(i)=y(i)+v*dt*(y(i+1)-y(i)/d;plot(x(i),y(i),.),hold onendend,计算程序:,To MATLAB(chase),案例六 在我方某前沿防守地域，敌人以一个炮排（含两门火炮）为单位对我方进行干扰和破坏为躲避我方打击，敌方对其阵地进行了伪装并经常变换射击地点,经过长期观察发现，我方指挥所对敌方目标的指示有50是准确的，而我方火力单位，在指示正确时，有1/3的射击效果能毁伤敌人一门火炮，有1/6的射击效果能全部消灭敌人,现在希望能用某种方式把我方将要对敌人实施的20次打击结果显现出来，确定有效射击的比率及毁伤敌方火炮的平均值,分析：这是一个概率问
16、题，可以通过理论计算得到相应的概率和期望值.但这样只能给出作战行动的最终静态结果，而显示不出作战行动的动态过程.,为了能显示我方20次射击的过程，现采用模拟的方式,需要模拟出以下两件事：,1. 问题分析,2 当指示正确时，我方火力单位的射击结果情况,1 观察所对目标的指示正确与否,模拟试验有两种结果，每种结果出现的概率都是1/2,因此，可用投掷1枚硬币的方式予以确定，当硬币出现正面时为指示正确，反之为不正确,模拟试验有三种结果：毁伤1门火炮的可能性为1/3(即2/6)，毁伤两门的可能性为1/6，没能毁伤敌火炮的可能性为1/2(即3/6),这时可用投掷骰子的方法来确定： 如果出现的是、点：则认为
17、没能击中敌人； 如果出现的是、点：则认为毁伤敌人一门火炮； 若出现的是点：则认为毁伤敌人两门火炮,2. 符号假设,i：要模拟的打击次数； k1：没击中敌人火炮的射击总数； k2：击中敌人一门火炮的射击总数；k3：击中敌人两门火炮的射击总数 E：有效射击比率； E1：20次射击平均每次毁伤敌人的火炮数,3. 模拟框图,4. 模拟结果,6. 结果比较,虽然模拟结果与理论计算不完全一致，但它却能更加真实地表达实际战斗动态过程,用蒙特卡罗方法进行计算机模拟的步骤：,1 设计一个逻辑框图，即模拟模型这个框图要正确反映系统各部分运行时的逻辑关系 2 模拟随机现象可通过具有各种概率分布的模拟随机数来模拟随机
18、现象,clear k1=0;k2=0;k3=0; for i=1:20R1=rand if R15/6k3=k3+1;elsek2=k2+1;endend elsek1=k1+1; end end i E=(k2+k3)/20 E1=(0*k1+1*k2+2*k3)/20,离散系统模拟实例: 排队问题,排队论主要研究随机服务系统的工作过程,在排队系统中，服务对象的到达时间和服务时间都是随机的排队论通过对每个个别的随机服务现象的统计研究，找出反映这些随机现象平均特性的规律，从而为设计新的服务系统和改进现有服务系统的工作提供依据,对于排队服务系统, 顾客常常注意排队的人是否太多, 等候的时间是否长
19、, 而服务员则关心他空闲的时间是否太短. 于是人们常用排队的长度、等待的时间及服务利用率等指标来衡量系统的性能.,1 系统的假设： （1） 顾客源是无穷的；（2） 排队的长度没有限制； （ 3） 到达系统的顾客按先后顺序依次进入服务， 即“先到先服务”,单服务员的排队模型：在某商店有一个售货员，顾客陆续来到，售货员逐个地接待顾客当到来的顾客较多时，一部分顾客便须排队等待，被接待后的顾客便离开商店设： 1顾客到来间隔时间服从参数为0.1的指数分布 对顾客的服务时间服从，上的均匀分布 排队按先到先服务规则，队长无限制,假定一个工作日为8小时，时间以分钟为单位 1模拟一个工作日内完成服务的个数及顾客
20、平均等待时间t 2模拟100个工作日，求出平均每日完成服务的个数及每日顾客的平均等待时间,2 符号说明 w：总等待时间；ci：第i个顾客的到达时刻；bi：第i个顾客开始服务时刻； ei：第i个顾客服务结束时刻xi：第i-1个顾客与第i个顾客到达之间的时间间隔yi：对第i个顾客的服务时间,c1,b1,c3,c4,c5,c2,e1,b2,e2,b3,e3,b4,e4,b5,ci=ci-1+ xi ei=bi+yi bi=max(ci,ei-1),t,1模拟1日To MATLAB(simu1),2模拟100日 To MATLAB(simu2),clear i=2; w=0; x(i)=exprnd(
21、10); c(i)=x(i); b(i)=x(i);while b(i)=480y(i)=unifrnd(4,15); e(i)=b(i)+y(i); w=w+b(i)-c(i); i=i+1; x(i)=exprnd(10); c(i)=c(i-1)+x(i); b(i)=max(c(i),e(i-1); endi=i-2; t=w/i m=i,作 业,1用matlab编一个福利彩票双色球电脑选号的程序并要求前六个数输出结果按升序排序,3. 某设备上安装有4只型号规格完全相同的电子管，已知电子管寿命服从100200h之间的均匀分布电子管损坏时有两种维修方案，一是每次更换损坏的那只；二是当其中1只损坏时4只同时更换已知更换时间为换1只时需1h，4只同时换为2h更换时机器因停止运转每小时的损失为20元，又每只电子管价格10元，试用模拟方法确定哪一个方案经济合理？,4. 设位于坐标原点的甲舰向位于x轴上点A(1, 0)处的乙舰发射导弹，导弹头始终对准乙舰.如果乙舰以最大的速度(常数)沿平行于y轴的直线行驶，导弹的速度为5. 模拟导弹运行的轨迹.乙舰行驶多远时，导弹将它击中？,
