非常重要和有趣的計(jì)算方法——蒙特卡洛方法

本文將向大家介紹一種在科學(xué)研究中非常重要和有趣的計(jì)算方法——蒙特卡洛方法，這種方法在數(shù)學(xué)、物理學(xué)、化學(xué)、工程、經(jīng)濟(jì)學(xué)、環(huán)境動(dòng)力學(xué)等多個(gè)領(lǐng)域都有廣泛的應(yīng)用。

到底什么是蒙特卡洛方法？我們可以先從它的名字開(kāi)始了解，蒙特卡洛（Monte Carlo）是摩納哥公國(guó)的一座城市，是世界著名的“賭城”。以“蒙特卡洛”來(lái)命名這種計(jì)算方法就是因?yàn)槠浔旧肀闶且环N概率算法，其核心思路是通過(guò)概率實(shí)驗(yàn)所求的概率來(lái)計(jì)算我們感興趣的一個(gè)量。

概率算法的1.0版

為了更好地理解蒙特卡洛方法，我們先簡(jiǎn)單了解一下“蒲豐投針問(wèn)題”，這個(gè)問(wèn)題的提出被認(rèn)為是蒙特卡洛方法的起源。 18世紀(jì)，法國(guó)數(shù)學(xué)家蒲豐提出了一種計(jì)算圓周率π的方法——隨機(jī)投針?lè)ǎ杭僭O(shè)我們有一個(gè)以平行且等距為a的木紋鋪成的地板，隨意拋一支長(zhǎng)度為l（比木紋之間距離小）的針，通過(guò)針和其中一條木紋相交的概率p，即可計(jì)算圓周率π。計(jì)算公式為：

，其中n是投針的總次數(shù)，m是針與平行直線交點(diǎn)的總數(shù)目。

圖 蒲豐投針問(wèn)題示意圖

這個(gè)方法的原理可以通過(guò)概率學(xué)的推導(dǎo)計(jì)算來(lái)進(jìn)行證明。由于投針掉落的位置與方向都是隨機(jī)且獨(dú)立的。我們假定落地后針的中心距最近的地板條紋的距離為X，那么X在[0，a/2]上均勻分布；針與地板條紋的夾角為Y，則Y在[0, π/2]之間均勻分布。當(dāng)

時(shí)，針與木紋相交，因此（X,Y）的概率密度函數(shù)和相交的概率P分別為：

相信大家都聽(tīng)懂了，下面我們可以…… 好吧，這里還有一種雖然不夠嚴(yán)謹(jǐn)、但易于理解的解釋: 我們想象一個(gè)長(zhǎng)度為πa的鐵絲，被繞成了一個(gè)直徑為a的圓環(huán)。那么無(wú)論我們?cè)趺慈舆@個(gè)鐵環(huán)，它與條紋的交點(diǎn)恒為兩個(gè)，因此當(dāng)投針n次后，相交的次數(shù)恒為2n。如果我們把鐵絲拉直再扔，這樣的鐵絲扔下時(shí)與平行線相交的情形要比圓圈復(fù)雜些，可能有4個(gè)交點(diǎn)、3個(gè)交點(diǎn)、2個(gè)交點(diǎn)、1個(gè)交點(diǎn)，甚至于都不相交。由于圓圈和直桿的長(zhǎng)度同為πa，根據(jù)機(jī)會(huì)均等的原理，投擲n次，直桿與平行線組交點(diǎn)的總數(shù)期望也是2n。同時(shí)還有一個(gè)規(guī)律，當(dāng)投針次數(shù)n固定時(shí)，鐵絲的長(zhǎng)度l與交點(diǎn)總數(shù)m應(yīng)為正比關(guān)系，即m=kl。考慮到l=πa時(shí)，m=2n，將

代入前式可得：

。 投針試驗(yàn)既然是依靠概率的算法，那么隨著投針次數(shù)越來(lái)越多，計(jì)算求得的π值也會(huì)越來(lái)越接近于真實(shí)值。下表給出了一些比較出名的投針試驗(yàn)得到的圓周率估計(jì)值，可以看到在投擲數(shù)千次后，計(jì)算得到的圓周率與我們所熟知的π值的誤差仍較大。

表 一些投針試驗(yàn)的計(jì)算結(jié)果

1995年，馬修斯發(fā)表了他如何通過(guò)觀察天空中亮星的分布計(jì)算圓周率。他的試驗(yàn)方法基于一個(gè)基本的原理：任意兩個(gè)自然數(shù)互質(zhì)的概率為。他從眾多星星中選擇100個(gè)亮星，將這些亮星兩個(gè)分成一對(duì)，然后計(jì)算每對(duì)星之間的角距，得出一堆數(shù)據(jù)，然后檢查這些數(shù)據(jù)的因子情況，從中計(jì)算出π值約為3.12772。

從上述兩個(gè)例子來(lái)看，依靠重復(fù)的物理、觀測(cè)等試驗(yàn)行為來(lái)獲取隨機(jī)性數(shù)據(jù)的方法往往很難得到令人滿意的計(jì)算結(jié)果，這主要是受到了樣本數(shù)量的限制。而早在魏晉時(shí)，我國(guó)的劉徽便通過(guò)割圓術(shù)求得了π的近似值3.1416。

“史詩(shī)級(jí)加強(qiáng)”X.0版

20世紀(jì)40年代，美國(guó)“曼哈頓計(jì)劃”的成員S.M.烏拉姆和J·馮·諾伊曼第一次把這種通過(guò)概率事件來(lái)計(jì)算關(guān)注的確定值的方法命名為“蒙特卡洛方法”。隨著電子計(jì)算機(jī)的發(fā)明和科學(xué)技術(shù)的發(fā)展，蒙特卡洛方法得到了“史詩(shī)級(jí)加強(qiáng)”。

圖 S.M.烏拉姆（左）和J·馮·諾伊曼（右）

計(jì)算機(jī)在進(jìn)行蒙特卡洛模擬的過(guò)程中獲取隨機(jī)性最根本的方法是通過(guò)固定算法得到符合[0,1]均勻分布的“偽隨機(jī)數(shù)”，它并不真正的隨機(jī)，但具有類似于隨機(jī)數(shù)的統(tǒng)計(jì)特征，如均勻性、獨(dú)立性等。

這里介紹另一種計(jì)算π值的蒙特卡洛方法——“撒豆法”。該方法假定有無(wú)數(shù)個(gè)豆子被均勻地撒在下圖所示的正方形中，那么落在圓內(nèi)的豆子數(shù)m與落在正方形內(nèi)的豆子總數(shù)n的比值的期望應(yīng)與它們面積的比值一致，即，這樣就可以計(jì)算得到π的值。

圖 “撒豆法”求解π值模型示意圖

利用計(jì)算機(jī)開(kāi)展上述計(jì)算，僅需不到一分鐘的時(shí)間，便可以完成十億次“撒豆”，并得到相應(yīng)的計(jì)算結(jié)果。計(jì)算的python代碼及運(yùn)行結(jié)果見(jiàn)下圖。

圖 計(jì)算所用python代碼及計(jì)算結(jié)果

計(jì)算機(jī)時(shí)代的蒙特卡洛模擬無(wú)疑具有超高的計(jì)算效率，且其計(jì)算效率隨著計(jì)算機(jī)技術(shù)的飛速發(fā)展而不斷提升。

核科學(xué)領(lǐng)域中的應(yīng)用

蒙特卡洛方法在核科學(xué)領(lǐng)域中有著廣泛的應(yīng)用。核物理領(lǐng)域的基本參數(shù)（如反應(yīng)截面、散射發(fā)射角度分布、能譜分布、衰變、衰減等）主要來(lái)自核物理實(shí)驗(yàn)及理論模型等，而蒙特卡洛模擬程序也在與核物理實(shí)驗(yàn)、理論模型等的相互參照、驗(yàn)證和迭代更新中不斷發(fā)展至今。

如今在核物理領(lǐng)域有許多廣泛使用的蒙特卡洛程序，如FLUKA、MCNP、PHITS、GEANT4等，這些蒙特卡洛程序?qū)宋锢韺W(xué)的發(fā)展至關(guān)重要，同時(shí)它們的運(yùn)用過(guò)程也非常的有趣和巧妙。

具體以加速器輻射防護(hù)領(lǐng)域?yàn)槔?，粒子加速器產(chǎn)生的高速運(yùn)動(dòng)的微觀粒子（一般每秒可達(dá)上億個(gè)粒子）在與其他物質(zhì)碰撞時(shí)，會(huì)通過(guò)核反應(yīng)產(chǎn)生帶電粒子、中子、γ射線等次級(jí)輻射，這些次級(jí)輻射又會(huì)繼續(xù)與材料發(fā)生核反應(yīng)，產(chǎn)生更多的次級(jí)輻射。如此周而復(fù)始，其反應(yīng)過(guò)程非常復(fù)雜，無(wú)法通過(guò)人力模擬計(jì)算。而蒙特卡洛方法，正好能夠解決這個(gè)難題。

圖 加速器產(chǎn)生電離輻射示意圖

為了使加速器產(chǎn)生的強(qiáng)電離輻射降低到可接受的低水平，科研人員需要對(duì)次級(jí)粒子的輸運(yùn)過(guò)程進(jìn)行仿真模擬，并根據(jù)模擬結(jié)果進(jìn)行屏蔽阻擋設(shè)計(jì)。

接下來(lái)我們以一個(gè)簡(jiǎn)單模型為例，介紹蒙特卡洛方法應(yīng)用于粒子輸運(yùn)模擬的基本思路。如下圖所示，假定在一個(gè)二維的矩形屏蔽中，充滿了物質(zhì)B，左、上、下側(cè)設(shè)置黑色隔板，右側(cè)設(shè)置紅色隔板。

圖 模型示意圖

假設(shè)粒子A具有以下特質(zhì)：

粒子A在物質(zhì)B中直線運(yùn)動(dòng)。但每直線前進(jìn)1m，就會(huì)停止前進(jìn)并與B發(fā)生反應(yīng)，反應(yīng)后繼續(xù)直線前進(jìn)。反應(yīng)有三種可能：

①：A忽略與B的反應(yīng)，繼續(xù)按原方向前進(jìn)1m，發(fā)生概率為1/3；

②：A與B發(fā)生正碰撞，向左偏轉(zhuǎn)45°前進(jìn)1m，發(fā)生概率為1/3；

③：A與B發(fā)生反碰撞，向右偏轉(zhuǎn)45°前進(jìn)1m，發(fā)生概率為1/3。

2. 當(dāng)A碰到黑色隔板后，會(huì)立刻消失；

3. 當(dāng)A碰到紅色隔板后，會(huì)進(jìn)入環(huán)境。 如果我們想知道當(dāng)大量粒子A從上圖中P位置向右水平射出，會(huì)有多少個(gè)粒子通過(guò)紅色隔板進(jìn)入環(huán)境，就可以利用蒙特卡洛程序進(jìn)行大量模擬。 對(duì)單個(gè)事例，當(dāng)粒子A從P點(diǎn)出發(fā)前進(jìn)1米后，抽取偽隨機(jī)數(shù)N為[0,1]的均勻分布，根據(jù)抽取的數(shù)字決定其下一步的運(yùn)動(dòng)軌跡：當(dāng)0

下圖給出了粒子A可能的兩種運(yùn)行軌跡，利用計(jì)算機(jī)程序可以高效地完成大量的粒子模擬，得到計(jì)算結(jié)果。根據(jù)統(tǒng)計(jì)學(xué)特性，計(jì)算的事例越多，結(jié)果也就越接近期望值。

圖 粒子A可能的兩種運(yùn)行軌跡

在理解了上述示例的計(jì)算思路后，不難想象在輻射防護(hù)領(lǐng)域的研究中，只要我們知道了粒子在運(yùn)行過(guò)程中每一步可能發(fā)生的反應(yīng)類型及其概率等基本的核物理參數(shù)，就可以通過(guò)計(jì)算機(jī)程序?qū)崿F(xiàn)各種不同情形的模擬計(jì)算。

比如FLUKA程序可被用于計(jì)算國(guó)內(nèi)醫(yī)用重離子加速器HIMM治療室的輻射劑量率分布。如下圖所示，可以看出碳離子集中損失的位置輻射劑量率最高，約為106μSv/h；而經(jīng)過(guò)混凝土屏蔽后，屏蔽外的輻射劑量率衰減到了2.5μSv/h以下。

圖 HIMM治療室劑量率分布圖。HIMM裝置每秒鐘能夠產(chǎn)生1億個(gè)最高能量為400MeV/u的碳離子束用于治療，束流照射人體時(shí)幾乎全部損失，從而產(chǎn)生各種次級(jí)輻射。

利用FLUKA程序，還可以開(kāi)展更為復(fù)雜的輻射防護(hù)模擬研究，如更加多樣的束流損失模式、更多的粒子種類與能量、更加復(fù)雜的建筑結(jié)構(gòu)等。

圖 各類加速器輻射劑量分布圖示例

讀到這里，相信大家對(duì)蒙特卡洛方法已經(jīng)有了一個(gè)基本的了解。

審核編輯 ：李倩