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

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

2、dTime(1)=unifrnd(45,90,1,1);%船只卸貨時間隨機化船只卸貨時間隨機化(均勻分布均勻分布)通用程序通用程序haibor1.m可實現(xiàn)多次模擬，結(jié)果保存到可實現(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并按列取平均值并按列取平均值第六講第六講 蒙特卡羅方法蒙特卡羅方法與計算機模擬與計算機模擬3海港系統(tǒng)的卸載貨物-問題2問題問題2 海港系統(tǒng)的卸載貨物海港系統(tǒng)的卸載貨物 (p124)程序片段程序片段(船只到港時間

3、船只到港時間指數(shù)指數(shù)分布分布,船只卸貨時間船只卸貨時間均勻分布均勻分布)ShipBetweenTime(1)=exprnd(60,1,1);%船只到港間隔時間隨機化船只到港間隔時間隨機化(指數(shù)分布指數(shù)分布)ShipUnloadTime(1)=unifrnd(45,90,1,1);%船只卸貨時間隨機化船只卸貨時間隨機化(均勻分布均勻分布)通用程序通用程序haibor2.m可實現(xiàn)多次模擬，結(jié)果保存到可實現(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并按列取平均值并按列取平均值第六講第六講 蒙特卡羅方法蒙特卡羅方法與計算機模擬與計算機模擬4海港系統(tǒng)的卸載貨物-問題3-1問題問題3 海港系統(tǒng)的卸載貨物海港系統(tǒng)的卸載貨物 (p124-129)程序片段程序片段(船只到港時間船只到港時間離散分布離散分布,船只卸貨時間船只卸貨時間離散分布離散分布)1 編寫編寫船只到港間隔船只到港間隔離散累積分布函數(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;第六講第六講 蒙特卡羅方法蒙特卡羅方法與計算機模擬與計算機模擬5海港系統(tǒng)的卸載貨物-問題3-2問題問題3 海港系統(tǒng)的卸載貨物海港系統(tǒng)的卸載貨物 (p124-129)程序片段程序片段(船只到港時間船只到港時間離散分布離散

6、分布,船只卸貨時間船只卸貨時間離散分布離散分布)2 編寫編寫船只到港間隔船只到港間隔離散累積分布反函數(shù)并作線性插值：離散累積分布反函數(shù)并作線性插值：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%以上程序已編寫通用以上程序已編寫通用M函數(shù)文件函數(shù)文件 harborrnd(xs,px,n)%即給出即給出n個滿足離

7、散分布個滿足離散分布(xs,px)的的船只到港間隔船只到港間隔隨機數(shù)隨機數(shù)第六講第六講 蒙特卡羅方法蒙特卡羅方法與計算機模擬與計算機模擬6海港系統(tǒng)的卸載貨物-問題3-3問題問題3 海港系統(tǒng)的卸載貨物海港系統(tǒng)的卸載貨物 (p124-129)程序片段程序片段(船只到港時間船只到港時間離散分布離散分布,船只卸貨時間船只卸貨時間離散分布離散分布)3 編寫編寫船只卸貨時間船只卸貨時間離散累積分布函數(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第六講第六講 蒙特卡羅方法蒙特卡羅方法與計算機模擬與計算機模擬7海港系

9、統(tǒng)的卸載貨物-問題3-4問題問題3 海港系統(tǒng)的卸載貨物海港系統(tǒng)的卸載貨物 (p124-129)程序片段程序片段(船只到港時間船只到港時間離散分布離散分布,船只卸貨時間船只卸貨時間離散分布離散分布)4 編寫編寫船只卸貨時間船只卸貨時間離散累積分布反函數(shù)并作線性插值：離散累積分布反函數(shù)并作線性插值：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%以上程序已編寫通用以上程序已編寫通用M函數(shù)文件函數(shù)文件 harborrnd(xs,px,n)%即給出即給出n個滿足離散分布個滿足離散分布(xs,px)的的船只卸貨時間船只卸貨時間隨機數(shù)隨機數(shù)第六講第六講 蒙特卡羅方法蒙特卡羅方法與計算機模擬與計算機模擬8海港系統(tǒng)的卸載貨物-問題3-5問題問題3 海港系統(tǒng)的卸載貨物海港系統(tǒng)的卸載貨物 (p124-129)程序片段程序片段(船只到港時間船只到港時間指數(shù)指數(shù)分布分布,船只卸貨時間船只卸貨時間均勻分布均勻分布)5 模擬模擬船只到港間隔船只到港間隔 / 卸貨時間卸貨時間均為離散分布的海港系統(tǒng)均為離散分布的海港系統(tǒng)ShipB

11、etweenTime(1)=harborrnd(sbtxs,sbtpx,1);%船只到港間隔時間隨機化船只到港間隔時間隨機化(離散分布離散分布)ShipUnloadTime(1)=harborrnd(sutxs,sutpx,1);%船只卸貨時間隨機化船只卸貨時間隨機化(離散分布離散分布)通用程序通用程序haibor3.m可實現(xiàn)多次模擬，結(jié)果保存到可實現(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); 第六講第六講 蒙特卡羅方法蒙特卡羅方法與計算機模擬與計算機模擬9離散系統(tǒng)的計算機模擬實例二范例范例 報童售報問題報童售報問題:一報童每天清晨從郵局訂購報一報童每天清晨從郵局訂購報紙后零售，每份報紙進價紙后零售，每份報紙進價0.35元，售價元，售價0.5元，郵元，郵局要求最低訂購數(shù)量為局要求最低訂購數(shù)量為60份，根據(jù)過去經(jīng)驗一個份，根據(jù)過去經(jīng)驗一個報童一天平均售出報紙報童一天平均售出報紙120份份(且滿足泊松分布且滿足泊松分布)，未售出的報紙只要沒有破損可退給郵局，試求報未售出的報紙只要沒有破損可退給郵局，試求報童每天清晨訂購多少份報紙可獲最大利潤？童每天清晨訂購多少份報紙

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

14、軍導(dǎo)彈基地發(fā)現(xiàn)正北方向120km處有一艘處有一艘敵艦以敵艦以90km/h的速度向正東方向行駛，該基地即的速度向正東方向行駛，該基地即刻發(fā)射導(dǎo)彈進行攔擊，導(dǎo)彈速率刻發(fā)射導(dǎo)彈進行攔擊，導(dǎo)彈速率450km/h，制導(dǎo)系，制導(dǎo)系統(tǒng)確保在任一時刻導(dǎo)彈都能對準(zhǔn)敵艦統(tǒng)確保在任一時刻導(dǎo)彈都能對準(zhǔn)敵艦問題問題1 試問導(dǎo)彈何時何處擊中敵艦試問導(dǎo)彈何時何處擊中敵艦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);第六講第六講 蒙特卡羅方法蒙特卡羅方法與計算機模擬與計算機模擬12導(dǎo)彈追擊問題-問題1-11 將隨等距時間將隨等距時間連續(xù)變化的狀態(tài)變量軌跡連續(xù)變化的狀態(tài)變量軌跡x(t),y(t)用用歐拉法歐拉法離散化：離散化：12212290450(90)(120)120450(90)(120)kkkkkkkkkkkxxxkxyyyykxy1111(,)()kkkkkkkkkkkkkkyyf tyttdxxxdtdyyydt向前歐拉法第六講第六講

16、 蒙特卡羅方法蒙特卡羅方法與計算機模擬與計算機模擬13導(dǎo)彈追擊問題-問題1-22 編寫程序模擬導(dǎo)彈攔擊敵艦過程編寫程序模擬導(dǎo)彈攔擊敵艦過程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%微分方程求解和計算機模擬過程已整合進微分方程求解和計算機模擬過程已整合進daodan1.m第六講第六講 蒙特卡羅方法蒙特卡羅方法與計算機模擬與計算機模擬14導(dǎo)彈追擊問題-問題1-33 導(dǎo)彈攔擊敵艦過程模擬動畫導(dǎo)彈攔擊敵艦過程模擬動畫(AVI視頻視頻)warning offmovie2avi(frame,daodan1,compression,Indeo5,quality,100,fps,6);

18、%將模擬過程編碼成將模擬過程編碼成AVI視頻視頻導(dǎo)彈于導(dǎo)彈于0.2770小時在小時在(25.0018千米千米, 119.9306千米千米)擊中敵艦擊中敵艦第六講第六講 蒙特卡羅方法蒙特卡羅方法與計算機模擬與計算機模擬15導(dǎo)彈追擊問題-問題2-1問題問題2 如果敵艦即刻發(fā)現(xiàn)導(dǎo)彈，并以垂直導(dǎo)彈方如果敵艦即刻發(fā)現(xiàn)導(dǎo)彈，并以垂直導(dǎo)彈方向向135km/h速度逃逸，試問導(dǎo)彈何時何處擊中敵艦速度逃逸，試問導(dǎo)彈何時何處擊中敵艦1 建立模型建立模型12112122122122112222tantan()()450()()135dyyydxxxdxyydyxxdydxdtdtdydxdtdto(x1,y1)x2-x1y2-y1(x2,y2)第六講第六講 蒙特卡羅方法蒙特卡羅方法與計算機模擬與計算機模擬16導(dǎo)彈追擊問題-問題2-22 求解并模擬模型求解并模擬模型 (下為片段下為片段,完整見完整見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