計算機(jī)模擬方法

1、 計算機(jī)模擬方法一四人追逐實(shí)驗(yàn)?zāi)M如圖1,在正方形ABCD的四個頂點(diǎn)各有一個人。設(shè)在初始時刻時，四人同時出發(fā)勻速以沿順時針走向下一個人。如果他們始終對準(zhǔn)下一個人為目標(biāo)行進(jìn)，最終結(jié)果會如何。作出各自的運(yùn)動軌跡。 解：該問題可以通過計算機(jī)模擬來實(shí)現(xiàn)。這需要將時間離散化。設(shè)時間間隔為，時刻表示時間 設(shè)第個人時刻的位置坐標(biāo)為： 對前面3個人表達(dá)式為： 其中 對第4個人表達(dá)式為： 其中 Matlab實(shí)現(xiàn)程序run.m如下：%模擬運(yùn)動n=240;x=zeros(4,n);y=zeros(4,n);dt=0.05; %時間間隔v=10; %速度x(1,1)=100; y(1,1)=0; %第1個人初始坐標(biāo)x

2、(2,1)=0; y(2,1)=0; %第2個人初始坐標(biāo)x(3,1)=0; y(3,1)=100; %第3個人初始坐標(biāo)x(4,1)=100; y(4,1)=100; %第4個人初始坐標(biāo)for j=1:n-1 for i=1:3 d=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坐標(biāo) y(i,j+1)=y(i,j)+v*dt*sinx; %求

3、新y坐標(biāo) end %考慮第1,2,3人運(yùn)動一步 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點(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)for j=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)圖ho

4、ld on %保持每次作圖,實(shí)現(xiàn)各次圖行迭加end執(zhí)行結(jié)果見圖1 圖1 模擬結(jié)果圖形二、電梯問題隨機(jī)模擬 設(shè)有個人在一樓進(jìn)入電梯，樓上共有層。設(shè)每個乘客在任何一層樓出電梯的可能性相同，求直到電梯中的人下完為止，電梯需停次數(shù)的數(shù)學(xué)期望。并對進(jìn)行計算機(jī)模擬驗(yàn)證。 計算機(jī)模擬算法思想： 樓上層的序號記為，定義整型數(shù)組,并將其個元素的初始值置0。計算機(jī)隨機(jī)產(chǎn)生1到的個隨機(jī)整數(shù)，若該數(shù)為，則令，表示第層電梯有人要下，從而電梯在該層要停；當(dāng)有多人在第層要下時，仍用表示。統(tǒng)計數(shù)組中為1的元素個數(shù)，即為該次模擬實(shí)驗(yàn)電梯要停的次數(shù)?？偣材M實(shí)驗(yàn)次，求平均值就是電梯所停次數(shù)的模擬值。 Matlab模擬程序dian

5、ti.mN=5000; %模擬次數(shù)n=28; %電梯層數(shù)r=14; %電梯開始進(jìn)的人數(shù)s=0;x=zeros(n,1); for k=1:N %模擬N次 s1=0; for i=1:n x(i)=0; end for j=1:r %對每個人是否下電梯進(jìn)行模擬 i=1+floor(rand(1,1)*n); x(i)=1; %第i層有人下 end s1=sum(x); %該次模擬中總共要下的人數(shù) s=s+s1; %統(tǒng)計各次模擬中下的人數(shù) end eq=s/N; %模擬平均值 ei=n*(1-(1-1/n)r); fprintf(模擬電梯要停次數(shù)的理論值%6.2f,模擬值為%6.2fn,ei,eq

6、);某兩次運(yùn)行結(jié)果如下：11.15,11.18。與理論數(shù)值11.17很接近。三、 礦工選門問題模擬實(shí)驗(yàn) 一名礦工陷入一個有三扇門的礦井中，第一扇門通到一個隧道，走2小時可以到達(dá)安全區(qū)；第二扇門通到另一個隧道，走3小時回到礦井中；第三扇門通到又一個隧道，走5小時回到礦井中。假定礦工總是等可能的選擇三扇門的任何一扇門，求到達(dá)安全區(qū)的平均時間計算機(jī)模擬的Matlab程序： N=50000; sum=0;for index=1:N flag=1; Time=0;while flag=1 %一次試驗(yàn)直到走出為止door=1+floor(3*rand(1,1); %隨機(jī)選擇1個門if(door=1) Ti

7、me=Time+2; flag=0;%選擇1號門elseif(door=2) Time=Time+3; %選擇2號門elseif(door=3) Time=Time+5; %選擇3號門 end end sum=sum+Time; end %end for index aver=1.0*sum/N; fprintf(平均時間為%5.2f小時.n,aver);某三次模擬結(jié)果為：10.04，10.14，9.91。與理論值10小時很接近。 基于時間步長的計算機(jī)仿真，就是按照時間流逝的順序，一步一步地對系統(tǒng)進(jìn)行仿真，在整個仿真過程中，通常將時間步長固定不變。 在進(jìn)行系統(tǒng)仿真過程中，可以把整個仿真過程分為

8、許多相等的時間間隔，時間步長的長度可根據(jù)實(shí)際問題分別取為秒、分、小時等。程序中按此步長前進(jìn)的時鐘就是仿真的時鐘。四 容器水流問題仿真實(shí)驗(yàn) 某容器內(nèi)盛入鹽水池有3000m3水，其中含鹽3kg，以每分鐘10 m3的速率向水池中注入含鹽率為0.5kg/ m3的鹽水。同時以每分鐘7 m3的速率從水池流出攪拌均勻的鹽水。求1小時后水池中的含鹽量。并做出水池含鹽量的變化曲線。分析：1 理論求解設(shè)初始時刻為0，t時刻水池含鹽量為kg。t時刻水池的水量為 m3經(jīng)過時間流入水池的鹽量為 t時刻水池的含鹽率為：經(jīng)過時間流出水池的鹽量為 則水池含鹽量滿足關(guān)系： 計算機(jī)仿真設(shè)仿真的時間步長為分鐘，時刻時間為，水池含鹽

9、量為kg由前面分析可得遞推關(guān)系為： 終止時間為1小時=60分鐘，對應(yīng)步數(shù) Matlab程序?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;for i=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.3f kgn,t(n+1),x(n+1);plot(t,x);grid on;輸出結(jié)果為：時間 60.0分鐘,水池含鹽量為283.306 kg.仿真值和理論值十分接近。圖形

10、為： 圖2.16 水池含鹽量曲線圖五 庫存問題仿真實(shí)驗(yàn) 在物資供應(yīng)過程中，由于到貨與銷售不可能做到同步同量，因此總要保持一定的庫存量。如果庫存過多，就會造成積壓浪費(fèi)和保管費(fèi)的增大；如果庫存過少，會造成缺貨。如何選擇庫存和缺貨策略，是一個值得研究的問題。庫存問題有多種類型，通常比較復(fù)雜，這里我們通過計算機(jī)仿真求解一種較為簡單的情況。某自行車商店的倉庫管理人員采取一種簡單定貨策略，當(dāng)庫存量降低到P輛自行車就向廠家訂貨，每次訂貨Q量，如果某一天的需求量超過了庫存量，商店就有缺貨損失和信譽(yù)損失；但如果庫存量過多，將會導(dǎo)致資金積壓和保管費(fèi)增加。若現(xiàn)在有六種策略見表2.12,試選擇一種總費(fèi)用最少的策略。

11、表2.12 訂貨策略方案編號訂貨點(diǎn)P量訂貨量Q輛112515021252503150250417525051753006200300已知條件為：（1） 從發(fā)出貨到收到貨需要隔3天。如周一訂貨則周四才能到貨。（2） 每輛自行車每天保管費(fèi)為0.75元，每輛自行車每天缺貨損失為1.80元，每次訂貨費(fèi)為75元。（3） 每天自行車的市場需求量服從參數(shù)為50輛的Poisson分布。（4） 原始庫存為120輛，并假設(shè)第一天沒有到貨也沒有發(fā)出訂貨。 該問題采用理論求解比較麻煩，我們采用計算機(jī)仿真，計算每天的庫存量，銷售量，計算一段時間如n=200天的總費(fèi)用及平均費(fèi)用。然后選取總費(fèi)用最少的策略。 計算機(jī)仿真分析

12、： 設(shè)第天早上（還沒有營業(yè)前）商店擁有量為輛，白天銷售后到晚上（停止?fàn)I業(yè)后）擁有量為輛(可為負(fù)，負(fù)表示該天缺貨量)，市場需求量為輛，每天到貨量為輛（沒有到貨則為0）。第天缺貨量為輛，庫存量為輛。保管費(fèi)用元，缺貨損失費(fèi)元，訂貨費(fèi)。 則有如下關(guān)系： 1）初始值第一天早上商店擁有量，到貨量. 2）隨機(jī)產(chǎn)生每天服從Poisson分布的市場需求量, 3）白天銷售后到晚上（停止?fàn)I業(yè)后）擁有量為滿足： 4）缺貨量計算式為： 5）庫存量計算式為： 6）到貨量及訂貨費(fèi)計算： 初始置所有， 計算當(dāng)庫存量，則,7）第天早上商店擁有量為滿足： 8) 第天保管費(fèi)計算式為： 9) 第天缺貨損失費(fèi)計算式為： 10)總費(fèi)用計

13、算11)平均每天費(fèi)用 Matlab程序：n=200; %總天數(shù)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); %保管費(fèi)r=zeros(1,n+3); %缺貨損失費(fèi)v=zeros(1,n+3); %訂貨費(fèi)P=120; %訂貨點(diǎn)Q=150; %訂貨量x(1)=120; %初始早上擁有量a(1)=0; %初始到貨量d=poissrnd(50,1,n); %產(chǎn)生每天的市場需

14、求for i=1:n y(i)=x(i)-d(i)+a(i); %第i天晚上擁有量 if(y(i)0) s(i)=y(i); end %當(dāng)天庫存量 if(s(i)P) a(i+3)=Q; %后3天將有訂貨到達(dá) v(i)=75 ; %當(dāng)天訂貨費(fèi) end x(i+1)=s(i); %第二天早上擁有量為前一晚庫存量 c(i)=0.75*s(i); %當(dāng)天庫存費(fèi) r(i)=1.8*u(i) ;%當(dāng)天缺貨損失費(fèi) endCost=sum(c(1:n)+sum(r(1:n)+sum(v(1:n); %總費(fèi)用Aver=Cost/n; %平均費(fèi)用fprintf(當(dāng)P=%2d,Q=%2d,n=%3d天,總費(fèi)用%6.2f,平均費(fèi)用%6.2fn,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=125,Q=250,n=200天,總費(fèi)用48503.70,平均費(fèi)用242.52當(dāng)P=150