计算机模拟方法

计算机模拟方法一．四人追逐实验模拟如图1,在正方形ABCD的四个顶点各有一个人。设在初始时刻0t时，四人同时出发匀速以v沿顺时针走向下一个人。如果他们始终对准下一个人为目标行进，最终结果会如何。作出各自的运动轨迹。解：该问题可以通过计算机模拟来实现。这需要将时间离散化。设时间间隔为t，j时刻表示时间.tjt设第i个人j时刻的位置坐标为：(,),(1,2,3,4;1,2,3,)ijijxyij对前面3个人表达式为：,1,,1,..cos(1,2,3;1,2,3,...,1)..sinijijijijxxvtxijnyyvtx其中1,,221,,1,,cos()()ijijijijijijxxxxxyy(1,2,3)i1,,221,,1,,sin()()ijijijijijijyyxxxyy(1,2,3)i对第4个人表达式为：4,14,4,14,..cos..sinjjjjxxvtxyyvtx其中1,4,221,4,1,4,cos()()jjjjjjxxxxxyy1,4,221,4,1,4,sin()()jjjjjjyyxxxyyMatlab实现程序run.m如下：%模拟运动n=240;x=zeros(4,n);y=zeros(4,n);dt=0.05;%时间间隔v=10;%速度x(1,1)=100;y(1,1)=0;%第1个人初始坐标x(2,1)=0;y(2,1)=0;%第2个人初始坐标x(3,1)=0;y(3,1)=100;%第3个人初始坐标x(4,1)=100;y(4,1)=100;%第4个人初始坐标forj=1:n-1fori=1:3d=sqrt((x(i+1,j)-x(i,j))^2+(y(i+1,j)-y(i,j))^2);%第i个人和第i+1个人距离cosx=(x(i+1,j)-x(i,j))/d;%求cos值sinx=(y(i+1,j)-y(i,j))/d;%求sin值x(i,j+1)=x(i,j)+v*dt*cosx;%求新x坐标y(i,j+1)=y(i,j)+v*dt*sinx;%求新y坐标end%考虑第1,2,3人运动一步d=sqrt((x(1,j)-x(4,j))^2+(y(1,j)-y(4,j))^2);%第4个人和第1个人距离cosx=(x(1,j)-x(4,j))/d;%求cos值sinx=(y(1,j)-y(4,j))/d;%求sin值x(4,j+1)=x(4,j)+v*dt*cosx;%求第4点新x坐标y(4,j+1)=y(4,j)+v*dt*sinx;%求第4点新y坐标end%plot(x,y)forj=1:nplot(x(1,j),y(1,j),x(2,j),y(2,j),x(3,j),y(3,j),x(4,j),y(4,j))%作点图holdon%保持每次作图,实现各次图行迭加end执行结果见图1图1模拟结果图形二、电梯问题随机模拟设有r个人在一楼进入电梯，楼上共有n层。设每个乘客在任何一层楼出电梯的可能性相同，求直到电梯中的人下完为止，电梯需停次数的数学期望。并对14,28rn进行计算机模拟验证。计算机模拟算法思想：楼上n层的序号记为1,2,,n，定义整型数组()xn,并将其n个元素的初始值置0。计算机随机产生1到n的r个随机整数，若该数为i，则令()1xi，表示第i层电梯有人要下，从而电梯在该层要停；当有多人在第i层要下时，仍用()1xi表示。统计数组()xn中为1的元素个数，即为该次模拟实验电梯要停的次数。总共模拟实验N次，求平均值就是电梯所停次数的模拟值。Matlab模拟程序dianti.mN=5000;%模拟次数n=28;%电梯层数r=14;%电梯开始进的人数s=0;x=zeros(n,1);fork=1:N%模拟N次s1=0;fori=1:nx(i)=0;endforj=1:r%对每个人是否下电梯进行模拟i=1+floor(rand(1,1)*n);x(i)=1;%第i层有人下ends1=sum(x);%该次模拟中总共要下的人数s=s+s1;%统计各次模拟中下的人数endeq=s/N;%模拟平均值ei=n*(1-(1-1/n)^r);fprintf('模拟电梯要停次数的理论值%6.2f,模拟值为%6.2f\n',ei,eq);某两次运行结果如下：11.15,11.18。与理论数值11.17很接近。三、矿工选门问题模拟实验一名矿工陷入一个有三扇门的矿井中，第一扇门通到一个隧道，走2小时可以到达安全区；第二扇门通到另一个隧道，走3小时回到矿井中；第三扇门通到又一个隧道，走5小时回到矿井中。假定矿工总是等可能的选择三扇门的任何一扇门，求到达安全区的平均时间计算机模拟的Matlab程序：N=50000;sum=0;forindex=1:Nflag=1;Time=0;whileflag==1%一次试验直到走出为止door=1+floor(3*rand(1,1));%随机选择1个门if(door==1)Time=Time+2;flag=0;%选择1号门elseif(door==2)Time=Time+3;%选择2号门elseif(door==3)Time=Time+5;%选择3号门endendsum=sum+Time;end%endforindexaver=1.0*sum/N;fprintf('平均时间为%5.2f小时.\n',aver);某三次模拟结果为：10.04，10.14，9.91。与理论值10小时很接近。基于时间步长的计算机仿真，就是按照时间流逝的顺序，一步一步地对系统进行仿真，在整个仿真过程中，通常将时间步长固定不变。在进行系统仿真过程中，可以把整个仿真过程分为许多相等的时间间隔，时间步长的长度可根据实际问题分别取为秒、分、小时等。程序中按此步长前进的时钟就是仿真的时钟。四容器水流问题仿真实验某容器内盛入盐水池有3000m3水，其中含盐3kg，以每分钟10m3的速率向水池中注入含盐率为0.5kg/m3的盐水。同时以每分钟7m3的速率从水池流出搅拌均匀的盐水。求1小时后水池中的含盐量。并做出水池含盐量的变化曲线。分析：1．理论求解设初始时刻为0，t时刻水池含盐量为()Xtkg。t时刻水池的水量为3000(107)30003ttm3经过t时间流入水池的盐量为100.55.ttt时刻水池的含盐率为：()30003Xtt经过t时间流出水池的盐量为().7.30003Xttt则水池含盐量满足关系：().7.()()5.30003XttXttXttt计算机仿真设仿真的时间步长为t分钟，n时刻时间为.nt，水池含盐量为()Xnkg由前面分析可得递推关系为：().7.(1)()5.30003..(0)3XntXnXntntX终止时间为1小时=60分钟，对应步数60ntMatlab程序为：dt=0.01;%分钟n=60/dt;%终止步数t=zeros(1,n+1);x=zeros(1,n+1);x(1)=3;t(1)=0;fori=1:nt(i+1)=t(i)+dt;x(i+1)=x(i)+5*dt-x(i)*7*dt/(3000+3*i*dt);endfprintf('时间%5.1f分钟,水池含盐量为%6.3fkg\n',t(n+1),x(n+1));plot(t,x);gridon;输出结果为：时间60.0分钟,水池含盐量为283.306kg.仿真值和理论值十分接近。图形为：0102030405060050100150200250300图2.16水池含盐量曲线图五库存问题仿真实验在物资供应过程中，由于到货与销售不可能做到同步同量，因此总要保持一定的库存量。如果库存过多，就会造成积压浪费和保管费的增大；如果库存过少，会造成缺货。如何选择库存和缺货策略，是一个值得研究的问题。库存问题有多种类型，通常比较复杂，这里我们通过计算机仿真求解一种较为简单的情况。某自行车商店的仓库管理人员采取一种简单定货策略，当库存量降低到P辆自行车就向厂家订货，每次订货Q量，如果某一天的需求量超过了库存量，商店就有缺货损失和信誉损失；但如果库存量过多，将会导致资金积压和保管费增加。若现在有六种策略见表2.12,试选择一种总费用最少的策略。表2.12订货策略方案编号订货点P量订货量Q辆112515021252503150250417525051753006200300已知条件为：（1）从发出货到收到货需要隔3天。如周一订货则周四才能到货。（2）每辆自行车每天保管费为0.75元，每辆自行车每天缺货损失为1.80元，每次订货费为75元。（3）每天自行车的市场需求量服从参数为50辆的Poisson分布。（4）原始库存为120辆，并假设第一天没有到货也没有发出订货。该问题采用理论求解比较麻烦，我们采用计算机仿真，计算每天的库存量，销售量，计算一段时间如n=200天的总费用及平均费用。然后选取总费用最少的策略。计算机仿真分析：设第i天早上（还没有营业前）商店拥有量为ix辆，白天销售后到晚上（停止营业后）拥有量为iy辆(可为负，负表示该天缺货量)，市场需求量为id辆，每天到货量为ia辆（没有到货则为0）。第i天缺货量为iu辆，库存量为is辆。保管费用ic元，缺货损失费ir元，订货费iv。则有如下关系：1）初始值第一天早上商店拥有量1120x，到货量10a.2）随机产生每天服从Poisson分布的市场需求量id,1,2,,in3）白天销售后到晚上（停止营业后）拥有量为iy满足：iiiiyxda1,2,,in4）缺货量计算式为：000iiiiyuyy5）库存量计算式为：000iiiiysyy6）到货量及订货费计算：初始置所有0,0iiav，1,2,,in计算当库存量isP，则3iaQ,75iv7）第1i天早上商店拥有量为1ix满足：1iixs1,2,,1in8)第i天保管费计算式为：0.75iics9)第i天缺货损失费计算式为：1.80iiru10)总费用计算1niiiiCostcrv11)平均每天费用/AverCostnMatlab程序：n=200;%总天数x=zeros(1,n+3);%早上拥有量y=zeros(1,n+3);%晚上拥有量a=zeros(1,n+3);%到货量d=zeros(1,n+3);%市场需求量u=zeros(1,n+3);%缺货量s=zeros(1,n+3);%库存量c=zeros(1,n+3);%保管费r=zeros(1,n+3);%缺货损失费v=zeros(1,n+3);%订货费P=120;%订货点Q=150;%订货量x(1)=120;%初始早上拥有量a(1)=0;%初始到货量d=poissrnd(50,1,n);%产生每天的市场需求fori=1:ny(i)=x(i)-d(i)+a(i);%第i天晚上拥有量if(y(i)0)u(i)=-y(i);end%当天缺货量if(y(i)0)s(i)=y(i);end%当天库存量if(s(i)P)a(i+3)=Q;%后3天将有订货到达v(i)=75;%当天订货费endx(i+1)=s(i);%第二天早上拥有量为前一晚库存量c(i)=0.75*s(i);%当天库存费r(i)=1.8*u(i);%当天缺货损失费endCost=sum(c(1:n))+sum(r(1:n))+sum(v(1:n));%总费用Aver=Cost/n;%平均费用fprintf('当P=%2d,Q=%2d,n=%3d天,总费用%6.2f,平均费用%6.2f\n',P,Q,n,Cost,Aver);某次仿真实验结果为：当P=125,Q=150,n=200天,总费用35630.10,平均费用178.15当P=125,Q=250,n=200天,总费用48503.70,平均费用242.5