動力學(xué)蒙特卡洛方法在計(jì)算物理教學(xué)中的引入

陳含爽

（安徽大學(xué)物理與材料科學(xué)學(xué)院，安徽合肥230601）

隨著計(jì)算機(jī)技術(shù)的快速發(fā)展，計(jì)算物理作為一門獨(dú)立的學(xué)科嶄露頭角，成為聯(lián)系理論物理和實(shí)驗(yàn)物理的橋梁。計(jì)算物理也是大多數(shù)本科高等院校物理以及相近專業(yè)開設(shè)的一門本科生專業(yè)課。大多數(shù)計(jì)算物理教材都會涉及蒙特卡洛方法的介紹，如圓周率和定積分的計(jì)算、任意分布隨機(jī)數(shù)的產(chǎn)生等[1]。蒙特卡洛方法是發(fā)展最為成熟的計(jì)算機(jī)模擬方法之一，最早是在1957年由Metropolis和Ulam等針對中子輸運(yùn)問題時提出的[2]。蒙特卡洛方法基本原理是基于多次隨機(jī)采樣來得到數(shù)值結(jié)果的計(jì)算機(jī)模型方法，已被廣泛地應(yīng)用于數(shù)值優(yōu)化和數(shù)值積分等問題。蒙特卡洛方法在許多領(lǐng)域均有廣泛的應(yīng)用，如統(tǒng)計(jì)物理、計(jì)算生物學(xué)、金融學(xué)以及人工智能等領(lǐng)域[3]。

動力學(xué)蒙特卡洛方法是模擬一些實(shí)際過程時間演化的有效的蒙特卡羅方法，其應(yīng)用范圍十分廣泛，包括晶體生長、表面反應(yīng)和擴(kuò)散等[4]。動力學(xué)蒙特卡洛方法的提出是為了解決常規(guī)蒙特卡洛方法效率不高的問題。例如，在模擬自旋模型的Metropolis算法中，自旋翻轉(zhuǎn)的嘗試可能不被執(zhí)行，這些嘗試被稱為空事件，尤其在低溫情況下，空事件出現(xiàn)的可能性很大，這時蒙特卡洛模擬的效率就不高，動力學(xué)蒙特卡洛方法可以避免空事件的發(fā)生[5]。動力學(xué)蒙特卡洛方法和Gillespie算法[6]本質(zhì)上是相同的，其算法實(shí)現(xiàn)是要回答兩個問題：一是下一個事件是什么時候發(fā)生，二是下一次發(fā)生的是哪一個事件。然而，大多數(shù)計(jì)算物理教科書都沒有涉及動力學(xué)蒙特卡洛方法的介紹。為此，本文嘗試通過簡潔的數(shù)學(xué)語言介紹動力學(xué)蒙特卡洛方法的基本原理以及實(shí)現(xiàn)方法。最后，通過一個簡單的例子——Susceptible-Infected-Susceptible：SIS傳播模型[7]，來闡述動力學(xué)蒙特卡洛方法實(shí)現(xiàn)的過程，以及如何進(jìn)行結(jié)果分析和討論。

1 動力學(xué)蒙特卡洛方法的基本原理

現(xiàn)有n個獨(dú)立無關(guān)的隨機(jī)過程，假設(shè)過程本身發(fā)生的時間忽略不計(jì)，相繼兩次過程之間的等待時間τ滿足指數(shù)分布（泊松過程），即

其中λi是第i個過程的發(fā)生速率。定義存活概率函數(shù)其含義是第i個過程的等待時間大于τ的概率，那么所有過程的存活概率為因此，整個過程的等待時間滿足概率密度函數(shù)是可以看出整個過程仍是泊松過程，任意一個過程發(fā)生的速率為所有過程發(fā)生速率之和，即因?yàn)榇婊罡怕适?到1之間的隨機(jī)數(shù)，所以等待時間可以寫做

其中r是0到1之間均勻分布的隨機(jī)數(shù)。第i個隨機(jī)過程發(fā)生的概率密度為

所以第i個隨機(jī)過程發(fā)生的概率為

基于上述理論，動力學(xué)蒙特卡洛方法可以總結(jié)為以下兩步：

（1）產(chǎn)生0到1之間的均勻分布的隨機(jī)數(shù)r，計(jì)算下一個過程發(fā)生的時間間隔更新時間

（2）產(chǎn)生另一個0到1之間的均勻分布的隨機(jī)數(shù)u，若滿足不等式則發(fā)生第j個隨機(jī)過程，更新粒子數(shù)和反應(yīng)速率。

2 動力學(xué)蒙特卡洛方法的應(yīng)用

作為一個具體例子，下面將動力學(xué)蒙特卡洛方法應(yīng)用到SIS模型。給定一個大小為N的網(wǎng)絡(luò)，網(wǎng)絡(luò)上每個節(jié)點(diǎn)放置一個粒子，粒子的狀態(tài)要么是易感態(tài)（susceptible），用S表示，要么是感染態(tài)（infected），用I表示。一方面，一個S粒子與一個I粒子接觸，即這兩個粒子是鄰居節(jié)點(diǎn)，那么S粒子被感染變成I粒子的速率為λ；另一方面，一個I粒子自發(fā)地恢復(fù)成S粒子的速率為μ。這兩個過程可以表示成下列兩個反應(yīng)：

SIS模型的動力學(xué)蒙特卡洛方法的實(shí)施如下：

（1）初始條件設(shè)置：隨機(jī)選擇一部分粒子作為感染種子節(jié)點(diǎn)（例如10%的節(jié)點(diǎn)），即這些種子節(jié)點(diǎn)的初始狀態(tài)設(shè)為I，剩下節(jié)點(diǎn)的初始狀態(tài)設(shè)為S。設(shè)定初始時間t=0，產(chǎn)生初始I粒子和SI邊的列表，記NI和NSI為I粒子和SI邊的數(shù)目。

（2）產(chǎn)生均勻分布的隨機(jī)數(shù)r∈[0,1]，計(jì)算下一個過程發(fā)生的時間間隔更新時間

（4）更新I粒子和SI邊的列表?；氐剑?）。

3 結(jié)果分析

定義有效感染系數(shù)β=λ μ。圖1給出了ER隨機(jī)網(wǎng)絡(luò)[8]上SIS模型的感染范圍ρ(t)=NI(t)N的3條時間序列，對應(yīng)于3個不同的β值。當(dāng)β較小時，最終的感染范圍變?yōu)榱?，稱之為健康態(tài)，該狀態(tài)是吸收態(tài)。吸收態(tài)的意思是一旦體系達(dá)到這個狀態(tài)，會永遠(yuǎn)停留在此態(tài)[9]。對SIS模型來說，一旦所有節(jié)點(diǎn)的狀態(tài)都是S態(tài)，(5)式兩個反應(yīng)都無法進(jìn)行。當(dāng)β較大時，最終的感染范圍會維持在某個值附近。由于有限大小網(wǎng)絡(luò)原因，感染范圍會在穩(wěn)態(tài)值附近漲落，網(wǎng)絡(luò)規(guī)模越小，漲落越大。

圖1 3個不同感染速率下的ER網(wǎng)絡(luò)上感染范圍ρ(t)的時間序列。網(wǎng)絡(luò)大小N=1 000，平均度 k=10

圖2 給出了β=0.15時兩次相繼反應(yīng)之間的等待時間Δt的分布P( )Δt，如果將縱坐標(biāo)設(shè)置成對數(shù)值，如圖2中的插圖所示，可以看出滿足指數(shù)分布，即這與之前泊松過程假設(shè)是一致的。

圖2 β=0.15時兩次相繼反應(yīng)之間的等待時間Δt的分布P( )Δt

圖3 給出了感染范圍的穩(wěn)態(tài)平均值 ρ隨有效感染系數(shù)β的變化圖，其中是采樣時間長度。當(dāng)β＜βc時，ρ=0；當(dāng)β＞βc時，ρ＞0。體系在傳播速率閾值β=βc時發(fā)生了相變。然而，當(dāng)β稍稍大于βc時，感染范圍并不大，由于有限尺度漲落原因，體系很容易達(dá)到吸收態(tài)。因此，很難通過圖3來準(zhǔn)確地定出傳播速率閾值，為此，采取準(zhǔn)靜態(tài)方法來確定βc的值。準(zhǔn)靜態(tài)方法的基本思想是：一旦體系達(dá)到吸收態(tài)，隨機(jī)地賦予體系的一個采用過的非吸收態(tài)構(gòu)型。通過計(jì)算ρ的漲落有關(guān)的量，定義為極大值的位置給出傳播速率閾值βc的值，如圖4所示，βc=0.106。

圖3 感染范圍的穩(wěn)態(tài)平均值 ρ隨有效感染系數(shù)β的變化圖

圖4 χ隨有效感染系數(shù)β的變化圖

4 結(jié)論

本文介紹了動力學(xué)蒙特卡洛方法的基本原理以及程序?qū)崿F(xiàn)的基本步驟，并以一個簡單的例子——網(wǎng)絡(luò)上SIS疾病傳播過程來說明動力學(xué)蒙特卡洛方法的實(shí)際應(yīng)用。通過計(jì)算傳播范圍與傳播速率之間的關(guān)系來闡述該模型相變的性質(zhì)，并通過計(jì)算漲落相關(guān)的物理量來確定相變的位置，即傳播速率的閾值。本文對動力學(xué)蒙特卡洛方法的基本原理的闡述是考慮了本科生對數(shù)學(xué)和物理知識掌握的實(shí)際程度，相信只要本科生掌握了一些微積分和概率論的知識，本文的推導(dǎo)是完全可以看懂的。此外，本文討論的是泊松過程，若隨機(jī)過程的等待時間不滿足指數(shù)分布，即非泊松過程，也可以使用類似的方法推導(dǎo)[11]。