蒙特卡羅模擬

1、 蒙特卡羅蒙特卡羅( (MonteCarlo)模擬，又稱蒙特卡模擬，又稱蒙特卡羅方法、統(tǒng)計試驗法等羅方法、統(tǒng)計試驗法等. . MC模擬是靜態(tài)模擬，描述特定時間點上模擬是靜態(tài)模擬，描述特定時間點上的系統(tǒng)行為的系統(tǒng)行為.模擬過程中模擬過程中不出現(xiàn)時間不出現(xiàn)時間參數(shù)。參數(shù)?；舅枷牖舅枷? :把隨機(jī)事件把隨機(jī)事件（變量）的概率特征與（變量）的概率特征與數(shù)學(xué)分析的解聯(lián)系起來數(shù)學(xué)分析的解聯(lián)系起來. . 概率特征：概率特征：隨機(jī)事件的概率和隨機(jī)變量的隨機(jī)事件的概率和隨機(jī)變量的數(shù)學(xué)期望等數(shù)學(xué)期望等. . 用試驗方法確定用試驗方法確定 例例7.3.1 用用MC 模擬求圓周率模擬求圓周率的估計值的估計值. 1

2、10 設(shè)二維隨機(jī)變量設(shè)二維隨機(jī)變量(X, Y)在正方形內(nèi)在正方形內(nèi)服從均勻分布服從均勻分布. . (X, Y)落在圓內(nèi)的概率為落在圓內(nèi)的概率為: : 4122 YXP計算機(jī)上做計算機(jī)上做n 次擲點試驗：次擲點試驗： 產(chǎn)生產(chǎn)生n 對二維隨機(jī)點對二維隨機(jī)點(xi，yi) ，i1 ,2, , n .xi 和和yi 是是RND 隨機(jī)數(shù)對隨機(jī)數(shù)對. . 檢查每對隨機(jī)數(shù)是否滿足檢查每對隨機(jī)數(shù)是否滿足:122 iiyx相當(dāng)于第相當(dāng)于第i個隨機(jī)點落個隨機(jī)點落在在1/4圓內(nèi)圓內(nèi).若有若有k 個點落在個點落在l4圓內(nèi)圓內(nèi) 隨機(jī)事件隨機(jī)事件“點落入點落入1/4圓內(nèi)圓內(nèi)”的的頻率為頻率為 k/n 根據(jù)概率論中的大數(shù)定

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

4、值法平均值法的算法如下：的算法如下：(1)(1)產(chǎn)生產(chǎn)生RND 隨機(jī)數(shù)：隨機(jī)數(shù)：r1，r2，rn；（2）令令 ui=a(ba)ri，i=1，2，n； （3）計算計算 作為作為I 的估計值的估計值. . 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 足夠大時，得近似公式：足夠大時，得近似公式： niinbafabdxxfI11)()

5、()(注：注： 平均值法本質(zhì)上是用樣本平均值作為平均值法本質(zhì)上是用樣本平均值作為總體教學(xué)期望的估計?？傮w教學(xué)期望的估計。 MC 模擬是一種試驗近似方法模擬是一種試驗近似方法 , 試驗次數(shù)試驗次數(shù)如何確定？如何確定？希望：模擬次數(shù)較少、希望：模擬次數(shù)較少、 模擬精度較高模擬精度較高頻率法的討論頻率法的討論用事件用事件A出現(xiàn)的頻率作為概率出現(xiàn)的頻率作為概率p 的估計的估計: nkpn 問題：問題：試驗次數(shù)試驗次數(shù) n 多大時，對給定的置信度多大時，對給定的置信度1（01），），估計精度達(dá)到估計精度達(dá)到. 即問：取多大的即問：取多大的n 使使 1pnkPppPn成立？成立？證明證明頻率法是事件頻率法

6、是事件A出現(xiàn)的頻率作為概率出現(xiàn)的頻率作為概率p的估計的估計 nkpn 答案：答案：22)1( zppn其中其中, z是正態(tài)分布的臨界值是正態(tài)分布的臨界值.n次獨(dú)立試驗中次獨(dú)立試驗中A出現(xiàn)的次數(shù)出現(xiàn)的次數(shù)knB(n, p).由中由中心極限定理知心極限定理知 )()( pnkpnPppPn )1()1()1(pnpnpnpnpkpnpnPn)1()1(pnpnpnpn 1)1(2 pnpn 11)1(2pnpn令令平均值法平均值法 在給定在給定和和下所需的試驗次數(shù)下所需的試驗次數(shù)的估計式為的估計式為 01202222)(11/niixxnSSzn查得正態(tài)分布的臨界值查得正態(tài)分布的臨界值z, ,可解

7、得可解得 22)1( zppn試驗次數(shù)估計式的分析試驗次數(shù)估計式的分析 01202222)(11/niixxnSSzn為估計概率為估計概率p做做模擬模擬, ,卻又需要卻又需要用用p去估計模擬去估計模擬次數(shù)次數(shù)n. 如何計算如何計算S2 ？ 解決方法解決方法：先做：先做n0 次模擬（稱為學(xué)習(xí)樣本次模擬（稱為學(xué)習(xí)樣本) )，根，根據(jù)學(xué)習(xí)樣本據(jù)學(xué)習(xí)樣本. .（1）先求出先求出p的估計，再估計模擬次數(shù)的估計，再估計模擬次數(shù)n : : 22)1(.1 zppn22)1( zppn（2）計算出的樣本方差計算出的樣本方差S2 , ,用來估計用來估計n. 2. M C模擬的估計精度模擬的估計精度與試驗次數(shù)與試

8、驗次數(shù)n的平的平方根成反比方根成反比, 若精度若精度提高提高10倍倍, ,則試驗次數(shù)則試驗次數(shù)n要增大要增大100倍倍. . P197表表8.2中列出了置信度為中列出了置信度為0.95 時時, , 在不同在不同精度精度及概率及概率p條件下頻率法所需試驗次數(shù)。條件下頻率法所需試驗次數(shù)。 對該表進(jìn)行分析，能得到什么結(jié)論？對該表進(jìn)行分析，能得到什么結(jié)論？1. 精度提高，試驗次數(shù)大幅提高；精度提高，試驗次數(shù)大幅提高；2. 事件發(fā)生概率越接近事件發(fā)生概率越接近0.5，試驗次數(shù)越高；，試驗次數(shù)越高； 核反應(yīng)堆屏蔽層是用一定厚度的鉛包圍反應(yīng)核反應(yīng)堆屏蔽層是用一定厚度的鉛包圍反應(yīng)堆，用以阻擋或減弱反應(yīng)堆發(fā)出的

9、各種射線堆，用以阻擋或減弱反應(yīng)堆發(fā)出的各種射線.在各種射線中在各種射線中, 中子對人體傷害極大，因此，中子對人體傷害極大，因此，在屏蔽層的設(shè)計中在屏蔽層的設(shè)計中, 了解中子穿透屏蔽層的概了解中子穿透屏蔽層的概率對反應(yīng)堆的安全運(yùn)行至關(guān)重要率對反應(yīng)堆的安全運(yùn)行至關(guān)重要.例例7.3.2 核反應(yīng)堆屏蔽層設(shè)計問題核反應(yīng)堆屏蔽層設(shè)計問題1.問題背景問題背景假定屏蔽層是理想的均勻平板假定屏蔽層是理想的均勻平板 一個中子進(jìn)入屏蔽層后運(yùn)動的物理過程一個中子進(jìn)入屏蔽層后運(yùn)動的物理過程: :中中子以初速度子以初速度v0和方向角和方向角射入屏蔽層射入屏蔽層, ,運(yùn)動一運(yùn)動一段距離后與鉛核發(fā)生碰撞，中子獲得新的速度段距

10、離后與鉛核發(fā)生碰撞，中子獲得新的速度及方向及方向(v1,1). 再游動一段距離后再游動一段距離后, ,與鉛核發(fā)生與鉛核發(fā)生第二次碰撞第二次碰撞, ,并獲得新的狀態(tài)并獲得新的狀態(tài)(v2,2),如此等等如此等等, ,經(jīng)過若干次碰撞后經(jīng)過若干次碰撞后, ,出現(xiàn)下述情況之一時中子出現(xiàn)下述情況之一時中子終止運(yùn)動過程終止運(yùn)動過程 1）中子被彈回反應(yīng)堆；中子被彈回反應(yīng)堆；2）中子穿透屏蔽層；）中子穿透屏蔽層；3）第第n次碰撞后次碰撞后, ,中子被屏蔽層吸收中子被屏蔽層吸收. .D返回返回穿透穿透吸收吸收三種三種狀態(tài)狀態(tài)為使屏蔽層的厚度達(dá)到安全設(shè)計要求為使屏蔽層的厚度達(dá)到安全設(shè)計要求, ,在計在計算機(jī)上對中子

11、在屏蔽層的運(yùn)動過程進(jìn)行模擬算機(jī)上對中子在屏蔽層的運(yùn)動過程進(jìn)行模擬 闡述中子的運(yùn)動，闡述中子的運(yùn)動，為模擬做理論準(zhǔn)備為模擬做理論準(zhǔn)備2. 簡化假設(shè)：簡化假設(shè)： *1 假定屏蔽層平行板厚度為假定屏蔽層平行板厚度為D=3d，其中其中d 為兩次碰撞之間中子的平均游動距離；為兩次碰撞之間中子的平均游動距離； *2 假設(shè)在第假設(shè)在第10 次碰撞以后，中子速度下降次碰撞以后，中子速度下降到為某一很小數(shù)值而終止運(yùn)動（被引收）到為某一很小數(shù)值而終止運(yùn)動（被引收）. . 因每次碰撞后因每次碰撞后, 中子因損失一部分能量而速中子因損失一部分能量而速度下降度下降. *3 假定中子在屏蔽層內(nèi)相繼兩次碰撞之間假定中子在屏

12、蔽層內(nèi)相繼兩次碰撞之間游動的距離服從指數(shù)分布游動的距離服從指數(shù)分布； *4 中子經(jīng)碰撞后的彈射角中子經(jīng)碰撞后的彈射角 U(0, 2).思考思考：請仔細(xì)分析以上假設(shè)的合理性：請仔細(xì)分析以上假設(shè)的合理性.3. 中子運(yùn)動的數(shù)學(xué)描述中子運(yùn)動的數(shù)學(xué)描述引進(jìn)變量引進(jìn)變量： 彈射角彈射角i 第第i 次碰撞后中子的運(yùn)動方向與次碰撞后中子的運(yùn)動方向與x 軸正向的夾角軸正向的夾角. .xi 第第i 次碰撞后中子所處位置與屏蔽層內(nèi)次碰撞后中子所處位置與屏蔽層內(nèi) 壁的距離壁的距離.R i 中子在第中子在第 i 次碰撞前后的游動距離次碰撞前后的游動距離.DD0 xxi i三個變量三個變量均為隨機(jī)均為隨機(jī)變量變量 中子在

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

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

15、數(shù)( (個個) ) 穿透穿透(%) 吸收吸收(%) 返回返回( (%) 100100030005000 30.026.026.526.3 28.023.421.8 22.0 42.050.651.7 51.7 中子穿透屏蔽層的百分比超過了中子穿透屏蔽層的百分比超過了1/4, 模擬模擬結(jié)果表明屏蔽層厚度結(jié)果表明屏蔽層厚度D=3d不合適不合適. 多厚的屏蔽層才能使穿透的概率多厚的屏蔽層才能使穿透的概率W106？問問 題題： 如何解決這個問題？如何解決這個問題？思路？思路？1. 計算機(jī)收索法計算機(jī)收索法 增大屏蔽層的厚度，如增大屏蔽層的厚度，如D=6d、12d、24d、36d，交由計算機(jī)進(jìn)行模擬交由

16、計算機(jī)進(jìn)行模擬, 并搜索到所并搜索到所求解求解.2. 分析法分析法 12mDDD 設(shè)計屏蔽層的厚度：設(shè)計屏蔽層的厚度： x=mD 將屏蔽層視為將屏蔽層視為m層厚度均為層厚度均為D的平行板的平行板. . 由于碰撞的能量損失由于碰撞的能量損失, ,中子穿過屏蔽層的平中子穿過屏蔽層的平均速度會逐層下降均速度會逐層下降. . 設(shè)設(shè)WD 是中子穿過厚度為是中子穿過厚度為D 屏蔽層的概率屏蔽層的概率,則則穿過整個屏蔽層的概率穿過整個屏蔽層的概率W 滿足滿足 mDWW 利用模擬結(jié)果利用模擬結(jié)果：當(dāng)當(dāng)D=3d，WD0.25，令令(WD)m106， 或或 (m)410 6967. 9602. 064log10log6 m取屏蔽層的厚度取屏蔽層的厚度 x=10D=30d，可使穿透屏蔽層的概率可使穿透屏蔽層的概率w106 注注：模擬：模擬5000個中子個中子的運(yùn)動的運(yùn)動, ,用穿透屏蔽的用穿透屏蔽的頻率估計穿透概率頻率估計穿透概率, ,由表由表8.2可知精度大約只有可知精度大約只有1 1, ,模擬精度太低模擬精度太低, ,應(yīng)適當(dāng)增大模擬次數(shù)應(yīng)適當(dāng)增大模擬次數(shù). .總結(jié)：總結(jié)：模擬的意義？模擬的意義？ 1.模擬方法本質(zhì)上是模擬方法本質(zhì)上是試驗性試驗性的的, ,模擬系統(tǒng)是模擬系統(tǒng)是 現(xiàn)實系