計算物理學蒙特卡羅方法

1、第八講 蒙特卡羅方法蒙特卡羅（Monte Carlo簡稱MC）方法又稱隨機抽樣法（Random Sampling）、隨機模擬（Random Simulation）或統(tǒng)計試驗法(Statistic Testing）。這個方法的起源可以追溯到十七世紀或更早的年代。Monte Carlo 是摩納哥（Monaco）的一個著名城市，位于地中海之濱，以旅游賭博聞名。Von Neumann 等人把計算機隨機模擬方法定名為Monte Carlo方法，顯然反映了這種方法帶有隨機的性質。簡單地說，MC方法是一種利用隨機統(tǒng)計規(guī)律，進行計算和模擬的方法。它可用于數值計算，也可用于數字仿真。在數值計算方面，可用于多重積

2、分、線性代數求解、矩陣求逆以及用于方程求解，包括常微分方程、偏微分方程、本征方程、非齊次線性積分方程和非線性方程等。在數字仿真方面，常用于核系統(tǒng)臨界條件模擬、反應堆模擬以及實驗核物理、高能物理、統(tǒng)計物理、真空、地震、生物物理和信息物理等領域。§8.l蒙特卡羅方法的基礎知識8.1.1 基本概念為了對MC方法有一點初步的認識，請先看使用MC方法的幾個例子。蒲豐投針問題：蒲豐（Buffon-法國著名數學家）在1777年發(fā)現隨機投針的概率與無理數之間的關系這個問題是說，若在平面上畫有距離為a的平行線束，向平面上投擲長為的針，試求針與一平行線相交的概率。這個問題的解法如下：以M表示落下后針的中

3、點，表M與最近一平行線的距離，表針與此線的交角，見上圖?？梢?，這兩式決定平面上一矩形R；為了使針與一平行線（這線必定是與針中點M最近的平行線）相交，充分而且必要條件是這個不等式決定R中一個子集G。因此，我們的問題等價于向R中均勻分布地擲點而求點落于G中的概率P.根據概率的幾何意義，得此式提供了求值的一個方法：可以通過投針事件求得針與平行線相交概率P，求得值： （8.1）若投擲次數為m，針與平行線相交的次數為n，那么即 于是，可用投針試驗來求無理數的近似值下表列舉了歷史上若干學者投針試驗計算值的結果：試驗者年份投針次數的計算值Wolf185050003.1596Smith185532043.15

4、53Fox189411203.1419Lazzarini190134083.1415929射擊問題（打靶游戲）：設表示射擊運動員的彈著點到靶心的距離，表示擊中處相應的得分數（環(huán)數），分布密度函數表示該運動員的彈著點分布，它反映運動員射擊水平。積分 （8.2）表示這個運動員的射擊成績。用概率語言說，就是隨機變量的數學期望值，記為?，F在，假設這個射擊運動員射擊N次，彈著點依次是環(huán)數分別為，則自然地認為N次射擊得分的平均值 （8.3）這個平均值相當好地代表了這個射擊運動員的成績。換句話說，是積分（8.2）式的一個估計值（或近似值）。這個例子通常稱為打靶游戲，它直觀地說明了蒙特卡羅方法計算定積分的基本

5、思想。為進一步闡明這個思想，我們再舉個例子：計算積分 （8.4）直觀上，就是在邊長為1的正方形里隨機投點，當點落在曲線下面，對積分值有“貢獻”，否則對積分值無“貢獻”。為此，假設向這個邊長為1的正方形里隨機投點N次，點落在曲線下面n次，則（8.4）式積分值近似為來近似。 從上述幾個例子可以看到，當所要求解的問題是某種事件出現的概率，或者是某個隨機變量的期望值時，它們可以通過某種“試驗”的方法，得到這種事件出現的頻率，或者這個隨機變量的平均值，并用它們作為問題的解。這就是蒙特卡羅方法的基本思想。所以用蒙特卡羅方法求解問題時，首先要建立一個隨機模型，然后要構造一系列的隨機數用以模擬這個過程，最后要

6、作統(tǒng)計性的處理。關于建立隨機模型，因問題而異。8.1.2 隨機變量及其分布函數 在一定條件下發(fā)生的事件分為必然事件（必然發(fā)生）、不可能事件（恒不發(fā)生）和隨機事件（可能發(fā)生也可能不發(fā)生），事件發(fā)生的可能性大小用概率表示。必然事件的發(fā)生概率為1，不可能事件的概率為零；隨機事件發(fā)生的概率。由于測量的隨機誤差和物理現象本身的隨機性。一次測量得到的某個值是隨機的，因此，實驗觀測的物理量是隨機變量，被研究的物理問題是一個隨機事件通常，描寫隨機事件A發(fā)生的概率用來表示，顯然， 經常碰到的隨機變量有兩類：一類是離散型隨機變量，這種隨機變量只能取有限個數值，能夠一列舉出來；另一類是連續(xù)型隨機變量，這種隨機變量的

7、可能值連續(xù)的分布某個區(qū)間離散型隨機變量可用分布列： （8.5）來描寫，它表示取值的概率為（）即 （8.6）分布列描寫了離散型隨機變量的概率分布 連續(xù)型隨機變量可能取的值是不可列的，因此不能象離散型隨機變量那樣，用分布列來描寫它的概率分布，而要用概率分布密度函數來描寫考慮連續(xù)型隨機變量落在區(qū)間內的概率，如果極限 （8.7）存在，則函數描寫了在點的概率密度，把叫做隨機變量的概率分布密度，簡稱分布密度或密度函數于是，隨機變量落在內的概率可寫為 （8.8）顯然，上式只有當可積時才有意義。 此外，還可定義隨機變量的分布函數來描寫隨機變量的概率分布規(guī)律分布函數一定義為 （8.9）即分布函數在處的值等于隨機

8、變量取值小于或等于這樣一個隨機事件的概率顯然， 對于離散型隨機變量，分布函數為階梯函數 （8.10）對于連續(xù)型隨機變量，分布函數與分布密度滿足如下關系：所以分布密度函數和分布函數的性質： 對于任意實數： 若在點連續(xù)，則有為了描述隨機變量的數學特征還需引進數學期望和方差的概念對于離散型隨機變量其數學期望定義為 （8.11）方差定義為 （8.12）對連續(xù)型隨機變量，若已知分布密度，則數學期望是 （8.13）方差 （8.14）在實際問題中，對于一組N個實驗觀測數據，相應于某一個隨機變量的一個樣本，可以用直方圖來形象地表示樣本中數據分布的規(guī)律性先將隨機變量的取值范圍劃分為若干個區(qū)間，將落入每一區(qū)間的數

9、據的個數m（稱為頻數）與隨機變量的取值區(qū)間之間的關系畫成階躍曲線，即構成了直方圖數 mN是觀測數據落入這個小區(qū)間的頻率直方圖的根坐標為隨機變量的取值范圍，縱坐標為頻數圖表示的是某一個隨機事件的直方圖。8.1.3 大數定理和中心極限定理 概率論中的大數定理和中心極限定是蒙特卡羅方法的數學基礎。大數定理： 設為一隨機變量序列，獨立同分布，數學期望值存在，則對任意，有 （8.15）是算術平均。大數定理指出，當時，算水平均收斂到數學期望（或統(tǒng)計平均）也就是說，可以用算術平均代替數學期望。對于收斂的程度，和誤差估計，要用到下面的中心極限定理中心極限定理：設為一隨機變量序列，獨立同分布，數學期望為，方差，

10、則當時， （8.16）利用中心極限定理，當很大時，可得到誤差 （8.17）成立的概率為。通常將稱為置信度，為置信水平。的定義為 （8.18）和的關系可以計算得到，下面給出常用的幾組和的值：0.50.050.020.010.67451.96002038632.5758從（8.17）式可以看到，算術平均值收斂到的階為，可見，蒙特卡羅方法收斂的階很低，收斂速度很慢。當，誤差稱為概率誤差。誤差由和決定在固定的情況下，要提高 1位精確度，就要增加 100倍試驗次數相反，若減小10倍，就可以減少100倍工作量因此，控制方差是應用蒙特卡羅方法中很重要的一點如果要求置信水平為0.95（95%），計算得，則說明

11、誤差估計的概率是0.95。§8.2隨機數和隨機抽樣在計算機上用蒙特卡羅方法模擬某過程時，需要產生具有各種概率分布的隨機變量最簡單和最基本的隨機變量就是區(qū)間上均勻分布的隨機變量這些隨機變量的抽樣值就稱為隨機數以后談到隨機數，如果不加特別說明，專指區(qū)間上均勻分布的隨機數各種其他分布的隨機變量的抽樣值都可借助均勻分布的隨機數得到 真正的隨機數如同擲骰子那樣，產生1-6范圍內的隨機數整數。抽獎用的搖號碼機則可產生09范圍內的隨機整數。這些真正的隨機數除統(tǒng)計規(guī)律外無任何其它規(guī)律可循。偽隨機數，或稱贗隨機數，是指按照某種算法可以給出的似乎隨機地出現的數。既然數的給出是按某種算法，也就是按某種規(guī)律

12、，那這種隨機數就必然具有一定的周期。設其周期為n，則第nl個數就等于第一個數，此后均依次重復出現。當然，如果周期n足夠大，可使在整個使用過程中不至于表現出其周期性，偽隨機數也是實用的。例如，計算機中的偽隨機數發(fā)出器要求其周期大于計算機的記憶單元數。此外偽隨機數的統(tǒng)計性質是表征隨機數品質的又一重要指標。均勻分布的隨機數，既要求數的出現是隨機的，又要求數的分布是均勻的。至于如何評估隨機數分布的均勻程度，將在本節(jié)稍后討論。8.2.1 均勻分布隨機數的產生 均勻分布的隨機數的產生有許多方法，通常采用乘同余法。例如，可按下列公式產生隨機數 （8.19）或寫成其中為給定常數。給出后，就可以用式（8.19）

13、依次給出等一系列隨機數。如何確定常數和，這是個十分關鍵的問題，是人們仍在不斷研究探索的課題。下面僅給出確定和的一般原則。關于的取值一般取，其中 m為計算機中二進制數的字長，則為計算機所能表示的最大整數。例如字長16位時，可取 等。關于的取值，一般取，其中 M為任一正整數。例如有取 等。建議取，這樣統(tǒng)計性較好。關于的取值，一般取為奇數?？梢则炞C，當為奇數時，周期是T，其它參數不變，當為偶數時，周期則為T/2.例如設，得隨機數序列為：周期為16。按（8.19）產生的隨機數其值域為。如果要產生之間的隨機數，只需將原產生的每個隨機數再除以即可。用類似的方法，經過簡單的變換，就可產生任何所需值域的贗隨機

14、數序列。例如在16位微機上，可用下式（8.20），并令，來產生之間的贗隨機數序列 (8.20)見程序：rand1.f90下面給出另一個產生之間均勻隨機數程序：rand.for,plot_rand.mimplicit double precision(a-h,o-z)real Ropen(2,file='rand.dat')ix=32765do i=1,100 call rand(ix,r) write(*,*) 'r=',r write(2,*) i,rend doENDsubroutine rand(ix,r)i=ix*259ix=i-i/32768*3276

15、8r=float(ix)/32768returnend8.2.2 隨機性統(tǒng)計檢驗 一個好的隨機數發(fā)生器或一個好的隨機數生成程序必須滿足兩個條件：第一，所生成的隨機數序列應當具有足夠長的周期；第二，所生成的隨機數序列應當具有真正隨機數序列所具有的統(tǒng)計性質。其周期的長短比較容易測試和判斷，而統(tǒng)計性質的優(yōu)劣則不那么簡單。通常對統(tǒng)計性質的檢驗方法是采用頻數分布檢驗：對于一個均勻分布的隨機數發(fā)生器，設所產生隨機數序列的值域為，則所產生的隨機數字應與從之間均勻的頻數分布相一致。為了檢驗頻數分布情況，可按畫統(tǒng)計直方圖的方法，將整個值域分成M個寬度相等的子區(qū)間，設是第區(qū)間內出現的隨機數的個數，即第i子區(qū)間的頻

16、數，則所有子區(qū)間中隨機數個數的平均頻數第i子區(qū)間頻數的偏差則頻數的方均根偏差如果所產生的N個隨機數均勻分布于整個值域，則且在任一子區(qū)間內出現個數字的幾率服從高斯分布規(guī)律，即其中為標準偏差?？梢?，一個均勻分布的隨機數序列，其最可幾頻數應為，而其頻數在范圍內的幾率應為0.68，其頻數的方均根偏差應接近于標準偏差。通常檢查均勻分布隨機數分布的均勻程度就可從這兩個方面考核，即各頻數是否接近于平均頻數和頻數的方均根偏差是否接近于標準偏差。8.2.3 產生給定分布的隨機變量-隨機抽樣隨機抽樣的方法很多，在計算機上實現時要考慮運算量的大小，也就是所謂“抽樣費用”因為應用蒙特卡羅方法求解一個物理問題時，大量的

17、計算時間將用于隨機抽樣，所以隨機抽樣方法的選取往往決定算題的費用但對不同問題和不同機器，不同的方法也可以有不同的評價下面介紹幾種常用的隨機抽樣方法，供讀者選擇1連續(xù)型分布的直接抽樣法利用0,1區(qū)間上的均勻分布隨機數可以產生具有給定分布的隨機變量數列我們知道，若隨機變量具有分布密度，則其分布函數 （8.21）可知，為單調非減函數，而且存在反函數?？梢宰C明，這個分布就是0,1區(qū)間上滿足均勻分布的隨機變量。令均勻分布的隨機數等于分布函數，則得到：就是具有分布密度的隨機變量的抽樣。 事實上：由分布函數定義 如果是0,1區(qū)間均勻分布的隨機變量其中用到了0,1區(qū)間均勻分布的隨機變量的分布密度函數抽樣的步驟

18、為： 給定分布密度； 計算； 產生隨機數； 計算令 重復，。=例題8.1 求經常用來描述電子元件的穩(wěn)定時間，系統(tǒng)的可靠性和粒子游動的自由程等的指數密度分布的隨機抽樣。解：，令，則由于與同分布，故可得指數分布的抽樣公式=例題8.2 求均勻密度函數在區(qū)間a,b上均勻分布的隨機數解：令0,1區(qū)間的均勻分布的隨機數得區(qū)間a,b上均勻分布的隨機數=2離散型隨機變量的抽樣法設是離散型隨機變量，分別以概率取值，可按下面步驟抽樣： 選取隨機數； 確定滿足不等式的值（這里約定）； 對應的就是所求的抽樣值，即； 事實上，如果設：，則由=例題8.3 按學生的成績抽樣，某班26名學生成績分布如下：分數擋1-1011-

19、2021-3031-4041-5051-6061-7071-8081-9091-100等價5152535455565758595人數0000258731概率00000.080.190.310.270.110.0400000.080.270.580.850.961.012345678910上面，解：這個分布抽樣方法如下：（1） 產生隨機數（由程序計算），假設是（2） 由不等式確定,（3） 由，確定抽樣值為分計算程序：scores.f90,plot_scores.m=例題8.4 設能量為的光子與物質相互作用的產生光電效應的截面為，產生康普頓散射截面，電子對效應截面為，確定相互作用類型抽樣解：設總截

20、面：，抽樣方法如下：（1） 產生隨機數；（2） 計算，如果，則發(fā)生光電效應；否則計算，如果，則發(fā)生康普頓散射；否則，發(fā)生電子對效應=8.2.4 蒙特卡羅方法基本步驟和基本思想 用蒙特卡羅方法處理的問題可以分為兩類：一類是隨機性問題，例如中子在介質內的傳播問題和后面要介紹的原子核裂變問題等對于這一類問題，通常采用直接模擬方法首先，必須根據物理問題的規(guī)律，建立一個概率模型（隨機向量或隨機過程），然后用計算機進行抽樣試驗，從而得出對應于這一物理問題的隨機變量的分布。假定隨機變量y是我們的研究對象，它是m個互相獨立的隨機變量的函數，如果概率分布密度為，則用蒙特卡羅方法計算的基本方法是：在計算機上用隨機

21、抽樣的方法根據概率分布密度抽樣，產生N組隨機變量，計算的值，用這樣的樣本分布來近似y的函數，由此可計算出這些量的統(tǒng)計值 由上可知，蒙特卡羅方法的計算過程就是用數學方法模擬實際的物理過程，它主要是在計算機上產生已知分布的隨機變量樣本以代替昂貴的甚至難以實現的實驗蒙特卡羅方法又被看作是用計算機來完成物理實驗的一種方法另一類是確定性問題在解決確定性問題時，首先要建立一個有關的概率統(tǒng)計模型。使所求的解就是這個模型的概率分布或數學期望，然后對這個模型作隨機抽樣，最后用其算術平均值作為所求解的近似值。§8.3 M-C方法的應用8.3.1 方程求根的M-C方法考慮方程 （8.3.1）其中是定義在a

22、,b區(qū)間實的連續(xù)函數。如果已知則方程（8.3.1）在a,b區(qū)間內至少有一個根。方法一：用頻率近似概率先討論在a,b區(qū)間內只有一個根的情況。（1）當，則a就是所求的根；（2）否則，設是a,b區(qū)間內的一個根，是a,b區(qū)間內均勻分布的隨機數，則落在內的概率是因此，在區(qū)間a,b內均勻投N個點，設其中有M個落在根的左側，于是有 （8.3.2）這樣，用M-C 求單根的步驟如下： 定義隨機變量，當第i次試驗（投點）時；，其他情況；其中是(0,1)區(qū)間上均勻分布的隨機數。 令： 當時， 當時， 對于區(qū)間a,b內有多個根時，可用分析的方法確定每個單根的界限，在每個界限內應用上述方法?；蛘邚腶點出發(fā)，以h為基本步

23、長向前跨長為h的小區(qū)間（h要選的合適使每個區(qū)間至少有一個根。太大容易丟根，太小浪費時間），當跨的小區(qū)間兩端的函數同號時，繼續(xù)向前跨；異號時，用上述方法求出小區(qū)間中的根。=例題8.5-1 求方程的全部實根解：見計算程序:root_MC_1.f90=方法二：取離根左右最近的兩值平均即?。簞t：，是根的一個近似估計。=例題8.5-2：求方程在區(qū)間0，內的一個實根。解：見計算程序:root_MC_2.f90=8.3.2 用M-C方法計算定積分 設區(qū)間a,b的中的隨機變量的分布由密度函數給定是區(qū)間a,b上的連續(xù)函數，數學期望 （8.3.3）存在，則積分（8.3.3）式可用如下方法計算近似值設隨機變量按密度

24、分布抽樣的一系列值為相應產生的隨機抽樣值，則根據大數定理有可見，只要n充分大，積分（8.3.3）式有近似值 (8.3.4)現在來討論積分的值為此，選擇某種密度分布函數滿足且能很方便地生成具有密度函數為的隨機抽樣同時將積分于是歸結為積分（8.3.1）的形式。在多數情況下，往往取為a,b上均勻分布的概率密度，現在，在區(qū)間（a,b）上均勻分布的隨機數中抽取，計算，然后計算平均值 （1）方法一在0,1之間取均勻分布的隨機數序列，并令得到區(qū)間均勻分布的隨機數，只要足夠大，則有=例題8.6用蒙特卡羅方法計算定積分解：計算程序MC_1.for及結果如下：=對于多重積分等復雜問題，基本方法都是相同的。例如計算

25、多重定積分取0,1區(qū)間均勻分布的隨機數序列：計算： 只要足夠大，則有（2） 方法二現在仍然討論積分假設函數在區(qū)間a,b上有最大值，的值。為此，在下圖作矩形ABCD，其寬度為，高為，則矩形面積 ,給出兩個隨機數，計算得如果滿足成立，隨機點落在矩形陰影區(qū)域內。設總共產生的隨機點數為N, 落在矩形陰影區(qū)域內的點數為M,當N 足夠大時，定積分=【例題8.7】：計算半徑為1圓的面積(計算值)解：見計算程序：pi_1.f, pi_2.f, pi_3.f90=【例題8.8】計算多重積分（這個積分的結果是歐拉常數）解：（1）在單位超立方體區(qū)域上構造一個聯合概率分布密度：其中當，其它：則積分可以寫為（2）用算術

26、平均值來近似數學期望?，F從分布密度抽取隨機向量的N個樣本，每個由s個0,1區(qū)間的均勻分布的隨機數組成，即N個樣本為:（3） 就是積分的一個近似估計。（4） 下面是模擬程序：Euler.f90gama( 5)= -0.57557088, gama( 6)= -0.57453740 gama( 7)= -0.55329597, gama( 8)= -0.56524354 gama( 9)= -0.57235909可見積分結果與積分重數沒有多大關系。 =【例題8.9】一球體半徑，球上有一半徑的圓柱形空洞，其軸線與球的直徑重合。試用M-C方法求實體的體積。解：如右圖所示，令球體中心位于坐標系的原點O處

27、，作邊長為的正方體，其中心與球心重合，則正方體的體積;產生一組三個隨機數（，它們的值域均為;然后判斷該隨機點是否位于實體內，其判據是:且若共產生了N組隨機數，而滿足上述判據的有M組，則球體的體積為了提高測量精度，可重復上述計算過程，例如用完全相同的方法再做兩遍，取三次結果的平均值作為最終結果。=對于一般（變積分限）的多重定積分=【例題8.10】計算函數積分限的多重積分解：該積分的區(qū)域見圖92。由圖92可見，該積分x軸積分限從，y軸積分限是不規(guī)則的，由，在常數區(qū)間產生一組隨機數序列在函數區(qū)間產生一組隨機數序列利用計算多（s）重積分的平均值方法：（對于常數積分限）是s 維空間積分區(qū)域的超體積（見下

28、面說明）。但是要注意，對于，積分限變化的情況，例如二維情況，對于中關于y是隨x變化的部分要放到求和里面，因為，對于此時的，的變化區(qū)間變化了，在本文題中是，所以要先計算然后再乘上x軸區(qū)間寬度（）=1計算程序：/* mc2.c */，計算結果為：，積分值精確到小數點后第6位是589/90=6.544444?？梢奙C 方法計算函數積分限的二重積分，僅用15000個隨機數，結果就相當滿意了。=【例題8.11】計算程序：mc3.c, 計算結果為： MC方法求解Laplace方程利用蒙特卡羅的隨機游動方法求解二維Laplace方程：其差分格式為：設Q是邊界上的點，P是區(qū)域內的點，求的M-C方法為隨機游動方

29、法，步驟如下： 粒子的狀態(tài)參數為； 粒子由狀態(tài)出發(fā)； 假設粒子已游動 n 步達到，從此點出發(fā)，以1/4等概率到達四個鄰點的一個，記為；若是邊界點，則游動終止，并記下邊界處的函數值；若不是邊界點，則重復和，直至達到邊界。由狀態(tài)出發(fā)重復m 次，這樣產生粒子游動的m次歷史，得然后，從其它點出發(fā)，求其它點的值。這種方法對于復雜邊界情況特別適用，特別是求區(qū)域內任意點的值很有效。=【例題8.12】利用蒙特卡羅方法求解例題7.3的問題解：計算程序MC_laplace.for,plot_MC_laplace.m= 8.3.4 核鏈式反應的模擬 放射性物質的鏈式反應是一個隨機過程，可借助計算機用MC方法模擬和研

30、究。由原子核物理知識可知，的原子核本質上是不穩(wěn)定的，會自發(fā)地發(fā)生裂變。裂變的激烈程度可用放射性物質的半衰期來描述，半衰期是指大量核中有12的核發(fā)生裂變所需要的時間。的半衰期為 7億多年。因此任何時刻發(fā)生裂變的核只是相對很小部分，其釋放的能量只能使其本身微微溫熱。但是，在一定條件下，自發(fā)裂變放出的兩個中子轟擊其它核而被吸收，引起新的裂變而放出更多的中子，這更多的中子又引起新一輪更多的裂變，依次類推，可迅速釋放出大量能量，甚至引起爆炸，這就是鏈式反應。 設開始有 N個核發(fā)生裂變，每個核放出兩個中子，稱為第一代中子，共2N個。2N個中子又感生新一輪裂變，產生第2代中子，為4N個。如此進行下去，直至第

31、n次裂變，產生第n代中子為個。按此計算，30代可產生裂變的核數為億，即為第一次裂變核數的10億倍。現在的問題是，在什么條件下才能發(fā)生鏈式反應呢？其基本要求是裂變所產生的兩個中子中至少有一個能使第二個鈾核發(fā)生裂變。為此要求核材料中雜質的含量，包括的含量應足夠少，以避免中子被和其它雜質所吸收。另外，由于熱中子使裂變的機會很大，所以在鈾堆中還必須加入減速劑，如重水或石墨等，以使快中子減速到熱中子。最后，非常重要的條件是鈾堆的體積必須足夠大，以避免裂變所放出的中子過多地未與鈾核相遇而飛出鈾體外。這就涉及到臨界體積和臨界質量的概念。所謂臨界質量是指可裂變物質能發(fā)生自持鏈式反應的最小質量。由于鈾核體積很小

32、，一鈾核裂變放出的中子在和另一鈾核作用并使之發(fā)生裂變之前，平均地說要經過一定相對很長的距離，約為厘米數量級。因此，假定有個核發(fā)生自發(fā)裂變而放出2個中子，其中N個中子在鈾塊中引起另外的核發(fā)生裂變，其余的中子未與其它核碰撞而飛出鈾塊。為描述一次裂變能引起下一次裂變的可能程度，定義裂變過程的倍增系數：不難理解，維持自持鏈式反應的條件是：，即。倍增系數kl是臨界質量Mc的條件。k的值與前面論及的諸多因素有關，本節(jié)將只限于討論k與鈾塊的質量和形狀有關的問題，用計算機程序來模擬具有一定大小和形狀的鈾塊中大量隨機的裂變過程，統(tǒng)計算出相應的倍增系數k。設鈾塊為長方體，發(fā)生裂變的鈾核位于鈾塊內隨機點處，如下圖所

33、示。隨機點坐標的值域為：鏈式反應模擬示意圖該核子裂變反應產生兩個中子，其運動方向可以用兩個角坐標和來描述，見上圖（b）。釋放出的每一個中子按飛向各個方向的幾率均等來考慮，或者說中子飛行方向的幾率是按以為頂點的立體角均勻分布的。立體角元可表示為可見，按立體角均勻分布是按角均勻分布和按均勻分布，而并非按角均勻分布。因此，對應的兩個隨機數的值域為下面積算按極角正余弦分布和方位角分布的抽樣公式：各向同性散射角余弦分布和散射方位角分布按分布密度：和抽樣的方法就是采用直接抽樣： 平均地說，能否擊中另一個核只取決于中子在鈾塊體內飛行的距離。假設在0至1cm距離之間經過任何一段相同距離擊中另一鈾核的幾率均等。

34、或者說，中子在擊中另一鈾核之前飛行的距離為01之間均勻分布的隨機數。因此與飛行距離相應的隨機數為 由此可計算出被擊中的鈾核的位置最后，檢查計算出的碰撞點是否位于鈾塊體內。若在鈾塊體內，則累計引起新裂變中子數N。按照上述原則，歸納計算k的具體步驟為：1給定鈾塊質量M、鈾塊邊長比和用于計算k的隨機自發(fā)裂變核的個數，即舊裂變核個數，并設所選約化單位可使鈾塊的密度為1，體積為V，則得或 2產生九個之間的隨機數：3舊裂變核位置：4.舊裂變放出的兩個中子的方向5.中子的飛行距離6.可能發(fā)生新裂變的位置7. 檢查上述位置，是否在鈾塊體內。如果均滿足，則的值增加1；同樣，如果均滿足，則的值也增加1；8.重復步

35、驟2-7，共執(zhí)行次，然后計算： 若計算結果，則調整M和s的值后再進行上述步驟，直至，此時M即為臨界質量。顯然臨界質量與s有關，與核材料的形狀有關。計算程序：MC_2.f90和plot_k_m.m8.3.5 關于中子貫穿概率問題 原子反應堆的外壁是鉛板圍成的。中子從鉛壁的內側向外輻射，稱中子貫穿概率問題。可以用MC方法模擬中子貫穿鋁板的概率。為簡化問題，假定中子源用一個厚5個單位的鉛板圍住，中子貫穿方向為X方向，Y方向是無限高。實驗表明，中子在10次相撞后能量消耗殆盡即被鉛原子吸收。中子進入鉛板，走一定距離與鉛原子核碰撞（鉛原子直徑d），之后隨機改變方向。又走一定距離，與另一個原子核碰撞。如圖9

36、4所示，如此經過多次碰撞后，中子可能穿透鉛壁輻射到反應堆外，也可能將其能量耗盡被鉛壁吸收，還可能被反射回反應堆內。顯然，鉛壁設計得越厚，穿透的概率就越小，反應堆就越安全。由于每次碰撞后彈出的角度是隨機的，因此對一個中子來說是穿貫，還是被吸收或返回，也是隨機的。中子貫穿概率問題是由大量中子運動的統(tǒng)計規(guī)律決定的。要得到鉛壁厚和穿透率之間的關系以及貫穿概率，用解析方法是極為困難的，而用MC方法模擬中子貫穿鉛板問題的求解大大簡化。中子在壁內的運動與其每次與鉛原子核碰撞后的散射角有關，由三角學可知，取中子每次的自由程為斜邊1，直角三角形中直角邊是。注意隨機散射角為，當和時，中子在鉛壁厚度方向上前進單位距

37、離，當時，中子后退單位距離。在區(qū)間產生一個隨機數作為每次碰撞隨機散射角其中：rand()是0，1區(qū)間的隨機數，中子在x 方向每次運動的距離（初始s=0）模擬一個中子在鉛板中的運動得到一個結果，對大量中子的運動進行模擬的結果，便可以統(tǒng)計出穿透率PP%,吸收率PA%,和返回率PR%.下面是C語言編寫的模擬程序：/* mc1.c */模擬結果是：鉛板厚5個單位，中子數取10000個，得到穿透率PP=1.88%,吸收率PA=15.41%,和返回率PR82.71%.【例題8.13】直線加速武器的邊耦合腔是無限長圓柱筒，內半徑，空心，外半徑，在與之間的圓環(huán)均勻地布滿放射源銅，在內半徑內空心柱內等距離放入一

38、個個的銅片（厚度），求整個柱筒對距離筒軸處觀察點o的放射性貢獻，即計算積分式中的表示由點出發(fā)到點o間經過圓柱筒中銅的總長度。是衰減常數。積分區(qū)域是圓柱筒中銅所占的幾何區(qū)域。解：我們知道，在滿足一定的條件下，放射源與觀察點是可以進行倒易的。這就是通常的光學倒易定理。本文題與能量無關，倒易定理成立。根據光學倒易定理，可以用點源對圓柱筒中銅的通量貢獻來代替圓柱源對點通量的貢獻。即將O點看作放射性點源，將放射性銅看作是觀察者。，表示從點O沿方向經過圓柱筒中銅的總長度。表示各向同性散射角分布密度；表示散射方位角分布密度；計算程序：MC_3.for,計算結果：8.3.6 其他例子1.隨機行走問題定義： 設是一個獨立同分布的隨機變量序列，對于每一個正整數n,設表示和，序列：稱為隨機行走（a random walk）.醉漢行走問題 醉漢開始從一根電線桿的位置出發(fā)（其坐標為，坐標向右為正，向左為負），假定醉漢的步長為，他走的每一步的取向是隨機的，與前一步的方向無