数学建模之计算机仿真

数学建模之计算机仿真第1页，课件共61页，创作于2023年2月概述

计算机科学技术的迅猛发展，给许多学科带来了巨大的影响．计算机不但使问题的求解变得更加方便、快捷和精确，而且使得解决实际问题的领域更加广泛．计算机适合于解决那些规模大、难以解析化以及不确定的数学模型．例如对于一些带随机因素的复杂系统，用分析方法建模常常需要作许多简化假设，与面临的实际问题可能相差甚远，以致解答根本无法应用，这时仿真几乎成为人们的唯一选择．在历届的美国和中国大学生的数学建模竞赛（MCM）中，学生们经常用到计算机仿真方法去求解、检验等．计算机仿真(computersimulation)是建模过程中较为重要的一类方法．第2页，课件共61页，创作于2023年2月计算机仿真的基本概念

计算机仿真，是根据已知的信息和知识（如数学、物理规律），利用计算机模拟现实情况和系统的演变过程，发现新的知识或规律，从而解决问题的一种方法第3页，课件共61页，创作于2023年2月计算机仿真的基本概念

独立于理论研究与实验研究的认识世界的第三中方法第4页，课件共61页，创作于2023年2月计算机仿真的基本概念计算机仿真的特点

代价小、时间短、可重复、参数设置灵活第5页，课件共61页，创作于2023年2月计算机仿真可以解决以下5类问题:(1)难以用数学公式表示的系统,或者没有建立和求解的有效方法.(2)虽然可以用解析的方法解决问题,但数学的分析与计算过于复杂,这时计算机仿真可能提供简单可行的求解方法.(3)希望能在较短的时间内观察到系统发展的全过程,以估计某些参数对系统行为的影响.(4)难以在实际环境中进行试验和观察时,计算机仿真是唯一可行的方法,如太空飞行的研究.(5)需要对系统或过程进行长期运行比较,从大量方案中寻找最优方案.第6页，课件共61页，创作于2023年2月

计算机仿真的分类计算机仿真在计算机中运行实现,不怕破坏,易修改,可重用,安全经济,不受外界条件和场地空间的限制.仿真分为静态仿真(staticsimulation)和动态仿真(dynamicsimulation).数值积分中的蒙特卡洛方法（统计模拟方法）是典型的静态仿真．动态仿真又分为连续系统仿真和离散系统仿真．连续系统是指状态变量随着时间连续变化的系统，例如传染病的检测与预报系统.离散系统是指系统状态变量只在有限的时间点或可数的时间点上发生变化的系统,例如排队系统.第7页，课件共61页，创作于2023年2月

仿真系统，必须设置一个仿真时钟(simulateclock)，它能将时间从一个时刻向另一个时刻进行推进，并且能随时反映系统时间的当前值．其中，模拟时间推进方式有两种:时间步长法(均匀间隔时间推进法,连续系统常用)和事件步长法(下次事件推进法,离散系统常用)．第8页，课件共61页，创作于2023年2月主要内容一:准备知识:随机数的产生二:随机变量的模拟三:连续系统的模拟-时间步长法四:离散系统的模拟-事件步长法●五：蒙特卡洛方法第9页，课件共61页，创作于2023年2月一:准备知识:随机数的产生由于仿真研究的实际系统要受到多种随机因素的作用和影响,在仿真过程中必须处理大量的随机因素.要解决此问题的前提是确定随机变量的类型和选择合适的随机数产生的方法.对随机现象进行模拟,实质是要给出随机变量的模拟,也就是说要利用计算机随机产生一系列数值,使它们服从一定的概率分布,称这些数值为随机数.最基本,最常用的是(0,1)区间内均匀分布的随机数.其他分布的随机数均可利用它来产生.第10页，课件共61页，创作于2023年2月1:产生模拟随机数的计算机命令在MATLAB中,可以直接产生满足各种分布的随机数,命令如下:常见的分布函数MATLAB语句均匀分布U[0,1]R=rand(m,n)均匀分布U[a,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,n)泊松分布P(λ)R=poissrnd(λ,m,n)以上语句均产生m×n的矩阵．第11页，课件共61页，创作于2023年2月2:案例分析例1:unifrnd(2,3)unifrnd(1,32,1,4)normrnd(1,2)normrnd(1,2,2,3)rand(2,3)randn(2,3)第12页，课件共61页，创作于2023年2月2:案例分析ans=2.8132ans=1.30575.30567.28577.1604ans=0.2527ans=2.74290.02192.77592.27560.0992-0.9560ans=0.60380.19880.74680.27220.01530.4451ans=-0.0945-1.3089-0.2440-0.21410.8248-0.1778>>第13页，课件共61页，创作于2023年2月2:案例分析例2:敌空战部队对我方港口进行空袭,其到达规律服从泊松分布,平均每分钟到达4架飞机.(1)模拟敌机在3分钟内到达目标区域的数量,以及在第1,2,3分钟内各到达几架飞机;(2)模拟在3分钟内每架飞机的到达时刻.分析:(1)n1=poissrnd(4),n2=poissrnd(4),n3=poissrnd(4),n=n1+n2+n3(2)由排队论知识,敌机到达规律服从泊松分布等价于敌机到达港口的间隔时间服从参数为1/4的指数分布,故可由指数分布模拟每架飞机的到达时刻.第14页，课件共61页，创作于2023年2月2:案例分析cleart=0;j=0;%到达的飞机数whilet<3j=j+1t=t+exprnd(1/4)end第15页，课件共61页，创作于2023年2月二:随机变量的模拟利用均匀分布的随机数可以产生具有任意分布的随机变量的样本,从而可以对随机变量的取值情况进行模拟.1连续型随机变量的模拟具有给定分布的连续型随机变量可以利用在区间(0,1)上均匀分布的随机数来模拟,最常用的方法是逆变换法.结论:若随机变量Y有连续的分布函数F(y),

则Z与Y有相同的分布.第16页，课件共61页，创作于2023年2月1连续型随机变量的模拟若已知Y的概率密度为如果给定区间(0,1)上均匀分布的随机数,则具有给定分布Y的随机数可由方程解出.例:模拟服从参数为的指数分布时,由可得第17页，课件共61页，创作于2023年2月2离散型随机变量的模拟设随机变量X的分布律为:有相同的发生的概率.因此我们可以用随机变量R落在小区间内的情况来模拟离散的随机变量X的取值情况.第18页，课件共61页，创作于2023年2月2离散型随机变量的模拟例3:随机变量表示每分钟到达银行柜台的顾客数.X的分布列见下表,试模拟10分钟内顾客到达柜台的情况.

表110分钟内顾客到达柜台的情况

Xk012pk0.40.30.3分析:因为每分钟到达柜台的人数是随机的,所以可用计算机随机生成一组(0,1)的数据,由X的概率分布情况,可认为随机数在(0,0.4)范围内时没有顾客光顾,在[0.4,0.7)时,有一个顾客光顾,在[0.7,1)时,有两个顾客光顾.

从而有MATLAB程序:第19页，课件共61页，创作于2023年2月2离散型随机变量的模拟r=rand(1,10);fori=1:10;ifr(i)<0.4n(i)=0;elseif0.4<=r(i)&r(i)<0.7n(i)=1;elsen(i)=2;end;endrn第20页，课件共61页，创作于2023年2月三:连续系统的模拟-时间步长法对连续系统的计算机模拟是近似地获取系统状态在一些离散时刻点上的数值．在一定假设条件下，利用数学运算模拟系统的运行过程．连续系统模型一般是微分方程，它在数值模拟中最基本的算法是数值积分算法．例如有一系统可用微分方程来描述:已知输出量y的初始条件，现在要求出输出量y随时间变化的过程y(t)。第21页，课件共61页，创作于2023年2月1理论介绍最直观的想法是：首先将时间离散化，令，称为第k步的计算步距（一般是等间距的），然后按以下算法计算状态变量在各时刻上的近似值:

其中初始点按照这种作法即可求出整个的曲线．这种最简单的数值积分算法称为欧拉法．除此之外，还有其他一些算法．第22页，课件共61页，创作于2023年2月1理论介绍因此，连续系统模拟方法是：首先确定系统的连续状态变量，然后将它在时间上进行离散化处理，并由此模拟系统的运行状态．模拟过程分为许多相等的时间间隔,时间步长的长度可以根据实际问题分别取为秒,分,小时,天等.仿真时钟按时间步长等距推进,每次推进都要扫描系统中所有活动,按预定计划和目标进行分析,计算,记录系统状态的变化,直到满足某个终止条件为止.第23页，课件共61页，创作于2023年2月

例:追逐问题1．问题提出假设正方形ABCD的四个顶点处各站一人．在某一时刻，四人同时以匀速v沿顺时针方向追逐下一个人，并且在任意时刻他们始终保持追逐的方向是对准追逐目标，例如，A追逐B，任意时刻A始终向着B追．可以证明四人的运动轨迹将按螺旋曲线状汇合于中心O．怎样证明呢？有两种证明方法．一是分别求出四人的运动轨迹曲线解析式，求证四条曲线在某时刻相交于一点．另一方法则是用计算机模拟将四人的运动轨迹直观地表示在图形上.2应用举例第24页，课件共61页，创作于2023年2月2．建立模型及模拟方法模拟步骤：1）建立平面直角坐标系．

2）以时间间隔Δt进行采样，在每一时t计算每个人在下一时t+Δt时的坐标．

3）不妨设甲的追逐对象是乙，在时间t时，甲第25页，课件共61页，创作于2023年2月2．建立模型及模拟方法4)选取足够小的，模拟到时为止5）连接四人在各时刻的位置，就得到所求的轨迹．根据以上模拟步骤，编出MATLAB程序(simu2.m)如下：第26页，课件共61页，创作于2023年2月%取v=1,t=12,A,B,C,D点的坐标分另为（0，10），（10，10），（10，0），(0,0)v=1;dt=0.05;d=20;x=[000101010100];x(9)=x(1);x(10)=x(2);holdaxis('equal')axis([010010]);3MATLAB实现第27页，课件共61页，创作于2023年2月fork=1:2:7plot(x(k),x(k+1),'.')endwhile(d>0.1)fori=1:2:7d=sqrt((x(i)-x(i+1))^2+(x(i+1)-x(i+3))^2);x(i)=x(i)+v*dt*(x(i+2)-x(i))/d;x(i+1)=x(i+1)+v*dt*(x(i+3)-x(i+1))/d;plot(x(i),x(i+1),'.')endx(9)=x(1);x(10)=x(2);endhold3MATLAB实现第28页，课件共61页，创作于2023年2月3MATLAB实现第29页，课件共61页，创作于2023年2月例9.6水池含盐量问题某水池有2000m3水,其中含盐2kg,以6m3/min的速率向水池内注入含盐为0.5kg/m3的盐水,同时又以4m3/min的速率从水池流出搅拌均匀的盐水.试用计算机仿真该水池内盐水的变化过程,并每隔10min计算水池中水的体积,含盐量,含盐率.欲使池中盐水含盐率达到0.2kg/m3,需经过多长时间?分析:这是一个连续系统,首先要将系统离散化,在一些离散点上进行考察,这些离散点的间隔就是时间步长.可取步长为1min,即隔1min考察一次系统的状态,并相应地记录和分析.在注入和流出的作用下,池中水的体积与含盐量,含盐率均随时间变化,初始时刻含盐率为0.001kg/m3,以后每分钟注入含盐率为0.5kg/m3的水6m3,流出混合均匀的盐水为4m3,当池中水的含盐率达到0.2kg/m3时,仿真过程结束.第30页，课件共61页，创作于2023年2月例9.6水池含盐量问题记T时刻的体积为w(m3)，水的含盐量为s(kg)，水的含盐率为r=s/w(kg/m3)，每隔1min池水的动态变化过程如下：每分钟水的体积增加6-4=2(m3)；每分钟向池内注入盐6×0.5=3(kg)；每分钟向池外流出盐4×r(kg)；每分钟池内增加盐3-4×r(kg).本例还可以用微分方程建立数学模型，并求出它的解析解，这个解析解就是问题的精确解，有兴趣的读者可以按照这个思路求出该问题的精确解，考察相应时刻精确解与仿真解的差异，还可以进一步调整仿真过程的时间步长，通过与精确解的比较来研究时间步长的大小对仿真度的影响。第31页，课件共61页，创作于2023年2月MATLAB实现clearh=1;%时间步长为1s0=2;%初始含盐2kgw0=2000;%初始水池有水2000m3r0=s0/w0;%初始浓度s(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*h^2+h^3)/(1000+h)^2;fori=2:200t(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(i))^2;m=floor(i/10);第32页，课件共61页，创作于2023年2月MATLAB实现ifi/10-m<0.1tm(m)=m;wm(m)=w(i);sm(m)=s(i);rm(m)=r(i);endifr(i)>0.2%若第i步后的盐水浓度大于0.2t02=i*h;r02=r(i);breakendend[t02,r02][10*tm',sm',rm']%'表示逆subplot(1,2,1),plot(t,s,'blue');holdonsubplot(1,2,2),plot(t,y,'red');第33页，课件共61页，创作于2023年2月

四:离散系统的模拟-事件步长法

离散系统(discretesystem)是指系统状态只在有限的时间点或可数的时间点上有随机事件驱动的系统．例如排队系统（queuesystem），显然状态量的变化只是在离散的随机时间点上发生．假设离散系统状态的变化是在一个时间点上瞬间完成的．常用的是事件步长法（下次事件推进法）．其过程是：置模拟时钟的初值为0，跳到第一个事件发生的时刻，计算系统的状态，产生未来事件并加入到队列中去，跳到下一事件，计算系统状态，……，重复这一过程直到满足某个终止条件为止．第34页，课件共61页，创作于2023年2月例9.7收款台前的排队过程假设:

(1)顾客到达收款台是随机的,平均时间间隔为0.5min,即间隔时间服从lamda=2的指数分布.(2)对不同的顾客收款和装袋的时间服从正态分布N(1,1/3).试模拟20位顾客到收款台前的排队情况,我们关心的问题是每个顾客的平均等待时间,队长及服务员的工作效率第35页，课件共61页，创作于2023年2月例9.7收款台前的排队过程分析:单服务台结构的排队系统有两类原发事件即到来和离去,顾客到来的后继事件是顾客接受服务,顾客离去的后继事件是服务台寻找服务,这四类事件各自的子程序框图如图1所示.假设:t(i)为第i位顾客到达时刻;t2(i)为第i位顾客受到服务的时间(随机变量);T(i)为第i位顾客离去时刻.将第i位顾客到达作为件事发生:t(i+1)-t(i)=r(i)(随机变量);平衡关系:当t(i+1)≥T(i)时,T(i+1)=t(i+1)+t2(i+1);否则,T(i+1)=T(i)+t2(i+1).第36页，课件共61页，创作于2023年2月图1到来事件子程序

系统人数+1产生下一个顾客到来时刻调用接受服务事件的程序第37页，课件共61页，创作于2023年2月图2离去事件子程序

系统人数-1置服务台”闲”已服务人数+1调用寻找服务事件子程序第38页，课件共61页，创作于2023年2月图3:接受服务事件子程序

服务台空置服务台”忙”产生服务结束时刻登记到事件表排队人数+1否是第39页，课件共61页，创作于2023年2月图4:寻找服务事件子程序

排队人数≥1置服务台”忙”产生服务结束时刻排队人数-1排队人数+1否是第40页，课件共61页，创作于2023年2月五：蒙特卡洛方法在用传统方法难以解决的问题中，有很大一部分可以用概率模型进行描述．由于这类模型含有不确定的随机因素，分析起来通常比确定性的模型困难．有的模型难以作定量分析，得不到解析的结果，或者是虽有解析结果，但计算代价太大以至不能使用．在这种情况下，可以考虑采用MonteCarlo方法。MonteCarlo方法是计算机模拟的基础，它的名字来源于世界著名的赌城——摩纳哥的蒙特卡洛，其历史起源于1777年法国科学家蒲丰提出的一种计算圆周π的方法——随机投针法，即著名的蒲丰投针问题。第41页，课件共61页，创作于2023年2月五：蒙特卡洛方法MonteCarlo方法的基本思想是首先建立一个概率模型，使所求问题的解正好是该模型的参数或其他有关的特征量．然后通过模拟一统计试验，即多次随机抽样试验（确定m和n），统计出某事件发生的百分比．只要试验次数很大，该百分比便近似于事件发生的概率．这实际上就是概率的统计定义．利用建立的概率模型，求出要估计的参数．蒙特卡洛方法属于试验数学的一个分支．第42页，课件共61页，创作于2023年2月五：蒙特卡洛方法蒙特卡洛方法适用范围很广泛，它既能求解确定性的问题，也能求解随机性的问题以及科学研究中的理论问题．例如利用蒙特卡洛方法可以近似地计算定积分，即产生数值积分问题．第43页，课件共61页，创作于2023年2月蒙特卡洛法求圆周率clearn=50000X=rand(n,1);Y=rand(n,1);k=0;fori=1:n;ifX(i)^2+Y(i)^2<=1k=k+1;endend4*k/n第44页，课件共61页，创作于2023年2月普丰投针求圆周率的M程序functionpi_value=pinpi(k,a,l)%k投针次数%a线距%l针长（必须小于等于d）%pi_value返回pi值m=0;ifl>a;fprintf('error：针长必须小于等于%d\n',a);fprintf('请重新调用函数pinpi(k,d,l)\n');pi_value=0;elsefori=1:kifa*rand(1)<=l*sin(pi*rand(1))m=m+1;endendp=m/k;pi_value=2*l/(a*p);fprintf('投针法求得pi=%d\n',pi_value);end第45页，课件共61页，创作于2023年2月formatlongg>>pinpi(100000,4,3)投针法求得pi=3.143797e+000ans=3.14379728795087第46页，课件共61页，创作于2023年2月任意曲边梯形面积的近似计算

一个古老的问题：用一堆石头测量一个水塘的面积．应该怎样做呢？测量方法如下：假定水塘位于一块面积已知的矩形农田之中．如图8．2所示．随机地向这块农田扔石头使得它们都落在农田内．被扔到农田中的石头可能溅上了水，也可能没有溅上水，估计被“溅上水的”石头量占总的石头量的百分比．试想如何利用这估计的百分比去近似计算该水塘面积？

第47页，课件共61页，创作于2023年2月任意曲边梯形面积的近似计算

第48页，课件共61页，创作于2023年2月任意曲边梯形面积的近似计算

结合图8．2中的图形(1)分析，只要已知各种参数及函数（a，b，H，f(x)），有以下两种方法可近似计算水塘面积．1．随机投点法1）赋初值：试验次数n=0，成功次数m=0；规定投点试验的总次数N；2）随机选择m个数对xi,yi,1<i<m，其中ya<xi<b,

0<yi<H，置n＝n＋l；3）判断n≤N，若是，转4，否则停止计算；4）判断条件(表示一块溅水的石头)是否成立，若成立则置m=m+1，转2，否则转2；5）计算水塘面积的近似值S=H×(b-a)×m/N．第49页，课件共61页，创作于2023年2月任意曲边梯形面积的近似计算

2．平均值估计法1）产生[a,b]区间的均匀随机数x,i=1,2…,N2)计算f(xi），i=1,2…,N3）计算该方法的特点是估计函数f(x)在[a,b]上的平均值，面积近似等于该平均值乘以(b-a)．第50页，课件共61页，创作于2023年2月例9库存问题在物资的供应过程中,由于到货和销售不可能做到同步同量,故总要保持一定的库存储备.如果库存过多,就会造成积压浪费以及保管费的上升;如果库存过少,会造成缺货.如何选择库存和订货策略,就是一个需要研究的问题.库存问题有多种类型,一般都比较复杂,下面讨论一种简单的情形.

某电动车行的仓库管理人员采取一种简单的订货策略,当库存降低到P辆电动车时就向厂家第51页，课件共61页，创作于2023年2月例9库存问题订货,每次订货Q辆,如果某一天的需求量超过了库存量,商店就有销售损失和信誉损失,但如果库存量过多,就会导致资金积压和保管费增加.若现在有如下表所示的两种库存策略,试比较选择一种策略以使总费用最少.

表订货方案方案重新订货点P/辆重新订货量Q/辆方案1125150

方案2150250第52页，课件共61页，创作于2023年2月例9库存问题这个问题的已知条件是:(1)从发出订货到收到货物需隔3天.(2)每辆电动车保管费为0.50元/天,每辆电动车的缺货损失为1.60元/天,每次的订货费为75元.(3)每天电动车需求量是0到99之间均匀分布的随机数.(4)原始库存为110辆,假设第一天没有发出订货.第53页，课件共61页，创作于2023年2月例9库存问题分析:这一问题用解析法讨论比较麻烦,但用计算机按天仿真仓库货物的变动情况却很方便.我们以30天为例,依次对这两种方案进行仿真,最后比较各方案的总费用,从而就可以作出决策.计算机仿真时的工作流程是早上到货,全天销售,晚上订货,输入一些常数和初始数据后,以一天为时间步长进行仿真.首先检查这一天是否为预订到货日期,如果是,则原有库存量加Q,并把预定到货量清为0;如果不是,则库存量不变.接着用计算机中的随机函数仿真随机需求量,若库存量大于需求量,则新的库存量减去需求量;反之,则新库存量变为0,并且要在总费用上加缺货损失,然后检查实际第54页，课件共61页，创作于2023年2月例9库存问题库存量加上预定到货量是否小于重新订货点P,如果是,则需要重新订货,这时就加一次订货费.如此重复运行30天,即可得所需费用总值.由此比较这两种方案的总费用,即得最好方案.第55页，课件共61页，创作于2023年2月cleardays=30;P=[125,150];Q=[150,250];cost=[0,0];arrivalinterval=2;storagefee=0.5;lossfee=1.6;bookfee=75;storage0=110;booknumber=0;arrivedate=0;nr=rand(days,1);fori=1:2storage(1)=storage0;n=round(99*nr(1));sale=n;remain=storage(1)-n;ifremain<=P(i);booknumber=Q(i);arrivedate=4;orderfee=bookfee;elseorderfee=0;endstorage(1)=remain;cost(i)=cost(i)+remain*storagefee+orderfee;forj=2:daysdh=j;ifabs(dh-arrivedate)<0.01storage(j)=storage(j-1)+booknumber;booknumber=0;第56页，课件共61页，创作于2023年2月arrivedate=j;elsestorage(j)=storage(j-1);endn=round(99*nr(j));ifstorage(j)>=nsale=n;remain=storage(j)-n;shortagenumber=0;elsesale=storage(j);remain=0;shortagenumber=n-storage(j);endstorage(j)=remain;ifremain+booknumber<=P(i);booknumber=Q(i);