計算物理學蒙特卡羅方法

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

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

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

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

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

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

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

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

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

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

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

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

13、依次給出等一系列隨機數(shù)。如何確定常數(shù)和，這是個十分關(guān)鍵的問題，是人們?nèi)栽诓粩嘌芯刻剿鞯恼n題。下面僅給出確定和的一般原則。關(guān)于的取值一般取，其中 m為計算機中二進制數(shù)的字長，則為計算機所能表示的最大整數(shù)。例如字長16位時，可取 等。關(guān)于的取值，一般取，其中 M為任一正整數(shù)。例如有取 等。建議取，這樣統(tǒng)計性較好。關(guān)于的取值，一般取為奇數(shù)。可以驗證，當為奇數(shù)時，周期是T，其它參數(shù)不變，當為偶數(shù)時，周期則為T/2.例如設(shè)，得隨機數(shù)序列為：周期為16。按（8.19）產(chǎn)生的隨機數(shù)其值域為。如果要產(chǎn)生之間的隨機數(shù)，只需將原產(chǎn)生的每個隨機數(shù)再除以即可。用類似的方法，經(jīng)過簡單的變換，就可產(chǎn)生任何所需值域的贗隨機

14、數(shù)序列。例如在16位微機上，可用下式（8.20），并令，來產(chǎn)生之間的贗隨機數(shù)序列 (8.20)見程序：rand1.f90下面給出另一個產(chǎn)生之間均勻隨機數(shù)程序：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)計檢驗 一個好的隨機數(shù)發(fā)生器或一個好的隨機數(shù)生成程序必須滿足兩個條件：第一，所生成的隨機數(shù)序列應(yīng)當具有足夠長的周期；第二，所生成的隨機數(shù)序列應(yīng)當具有真正隨機數(shù)序列所具有的統(tǒng)計性質(zhì)。其周期的長短比較容易測試和判斷，而統(tǒng)計性質(zhì)的優(yōu)劣則不那么簡單。通常對統(tǒng)計性質(zhì)的檢驗方法是采用頻數(shù)分布檢驗：對于一個均勻分布的隨機數(shù)發(fā)生器，設(shè)所產(chǎn)生隨機數(shù)序列的值域為，則所產(chǎn)生的隨機數(shù)字應(yīng)與從之間均勻的頻數(shù)分布相一致。為了檢驗頻數(shù)分布情況，可按畫統(tǒng)計直方圖的方法，將整個值域分成M個寬度相等的子區(qū)間，設(shè)是第區(qū)間內(nèi)出現(xiàn)的隨機數(shù)的個數(shù)，即第i子區(qū)間的頻

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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