第6講蒙特卡羅方法與計(jì)算機(jī)模擬(第2次課)

1、第六講第六講 蒙特卡羅方法蒙特卡羅方法與計(jì)算機(jī)模擬與計(jì)算機(jī)模擬1離散系統(tǒng)的計(jì)算機(jī)模擬實(shí)例一范例范例 海港系統(tǒng)的卸載貨物問(wèn)題海港系統(tǒng)的卸載貨物問(wèn)題(p110-111,p119-129)第六講第六講 蒙特卡羅方法蒙特卡羅方法與計(jì)算機(jī)模擬與計(jì)算機(jī)模擬2海港系統(tǒng)的卸載貨物-問(wèn)題1 問(wèn)題問(wèn)題1 海港系統(tǒng)的卸載貨物海港系統(tǒng)的卸載貨物 (p122)程序片段程序片段(船只到港時(shí)間船只到港時(shí)間均勻分布均勻分布,船只卸貨時(shí)間船只卸貨時(shí)間均勻分布均勻分布)ShipBetweenTime(1)=unifrnd(15,145,1,1);%船只到港間隔時(shí)間隨機(jī)化船只到港間隔時(shí)間隨機(jī)化(均勻分布均勻分布)ShipUnloa

2、dTime(1)=unifrnd(45,90,1,1);%船只卸貨時(shí)間隨機(jī)化船只卸貨時(shí)間隨機(jī)化(均勻分布均勻分布)通用程序通用程序haibor1.m可實(shí)現(xiàn)多次模擬，結(jié)果保存到可實(shí)現(xiàn)多次模擬，結(jié)果保存到H1.txtdelete H1.txt %清除歷史數(shù)據(jù)清除歷史數(shù)據(jù)harbor1(100,15,145,45,90)load H1.txt; Hmean1=mean(H1); %導(dǎo)入導(dǎo)入H1并按列取平均值并按列取平均值第六講第六講 蒙特卡羅方法蒙特卡羅方法與計(jì)算機(jī)模擬與計(jì)算機(jī)模擬3海港系統(tǒng)的卸載貨物-問(wèn)題2問(wèn)題問(wèn)題2 海港系統(tǒng)的卸載貨物海港系統(tǒng)的卸載貨物 (p124)程序片段程序片段(船只到港時(shí)間

3、船只到港時(shí)間指數(shù)指數(shù)分布分布,船只卸貨時(shí)間船只卸貨時(shí)間均勻分布均勻分布)ShipBetweenTime(1)=exprnd(60,1,1);%船只到港間隔時(shí)間隨機(jī)化船只到港間隔時(shí)間隨機(jī)化(指數(shù)分布指數(shù)分布)ShipUnloadTime(1)=unifrnd(45,90,1,1);%船只卸貨時(shí)間隨機(jī)化船只卸貨時(shí)間隨機(jī)化(均勻分布均勻分布)通用程序通用程序haibor2.m可實(shí)現(xiàn)多次模擬，結(jié)果保存到可實(shí)現(xiàn)多次模擬，結(jié)果保存到H2.txtdelete H2.txt %清除歷史數(shù)據(jù)清除歷史數(shù)據(jù)harbor2(100,60,45,90)load H2.txt; Hmean2=mean(H2); %導(dǎo)入導(dǎo)

4、入H2并按列取平均值并按列取平均值第六講第六講 蒙特卡羅方法蒙特卡羅方法與計(jì)算機(jī)模擬與計(jì)算機(jī)模擬4海港系統(tǒng)的卸載貨物-問(wèn)題3-1問(wèn)題問(wèn)題3 海港系統(tǒng)的卸載貨物海港系統(tǒng)的卸載貨物 (p124-129)程序片段程序片段(船只到港時(shí)間船只到港時(shí)間離散分布離散分布,船只卸貨時(shí)間船只卸貨時(shí)間離散分布離散分布)1 編寫(xiě)編寫(xiě)船只到港間隔船只到港間隔離散累積分布函數(shù)并作階梯圖：離散累積分布函數(shù)并作階梯圖：xs=15:10:145; for i=1:length(xs)-1,x(i)=(xs(i)+xs(i+1)/2;endpx=0.009,0.029,0.035,0.051,0.090,0.161,0.200

5、,0.172,0.125,0.071,0.037,0.017,0.003; Fx=0; for i=1:length(px), Fx=Fx,sum(px(1:i); endplot(10,x,Fx,-rs); hold on; stairs(0,x-5,145,Fx,1);set(gca,xtick,0:5:145); set(gca,xgrid,on); axis tight;第六講第六講 蒙特卡羅方法蒙特卡羅方法與計(jì)算機(jī)模擬與計(jì)算機(jī)模擬5海港系統(tǒng)的卸載貨物-問(wèn)題3-2問(wèn)題問(wèn)題3 海港系統(tǒng)的卸載貨物海港系統(tǒng)的卸載貨物 (p124-129)程序片段程序片段(船只到港時(shí)間船只到港時(shí)間離散分布離散

6、分布,船只卸貨時(shí)間船只卸貨時(shí)間離散分布離散分布)2 編寫(xiě)編寫(xiě)船只到港間隔船只到港間隔離散累積分布反函數(shù)并作線(xiàn)性插值：離散累積分布反函數(shù)并作線(xiàn)性插值：Fxi=0:0.001:1-eps;xi=interp1(Fx,0,x,Fxi,linear); index=find(xixs(end);xi(index)=; Fxi(index)=; r=rand(1,n);rnd=;for i=1:n index=find(r(i)=Fxi); rnd=rnd,xi(index(1);end%以上程序已編寫(xiě)通用以上程序已編寫(xiě)通用M函數(shù)文件函數(shù)文件 harborrnd(xs,px,n)%即給出即給出n個(gè)滿(mǎn)足離

7、散分布個(gè)滿(mǎn)足離散分布(xs,px)的的船只到港間隔船只到港間隔隨機(jī)數(shù)隨機(jī)數(shù)第六講第六講 蒙特卡羅方法蒙特卡羅方法與計(jì)算機(jī)模擬與計(jì)算機(jī)模擬6海港系統(tǒng)的卸載貨物-問(wèn)題3-3問(wèn)題問(wèn)題3 海港系統(tǒng)的卸載貨物海港系統(tǒng)的卸載貨物 (p124-129)程序片段程序片段(船只到港時(shí)間船只到港時(shí)間離散分布離散分布,船只卸貨時(shí)間船只卸貨時(shí)間離散分布離散分布)3 編寫(xiě)編寫(xiě)船只卸貨時(shí)間船只卸貨時(shí)間離散累積分布函數(shù)并作階梯圖：離散累積分布函數(shù)并作階梯圖：xs=45:5:90; for i=1:length(xs)-1,x(i)=(xs(i)+xs(i+1)/2;endpx=0.017,0.045,0.095,0.086

8、,0.130,0.185,0.208,0.143,0.091; Fx=0; for i=1:length(px), Fx=Fx,sum(px(1:i); endplot(40,x,Fx,-rs); hold on; stairs(40,x-2.5,90,Fx,1);set(gca,xtick,40:2.5:90); set(gca,xgrid,on); axis tight;4042.54547.55052.55557.56062.56567.57072.57577.58082.58587.59000.20.40.60.81第六講第六講 蒙特卡羅方法蒙特卡羅方法與計(jì)算機(jī)模擬與計(jì)算機(jī)模擬7海港系

9、統(tǒng)的卸載貨物-問(wèn)題3-4問(wèn)題問(wèn)題3 海港系統(tǒng)的卸載貨物海港系統(tǒng)的卸載貨物 (p124-129)程序片段程序片段(船只到港時(shí)間船只到港時(shí)間離散分布離散分布,船只卸貨時(shí)間船只卸貨時(shí)間離散分布離散分布)4 編寫(xiě)編寫(xiě)船只卸貨時(shí)間船只卸貨時(shí)間離散累積分布反函數(shù)并作線(xiàn)性插值：離散累積分布反函數(shù)并作線(xiàn)性插值：Fxi=0:0.001:1-eps;xi=interp1(Fx,0,x,Fxi,linear); index=find(xixs(end);xi(index)=; Fxi(index)=; r=rand(1,n);rnd=;for i=1:n index=find(r(i)=Fxi); rnd=rnd,

10、xi(index(1);end%以上程序已編寫(xiě)通用以上程序已編寫(xiě)通用M函數(shù)文件函數(shù)文件 harborrnd(xs,px,n)%即給出即給出n個(gè)滿(mǎn)足離散分布個(gè)滿(mǎn)足離散分布(xs,px)的的船只卸貨時(shí)間船只卸貨時(shí)間隨機(jī)數(shù)隨機(jī)數(shù)第六講第六講 蒙特卡羅方法蒙特卡羅方法與計(jì)算機(jī)模擬與計(jì)算機(jī)模擬8海港系統(tǒng)的卸載貨物-問(wèn)題3-5問(wèn)題問(wèn)題3 海港系統(tǒng)的卸載貨物海港系統(tǒng)的卸載貨物 (p124-129)程序片段程序片段(船只到港時(shí)間船只到港時(shí)間指數(shù)指數(shù)分布分布,船只卸貨時(shí)間船只卸貨時(shí)間均勻分布均勻分布)5 模擬模擬船只到港間隔船只到港間隔 / 卸貨時(shí)間卸貨時(shí)間均為離散分布的海港系統(tǒng)均為離散分布的海港系統(tǒng)ShipB

11、etweenTime(1)=harborrnd(sbtxs,sbtpx,1);%船只到港間隔時(shí)間隨機(jī)化船只到港間隔時(shí)間隨機(jī)化(離散分布離散分布)ShipUnloadTime(1)=harborrnd(sutxs,sutpx,1);%船只卸貨時(shí)間隨機(jī)化船只卸貨時(shí)間隨機(jī)化(離散分布離散分布)通用程序通用程序haibor3.m可實(shí)現(xiàn)多次模擬，結(jié)果保存到可實(shí)現(xiàn)多次模擬，結(jié)果保存到H3.txtdelete H3.txt %清除歷史數(shù)據(jù)清除歷史數(shù)據(jù)load harbor.mat %載入數(shù)據(jù)載入數(shù)據(jù)harbor3(100,sbtxs,sbtpx,sutxs,sutpx)load H3.txt; Hmean3

12、=mean(H3); 第六講第六講 蒙特卡羅方法蒙特卡羅方法與計(jì)算機(jī)模擬與計(jì)算機(jī)模擬9離散系統(tǒng)的計(jì)算機(jī)模擬實(shí)例二范例范例 報(bào)童售報(bào)問(wèn)題報(bào)童售報(bào)問(wèn)題:一報(bào)童每天清晨從郵局訂購(gòu)報(bào)一報(bào)童每天清晨從郵局訂購(gòu)報(bào)紙后零售，每份報(bào)紙進(jìn)價(jià)紙后零售，每份報(bào)紙進(jìn)價(jià)0.35元，售價(jià)元，售價(jià)0.5元，郵元，郵局要求最低訂購(gòu)數(shù)量為局要求最低訂購(gòu)數(shù)量為60份，根據(jù)過(guò)去經(jīng)驗(yàn)一個(gè)份，根據(jù)過(guò)去經(jīng)驗(yàn)一個(gè)報(bào)童一天平均售出報(bào)紙報(bào)童一天平均售出報(bào)紙120份份(且滿(mǎn)足泊松分布且滿(mǎn)足泊松分布)，未售出的報(bào)紙只要沒(méi)有破損可退給郵局，試求報(bào)未售出的報(bào)紙只要沒(méi)有破損可退給郵局，試求報(bào)童每天清晨訂購(gòu)多少份報(bào)紙可獲最大利潤(rùn)？童每天清晨訂購(gòu)多少份報(bào)紙

13、可獲最大利潤(rùn)？1 數(shù)學(xué)建模數(shù)學(xué)建模()(,)() ()outoutininoutininoutinoutinoutinnpnpnnProfit nnnppnnPin=0.35元元 Pout=0.5元元Nin60,200區(qū)間區(qū)間Nout=poissrnd(120)隨機(jī)數(shù)隨機(jī)數(shù)第六講第六講 蒙特卡羅方法蒙特卡羅方法與計(jì)算機(jī)模擬與計(jì)算機(jī)模擬10報(bào)童售報(bào)問(wèn)題-22 計(jì)算機(jī)模擬計(jì)算機(jī)模擬當(dāng)訂購(gòu)當(dāng)訂購(gòu) 122 份報(bào)紙時(shí)，可獲得最大利潤(rùn)份報(bào)紙時(shí)，可獲得最大利潤(rùn) 17.3584 元元第六講第六講 蒙特卡羅方法蒙特卡羅方法與計(jì)算機(jī)模擬與計(jì)算機(jī)模擬11連續(xù)系統(tǒng)的計(jì)算機(jī)模擬實(shí)例一范例范例 某軍導(dǎo)彈基地發(fā)現(xiàn)正北方向某

14、軍導(dǎo)彈基地發(fā)現(xiàn)正北方向120km處有一艘處有一艘敵艦以敵艦以90km/h的速度向正東方向行駛，該基地即的速度向正東方向行駛，該基地即刻發(fā)射導(dǎo)彈進(jìn)行攔擊，導(dǎo)彈速率刻發(fā)射導(dǎo)彈進(jìn)行攔擊，導(dǎo)彈速率450km/h，制導(dǎo)系，制導(dǎo)系統(tǒng)確保在任一時(shí)刻導(dǎo)彈都能對(duì)準(zhǔn)敵艦統(tǒng)確保在任一時(shí)刻導(dǎo)彈都能對(duì)準(zhǔn)敵艦問(wèn)題問(wèn)題1 試問(wèn)導(dǎo)彈何時(shí)何處擊中敵艦試問(wèn)導(dǎo)彈何時(shí)何處擊中敵艦o(x,y)120-y90t-x22120tan90()()450dyydxtxdydxdtdtsyms x y t dydt dxdt; solve(dydt/dxdt=(120-y)/(90*t-x),dydt2+dxdt2=4502,dydt,dxdt

15、);ans=ans.dxdt,ans.dydt; dxdt=ans(2,1); dydt=ans(2,2); pretty(dxdt); pretty(dydt);第六講第六講 蒙特卡羅方法蒙特卡羅方法與計(jì)算機(jī)模擬與計(jì)算機(jī)模擬12導(dǎo)彈追擊問(wèn)題-問(wèn)題1-11 將隨等距時(shí)間將隨等距時(shí)間連續(xù)變化的狀態(tài)變量軌跡連續(xù)變化的狀態(tài)變量軌跡x(t),y(t)用用歐拉法歐拉法離散化：離散化：12212290450(90)(120)120450(90)(120)kkkkkkkkkkkxxxkxyyyykxy1111(,)()kkkkkkkkkkkkkkyyf tyttdxxxdtdyyydt向前歐拉法第六講第六講

16、 蒙特卡羅方法蒙特卡羅方法與計(jì)算機(jī)模擬與計(jì)算機(jī)模擬13導(dǎo)彈追擊問(wèn)題-問(wèn)題1-22 編寫(xiě)程序模擬導(dǎo)彈攔擊敵艦過(guò)程編寫(xiě)程序模擬導(dǎo)彈攔擊敵艦過(guò)程x1=0; y1=0; x2=0; y2=120; t=0.001;v1=450; v2=90; dis=120;axis(0,40,0,140); grid on;set(gca,nextplot,add); for k=1:1000 x1=x1+v1*t*(v2*k*t-x1)/sqrt(v2*k*t-x1)2+(dis-y1)2);y1=y1+v1*t*(dis-y1)/sqrt(v2*k*t-x1)2+(dis-y1)2);x2=x2+v2*t;y2

17、=y2;plot(x1,y1,ro,x2,y2,bs); frame(k)=getframe;if sqrt(x1-x2)2+(y1-y2)2)=0.1,break;endendT=k*t,x1,y1%微分方程求解和計(jì)算機(jī)模擬過(guò)程已整合進(jìn)微分方程求解和計(jì)算機(jī)模擬過(guò)程已整合進(jìn)daodan1.m第六講第六講 蒙特卡羅方法蒙特卡羅方法與計(jì)算機(jī)模擬與計(jì)算機(jī)模擬14導(dǎo)彈追擊問(wèn)題-問(wèn)題1-33 導(dǎo)彈攔擊敵艦過(guò)程模擬動(dòng)畫(huà)導(dǎo)彈攔擊敵艦過(guò)程模擬動(dòng)畫(huà)(AVI視頻視頻)warning offmovie2avi(frame,daodan1,compression,Indeo5,quality,100,fps,6);

18、%將模擬過(guò)程編碼成將模擬過(guò)程編碼成AVI視頻視頻導(dǎo)彈于導(dǎo)彈于0.2770小時(shí)在小時(shí)在(25.0018千米千米, 119.9306千米千米)擊中敵艦擊中敵艦第六講第六講 蒙特卡羅方法蒙特卡羅方法與計(jì)算機(jī)模擬與計(jì)算機(jī)模擬15導(dǎo)彈追擊問(wèn)題-問(wèn)題2-1問(wèn)題問(wèn)題2 如果敵艦即刻發(fā)現(xiàn)導(dǎo)彈，并以垂直導(dǎo)彈方如果敵艦即刻發(fā)現(xiàn)導(dǎo)彈，并以垂直導(dǎo)彈方向向135km/h速度逃逸，試問(wèn)導(dǎo)彈何時(shí)何處擊中敵艦速度逃逸，試問(wèn)導(dǎo)彈何時(shí)何處擊中敵艦1 建立模型建立模型12112122122122112222tantan()()450()()135dyyydxxxdxyydyxxdydxdtdtdydxdtdto(x1,y1)x2-x1y2-y1(x2,y2)第六講第六講 蒙特卡羅方法蒙特卡羅方法與計(jì)算機(jī)模擬與計(jì)算機(jī)模擬16導(dǎo)彈追擊問(wèn)題-問(wèn)題2-22 求解并模擬模型求解并模擬模型 (下為片段下為片段,完整見(jiàn)完整見(jiàn)daodan2.m)syms x1 y1 x2 y2 t dy1dt dx1dt dy2dt dx2dt; solve(dy1dt/dx1dt=(y2-y1)/(x2-x1),dy2dt/dx2dt=(x2-x1)/(y2-y1),dy1dt2+dx1dt2=45