蒙特卡羅模擬及例課件

1、 例例1 用用MC 模擬求圓周率模擬求圓周率的估計(jì)值的估計(jì)值. 設(shè)二維隨機(jī)變量設(shè)二維隨機(jī)變量(X, Y)在正方形內(nèi)服從均勻分布在正方形內(nèi)服從均勻分布. . (X, Y)落在圓內(nèi)的落在圓內(nèi)的概率為概率為: : 110PX2+Y21=4計(jì)算機(jī)上做計(jì)算機(jī)上做n次擲點(diǎn)試驗(yàn)：次擲點(diǎn)試驗(yàn)： 產(chǎn)生產(chǎn)生n 對(duì)二維隨機(jī)點(diǎn)對(duì)二維隨機(jī)點(diǎn)(xi，yi) ，i1 ,2, , n .若有若有k 個(gè)點(diǎn)落在個(gè)點(diǎn)落在l4圓內(nèi)圓內(nèi) 其中其中，xi 和和yi 是是RND 隨機(jī)數(shù)對(duì)隨機(jī)數(shù)對(duì). . 相當(dāng)于第相當(dāng)于第i個(gè)隨機(jī)點(diǎn)個(gè)隨機(jī)點(diǎn)落在落在1/4圓內(nèi)圓內(nèi).檢查每對(duì)數(shù)是否滿足檢查每對(duì)數(shù)是否滿足: 隨機(jī)事件隨機(jī)事件“點(diǎn)落入點(diǎn)落入1/4圓內(nèi)

2、圓內(nèi)”的的頻率為頻率為 k/n 根據(jù)概率論中的大數(shù)定律，事件發(fā)生的頻率依根據(jù)概率論中的大數(shù)定律，事件發(fā)生的頻率依概率收斂于事件發(fā)生的概率概率收斂于事件發(fā)生的概率p，即有即有1lim pPnkn122 iiyx得圓周率得圓周率的估計(jì)值為的估計(jì)值為=4k/n且當(dāng)且當(dāng)試驗(yàn)次數(shù)足夠大試驗(yàn)次數(shù)足夠大時(shí)，其精度也隨之提高時(shí)，其精度也隨之提高。 分析：分析：實(shí)際上概率值為實(shí)際上概率值為41102 dxx恰為恰為1/4圓的面積圓的面積 頻率法：頻率法： 利用隨機(jī)變量落進(jìn)指定區(qū)域內(nèi)的頻率來(lái)計(jì)算定積分。利用隨機(jī)變量落進(jìn)指定區(qū)域內(nèi)的頻率來(lái)計(jì)算定積分。 平均值法：平均值法： 利用隨機(jī)變量的平均值（數(shù)學(xué)期望利用隨機(jī)變量

3、的平均值（數(shù)學(xué)期望 )來(lái)計(jì)算定積分。來(lái)計(jì)算定積分。 badxxfI)(平均值法平均值法的算法如下：的算法如下：(1)(1) 產(chǎn)生產(chǎn)生RND 隨機(jī)數(shù)：隨機(jī)數(shù)：r1，r2，rn；（2）令令 ui=a(ba)ri，i=1，2，n； （3）計(jì)算計(jì)算 作為作為I的估計(jì)值的估計(jì)值。 niiufnab1)(原理分析：原理分析： 設(shè)隨機(jī)變量設(shè)隨機(jī)變量1，2，n相互獨(dú)立相互獨(dú)立，且且iU(0,1)f(i)，i=1，2，n 相互獨(dú)立同分布相互獨(dú)立同分布abIdxxfabfEbai )(1)( 由由(強(qiáng)強(qiáng))大數(shù)定律知大數(shù)定律知.)(1lim1eaabIfnniin 即以概率為即以概率為1 成立。成立。 當(dāng)當(dāng)n 足夠

4、大時(shí)，得近似公式：足夠大時(shí)，得近似公式： niinbafabdxxfI11)()()(注：注： 平均值法本質(zhì)上是用樣本平均值作為平均值法本質(zhì)上是用樣本平均值作為總體教學(xué)期望的估計(jì)?？傮w教學(xué)期望的估計(jì)。 MC 模擬是一種試驗(yàn)近似方法模擬是一種試驗(yàn)近似方法 , 試驗(yàn)次數(shù)如何確定？試驗(yàn)次數(shù)如何確定？希望：模擬次數(shù)較少、希望：模擬次數(shù)較少、 模擬精度較高模擬精度較高頻率法的討論頻率法的討論 頻率法是用事件頻率法是用事件A出現(xiàn)的頻率作為概率出現(xiàn)的頻率作為概率p 的估計(jì)的估計(jì): nkpn 問(wèn)題：?jiǎn)栴}：試驗(yàn)次數(shù)試驗(yàn)次數(shù) n 多大時(shí)，對(duì)給定的置信度多大時(shí)，對(duì)給定的置信度1（01），），估計(jì)精度達(dá)到估計(jì)精度達(dá)到

5、. 即問(wèn)：取多大的即問(wèn)：取多大的n 使使 1pnkPppPn成立？成立？答案：答案：22)1( zppn其中，其中， z是正態(tài)分布的臨界值。是正態(tài)分布的臨界值。證明證明平均值法平均值法 在給定在給定和和下所需的試驗(yàn)次數(shù)下所需的試驗(yàn)次數(shù)的估計(jì)式為的估計(jì)式為 01202222)(11/niixxnSSzn試驗(yàn)次數(shù)估計(jì)式的分析試驗(yàn)次數(shù)估計(jì)式的分析 01202222)(11/niixxnSSzn為估計(jì)概率為估計(jì)概率p做模擬，卻又需做模擬，卻又需要用要用p去估計(jì)模擬次數(shù)去估計(jì)模擬次數(shù)n. 如何計(jì)算如何計(jì)算S2 ？ 解決方法解決方法：先做：先做n0 次模擬（稱為學(xué)習(xí)樣本次模擬（稱為學(xué)習(xí)樣本) )，根據(jù)學(xué)習(xí)

6、樣本，根據(jù)學(xué)習(xí)樣本（1）先求出先求出p的估計(jì)，再估計(jì)模擬次數(shù)的估計(jì)，再估計(jì)模擬次數(shù)n : : 22)1( zppn1.22)1( zppn（2）計(jì)算出的樣本方差計(jì)算出的樣本方差S2 , ,用來(lái)估計(jì)用來(lái)估計(jì)n. 2. M C模擬的估計(jì)精度模擬的估計(jì)精度與試驗(yàn)次數(shù)與試驗(yàn)次數(shù)n的平的平方根成反比方根成反比. 核反應(yīng)堆屏蔽層是用一定厚度的鉛包圍反應(yīng)核反應(yīng)堆屏蔽層是用一定厚度的鉛包圍反應(yīng)堆，用以阻擋或減弱反應(yīng)堆發(fā)出的各種射線。在堆，用以阻擋或減弱反應(yīng)堆發(fā)出的各種射線。在各種射線中，中子對(duì)人體傷害極大，因此，在屏各種射線中，中子對(duì)人體傷害極大，因此，在屏蔽層的設(shè)計(jì)中，了解中子穿透屏蔽層的概率，對(duì)蔽層的設(shè)計(jì)

7、中，了解中子穿透屏蔽層的概率，對(duì)反應(yīng)堆的安全運(yùn)行至關(guān)重要反應(yīng)堆的安全運(yùn)行至關(guān)重要。例例2 核反應(yīng)堆屏蔽層設(shè)計(jì)問(wèn)題核反應(yīng)堆屏蔽層設(shè)計(jì)問(wèn)題1）問(wèn)題背景）問(wèn)題背景D返回返回穿透穿透吸收吸收三種狀態(tài)三種狀態(tài)闡述中子的運(yùn)動(dòng)，闡述中子的運(yùn)動(dòng)，為模擬做理論準(zhǔn)備為模擬做理論準(zhǔn)備2）簡(jiǎn)化假設(shè)：）簡(jiǎn)化假設(shè)： *1 假定屏蔽層平行板厚度為假定屏蔽層平行板厚度為D=3d，其中其中d 為兩次碰撞之間中子的平均游動(dòng)距離；為兩次碰撞之間中子的平均游動(dòng)距離； *2 假設(shè)在第假設(shè)在第10 次碰撞以后，中子速度下降到次碰撞以后，中子速度下降到為某一很小數(shù)值而終止運(yùn)動(dòng)（被引收）為某一很小數(shù)值而終止運(yùn)動(dòng)（被引收）. . *3 假定中

8、子在屏蔽層內(nèi)相繼兩次碰撞之間游假定中子在屏蔽層內(nèi)相繼兩次碰撞之間游動(dòng)的距離服從指數(shù)分布動(dòng)的距離服從指數(shù)分布； *4 中子經(jīng)碰撞后的彈射角中子經(jīng)碰撞后的彈射角 U(0, 2).3）中子運(yùn)動(dòng)的數(shù)學(xué)描述）中子運(yùn)動(dòng)的數(shù)學(xué)描述引進(jìn)變量引進(jìn)變量： 彈射角彈射角i 第第i 次碰撞后中子的運(yùn)動(dòng)方向與次碰撞后中子的運(yùn)動(dòng)方向與x 軸正向的夾角軸正向的夾角. . xi 第第i 次碰撞后中子所處位置與屏蔽層內(nèi)壁的距離次碰撞后中子所處位置與屏蔽層內(nèi)壁的距離.DD0 xxi iR i 中子在第中子在第 i 次碰撞前后的游動(dòng)距離次碰撞前后的游動(dòng)距離.以上三個(gè)變量均為以上三個(gè)變量均為隨機(jī)變量隨機(jī)變量（1） xi0，中子返回反

9、應(yīng)堆中子返回反應(yīng)堆；（2） xiD，中子穿透屏蔽層；中子穿透屏蔽層； 經(jīng)過(guò)第經(jīng)過(guò)第i 次碰撞，中子在屏蔽層內(nèi)的位置是次碰撞，中子在屏蔽層內(nèi)的位置是 xi=xi1+Ricosi ，i=1，2，10 ， 中子在屏蔽層里隨機(jī)游動(dòng)，第中子在屏蔽層里隨機(jī)游動(dòng)，第 i 次碰撞以后，按照它的位置坐標(biāo)次碰撞以后，按照它的位置坐標(biāo) xi，可能有以下三可能有以下三種情況發(fā)生：種情況發(fā)生：4）模擬過(guò)程）模擬過(guò)程(1) 產(chǎn)生產(chǎn)生RND隨機(jī)數(shù)對(duì)隨機(jī)數(shù)對(duì)(ri, ui )；（3） 0 xiD，若若i10，中子中子在屏蔽層內(nèi)繼續(xù)運(yùn)動(dòng)，在屏蔽層內(nèi)繼續(xù)運(yùn)動(dòng)，若若 i=10，中子被屏蔽層吸收。中子被屏蔽層吸收。中子三狀態(tài)中子三狀

10、態(tài)判別準(zhǔn)則判別準(zhǔn)則(2) 將將(ri, ui )代入公式計(jì)算代入公式計(jì)算 ,2,lniiiiurdR (i=1,2,3，10) (3) 將將(Ri, i) 代入公式代入公式 xi = xi1+Ricosi ，i=1，2，10 (4) 判斷中子是否穿透屏蔽層判斷中子是否穿透屏蔽層.5） 模擬結(jié)果分析模擬結(jié)果分析 要求穿透屏蔽層的概率數(shù)量級(jí)為要求穿透屏蔽層的概率數(shù)量級(jí)為1061010，按假設(shè)條件得到一次按假設(shè)條件得到一次模擬結(jié)果模擬結(jié)果如下如下： 第第i 次中子的移動(dòng)距離和彈射角。次中子的移動(dòng)距離和彈射角。計(jì)算出第計(jì)算出第i 次碰撞中子與內(nèi)壁的距離次碰撞中子與內(nèi)壁的距離xi . 中子數(shù)中子數(shù)( (

11、個(gè)個(gè)) ) 穿透穿透(%) 吸收吸收(%) 返回返回( (%) 100100030005000 30.026.026.526.3 28.023.421.8 22.0 42.050.651.7 51.7 中子穿透屏蔽層的百分比超過(guò)了中子穿透屏蔽層的百分比超過(guò)了1/4，模擬結(jié)模擬結(jié)果表明屏蔽層厚度果表明屏蔽層厚度D=3d不合適不合適. 多厚的屏蔽層才能使穿透的概率多厚的屏蔽層才能使穿透的概率W106？問(wèn)問(wèn) 題題： 如何解決這個(gè)問(wèn)題？如何解決這個(gè)問(wèn)題？思路？思路？1.計(jì)算機(jī)收索法計(jì)算機(jī)收索法 2. 分析法分析法 12mDDD 設(shè)計(jì)屏蔽層的厚度：設(shè)計(jì)屏蔽層的厚度： x=mD 將屏蔽層視為將屏蔽層視為m層厚度均為層厚度均為D的平行板的平行板. . 增大屏蔽層的厚度，如增大屏蔽層的厚度，如D=6d、12d、24d、36d，交由計(jì)算機(jī)進(jìn)行模擬，并搜索到所交由計(jì)算機(jī)進(jìn)行模擬，并搜索到所求解求解. 由于碰撞的能量損失，中子穿過(guò)屏蔽層的平由于碰撞的能量損失，中子穿過(guò)屏蔽層的平均速度會(huì)逐層下降均速度會(huì)逐層下降. . 設(shè)設(shè)WD 是中子穿過(guò)厚度為是中子穿過(guò)厚度為D 屏蔽層的概率，則屏蔽層的概率，則穿過(guò)整個(gè)屏蔽層的概率穿過(guò)整個(gè)屏蔽層的概率W滿足滿足 mDWW 利用模擬結(jié)果利用模擬結(jié)果：當(dāng)當(dāng)D=3d，WD0.25，令令(WD)m106， 或或 (m