蒙特卡羅方法與計(jì)算機(jī)模擬資料

1、蒙特卡羅與計(jì)算機(jī)模擬,內(nèi)容：計(jì)算機(jī)模擬(或稱仿真)是一種廣義數(shù)值計(jì)算 方法，適合解決一些規(guī)模大、難以解析化 以及受隨機(jī)因素影響的不確定數(shù)學(xué)模型 目的：了解蒙特卡羅方法的基本思想，掌握利用 Matlab對離散/連續(xù)系統(tǒng)進(jìn)行模擬的方法 要求：掌握Matlab隨機(jī)數(shù)函數(shù)，處理應(yīng)用問題 了解蒙特卡羅方法的起源和基本思想 掌握離散系統(tǒng)和連續(xù)系統(tǒng)計(jì)算機(jī)模擬實(shí)例 掌握隨機(jī)函數(shù) rand unifrnd normrnd exprnd 了解Matlab仿真模塊 Simulink,蒙特卡羅方法的起源和基本思想,蒙特卡羅方法(Monte Carlo method)，或稱計(jì)算機(jī)隨機(jī)模擬方法，是一種基于“隨機(jī)數(shù)”的計(jì)算

2、方法。源于美國在第二次世界大戰(zhàn)研制原子彈的“曼哈頓計(jì)劃”，該計(jì)劃的主持人之一數(shù)學(xué)家馮諾伊曼用馳名世界的賭城摩納哥的Monte Carlo來命名這種方法，為它蒙上了一層神秘色彩。 蒙特卡羅方法的基本思想很早以前就被人們所發(fā)現(xiàn)和利用。早在17世紀(jì)，人們就知道用事件發(fā)生的“頻率”來決定事件的“概率”。19世紀(jì)人們用蒲豐投針的方法來計(jì)算圓周率，上世紀(jì)40年代電子計(jì)算機(jī)的出現(xiàn)，特別是近年來高速電子計(jì)算機(jī)的出現(xiàn)，使得用數(shù)學(xué)方法在計(jì)算機(jī)上大量、快速地模擬這樣的試驗(yàn)成為可能。,蒲豐投針實(shí)驗(yàn)近似計(jì)算圓周率,蒲豐投針實(shí)驗(yàn)： 法國科學(xué)家蒲豐(Buffon)在1777年提 出的蒲豐投針實(shí)驗(yàn)是早期幾何概率一個(gè) 非常著名

3、的例子。蒲豐投針實(shí)驗(yàn)的重要 性并非是為了求得比其它方法更精確的 值，而是它開創(chuàng)了使用隨機(jī)數(shù)處理確定性數(shù)學(xué)問題的先河，是用偶然性方法去解決確定性計(jì)算的前導(dǎo)，由此可以領(lǐng)略到從“概率土壤”上開出的一朵瑰麗的鮮花蒙特卡羅方法(MC) 蒲豐投針實(shí)驗(yàn)可歸結(jié)為下面的數(shù)學(xué)問題：平面上畫有距離為a的一些平行線，向平面上任意投一根長為l(la)的針，假設(shè)針落在任意位置的可能性相同，試求針與平行線相交的概率P(從而求),蒲豐投針實(shí)驗(yàn)近似計(jì)算圓周率,蒲豐投針實(shí)驗(yàn)： 如右圖所示，以M 表示針落下后的中點(diǎn)， 以x表示M到最近一條 平行線的距離，以表示針與此線的交角： 針落地的所有可能結(jié)果滿足： 其樣本空間視作矩形區(qū)域,

4、面積是: 針與平行線相交的條件： 它是樣本空間子集A，面積是： syms l phi; int(l/2*sin(phi),phi,0,pi); %ans=l 因此，針與平行線相交的概率為： 從而有： 特別當(dāng) 時(shí),蒲豐投針實(shí)驗(yàn)近似計(jì)算圓周率,蒲豐投針實(shí)驗(yàn)的計(jì)算機(jī)模擬： format long; %設(shè)置15位顯示精度 a=1; l=0.6; %兩平行線間的寬度和針長 figure; axis(0,pi,0,a/2); %初始化繪圖板 set(gca,nextplot,add); %初始化繪圖方式為疊加 counter=0; n=2010; %初始化計(jì)數(shù)器和設(shè)定投針次數(shù) x=unifrnd(0,a/

5、2,1,n); phi=unifrnd(0,pi,1,n); %樣本空間 for i=1:n if x(i)l*sin(phi(i)/2 %滿足此條件表示針與線的相交 plot(phi(i),x(i),r.); counter=counter+1; %統(tǒng)計(jì)針與線相交的次數(shù) frame(counter)=getframe; %描點(diǎn)并取幀 end end fren=counter/n; pihat=2*l/(a*fren) %用頻率近似計(jì)算 figure(2) movie(frame,1) %播放幀動畫1次,蒲豐投針實(shí)驗(yàn)近似計(jì)算圓周率,蒲豐投針實(shí)驗(yàn)的計(jì)算機(jī)模擬：,意大利數(shù)學(xué)家拉澤里尼得到了準(zhǔn)確到6

6、位小數(shù)的值，不過他的實(shí)驗(yàn)因?yàn)樘珳?zhǔn)確而受到了質(zhì)疑,蒲豐投針實(shí)驗(yàn)計(jì)算圓周率,蒙特卡羅投點(diǎn)法是蒲豐投針實(shí)驗(yàn)的推廣： 在一個(gè)邊長為a的正方形內(nèi)隨機(jī)投點(diǎn)，該點(diǎn)落在此正方形的內(nèi)切圓中的概率應(yīng)為該內(nèi)切圓與正方形的面積比值，即 n=10000; a=2; m=0; for i=1:n x=rand(1)*a; y=rand(1)*a; if ( (x-a/2)2+(y-a/2)2 = (a/2)2 ) m=m+1; end end disp(投點(diǎn)法近似計(jì)算的為: ,num2str(4*m/n);,常見分布的隨機(jī)數(shù)產(chǎn)生語句,蒙特卡羅方法的關(guān)鍵步驟在于隨機(jī)數(shù)的產(chǎn)生，計(jì)算機(jī)產(chǎn)生的隨機(jī)數(shù)都不是真正的隨機(jī)數(shù)(由算法確定

7、的緣故)，如果偽隨機(jī)數(shù)能夠通過一系列統(tǒng)計(jì)檢驗(yàn)，我們也可以將其當(dāng)作真正的隨機(jī)數(shù)使用：,MATLAB隨機(jī)數(shù)的“重置”問題,Matlab的隨機(jī)數(shù)是偽隨機(jī)數(shù)，但在一定的信度之下可以看作真正的隨機(jī)數(shù)。問題是rand函數(shù)產(chǎn)生的隨機(jī)數(shù)從一個(gè)隨機(jī)數(shù)序列中取出來，而每次啟動Matlab時(shí)，rand的狀態(tài)都會被重置(相當(dāng)于把序列的指針移到了隨機(jī)數(shù)序列的開始)，換言之第一次啟動Matlab調(diào)用的第n次rand函數(shù)與下一次啟動調(diào)用的第n個(gè)rand函數(shù)產(chǎn)生相同的數(shù)值。 如果想打亂這種狀態(tài)，可以為rand指定一個(gè)與當(dāng)前時(shí)間相關(guān)的初始狀態(tài)，而不用默認(rèn)狀態(tài)： rand(state,sum(100*clock); 或者 ran

8、d(state,sum(100*clock)*rand); 若使每次運(yùn)行程序產(chǎn)生的值是相同的 rand(seed,0.1); rand(1),非常見分布的隨機(jī)數(shù)的產(chǎn)生,1.逆變換方法,非常見分布的隨機(jī)數(shù)的產(chǎn)生,舉例說明： 1 連續(xù)型隨機(jī)變量(以指數(shù)分布為例)：,syms t x lambda; Fx=int(lambda*exp(-lambda*t),t,0,x) %分布函數(shù) syms r; Fxinv=finverse(Fx,x); %求反函數(shù) Fxinv=subs(Fxinv,x,r) %替換反函數(shù)變量x為r Fxinv=inline(Fxinv) x=Fxinv(3,rand) %產(chǎn)生參

9、數(shù) lambda=3 指數(shù)分布的隨機(jī)數(shù) %指數(shù)分布隨機(jī)數(shù)產(chǎn)生函數(shù)已經(jīng)提供 exprnd(1/3,1,1),非常見分布的隨機(jī)數(shù)的產(chǎn)生,2 離散型隨機(jī)變量(以離散分布為例)：,x=2,4,6,8; px=0.1,0.4,0.3,0.2; %以下為程序片段 Fx=0; for n=1:length(px),Fx=Fx,sum(px(1:n);end r=rand; index=find(rFx); x(index(1)-1) %已編寫通用離散分布隨機(jī)數(shù)產(chǎn)生程序 scatrnd(x,px,n),2.Acceptance-Rejection 方法 最早由 Von Neumann提出，現(xiàn)在已經(jīng) 廣泛應(yīng)用于

10、各種隨機(jī)數(shù)的生成。,離散系統(tǒng)的計(jì)算機(jī)模擬實(shí)例,范例 報(bào)童售報(bào)問題:一報(bào)童每天清晨從郵局訂購報(bào)紙后零售，每份報(bào)紙進(jìn)價(jià)0.35元，售價(jià)0.5元，郵局要求最低訂購數(shù)量為60份，根據(jù)過去經(jīng)驗(yàn)一個(gè)報(bào)童一天平均售出報(bào)紙120份(且滿足泊松分布)，未售出的報(bào)紙只要沒有破損可退給郵局，試求報(bào)童每天清晨訂購多少份報(bào)紙可獲最大利潤？ 1 數(shù)學(xué)建模,Pin=0.35元 Pout=0.5元 Nin60,200區(qū)間 Nout=poissrnd(120)隨機(jī)數(shù),離散系統(tǒng)的計(jì)算機(jī)模擬實(shí)例,2 計(jì)算機(jī)模擬 當(dāng)訂購 122 份報(bào)紙時(shí)，可獲得最大利潤 17.3584 元,離散系統(tǒng)的計(jì)算機(jī)模擬實(shí)例,范例 海港系統(tǒng)的卸載貨物問題,離

11、散系統(tǒng)的計(jì)算機(jī)模擬實(shí)例,范例 海港系統(tǒng)的卸載貨物問題1 程序片段(船只到港時(shí)間均勻分布,船只卸貨時(shí)間均勻分布) ShipBetweenTime(1)=unifrnd(15,145,1,1); %船只到港間隔時(shí)間隨機(jī)化(均勻分布) ShipUnloadTime(1)=unifrnd(45,90,1,1); %船只卸貨時(shí)間隨機(jī)化(均勻分布) 通用程序haibor.m可實(shí)現(xiàn)多次模擬，并將結(jié)果保存到H.txt delete H.txt %清除歷史數(shù)據(jù) harbor(100,15,145,45,90) load H.txt;Hmean=mean(H); %導(dǎo)入H并按列取平均值,離散系統(tǒng)的計(jì)算機(jī)模擬實(shí)例,

12、范例 海港系統(tǒng)的卸載貨物問題2 程序片段(船只到港時(shí)間指數(shù)分布,船只卸貨時(shí)間均勻分布) ShipBetweenTime(1)=exprnd(60,1,1); %船只到港間隔時(shí)間隨機(jī)化(指數(shù)分布) ShipUnloadTime(1)=unifrnd(45,90,1,1); %船只卸貨時(shí)間隨機(jī)化(均勻分布) 通用程序haibor2.m可實(shí)現(xiàn)多次模擬，結(jié)果保存到H2.txt delete H2.txt %清除歷史數(shù)據(jù) harbor2(100,60,45,90) load H2.txt;Hmean2=mean(H2); %導(dǎo)入H2并按列取平均值,離散系統(tǒng)的計(jì)算機(jī)模擬實(shí)例,范例 海港系統(tǒng)的卸載貨物問題3

13、 程序片段(船只到港時(shí)間離散分布,船只卸貨時(shí)間離散分布) 1 編寫船只到港間隔離散累積分布函數(shù)并作階梯圖： xs=15:10:145; for i=1:length(xs)-1,x(i)=(xs(i)+xs(i+1)/2;end px=0.009,0.029,0.035,0.051,0.090,0.161,0.200,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);end plot(10,x,Fx,-rs); hold on; stairs(0,x-5,145,Fx,1); set(

14、gca,xtick,0:5:145); set(gca,xgrid,on); axis tight;,離散系統(tǒng)的計(jì)算機(jī)模擬實(shí)例,范例 海港系統(tǒng)的卸載貨物問題3 程序片段(船只到港時(shí)間離散分布,船只卸貨時(shí)間離散分布) 2 編寫船只到港間隔離散累積分布反函數(shù)并作線性插值： Fxi=0:0.001:1-eps; xi=interp1(Fx,0,x,Fxi,linear); r=rand(1,n); rnd=; for i=1:n index=find(r(i)=Fxi); if xs(1)=xi(index(1)-1)=xs(length(xs) rnd=rnd,xi(index(1)-1); en

15、d end %以上程序已編寫通用M函數(shù)文件 harborrnd(xs,px,n) %即給出n個(gè)滿足離散分布(x,px)的船只到港間隔隨機(jī)數(shù),離散系統(tǒng)的計(jì)算機(jī)模擬實(shí)例,范例 海港系統(tǒng)的卸載貨物問題3 程序片段(船只到港時(shí)間離散分布,船只卸貨時(shí)間離散分布) 3 編寫船只卸貨時(shí)間離散累積分布函數(shù)并作階梯圖： xs=45:5:90; for i=1:length(xs)-1,x(i)=(xs(i)+xs(i+1)/2;end px=0.017,0.045,0.095,0.086,0.130,0.185,0.208,0.143,0.091; Fx=0; for i=1:length(px),Fx=Fx,

16、sum(px(1:i);end plot(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;,離散系統(tǒng)的計(jì)算機(jī)模擬實(shí)例,范例 海港系統(tǒng)的卸載貨物問題3 程序片段(船只到港時(shí)間離散分布,船只卸貨時(shí)間離散分布) 4 編寫船只卸貨時(shí)間離散累積分布反函數(shù)并作線性插值： Fxi=0:0.001:1-eps; xi=interp1(Fx,0,x,Fxi,linear); r=rand(1,n); rnd=; for i=1:n index=find

17、(r(i)=Fxi); if xs(1)=xi(index(1)-1)=xs(length(xs) rnd=rnd,xi(index(1)-1); end end %以上程序已編寫通用M函數(shù)文件 harborrnd(xs,px,n) %即給出n個(gè)滿足離散分布(x,px)的船只卸貨時(shí)間隨機(jī)數(shù),離散系統(tǒng)的計(jì)算機(jī)模擬實(shí)例,范例 海港系統(tǒng)的卸載貨物問題3 程序片段(船只到港時(shí)間指數(shù)分布,船只卸貨時(shí)間均勻分布) 5 模擬船只到港間隔 / 卸貨時(shí)間均為離散分布的海港系統(tǒng) ShipBetweenTime(1)=harborrnd(sbtxs,sbtpx,1); %船只到港間隔時(shí)間隨機(jī)化(離散分布) Ship

18、UnloadTime(1)=harborrnd(sutxs,sutpx,1); %船只卸貨時(shí)間隨機(jī)化(離散分布) 通用程序haibor3.m可實(shí)現(xiàn)多次模擬，結(jié)果保存到H3.txt delete H3.txt %清除歷史數(shù)據(jù) load harbor.mat %載入數(shù)據(jù) harbor3(100,sbtxs,sbtpx,sutxs,sutpx) load H3.txt;Hmean3=mean(H3); %導(dǎo)入H3并按列取平均值,連續(xù)系統(tǒng)的計(jì)算機(jī)模擬實(shí)例,范例 某軍導(dǎo)彈基地發(fā)現(xiàn)正北方向120km處有一艘敵艦以90km/h的速度向正東方向行駛，該基地即刻發(fā)射導(dǎo)彈進(jìn)行攔擊，導(dǎo)彈速率450km/h，制導(dǎo)系統(tǒng)

19、確保在任一時(shí)刻導(dǎo)彈都能對準(zhǔn)敵艦 (問題1) 試問導(dǎo)彈何時(shí)何處擊中敵艦,syms x y t dydt dxdt; solve(dydt/dxdt=(120-y)/(90*t-x),dydt2+dxdt2=4502,dydt,dxdt); ans=ans.dxdt,ans.dydt; dxdt=ans(2,1); dydt=ans(2,2); pretty(dxdt); pretty(dydt);,連續(xù)系統(tǒng)的計(jì)算機(jī)模擬實(shí)例,1 將隨等距時(shí)間連續(xù)變化的狀態(tài)變量軌跡x(t),y(t)用歐拉法離散化：,連續(xù)系統(tǒng)的計(jì)算機(jī)模擬實(shí)例,2 編寫程序模擬導(dǎo)彈攔擊敵艦過程 x1=0; y1=0; x2=0; y2

20、=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=y2; plot(x1,y1,ro,x2,y2,bs); frame(k)=getframe; if sqrt(x1-x2)2+(y1-y2)2)=0.1,brea

21、k;end end T=k*t,x1,y1 %微分方程求解和計(jì)算機(jī)模擬過程已整合進(jìn)daodan1.m,連續(xù)系統(tǒng)的計(jì)算機(jī)模擬實(shí)例,3 導(dǎo)彈攔擊敵艦過程模擬動畫(AVI視頻) warning off movie2avi(frame,daodan1,compression,Indeo5,quality,100,fps,6); %將模擬過程編碼成AVI視頻 導(dǎo)彈于0.2770小時(shí)在(25.0018千米, 119.9306千米)擊中敵艦,連續(xù)系統(tǒng)的計(jì)算機(jī)模擬實(shí)例,(問題2) 如果敵艦即刻發(fā)現(xiàn)導(dǎo)彈，并以垂直導(dǎo)彈方向135km/h速度逃逸，試問導(dǎo)彈何時(shí)何處擊中敵艦 1 建立模型,連續(xù)系統(tǒng)的計(jì)算機(jī)模擬實(shí)例,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=4502,dy2dt2+dx2dt2=1352,dy1dt,dx1dt,dy2dt,dx2dt); for k=1:1000 x1=x1+dx1dt(x1,x2,y1,y2)*t; y1=y1+dy1dt(x1,x2,y1,y2)*t; x2=x2+dx2dt(x1,x2,y1,y2)*t; y2=y

23、2-dy2dt(x1,x2,y1,y2)*t; plot(x1,y1,ro,x2,y2,bs); frame(k)=getframe; if sqrt(x1-x2)2+(y1-y2)2)=0.3,break;end end,連續(xù)系統(tǒng)的計(jì)算機(jī)模擬實(shí)例,3 過程動畫 (AVI視頻) warning off movie2avi(frame,daodan2,compression,Indeo5,quality,100,fps,6); %將模擬過程編碼成AVI視頻 導(dǎo)彈于0.2660小時(shí)在(32.8916千米, 109.8963千米)擊中敵艦,排隊(duì)問題隨機(jī)模擬,排隊(duì)論主要研究隨機(jī)服務(wù)系統(tǒng)的工作過程。,在

24、排隊(duì)系統(tǒng)中，服務(wù)對象的到達(dá)時(shí)間和服務(wù)時(shí)間都是隨機(jī)的。排隊(duì)論通過對隨機(jī)服務(wù)現(xiàn)象的統(tǒng)計(jì)研究，找出反映這些隨機(jī)現(xiàn)象平均特性的規(guī)律指標(biāo)，如排隊(duì)的長度、等待的時(shí)間及服務(wù)利用率。,1 系統(tǒng)的假設(shè)： （1） 顧客源是無窮的； （2） 排隊(duì)的長度沒有限制； （ 3）到達(dá)系統(tǒng)的顧客按先后順序依次進(jìn)入服務(wù)，“先到先服務(wù)”。,在某商店有一個(gè)售貨員，顧客陸續(xù)來到，售貨員逐個(gè)地接待顧客當(dāng)?shù)絹淼念櫩洼^多時(shí)，一部分顧客便須排隊(duì)等待，被接待后的顧客便離開商店設(shè)： 1顧客到來間隔時(shí)間服從參數(shù)為0.1的指數(shù)分布 對顧客的服務(wù)時(shí)間服從，上的均勻分布 排隊(duì)按先到先服務(wù)規(guī)則，隊(duì)長無限制,假定一個(gè)工作日為8小時(shí)，時(shí)間以分鐘為單位。 1模擬一個(gè)工作日內(nèi)完成服務(wù)的顧客個(gè)數(shù)及顧客平均等待時(shí)間t 2模擬100個(gè)工作日，求出平均每日完成服務(wù)的顧客個(gè)數(shù)及每日顧客的平均等待時(shí)間。,單服務(wù)員的排隊(duì)模型模擬,w：總等待時(shí)間； ci：第i個(gè)顧客的到達(dá)時(shí)刻； bi：第i個(gè)顧客開始服務(wù)時(shí)刻； ei：第i個(gè)顧客服務(wù)結(jié)束時(shí)刻； xi：第i-1個(gè)顧客與第i個(gè)顧客之間到達(dá)的間隔時(shí)間； yi：對第i個(gè)顧客的服務(wù)時(shí)間。,符號說明,c1,b1,c3,c4,c5,c2,e1,b2,e2,b3,e3,b4,e4,b5,ci=ci-1+ xi ei=bi+yi bi=max(ci,ei-1),