計算物理學(xué)-蒙特卡羅方法

/第八講蒙特卡羅方法蒙特卡羅（MonteCarlo簡稱MC）方法又稱隨機(jī)抽樣法（RandomSampling）、隨機(jī)模擬（RandomSimulation）或統(tǒng)計試驗(yàn)法(StatisticTesting）。這個方法的起源可以追溯到十七世紀(jì)或更早的年代。MonteCarlo是摩納哥（Monaco）的一個著名城市，位于地中海之濱，以旅游賭博聞名。VonNeumann等人把計算機(jī)隨機(jī)模擬方法定名為MonteCarlo方法，顯然反映了這種方法帶有隨機(jī)的性質(zhì)。簡單地說，MC方法是一種利用隨機(jī)統(tǒng)計規(guī)律，進(jìn)行計算和模擬的方法。它可用于數(shù)值計算，也可用于數(shù)字仿真。在數(shù)值計算方面，可用于多重積分、線性代數(shù)求解、矩陣求逆以及用于方程求解，包括常微分方程、偏微分方程、本征方程、非齊次線性積分方程和非線性方程等。在數(shù)字仿真方面，常用于核系統(tǒng)臨界條件模擬、反應(yīng)堆模擬以及實(shí)驗(yàn)核物理、高能物理、統(tǒng)計物理、真空、地震、生物物理和信息物理等領(lǐng)域?！?.l蒙特卡羅方法的基礎(chǔ)知識8.1.1為了對MC方法有一點(diǎn)初步的認(rèn)識，請先看使用MC方法的幾個例子。蒲豐投針問題：蒲豐（Buffon-法國著名數(shù)學(xué)家）在1777年發(fā)現(xiàn)隨機(jī)投針的概率與無理數(shù)之間的關(guān)系．這個問題是說，若在平面上畫有距離為a的平行線束，向平面上投擲長為的針，試求針與一平行線相交的概率。這個問題的解法如下：以M表示落下后針的中點(diǎn)，表M與最近一平行線的距離，表針與此線的交角，見上圖?？梢?，這兩式?jīng)Q定平面上一矩形R；為了使針與一平行線（這線必定是與針中點(diǎn)M最近的平行線）相交，充分而且必要條件是這個不等式?jīng)Q定R中一個子集G。因此，我們的問題等價于向R中均勻分布地擲點(diǎn)而求點(diǎn)落于G中的概率P.根據(jù)概率的幾何意義，得此式提供了求值的一個方法：可以通過投針事件求得針與平行線相交概率P，求得值：（8.1）若投擲次數(shù)為m，針與平行線相交的次數(shù)為n，那么即于是，可用投針試驗(yàn)來求無理數(shù)的近似值．下表列舉了歷史上若干學(xué)者投針試驗(yàn)計算值的結(jié)果：試驗(yàn)者年份投針次數(shù)的計算值Wolf185050003.1596Smith185532043.1553Fox189411203.1419Lazzarini190134083.1415929射擊問題（打靶游戲）：設(shè)表示射擊運(yùn)動員的彈著點(diǎn)到靶心的距離，表示擊中處相應(yīng)的得分?jǐn)?shù)（環(huán)數(shù)），分布密度函數(shù)表示該運(yùn)動員的彈著點(diǎn)分布，它反映運(yùn)動員射擊水平。積分（8.2）表示這個運(yùn)動員的射擊成績。用概率語言說，就是隨機(jī)變量的數(shù)學(xué)期望值，記為?，F(xiàn)在，假設(shè)這個射擊運(yùn)動員射擊N次，彈著點(diǎn)依次是環(huán)數(shù)分別為，則自然地認(rèn)為N次射擊得分的平均值（8.3）這個平均值相當(dāng)好地代表了這個射擊運(yùn)動員的成績。換句話說，是積分（8.2）式的一個估計值（或近似值）。這個例子通常稱為打靶游戲，它直觀地說明了蒙特卡羅方法計算定積分的基本思想。為進(jìn)一步闡明這個思想，我們再舉個例子：計算積分（8.4）直觀上，就是在邊長為1的正方形里隨機(jī)投點(diǎn)，當(dāng)點(diǎn)落在曲線下面，對積分值有“貢獻(xiàn)”，否則對積分值無“貢獻(xiàn)”。為此，假設(shè)向這個邊長為1的正方形里隨機(jī)投點(diǎn)N次，點(diǎn)落在曲線下面n次，則（8.4）式積分值近似為來近似。從上述幾個例子可以看到，當(dāng)所要求解的問題是某種事件出現(xiàn)的概率，或者是某個隨機(jī)變量的期望值時，它們可以通過某種“試驗(yàn)”的方法，得到這種事件出現(xiàn)的頻率，或者這個隨機(jī)變量的平均值，并用它們作為問題的解。這就是蒙特卡羅方法的基本思想。所以用蒙特卡羅方法求解問題時，首先要建立一個隨機(jī)模型，然后要構(gòu)造一系列的隨機(jī)數(shù)用以模擬這個過程，最后要作統(tǒng)計性的處理。關(guān)于建立隨機(jī)模型，因問題而異。8.1.2隨機(jī)變量及其分布在一定條件下發(fā)生的事件分為必然事件（必然發(fā)生）、不可能事件（恒不發(fā)生）和隨機(jī)事件（可能發(fā)生也可能不發(fā)生），事件發(fā)生的可能性大小用概率表示。必然事件的發(fā)生概率為1，不可能事件的概率為零；隨機(jī)事件發(fā)生的概率。由于測量的隨機(jī)誤差和物理現(xiàn)象本身的隨機(jī)性。一次測量得到的某個值是隨機(jī)的，因此，實(shí)驗(yàn)觀測的物理量是隨機(jī)變量，被研究的物理問題是一個隨機(jī)事件．通常，描寫隨機(jī)事件A發(fā)生的概率用來表示，顯然，．經(jīng)常碰到的隨機(jī)變量有兩類：一類是離散型隨機(jī)變量，這種隨機(jī)變量只能取有限個數(shù)值，能夠—一列舉出來；另一類是連續(xù)型隨機(jī)變量，這種隨機(jī)變量的可能值連續(xù)的分布某個區(qū)間．離散型隨機(jī)變量可用分布列：（8.5）來描寫，它表示取值的概率為（）．即（8.6）分布列描寫了離散型隨機(jī)變量的概率分布．連續(xù)型隨機(jī)變量可能取的值是不可列的，因此不能象離散型隨機(jī)變量那樣，用分布列來描寫它的概率分布，而要用概率分布密度函數(shù)來描寫．考慮連續(xù)型隨機(jī)變量落在區(qū)間內(nèi)的概率，如果極限（8.7）存在，則函數(shù)描寫了在點(diǎn)的概率密度，把叫做隨機(jī)變量的概率分布密度，簡稱分布密度或密度函數(shù)．于是，隨機(jī)變量落在內(nèi)的概率可寫為（8.8）顯然，上式只有當(dāng)可積時才有意義。此外，還可定義隨機(jī)變量的分布函數(shù)來描寫隨機(jī)變量的概率分布規(guī)律．分布函數(shù)一定義為（8.9）即分布函數(shù)在處的值等于隨機(jī)變量取值小于或等于這樣一個隨機(jī)事件的概率．顯然，．對于離散型隨機(jī)變量，分布函數(shù)為階梯函數(shù)（8.10）對于連續(xù)型隨機(jī)變量，分布函數(shù)與分布密度滿足如下關(guān)系：所以分布密度函數(shù)和分布函數(shù)的性質(zhì)：①②③對于任意實(shí)數(shù)：④若在點(diǎn)連續(xù)，則有為了描述隨機(jī)變量的數(shù)學(xué)特征還需引進(jìn)數(shù)學(xué)期望和方差的概念．對于離散型隨機(jī)變量其數(shù)學(xué)期望定義為（8.11）方差定義為（8.12）對連續(xù)型隨機(jī)變量，若已知分布密度，則數(shù)學(xué)期望是（8.13）方差（8.14）在實(shí)際問題中，對于一組N個實(shí)驗(yàn)觀測數(shù)據(jù)，相應(yīng)于某一個隨機(jī)變量的一個樣本，可以用直方圖來形象地表示樣本中數(shù)據(jù)分布的規(guī)律性．先將隨機(jī)變量的取值范圍劃分為若干個區(qū)間，將落入每一區(qū)間的數(shù)據(jù)的個數(shù)m（稱為頻數(shù)）與隨機(jī)變量的取值區(qū)間之間的關(guān)系畫成階躍曲線，即構(gòu)成了直方圖．?dāng)?shù)m／N是觀測數(shù)據(jù)落入這個小區(qū)間的頻率．直方圖的根坐標(biāo)為隨機(jī)變量的取值范圍，縱坐標(biāo)為頻數(shù)．圖表示的是某一個隨機(jī)事件的直方圖。8.1.3概率論中的大數(shù)定理和中心極限定是蒙特卡羅方法的數(shù)學(xué)基礎(chǔ)。大數(shù)定理：設(shè)為一隨機(jī)變量序列，獨(dú)立同分布，數(shù)學(xué)期望值存在，則對任意，有（8.15）是算術(shù)平均。大數(shù)定理指出，當(dāng)時，算水平均收斂到數(shù)學(xué)期望（或統(tǒng)計平均）．也就是說，可以用算術(shù)平均代替數(shù)學(xué)期望。對于收斂的程度，和誤差估計，要用到下面的中心極限定理．中心極限定理：設(shè)為一隨機(jī)變量序列，獨(dú)立同分布，數(shù)學(xué)期望為，方差，則當(dāng)時，（8.16）利用中心極限定理，當(dāng)很大時，可得到誤差（8.17）成立的概率為。通常將稱為置信度，為置信水平。的定義為（8.18）和的關(guān)系可以計算得到，下面給出常用的幾組和的值：0.50.050.020.010.67451.96002038632.5758從（8.17）式可以看到，算術(shù)平均值收斂到的階為，可見，蒙特卡羅方法收斂的階很低，收斂速度很慢。當(dāng)，誤差稱為概率誤差。誤差由和決定．在固定的情況下，要提高1位精確度，就要增加100倍試驗(yàn)次數(shù)．相反，若減小10倍，就可以減少100倍工作量．因此，控制方差是應(yīng)用蒙特卡羅方法中很重要的一點(diǎn)．如果要求置信水平為0.95（95%），計算得，則說明誤差估計的概率是0.95?！?.2隨機(jī)數(shù)和隨機(jī)抽樣在計算機(jī)上用蒙特卡羅方法模擬某過程時，需要產(chǎn)生具有各種概率分布的隨機(jī)變量．最簡單和最基本的隨機(jī)變量就是區(qū)間上均勻分布的隨機(jī)變量．這些隨機(jī)變量的抽樣值就稱為隨機(jī)數(shù)．以后談到隨機(jī)數(shù)，如果不加特別說明，專指區(qū)間上均勻分布的隨機(jī)數(shù)．各種其他分布的隨機(jī)變量的抽樣值都可借助均勻分布的隨機(jī)數(shù)得到．真正的隨機(jī)數(shù)如同擲骰子那樣，產(chǎn)生1-6范圍內(nèi)的隨機(jī)數(shù)整數(shù)。抽獎用的搖號碼機(jī)則可產(chǎn)生0－9范圍內(nèi)的隨機(jī)整數(shù)。這些真正的隨機(jī)數(shù)除統(tǒng)計規(guī)律外無任何其它規(guī)律可循。偽隨機(jī)數(shù)，或稱贗隨機(jī)數(shù)，是指按照某種算法可以給出的似乎隨機(jī)地出現(xiàn)的數(shù)。既然數(shù)的給出是按某種算法，也就是按某種規(guī)律，那這種隨機(jī)數(shù)就必然具有一定的周期。設(shè)其周期為n，則第n＋l個數(shù)就等于第一個數(shù)，此后均依次重復(fù)出現(xiàn)。當(dāng)然，如果周期n足夠大，可使在整個使用過程中不至于表現(xiàn)出其周期性，偽隨機(jī)數(shù)也是實(shí)用的。例如，計算機(jī)中的偽隨機(jī)數(shù)發(fā)出器要求其周期大于計算機(jī)的記憶單元數(shù)。此外偽隨機(jī)數(shù)的統(tǒng)計性質(zhì)是表征隨機(jī)數(shù)品質(zhì)的又一重要指標(biāo)。均勻分布的隨機(jī)數(shù)，既要求數(shù)的出現(xiàn)是隨機(jī)的，又要求數(shù)的分布是均勻的。至于如何評估隨機(jī)數(shù)分布的均勻程度，將在本節(jié)稍后討論。8.2.1均勻分布的隨機(jī)數(shù)的產(chǎn)生有許多方法，通常采用乘同余法。例如，可按下列公式產(chǎn)生隨機(jī)數(shù)（8.19）或?qū)懗善渲袨榻o定常數(shù)。給出后，就可以用式（8.19）依次給出等一系列隨機(jī)數(shù)。如何確定常數(shù)和，這是個十分關(guān)鍵的問題，是人們?nèi)栽诓粩嘌芯刻剿鞯恼n題。下面僅給出確定和的一般原則。關(guān)于的取值一般取，其中m為計算機(jī)中二進(jìn)制數(shù)的字長，則為計算機(jī)所能表示的最大整數(shù)。例如字長16位時，可取等。關(guān)于的取值，一般取，其中M為任一正整數(shù)。例如有取等。建議取，這樣統(tǒng)計性較好。關(guān)于的取值，一般取為奇數(shù)?？梢则?yàn)證，當(dāng)為奇數(shù)時，周期是T，其它參數(shù)不變，當(dāng)為偶數(shù)時，周期則為T/2.例如設(shè)，得隨機(jī)數(shù)序列為：周期為16。按（8.19）產(chǎn)生的隨機(jī)數(shù)其值域?yàn)?。如果要產(chǎn)生之間的隨機(jī)數(shù)，只需將原產(chǎn)生的每個隨機(jī)數(shù)再除以即可。用類似的方法，經(jīng)過簡單的變換，就可產(chǎn)生任何所需值域的贗隨機(jī)數(shù)序列。例如在16位微機(jī)上，可用下式（8.20），并令，來產(chǎn)生之間的贗隨機(jī)數(shù)序列(8.20)見程序：rand1.f90下面給出另一個產(chǎn)生之間均勻隨機(jī)數(shù)程序：rand.for,plot_rand.mimplicitdoubleprecision(a-h,o-zimplicitdoubleprecision(a-h,o-z) realR open(2,file='rand.dat') ix=32765 doi=1,100 callrand(ix,r) write(*,*)'r=',r write(2,*)i,r enddo ENDsubroutinerand(ix,r) i=ix*259 ix=i-i/32768*32768 r=float(ix)/32768 return end8.2.2一個好的隨機(jī)數(shù)發(fā)生器或一個好的隨機(jī)數(shù)生成程序必須滿足兩個條件：第一，所生成的隨機(jī)數(shù)序列應(yīng)當(dāng)具有足夠長的周期；第二，所生成的隨機(jī)數(shù)序列應(yīng)當(dāng)具有真正隨機(jī)數(shù)序列所具有的統(tǒng)計性質(zhì)。其周期的長短比較容易測試和判斷，而統(tǒng)計性質(zhì)的優(yōu)劣則不那么簡單。通常對統(tǒng)計性質(zhì)的檢驗(yàn)方法是采用頻數(shù)分布檢驗(yàn)：對于一個均勻分布的隨機(jī)數(shù)發(fā)生器，設(shè)所產(chǎn)生隨機(jī)數(shù)序列的值域?yàn)?，則所產(chǎn)生的隨機(jī)數(shù)字應(yīng)與從之間均勻的頻數(shù)分布相一致。為了檢驗(yàn)頻數(shù)分布情況，可按畫統(tǒng)計直方圖的方法，將整個值域分成M個寬度相等的子區(qū)間，設(shè)是第區(qū)間內(nèi)出現(xiàn)的隨機(jī)數(shù)的個數(shù)，即第i子區(qū)間的頻數(shù)，則所有子區(qū)間中隨機(jī)數(shù)個數(shù)的平均頻數(shù)第i子區(qū)間頻數(shù)的偏差則頻數(shù)的方均根偏差如果所產(chǎn)生的N個隨機(jī)數(shù)均勻分布于整個值域，則且在任一子區(qū)間內(nèi)出現(xiàn)個數(shù)字的幾率服從高斯分布規(guī)律，即其中為標(biāo)準(zhǔn)偏差?？梢姡粋€均勻分布的隨機(jī)數(shù)序列，其最可幾頻數(shù)應(yīng)為，而其頻數(shù)在范圍內(nèi)的幾率應(yīng)為0.68，其頻數(shù)的方均根偏差應(yīng)接近于標(biāo)準(zhǔn)偏差。通常檢查均勻分布隨機(jī)數(shù)分布的均勻程度就可從這兩個方面考核，即各頻數(shù)是否接近于平均頻數(shù)和頻數(shù)的方均根偏差是否接近于標(biāo)準(zhǔn)偏差。8.2.3隨機(jī)抽樣的方法很多，在計算機(jī)上實(shí)現(xiàn)時要考慮運(yùn)算量的大小，也就是所謂“抽樣費(fèi)用”．因?yàn)閼?yīng)用蒙特卡羅方法求解一個物理問題時，大量的計算時間將用于隨機(jī)抽樣，所以隨機(jī)抽樣方法的選取往往決定算題的費(fèi)用．但對不同問題和不同機(jī)器，不同的方法也可以有不同的評價．下面介紹幾種常用的隨機(jī)抽樣方法，供讀者選擇．1．連續(xù)型分布的直接抽樣法利用[0,1]區(qū)間上的均勻分布隨機(jī)數(shù)可以產(chǎn)生具有給定分布的隨機(jī)變量數(shù)列．我們知道，若隨機(jī)變量具有分布密度，則其分布函數(shù)（8.21）可知，為單調(diào)非減函數(shù)，而且存在反函數(shù)?？梢宰C明，這個分布就是[0,1]區(qū)間上滿足均勻分布的隨機(jī)變量。令均勻分布的隨機(jī)數(shù)等于分布函數(shù)，則得到：就是具有分布密度的隨機(jī)變量的抽樣。事實(shí)上：由分布函數(shù)定義如果是[0,1]區(qū)間均勻分布的隨機(jī)變量其中用到了[0,1]區(qū)間均勻分布的隨機(jī)變量的分布密度函數(shù)抽樣的步驟為：給定分布密度；計算；產(chǎn)生隨機(jī)數(shù)；④計算令重復(fù)③，④。※================================================================※例題8.1求經(jīng)常用來描述電子元件的穩(wěn)定時間，系統(tǒng)的可靠性和粒子游動的自由程等的指數(shù)密度分布的隨機(jī)抽樣。解：，令，則由于與同分布，故可得指數(shù)分布的抽樣公式※================================================================※例題8.2求均勻密度函數(shù)在區(qū)間[a,b]上均勻分布的隨機(jī)數(shù)解：令[0,1]區(qū)間的均勻分布的隨機(jī)數(shù)得區(qū)間[a,b]上均勻分布的隨機(jī)數(shù)※================================================================※2．離散型隨機(jī)變量的抽樣法設(shè)是離散型隨機(jī)變量，分別以概率取值，可按下面步驟抽樣：選取隨機(jī)數(shù)；確定滿足不等式的值（這里約定）；對應(yīng)的就是所求的抽樣值，即；事實(shí)上，如果設(shè)：，則由※================================================================※例題8.3按學(xué)生的成績抽樣，某班26名學(xué)生成績分布如下：分?jǐn)?shù)擋1-1011-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上面，解：這個分布抽樣方法如下：產(chǎn)生隨機(jī)數(shù)（由程序計算），假設(shè)是由不等式確定,由，確定抽樣值為分計算程序：scores.f90,plot_scores.m※================================================================※例題8.4設(shè)能量為的光子與物質(zhì)相互作用的產(chǎn)生光電效應(yīng)的截面為，產(chǎn)生康普頓散射截面，電子對效應(yīng)截面為，確定相互作用類型抽樣解：設(shè)總截面：，抽樣方法如下：產(chǎn)生隨機(jī)數(shù)；計算，如果，則發(fā)生光電效應(yīng)；否則計算，如果，則發(fā)生康普頓散射；否則，發(fā)生電子對效應(yīng)※================================================================※8.2.4用蒙特卡羅方法處理的問題可以分為兩類：一類是隨機(jī)性問題，例如中子在介質(zhì)內(nèi)的傳播問題和后面要介紹的原子核裂變問題等．對于這一類問題，通常采用直接模擬方法．首先，必須根據(jù)物理問題的規(guī)律，建立一個概率模型（隨機(jī)向量或隨機(jī)過程），然后用計算機(jī)進(jìn)行抽樣試驗(yàn)，從而得出對應(yīng)于這一物理問題的隨機(jī)變量的分布。假定隨機(jī)變量y是我們的研究對象，它是m個互相獨(dú)立的隨機(jī)變量的函數(shù)，如果概率分布密度為，則用蒙特卡羅方法計算的基本方法是：在計算機(jī)上用隨機(jī)抽樣的方法根據(jù)概率分布密度抽樣，產(chǎn)生N組隨機(jī)變量，計算的值，用這樣的樣本分布來近似y的函數(shù)，由此可計算出這些量的統(tǒng)計值．由上可知，蒙特卡羅方法的計算過程就是用數(shù)學(xué)方法模擬實(shí)際的物理過程，它主要是在計算機(jī)上產(chǎn)生已知分布的隨機(jī)變量樣本．以代替昂貴的甚至難以實(shí)現(xiàn)的實(shí)驗(yàn)．蒙特卡羅方法又被看作是用計算機(jī)來完成物理實(shí)驗(yàn)的一種方法．另一類是確定性問題．在解決確定性問題時，首先要建立一個有關(guān)的概率統(tǒng)計模型。使所求的解就是這個模型的概率分布或數(shù)學(xué)期望，然后對這個模型作隨機(jī)抽樣，最后用其算術(shù)平均值作為所求解的近似值。§8.3M8.3.1考慮方程（8.3.1）其中是定義在[a,b]區(qū)間實(shí)的連續(xù)函數(shù)。如果已知則方程（8.3.1方法一：用頻率近似概率先討論在[a,b]區(qū)間內(nèi)只有一個根的情況。（1）當(dāng)，則a就是所求的根；（2）否則，設(shè)是[a,b]區(qū)間內(nèi)的一個根，是[a,b]區(qū)間內(nèi)均勻分布的隨機(jī)數(shù)，則落在內(nèi)的概率是因此，在區(qū)間[a,b]內(nèi)均勻投N個點(diǎn)，設(shè)其中有M個落在根的左側(cè)，于是有（8.3.2）這樣，用M-C求單根的步驟如下：①定義隨機(jī)變量，當(dāng)?shù)趇次試驗(yàn)（投點(diǎn)）時；，其他情況；其中是(0,1)區(qū)間上均勻分布的隨機(jī)數(shù)。②令：③當(dāng)時，當(dāng)時，對于區(qū)間[a,b]內(nèi)有多個根時，可用分析的方法確定每個單根的界限，在每個界限內(nèi)應(yīng)用上述方法。或者從a點(diǎn)出發(fā)，以h為基本步長向前跨長為h的小區(qū)間（h要選的合適使每個區(qū)間至少有一個根。太大容易丟根，太小浪費(fèi)時間），當(dāng)跨的小區(qū)間兩端的函數(shù)同號時，繼續(xù)向前跨；異號時，用上述方法求出小區(qū)間中的根?！?===============================================================※例題8.5-1求方程的全部實(shí)根解：見計算程序:root_MC_1.f90※================================================================※方法二：取離根左右最近的兩值平均即?。簞t：，是根的一個近似估計?！?===============================================================※例題8.5-2：求方程在區(qū)間[0，]內(nèi)的一個實(shí)根。解：見計算程序:root_MC_2.f90※================================================================※8.3.設(shè)區(qū)間[a,b]的中的隨機(jī)變量的分布由密度函數(shù)給定．是區(qū)間[a,b]上的連續(xù)函數(shù)，數(shù)學(xué)期望（8.3.3存在，則積分（8.3.3設(shè)隨機(jī)變量按密度分布抽樣的一系列值為相應(yīng)產(chǎn)生的隨機(jī)抽樣值，則根據(jù)大數(shù)定理有可見，只要n充分大，積分（8.3.(8.3.4現(xiàn)在來討論積分的值．為此，選擇某種密度分布函數(shù)滿足且能很方便地生成具有密度函數(shù)為的隨機(jī)抽樣．同時將積分于是歸結(jié)為積分（8.3.1）的形式。在多數(shù)情況下，往往取為[a,b]上均勻分布的概率密度，現(xiàn)在，在區(qū)間（a,b）上均勻分布的隨機(jī)數(shù)中抽取，計算，然后計算平均值（1）方法一在[0,1]之間取均勻分布的隨機(jī)數(shù)序列，并令得到區(qū)間均勻分布的隨機(jī)數(shù)，只要足夠大，則有※================================================================※例題8.6用蒙特卡羅方法計算定積分解：計算程序MC_1.for及結(jié)果如下：※================================================================※對于多重積分等復(fù)雜問題，基本方法都是相同的。例如計算多重定積分取[0,1]區(qū)間均勻分布的隨機(jī)數(shù)序列：計算：只要足夠大，則有（2）方法二現(xiàn)在仍然討論積分假設(shè)函數(shù)在區(qū)間[a,b]上有最大值，的值。為此，在下圖作矩形ABCD，其寬度為，高為，則矩形面積,給出兩個隨機(jī)數(shù)，計算得如果滿足成立，隨機(jī)點(diǎn)落在矩形陰影區(qū)域內(nèi)。設(shè)總共產(chǎn)生的隨機(jī)點(diǎn)數(shù)為N,落在矩形陰影區(qū)域內(nèi)的點(diǎn)數(shù)為M,當(dāng)N足夠大時，定積分※================================================================※【例題8.7】：計算半徑為1圓的面積(計算值)解：見計算程序：pi_1.f,pi_2.f,pi_3.f90※================================================================※【例題8.8】計算多重積分（這個積分的結(jié)果是歐拉常數(shù)）解：（1）在單位超立方體區(qū)域上構(gòu)造一個聯(lián)合概率分布密度：其中當(dāng)，其它：則積分可以寫為（2）用算術(shù)平均值來近似數(shù)學(xué)期望。現(xiàn)從分布密度抽取隨機(jī)向量的N個樣本，每個由s個[0,1]區(qū)間的均勻分布的隨機(jī)數(shù)組成，即N個樣本為:就是積分的一個近似估計。下面是模擬程序：Euler.f90gama(5)=-0.57557088,gama(6)=-0.57453740gama(7)=-0.55329597,gama(8)=-0.56524354gama(9)=-0.57235909可見積分結(jié)果與積分重數(shù)沒有多大關(guān)系。===============================================================※【例題8.9】一球體半徑，球上有一半徑的圓柱形空洞，其軸線與球的直徑重合。試用M-C方法求實(shí)體的體積。解：如右圖所示，令球體中心位于坐標(biāo)系的原點(diǎn)O處，作邊長為的正方體，其中心與球心重合，則正方體的體積;產(chǎn)生一組三個隨機(jī)數(shù)（，它們的值域均為;然后判斷該隨機(jī)點(diǎn)是否位于實(shí)體內(nèi)，其判據(jù)是:且若共產(chǎn)生了N組隨機(jī)數(shù)，而滿足上述判據(jù)的有M組，則球體的體積為了提高測量精度，可重復(fù)上述計算過程，例如用完全相同的方法再做兩遍，取三次結(jié)果的平均值作為最終結(jié)果?！?===============================================================※對于一般（變積分限）的多重定積分※================================================================※【例題8.10】計算函數(shù)積分限的多重積分解：該積分的區(qū)域見圖9－2。由圖9－2可見，該積分x軸積分限從，y軸積分限是不規(guī)則的，由，在常數(shù)區(qū)間[]產(chǎn)生一組隨機(jī)數(shù)序列在函數(shù)區(qū)間[]產(chǎn)生一組隨機(jī)數(shù)序列利用計算多（s）重積分的平均值方法：（對于常數(shù)積分限）是s維空間積分區(qū)域的超體積（見下面說明）。但是要注意，對于，積分限變化的情況，例如二維情況，對于中關(guān)于y是隨x變化的部分要放到求和里面，因?yàn)?，對于此時的，的變化區(qū)間變化了，在本文題中是，所以要先計算然后再乘上x軸區(qū)間寬度（）=1計算程序：/*mc2.c*/，計算結(jié)果為：，積分值精確到小數(shù)點(diǎn)后第6位是589/90=6.544444?？梢奙C方法計算函數(shù)積分限的二重積分，僅用15000個隨機(jī)數(shù)，結(jié)果就相當(dāng)滿意了?！?===============================================================※【例題8.11】計算程序：mc3.c,計算結(jié)果為：8.3.3MC方法求解Laplace方程利用蒙特卡羅的隨機(jī)游動方法求解二維Laplace方程：其差分格式為：設(shè)Q是邊界上的點(diǎn)，P是區(qū)域內(nèi)的點(diǎn)，求的M-C方法為隨機(jī)游動方法，步驟如下：粒子的狀態(tài)參數(shù)為；粒子由狀態(tài)出發(fā)；假設(shè)粒子已游動n步達(dá)到，從此點(diǎn)出發(fā)，以1/4等概率到達(dá)四個鄰點(diǎn)的一個，記為；③若是邊界點(diǎn)，則游動終止，并記下邊界處的函數(shù)值④；若不是邊界點(diǎn)，則重復(fù)②和③，直至達(dá)到邊界。由狀態(tài)出發(fā)重復(fù)m次，這樣產(chǎn)生粒子游動的m次歷史，得然后，從其它點(diǎn)出發(fā)，求其它點(diǎn)的值。這種方法對于復(fù)雜邊界情況特別適用，特別是求區(qū)域內(nèi)任意點(diǎn)的值很有效。※================================================================※【例題8.12】利用蒙特卡羅方法求解例題7.3的問題解：計算程序MC_laplace.for,plot_MC_laplace.m※================================================================※8.3.4放射性物質(zhì)的鏈?zhǔn)椒磻?yīng)是一個隨機(jī)過程，可借助計算機(jī)用M－C方法模擬和研究。由原子核物理知識可知，的原子核本質(zhì)上是不穩(wěn)定的，會自發(fā)地發(fā)生裂變。裂變的激烈程度可用放射性物質(zhì)的半衰期來描述，半衰期是指大量核中有1／2的核發(fā)生裂變所需要的時間。的半衰期為7億多年。因此任何時刻發(fā)生裂變的核只是相對很小部分，其釋放的能量只能使其本身微微溫?zé)?。但是，在一定條件下，自發(fā)裂變放出的兩個中子轟擊其它核而被吸收，引起新的裂變而放出更多的中子，這更多的中子又引起新一輪更多的裂變，依次類推，可迅速釋放出大量能量，甚至引起爆炸，這就是鏈?zhǔn)椒磻?yīng)。設(shè)開始有N個核發(fā)生裂變，每個核放出兩個中子，稱為第一代中子，共2N個。2N個中子又感生新一輪裂變，產(chǎn)生第2代中子，為4N個。……如此進(jìn)行下去，直至第n次裂變，產(chǎn)生第n代中子為個。按此計算，30代可產(chǎn)生裂變的核數(shù)為億，即為第一次裂變核數(shù)的10億倍?，F(xiàn)在的問題是，在什么條件下才能發(fā)生鏈?zhǔn)椒磻?yīng)呢？其基本要求是裂變所產(chǎn)生的兩個中子中至少有一個能使第二個鈾核發(fā)生裂變。為此要求核材料中雜質(zhì)的含量，包括的含量應(yīng)足夠少，以避免中子被和其它雜質(zhì)所吸收。另外，由于熱中子使裂變的機(jī)會很大，所以在鈾堆中還必須加入減速劑，如重水或石墨等，以使快中子減速到熱中子。最后，非常重要的條件是鈾堆的體積必須足夠大，以避免裂變所放出的中子過多地未與鈾核相遇而飛出鈾體外。這就涉及到臨界體積和臨界質(zhì)量的概念。所謂臨界質(zhì)量是指可裂變物質(zhì)能發(fā)生自持鏈?zhǔn)椒磻?yīng)的最小質(zhì)量。由于鈾核體積很小，一鈾核裂變放出的中子在和另一鈾核作用并使之發(fā)生裂變之前，平均地說要經(jīng)過一定相對很長的距離，約為厘米數(shù)量級。因此，假定有個核發(fā)生自發(fā)裂變而放出2個中子，其中N個中子在鈾塊中引起另外的核發(fā)生裂變，其余的中子未與其它核碰撞而飛出鈾塊。為描述一次裂變能引起下一次裂變的可能程度，定義裂變過程的倍增系數(shù)：不難理解，維持自持鏈?zhǔn)椒磻?yīng)的條件是：，即。倍增系數(shù)k＝l是臨界質(zhì)量Mc的條件。k的值與前面論及的諸多因素有關(guān)，本節(jié)將只限于討論k與鈾塊的質(zhì)量和形狀有關(guān)的問題，用計算機(jī)程序來模擬具有一定大小和形狀的鈾塊中大量隨機(jī)的裂變過程，統(tǒng)計算出相應(yīng)的倍增系數(shù)k。設(shè)鈾塊為長方體，發(fā)生裂變的鈾核位于鈾塊內(nèi)隨機(jī)點(diǎn)處，如下圖所示。隨機(jī)點(diǎn)坐標(biāo)的值域?yàn)椋烘準(zhǔn)椒磻?yīng)模擬示意圖該核子裂變反應(yīng)產(chǎn)生兩個中子，其運(yùn)動方向可以用兩個角坐標(biāo)和來描述，見上圖（b）。釋放出的每一個中子按飛向各個方向的幾率均等來考慮，或者說中子飛行方向的幾率是按以為頂點(diǎn)的立體角均勻分布的。立體角元可表示為可見，按立體角均勻分布是按角均勻分布和按均勻分布，而并非按角均勻分布。因此，對應(yīng)的兩個隨機(jī)數(shù)的值域?yàn)橄旅娣e算按極角正余弦分布和方位角分布的抽樣公式：各向同性散射角余弦分布和散射方位角分布按分布密度：和抽樣的方法就是采用直接抽樣：平均地說，能否擊中另一個核只取決于中子在鈾塊體內(nèi)飛行的距離。假設(shè)在0至1cm距離之間經(jīng)過任何一段相同距離擊中另一鈾核的幾率均等?；蛘哒f，中子在擊中另一鈾核之前飛行的距離為0－1之間均勻分布的隨機(jī)數(shù)。因此與飛行距離相應(yīng)的隨機(jī)數(shù)為由此可計算出被擊中的鈾核的位置最后，檢查計算出的碰撞點(diǎn)是否位于鈾塊體內(nèi)。若在鈾塊體內(nèi)，則累計引起新裂變中子數(shù)N。按照上述原則，歸納計算k的具體步驟為：1．給定鈾塊質(zhì)量M、鈾塊邊長比和用于計算k的隨機(jī)自發(fā)裂變核的個數(shù)，即舊裂變核個數(shù)，并設(shè)所選約化單位可使鈾塊的密度為1，體積為V，則得或2．產(chǎn)生九個之間的隨機(jī)數(shù)：3．舊裂變核位置：4.舊裂變放出的兩個中子的方向5.中子的飛行距離6.可能發(fā)生新裂變的位置7.檢查上述位置，是否在鈾塊體內(nèi)。如果均滿足，則的值增加1；同樣，如果均滿足，則的值也增加1；8.重復(fù)步驟2-7，共執(zhí)行次，然后計算：若計算結(jié)果，則調(diào)整M和s的值后再進(jìn)行上述步驟，直至，此時M即為臨界質(zhì)量。顯然臨界質(zhì)量與s有關(guān)，與核材料的形狀有關(guān)。計算程序：MC_2.f90和plot_k_m.m8.3.5原子反應(yīng)堆的外壁是鉛板圍成的。中子從鉛壁的內(nèi)側(cè)向外輻射，稱中子貫穿概率問題。可以用MC方法模擬中子貫穿鋁板的概率。為簡化問題，假定中子源用一個厚5個單位的鉛板圍住，中子貫穿方向?yàn)閄方向，Y方向是無限高。實(shí)驗(yàn)表明，中子在10次相撞后能量消耗殆盡即被鉛原子吸收。中子進(jìn)入鉛板，走一定距離與鉛原子核碰撞（鉛原子直徑d），之后隨機(jī)改變方向。又走一定距離，與另一個原子核碰撞。如圖9－4所示，如此經(jīng)過多次碰撞后，中子可能穿透鉛壁輻射到反應(yīng)堆外，也可能將其能量耗盡被鉛壁吸收，還可能被反射回反應(yīng)堆內(nèi)。顯然，鉛壁設(shè)計得越厚，穿透的概率就越小，反應(yīng)堆就越安全。由于每次碰撞后彈出的角度是隨機(jī)的，因此對一個中子來說是穿貫，還是被吸收或返回，也是隨機(jī)的。中子貫穿概率問題是由大量中子運(yùn)動的統(tǒng)計規(guī)律決定的。要得到鉛壁厚和穿透率之間的關(guān)系以及貫穿概率，用解析方法是極為困難的，而用MC方法模擬中子貫穿鉛板問題的求解大大簡化。中子在壁內(nèi)的運(yùn)動與其每次與鉛原子核碰撞后的散射角有關(guān)，由三角學(xué)可知，取中子每次的自由程為斜邊1，直角三角形中直角邊是。注意隨機(jī)散射角為［］，當(dāng)和時，中子在鉛壁厚度方向上前進(jìn)單位距離，當(dāng)時，中子后退單位距離。在［］區(qū)間產(chǎn)生一個隨機(jī)數(shù)作為每次碰撞隨機(jī)散射角其中：rand()是[0，1]區(qū)間的隨機(jī)數(shù)，中子在x方向每次運(yùn)動的距離（初始s=0）模擬一個中子在鉛板中的運(yùn)動得到一個結(jié)果，對大量中子的運(yùn)動進(jìn)行模擬的結(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)等距離放入一個個的銅片（厚度），求整個柱筒對距離筒軸處觀察點(diǎn)o’的放射性貢獻(xiàn)，即計算積分式中的表示由點(diǎn)出發(fā)到點(diǎn)o’間經(jīng)過圓柱筒中銅的總長度。是衰減常數(shù)。積分區(qū)域是圓柱筒中銅所占的幾何區(qū)域。解：我們知道，在滿足一定的條件下，放射源與觀察點(diǎn)是可以進(jìn)行倒易的。這就是通常的光學(xué)倒易定理。本文題與能量無關(guān)，倒易定理成立。根據(jù)光學(xué)倒易定理，可以用點(diǎn)源對圓柱筒中銅的通量貢獻(xiàn)來代替圓柱源對點(diǎn)通量的貢獻(xiàn)。即將O’點(diǎn)看作放射性點(diǎn)源，將放射性銅看作是觀察者。，表示從點(diǎn)O’沿方向經(jīng)過圓柱筒中銅的總長度。表示各向同性散射角分布密度；表示散射方位角分布密度；計算程序：MC_3.for,計算結(jié)果：8.3.61.隨機(jī)行走問題定義：