計(jì)算機(jī)模擬方法

計(jì)算機(jī)模擬方法一．四人追逐實(shí)驗(yàn)?zāi)M如圖1,在正方形ABCD的四個(gè)頂點(diǎn)各有一個(gè)人。設(shè)在初始時(shí)刻t=0時(shí)，四人同時(shí)出發(fā)勻速以v沿順時(shí)針走向下一個(gè)人。如果他們始終對(duì)準(zhǔn)下一個(gè)人為目標(biāo)行進(jìn)，最終結(jié)果會(huì)如何。作出各自的運(yùn)動(dòng)軌跡。解：該問題可以通過計(jì)算機(jī)模擬來(lái)實(shí)現(xiàn)。這需要將時(shí)間離散化。設(shè)時(shí)間間隔為At，j時(shí)刻表示時(shí)間t=j.At設(shè)第i個(gè)人j時(shí)刻的位置坐標(biāo)為：(x,y),(i=1,2,3,4；j=1,2,3,)ijijTOC\o"1-5"\h\z對(duì)前面3個(gè)人表達(dá)式為： …x=x+vA.cosx-齊1 -j人. (i=1,2,3= 1,2,n, 1)y=y+vAtsinxi,A1 i,j其中cosx(i=1,2,3sinx其中cosx(i=1,2,3sinx=-x)2+(y-y)2i,j i+1,j i,j-yy<(x-x)2+(y-y)2理 i+1,j i,j i+1,j i,j(i=1,2,3對(duì)第4個(gè)人表達(dá)式為：x =x +v.At.cexs4j+1j4,y =y +v.At.sixn其中4j,+1j4,其中TOC\o"1-5"\h\zx-xcosx= 斗—(x-x)2+(y-y)21,j 4,j 1,j 4,jsinx=ysinx=9 4,j —(x-x)2+(y-y)21,j4,j 1,j4,jMatlab實(shí)現(xiàn)程序run.m如下:%模擬運(yùn)動(dòng)n=240;x=zeros(4,n);y=zeros(4,n);dt=0.05;%時(shí)間間隔v=10;%速度x(1,1)=100;y(1,1)=0;%第1個(gè)人初始坐標(biāo)x(2,1)=0;y(2,1)=0;%第2個(gè)人初始坐標(biāo)x(3,1)=0;y(3,1)=100;%第3個(gè)人初始坐標(biāo)x(4,1)=100;y(4,1)=100;%第4個(gè)人初始坐標(biāo)forj=1:n-1fori=1:3d=sqrt((x(i+l,j)-x(i,j))人2+(y(i+l,j)-y(i,j))人2);%第i個(gè)人和第i+1個(gè)人距離cosx=(x(i+l,j)-x(i,j))/d;%求cos值sinx=(y(i+l,j)-y(i,j))/d;%求sin值x(i,j+1)=x(i,j)+v*dt*cosx;%求新x坐標(biāo)y(i,j+1)=y(i,j)+v*dt*sinx;%求新y坐標(biāo)end%考慮第1,2,3人運(yùn)動(dòng)一步d=sqrt((x(1,j)-x(4,j))A2+(y(1,j)-y(4,j))A2);%第4個(gè)人和第1個(gè)人距離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點(diǎn)新x坐標(biāo)y(4,j+1)=y(4,j)+v*dt*sinx; %求第4點(diǎn)新y坐標(biāo)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))%作點(diǎn)圖holdon%保持每次作圖,實(shí)現(xiàn)各次圖行迭加end執(zhí)行結(jié)果見圖1

二、電梯問題隨機(jī)模擬設(shè)有r個(gè)人在一樓進(jìn)入電梯，樓上共有n層。設(shè)每個(gè)乘客在任何一層樓出電梯的可能性相同，求直到電梯中的人下完為止，電梯需停次數(shù)的數(shù)學(xué)期望。并對(duì)r=14,n=28進(jìn)行計(jì)算機(jī)模擬驗(yàn)證。計(jì)算機(jī)模擬算法思想：樓上n層的序號(hào)記為1,2, ,n，定義整型數(shù)組x(n),并將其n個(gè)元素的初始值置0。計(jì)算機(jī)隨機(jī)產(chǎn)生1到n的r個(gè)隨機(jī)?整數(shù)，若該數(shù)為i，則令x(i)=1，表示第i層電梯有人要下，從而電梯在該層要停；當(dāng)有多人在第i層要下時(shí)，仍用x(i)=1表示。統(tǒng)計(jì)數(shù)組x(n)中為1的元素個(gè)數(shù)，即為該次模擬實(shí)驗(yàn)電梯要停的次數(shù)?？偣材M實(shí)驗(yàn)N次，求平均值就是電梯所停次數(shù)的模擬值。Matlab模擬程序dianti.mN=5000;%模擬次數(shù)n=28;%電梯層數(shù)r=14;%電梯開始進(jìn)的人數(shù)s=0;x=zeros(n,1);fork=l:N%模擬N次s1=0;fori=l:nx(i)=0;endforj=l:r%對(duì)每個(gè)人是否下電梯進(jìn)行模擬i=l+floor(rand(l,l)*n);x(i)=1;%第：層有人下ends1=sum(x);%該次模擬中總共要下的人數(shù)s=s+s1;%統(tǒng)計(jì)各次模擬中下的人數(shù)endeq=s/N;%模擬平均值ei=n*(1-(1-1/n)Ar);fprintf(模擬電梯要停次數(shù)的理論值％6.2f,模擬值為％6.2f\n',ei,eq);某兩次運(yùn)行結(jié)果如下：11.15,11.18。與理論數(shù)值11.17很接近。三、礦工選門問題模擬實(shí)驗(yàn)一名礦工陷入一個(gè)有三扇門的礦井中，第一扇門通到一個(gè)隧道，走2小時(shí)可以到達(dá)安全區(qū)；第二扇門通到另一個(gè)隧道，走3小時(shí)回到礦井中；第三扇門通到又一個(gè)隧道，走5小時(shí)回到礦井中。假定礦工總是等可能的選擇三扇門的任何一扇門，求到達(dá)安全區(qū)的平均時(shí)間計(jì)算機(jī)模擬的Matlab程序：N=50000;sum=0;forindex=1:Nflag=1;Time=0;whileflag==1%一次試驗(yàn)直到走出為止door=1+floor(3*rand(1,1));%隨機(jī)選擇1個(gè)門if(door==1)Time=Time+2;flag=0;%選擇1號(hào)門elseif(door==2)Time=Time+3;%選擇2號(hào)門

elseif(door==3)Time=Time+5;%選擇3號(hào)門endendsum=sum+Time;end%endforindexaver=1.0*sum/N;fprintf(平均時(shí)間為%5.2f小時(shí).\n',aver);某三次模擬結(jié)果為：10.04，10.14，9.91。與理論值10小時(shí)很接近。基于時(shí)間步長(zhǎng)的計(jì)算機(jī)仿真，就是按照時(shí)間流逝的順序，一步一步地對(duì)系統(tǒng)進(jìn)行仿真，在整個(gè)仿真過程中，通常將時(shí)間步長(zhǎng)固定不變。在進(jìn)行系統(tǒng)仿真過程中，可以把整個(gè)仿真過程分為許多相等的時(shí)間間隔，時(shí)間步長(zhǎng)的長(zhǎng)度可根據(jù)實(shí)際問題分別取為秒、分、小時(shí)等。程序中按此步長(zhǎng)前進(jìn)的時(shí)鐘就是仿真的時(shí)鐘。四容器水流問題仿真實(shí)驗(yàn)?zāi)橙萜鲀?nèi)盛入鹽水池有3000m3水，其中含鹽3kg，以每分鐘10m3的速率向水池中注入含鹽率為0.5kg/m3的鹽水。同時(shí)以每分鐘7m3的速率從水池流出攪拌均勻的鹽水。求1小時(shí)后水池中的含鹽量。并做出水池含鹽量的變化曲線。分析：1．理論求解設(shè)初始時(shí)刻為0,t時(shí)刻水池含鹽量為X(t)kg。t時(shí)刻水池的水量為3000+(10—7)t=3000+3tm3經(jīng)過At時(shí)間流入水池的鹽量為10x0.5xAt=5.Att時(shí)刻水池的含鹽率為:t時(shí)刻水池的含鹽率為:X(t)

3000+3?經(jīng)過At時(shí)間流出水池的鹽量為則水池含鹽量滿足關(guān)系:x(t+A酉=X(酉+沁品3計(jì)算機(jī)仿真設(shè)仿真的時(shí)間步長(zhǎng)為At分鐘，n時(shí)刻時(shí)間為n.At，水池含鹽量為X(n)kg由前面分析可得遞推關(guān)系為:1、 、uA X(n).7.AtX(n+1)=X(n)+5.At— < 3000+3.n.At、X(0)=360終止時(shí)間為1小時(shí)=60分鐘，對(duì)應(yīng)步數(shù)n=-AtMatlab程序?yàn)椋篸t=0.01;%分鐘n=60/dt;%終止步數(shù)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(時(shí)間％5.1f分鐘,水池含鹽量為％6.3fkg\n',t(n+1),x(n+1));plot(t,x);gridon;輸出結(jié)果為：時(shí)間60.0分鐘,水池含鹽量為283.306kg.仿真值和理論值十分接近。圖形為：

圖2.16水池含鹽量曲線圖五庫(kù)存問題仿真實(shí)驗(yàn)在物資供應(yīng)過程中，由于到貨與銷售不可能做到同步同量，因此總要保持一定的庫(kù)存量。如果庫(kù)存過多，就會(huì)造成積壓浪費(fèi)和保管費(fèi)的增大；如果庫(kù)存過少，會(huì)造成缺貨。如何選擇庫(kù)存和缺貨策略，是一個(gè)值得研究的問題。庫(kù)存問題有多種類型，通常比較復(fù)雜，這里我們通過計(jì)算機(jī)仿真求解一種較為簡(jiǎn)單的情況。某自行車商店的倉(cāng)庫(kù)管理人員采取一種簡(jiǎn)單定貨策略，當(dāng)庫(kù)存量降低到P輛自行車就向廠家訂貨，每次訂貨Q量，如果某一天的需求量超過了庫(kù)存量，商店就有缺貨損失和信譽(yù)損失；但如果庫(kù)存量過多，將會(huì)導(dǎo)致資金積壓和保管費(fèi)增加。若現(xiàn)在有六種策略見表2.12,試選擇一種總費(fèi)用最少的策略。表2.12訂貨策略方案編號(hào)訂貨點(diǎn)P量訂貨量Q輛112515021252503150250417525051753006200300已知條件為：（1） 從發(fā)出貨到收到貨需要隔3天。如周一訂貨則周四才能到貨。（2） 每輛自行車每天保管費(fèi)為0.75元，每輛自行車每天缺貨損失為1.80元，每次訂貨費(fèi)為75元。（3） 每天自行車的市場(chǎng)需求量服從參數(shù)為50輛的Poisson分布。（4）原始庫(kù)存為120輛，并假設(shè)第一天沒有到貨也沒有發(fā)出訂貨。該問題采用理論求解比較麻煩，我們采用計(jì)算機(jī)仿真，計(jì)算每天的庫(kù)存量，銷售量，計(jì)算一段時(shí)間如n=200天的總費(fèi)用及平均費(fèi)用。然后選取總費(fèi)用最少的策略。計(jì)算機(jī)仿真分析：設(shè)第i天早上（還沒有營(yíng)業(yè)前）商店擁有量為x輛，白天銷售后到晚上（停止?fàn)I業(yè)后）i擁有量為y輛（可為負(fù)，負(fù)表示該天缺貨量），市場(chǎng)需求量為d輛，每天到貨量為a輛（沒i i i有到貨則為0）。第i天缺貨量為u輛，庫(kù)存量為s輛。保管費(fèi)用C元，缺貨損失費(fèi)r元，訂i i i i貨費(fèi)v。i則有如下關(guān)系：1） 初始值第一天早上商店擁有量x=120，到貨量a=0.112） 隨機(jī)產(chǎn)生每天服從Poisson分布的市場(chǎng)需求量d,i二1,2, ,ni3） 白天銷售后到晚上（停止?fàn)I業(yè)后）擁有量為y滿足：…iy=x一d+ai=1,2,n,iiiiTOC\o"1-5"\h\z4）缺貨量計(jì)算式為： ???‘0 y>0U=< ii［一y y<0ii5） 庫(kù)存量計(jì)算式為：‘0 y<0S=< iiy y>0ii6） 到貨量及訂貨費(fèi)計(jì)算：初始置所有a=0,v=0，i=1,2, ,nii計(jì)算當(dāng)庫(kù)存量s<P，則a=Q,?v=75i i+3 i7） 第i+1天早上商店擁有量為x滿足：i+1x=si=1,2,n,一1i+1 i8） 第i天保管費(fèi)計(jì)算式為：c=0.75xsii9） 第i天缺貨損失費(fèi)計(jì)算式為：r=1.8x0uii10)總費(fèi)用計(jì)算Cost=￡c+r+viiii=111)平均每天費(fèi)用Aver=Cost/nMatlab程序：n=200;%總天數(shù)x=zeros(l,n+3);%早上擁有量y=zeros(l,n+3);%晚上擁有量a=zeros(1,n+3);%到貨量d=zeros(l,n+3);%市場(chǎng)需求量u=zeros(l,n+3);%缺貨量s=zeros(l,n+3);%庫(kù)存量c=zeros(l,n+3);%保管費(fèi)r=zeros(l,n+3);%缺貨損失費(fèi)v=zeros(l,n+3);%訂貨費(fèi)P=l20;%訂貨點(diǎn)Q=l50;%訂貨量x(l)=l20;%初始早上擁有量a(l)=0;%初始到貨量d=poissrnd(50,l,n);%產(chǎn)生每天的市場(chǎng)需求fori=l:ny(i)=x(i)-d(i)+a(i);%第i天晚上擁有量if(y(i)<0)u(i)=-y(i);end%當(dāng)天缺貨量if(y(i)>0)s(i)=y(i);end %當(dāng)天庫(kù)存量if(s(i)<P)a(i+3)=Q;%后3天將有訂貨到達(dá)v(i)=75;%當(dāng)天訂貨費(fèi)endx(i+l)=s(i);%第二天早上擁有量為前一晚庫(kù)存量c(i)=0.75*s(i%當(dāng)天庫(kù)存費(fèi)r(i)=1.8*u(i)當(dāng)天缺貨損失費(fèi)endCost=sum(c(l:n))+sum(r(l:n))+sum(v(l:n總費(fèi)用Aver=Cost/n;平%均費(fèi)用fprint當(dāng)'P=%2d,Q=%2d,n=%3d天,總費(fèi)用％6.2f,平均費(fèi)用％6.2f\n',P,Q,n,Cost,Aver);某次仿真實(shí)驗(yàn)結(jié)果為：當(dāng)P=125,Q=150,n=200天,總費(fèi)用35630.10平均費(fèi)用178.15當(dāng)P=12