王玉輝

（湖南第一師范學(xué)院 信息科學(xué)與工程學(xué)院，湖南 長(zhǎng)沙 410205）

漂移—擴(kuò)散模型的蒙特卡洛隨機(jī)仿真

王玉輝

（湖南第一師范學(xué)院 信息科學(xué)與工程學(xué)院，湖南 長(zhǎng)沙 410205）

為比較直接蒙特卡洛仿真與基于隨機(jī)游走的蒙特卡洛仿真的精度，使用數(shù)值仿真對(duì)粒子漂移和擴(kuò)散過(guò)程中的一階、二階和二階中心矩的相對(duì)誤差進(jìn)行了測(cè)試。實(shí)驗(yàn)結(jié)果表明，兩種方法的精度和性能大致相當(dāng)。因此，在漂移—擴(kuò)散模型的蒙特卡洛隨機(jī)仿真中，固定或變化的時(shí)間步長(zhǎng)可以根據(jù)需要和連續(xù)空間或離散網(wǎng)格靈活組合使用，并不會(huì)嚴(yán)重影響仿真的精度。

蒙特卡洛隨機(jī)仿真；漂移—擴(kuò)散；隨機(jī)游走；仿真精度

漂移—擴(kuò)散是各種物理、化學(xué)、生物系統(tǒng)中一種常見(jiàn)的現(xiàn)象。對(duì)于這類(lèi)漂移—擴(kuò)散或更復(fù)雜的反應(yīng)—漂移—擴(kuò)散模型，當(dāng)問(wèn)題的解析解難以求得或明確需要采用單個(gè)粒子追蹤等方法模擬粒子的微觀(guān)、隨機(jī)運(yùn)動(dòng)過(guò)程時(shí)，蒙特卡洛隨機(jī)仿真就成了求解這類(lèi)模型的一種重要方法［1］。

相對(duì)而言，直接蒙特卡洛仿真和基于隨機(jī)游走的蒙特卡洛仿真是兩種較為常見(jiàn)和經(jīng)典的隨機(jī)仿真方法。前者根據(jù)漂移-擴(kuò)散模型在自然邊界條件下解的形式，使用隨機(jī)的、服從正態(tài)分布的跳躍步長(zhǎng)和固定或變化的時(shí)間步長(zhǎng)，在連續(xù)空間中模擬粒子的運(yùn)動(dòng)。后者則使用固定網(wǎng)格步長(zhǎng)以及固定或服從一定概率分布的隨機(jī)時(shí)間步長(zhǎng)，根據(jù)一定的停留概率和到鄰近網(wǎng)格的轉(zhuǎn)移概率，在離散網(wǎng)格格點(diǎn)上模擬粒子運(yùn)動(dòng)［2］。

一般而言，蒙特卡洛隨機(jī)仿真的精度主要取決于算法自身帶來(lái)的誤差和偽隨機(jī)數(shù)產(chǎn)生器等因素。但是，目前有關(guān)蒙特卡洛隨機(jī)仿真的算法精度的研究相對(duì)比較少，更難以從理論上進(jìn)行分析。對(duì)于基于隨機(jī)游走的蒙特卡洛仿真，若網(wǎng)格步長(zhǎng)足夠精細(xì)，其精度原則上可以無(wú)限逼近連續(xù)系統(tǒng)模型對(duì)應(yīng)的解析解［3］。為了定量分析和比較以上兩種蒙特卡洛仿真的精度，為具體應(yīng)用實(shí)踐中的算法選擇和時(shí)間步長(zhǎng)、網(wǎng)格步長(zhǎng)等仿真參數(shù)選擇提供參考，可以使用粒子漂移-擴(kuò)散過(guò)程中的一階、二階和二階中心矩的相對(duì)誤差度量和評(píng)估仿真精度，通過(guò)數(shù)值仿真方法進(jìn)行實(shí)驗(yàn)測(cè)試和驗(yàn)證。

1蒙特卡洛隨機(jī)仿真算法

不失一般性，考慮一維空間具有恒定擴(kuò)散系數(shù)D、漂移速度q的漂移-擴(kuò)散模型［4］：

在自然邊界和初始條件下，其解析解為：

任意時(shí)刻，相應(yīng)的均值和方差為：

原則上，蒙特卡洛隨機(jī)仿真的目標(biāo)是精確模擬和復(fù)現(xiàn)上述公式（3）和（2）描述的粒子宏觀(guān)運(yùn)動(dòng)特征和微觀(guān)分布狀態(tài)。

2.1直接蒙特卡洛仿真算法

直接蒙特卡洛仿真的思想比較簡(jiǎn)單：根據(jù)公式（2）描述的粒子概率密度分布函數(shù)及給定的時(shí)間步長(zhǎng)，可推斷下一時(shí)刻粒子的位置分布服從正態(tài)分布。因此，可以直接使用正態(tài)分布進(jìn)行隨機(jī)采樣，得到一個(gè)隨機(jī)的跳躍步長(zhǎng)，再使用該步長(zhǎng)更新粒子的位置。對(duì)所有粒子重復(fù)上述過(guò)程，直至仿真結(jié)束，即可實(shí)現(xiàn)蒙特卡洛隨機(jī)仿真。

顯然，上述隨機(jī)跳躍步長(zhǎng)采樣將導(dǎo)致粒子在連續(xù)空間的位置更新。若設(shè)定一個(gè)網(wǎng)格尺寸，可以將粒子在連續(xù)空間的位置映射到網(wǎng)格格點(diǎn)中，但同時(shí)仍需保留在連續(xù)空間中的位置狀態(tài)變量。否則，可能影響和降低仿真的精度。

原則上，直接蒙特卡洛仿真可以使用固定或者變化的時(shí)間步長(zhǎng)。步長(zhǎng)的大小只影響仿真的時(shí)間分辨率，而不會(huì)對(duì)仿真的精度產(chǎn)生太大的影響。因而，該方法具有簡(jiǎn)單、直觀(guān)、靈活的優(yōu)點(diǎn)。

2.2基于隨機(jī)游走的網(wǎng)格蒙特卡洛仿真算法

采用網(wǎng)格蒙特卡洛仿真算法對(duì)漂移-擴(kuò)散模型進(jìn)行求解首先需要將問(wèn)題物理空間劃分為網(wǎng)格（通常為均勻網(wǎng)格）。然后，計(jì)算網(wǎng)格中的粒子在離散的時(shí)間步長(zhǎng)上在當(dāng)前網(wǎng)格格點(diǎn)中的停留概率及向鄰近網(wǎng)格格點(diǎn)的轉(zhuǎn)移概率，同時(shí)按照一定概率分布對(duì)時(shí)間步長(zhǎng)進(jìn)行采樣，由此追蹤粒子在不同網(wǎng)格格點(diǎn)間的漂移-擴(kuò)散過(guò)程。顯然，以上轉(zhuǎn)移概率、停留概率和時(shí)間步長(zhǎng)均與網(wǎng)格步長(zhǎng)有關(guān)。

對(duì)于公式（1）中的漂移—擴(kuò)散模型，假設(shè)網(wǎng)格步長(zhǎng)為λ，則一個(gè)連續(xù)時(shí)間隨機(jī)游走粒子的首次跳出當(dāng)前網(wǎng)格的時(shí)間概率密度分布可以用以下解析形式給出：

向左右兩側(cè)轉(zhuǎn)移的時(shí)間積分概率為：

原則上，根據(jù)公式（4）和（5）即可設(shè)計(jì)出一個(gè)基于連續(xù)時(shí)間隨機(jī)游走的蒙特卡洛仿真算法［5］：首先根據(jù)公式（4）進(jìn)行采樣，得到一個(gè)隨機(jī)的時(shí)間步長(zhǎng)，再隨機(jī)產(chǎn)生一個(gè)在［0，1］區(qū)間均勻分布的隨機(jī)數(shù)，根據(jù)公式（5）計(jì)算和確定該粒子跳躍和轉(zhuǎn)移到的網(wǎng)格格點(diǎn)。事實(shí)上，由于公式（4）涉及無(wú)窮項(xiàng)級(jí)數(shù)的求和運(yùn)算，而很多數(shù)值計(jì)算應(yīng)用需要固定時(shí)間步長(zhǎng)算法，因而需要對(duì)上述算法進(jìn)行改進(jìn)。其主要思路是通過(guò)一個(gè)等價(jià)的離散時(shí)間隨機(jī)游走模型代替連續(xù)時(shí)間隨機(jī)游走模型，同時(shí)保證正確地復(fù)現(xiàn)粒子分布的均值和均方位移等特性［2］。

為方便比較分析，公式（4）的Laplace變換為：

使用Taylor展開(kāi)公式，可以得到首次跳出時(shí)間的一階和二階矩：

假設(shè)使用一個(gè)固定時(shí)間步長(zhǎng)隨機(jī)游走模型，其首次跳出時(shí)間的概率分布密度函數(shù)為：

其中，τ為時(shí)間步長(zhǎng)，p0為停留概率，δ為Dirac函數(shù)。易知，公式（10）對(duì)應(yīng)的Laplace變換為：

同樣使用公式Taylor展開(kāi)，可得：

比較公式（12）和公式（7）中的系數(shù)，可得以下關(guān)系：

根據(jù)公式（13）和（14）即可構(gòu)造固定步長(zhǎng)的網(wǎng)格蒙特卡洛仿真算法：首先根據(jù)給定的網(wǎng)格步長(zhǎng)λ、飄移速度q和擴(kuò)散系數(shù)D，計(jì)算停留概率p0和時(shí)間步長(zhǎng)τ，在每個(gè)時(shí)間步長(zhǎng)上，對(duì)每個(gè)粒子隨機(jī)產(chǎn)生一個(gè)在［0，1］區(qū)間均勻分布的隨機(jī)數(shù)，根據(jù)公式（13）計(jì)算和確定該粒子是否跳躍；若是，再產(chǎn)生一個(gè)在［0，1］區(qū)間均勻分布的隨機(jī)數(shù)，根據(jù)公式（5）計(jì)算和確定轉(zhuǎn)移到的網(wǎng)格格點(diǎn)。重復(fù)上述過(guò)程，直至仿真結(jié)束。

實(shí)際上，為了提高效率，可以將公式（13）和公式（5）結(jié)合，每次只采樣產(chǎn)生一個(gè)隨機(jī)數(shù)，即可確定其下一時(shí)刻的位置。

對(duì)于固定步長(zhǎng)的網(wǎng)格蒙特卡洛仿真算法，當(dāng)網(wǎng)格步長(zhǎng)和飄移—擴(kuò)散參數(shù)給定時(shí)，即可唯一地確定時(shí)間步長(zhǎng)、停留概率和轉(zhuǎn)移概率。使用其他時(shí)間步長(zhǎng)或者變化時(shí)間步長(zhǎng)都將損害仿真的精度，甚至導(dǎo)致不正確的結(jié)果。

2.數(shù)值結(jié)果

為了分析比較以上兩種仿真算法的精度，根據(jù)公式（3）給出的理論值，可以使用粒子在漂移-擴(kuò)散過(guò)程中的一階、二階和二階中心矩的相對(duì)誤差對(duì)仿真精度進(jìn)行度量和評(píng)估：

根據(jù)第2節(jié)中的蒙特卡洛仿真算法，使用標(biāo)準(zhǔn)C++語(yǔ)言編程，利用boost軟件庫(kù)提供的隨機(jī)數(shù)生成器，不難實(shí)現(xiàn)上述的直接蒙特卡洛仿真算法和固定時(shí)間步長(zhǎng)的網(wǎng)格蒙特卡洛仿真算法。

為了進(jìn)行數(shù)值仿真和便于比較分析，實(shí)驗(yàn)中選擇使用一個(gè)固定的擴(kuò)散系數(shù)和相對(duì)變化的其他參數(shù)，具體參數(shù)設(shè)置如表1中所示。

表1　飄移—擴(kuò)散模型仿真參數(shù)設(shè)置

，其取值在0.025到10之間，從而使仿真實(shí)驗(yàn)具有較好的代表性。數(shù)值仿真使用Windows8 64位操作系統(tǒng)，直接蒙特卡洛仿真算法采用和網(wǎng)格蒙特卡洛仿真算法相同的固定時(shí)間步長(zhǎng)。服從均勻分布和正態(tài)分布的隨機(jī)數(shù)生成均使用Mersenne Twister算法，以盡可能降低偽隨機(jī)數(shù)產(chǎn)生器帶來(lái)的誤差。仿真過(guò)程中，在每個(gè)時(shí)間步長(zhǎng)上，計(jì)算粒子漂移-擴(kuò)散的一階、二階和二階中心矩的相對(duì)誤差，取平均記錄分析。其中，直接蒙特卡洛仿真算法得到的粒子位置需要將其映射到網(wǎng)格格點(diǎn)上之后，再計(jì)算相對(duì)誤差。

圖1　平均相對(duì)誤差（粒子數(shù)目N=104）

對(duì)于不同的仿真配置和Péclet數(shù)，得到的一階、二階和二階中心矩的平均相對(duì)誤差分別在圖1（粒子數(shù)目N=104）和圖2（粒子數(shù)目N=105）中繪出。其中，實(shí)線(xiàn)代表直接蒙特卡洛仿真結(jié)果，虛線(xiàn)代表基于隨機(jī)游走的固定步長(zhǎng)網(wǎng)格蒙特卡洛仿真結(jié)果。為了清楚顯示，圖中的坐標(biāo)軸采用雙對(duì)數(shù)刻度方式。注意圖1和圖2中縱坐標(biāo)刻度的不同。

圖2　平均相對(duì)誤差（粒子數(shù)目N=10^5）

根據(jù)以上結(jié)果，易見(jiàn)：（1）不管粒子數(shù)目的多少，就一階、二階和二階中心矩的平均相對(duì)誤差而言，兩種仿真方法的數(shù)值精度基本相當(dāng)；（2）對(duì)于不同配置參數(shù)和Péclet數(shù)，相對(duì)于一階和二階矩而言，二階中心矩的平均相對(duì)誤差的變化較?。唬?）增加粒子數(shù)目或仿真次數(shù)，可以顯著降低兩種仿真方法的平均誤差，仿真精度提高的程度符合蒙特卡洛仿真的規(guī)律；（4）對(duì)于相同的漂移—擴(kuò)散模型參數(shù)，網(wǎng)格步長(zhǎng)的大小對(duì)仿真精度的影響相對(duì)較小，這表明兩種算法均比較穩(wěn)定；（5）一般而言，漂移速度越大，平均相對(duì)誤差越小，這可能與一階、二階矩的絕對(duì)取值有關(guān)，而二階中心矩的理論取值僅與擴(kuò)散系數(shù)有關(guān)，因而其平均相對(duì)誤差變化較小。

為了進(jìn)一步得到兩種仿真方法的直觀(guān)結(jié)果，圖3中給出了在一種典型配置條件（擴(kuò)散系數(shù)、漂移速度、網(wǎng)格步長(zhǎng)均為1，粒子數(shù)目為10^5，仿真次數(shù)為10）和不同時(shí)刻的粒子分布的概率密度曲線(xiàn)。其中，實(shí)線(xiàn)和虛線(xiàn)代表通過(guò)公式（2）計(jì)算的不同時(shí)刻的理論值，實(shí)心和空心的圓圈或方塊分別代表直接蒙特卡洛仿真和基于隨機(jī)游走的固定時(shí)間步長(zhǎng)網(wǎng)格蒙特卡洛仿真結(jié)果。

圖3　不同時(shí)刻的粒子分布概率密度

由圖3可以看出，兩種不同仿真方法得到的不同時(shí)刻粒子分布概率密度與理論值的誤差都很小，仿真精度也基本相當(dāng)。

就仿真算法的時(shí)間性能而言，兩者也相差無(wú)幾。其原因是，除了隨機(jī)數(shù)采樣使用不同的分布函數(shù)之外，兩種仿真方法使用基本相同的算法流程。而兩種仿真方法均使用相同的Mersenne Twister核心算法，經(jīng)變換后得到服從正態(tài)或均勻分布的隨機(jī)數(shù)，因而，所需的計(jì)算時(shí)間也相差不大。因此，兩者的運(yùn)行時(shí)間性能基本相當(dāng)，這與實(shí)際測(cè)試得到的結(jié)果完全一致，因此不再給出具體數(shù)值結(jié)果。

值得注意的是，對(duì)于蒙特卡洛仿真的精度評(píng)估還有其他度量指標(biāo)和方法。但是對(duì)于分析比較兩種蒙特卡洛仿真的精度而言，上述基于一階、二階和二階中心矩的平均相對(duì)誤差以及不同時(shí)刻的粒子分布概率密度的仿真結(jié)果已經(jīng)能夠基本上證明兩者具有大致相當(dāng)?shù)木取?/p>

此外，在仿真實(shí)驗(yàn)中測(cè)試得到的誤差不僅包含算法自身帶來(lái)的誤差，還包括偽隨機(jī)數(shù)產(chǎn)生器的誤差。進(jìn)一步的研究還應(yīng)考慮不同的偽隨機(jī)數(shù)產(chǎn)生器對(duì)仿真結(jié)果和精度的影響，從而更準(zhǔn)確地評(píng)估不同仿真算法的真實(shí)性能。

結(jié)論

蒙特卡洛隨機(jī)仿真是求解漂移—擴(kuò)散模型的一類(lèi)重要方法，其精度主要取決于算法自身帶來(lái)的誤差和偽隨機(jī)數(shù)產(chǎn)生器等因素，而算法的精度直接決定著仿真結(jié)果的可信度。

針對(duì)用于漂移—擴(kuò)散模型的直接蒙特卡洛仿真方法和基于隨機(jī)游走的蒙特卡洛仿真方法，使用數(shù)值仿真對(duì)粒子漂移和擴(kuò)散過(guò)程中的一階、二階和二階中心矩的平均相對(duì)誤差進(jìn)行了測(cè)試，分析和比較了兩種不同方法的精度。實(shí)驗(yàn)結(jié)果證實(shí)兩種方法具有大致相當(dāng)?shù)姆抡婢群蜁r(shí)間性能。因此，在漂移-擴(kuò)散模型的蒙特卡洛隨機(jī)仿真中，固定或變化的時(shí)間步長(zhǎng)可以根據(jù)實(shí)際需要和連續(xù)空間或離散網(wǎng)格靈活組合使用，而不會(huì)嚴(yán)重影響仿真的精度。

進(jìn)一步的研究應(yīng)充分考慮偽隨機(jī)數(shù)產(chǎn)生器在蒙特卡洛仿真中帶來(lái)的誤差以及影響，因?yàn)橐坏┱_設(shè)計(jì)和實(shí)現(xiàn)了蒙特卡洛仿真算法，其精度將主要取決于偽隨機(jī)數(shù)產(chǎn)生器的好壞。而在蒙特卡洛仿真中，通過(guò)嚴(yán)格隨機(jī)性指標(biāo)測(cè)試的偽隨機(jī)數(shù)產(chǎn)生器并不一定意味著帶來(lái)較高的仿真精度。除了采用粒子漂移和擴(kuò)散過(guò)程中的一階、二階和二階中心矩的平均相對(duì)誤差進(jìn)行仿真精度分析評(píng)估外，還可以考慮其他可能的度量指標(biāo)和評(píng)估方法，對(duì)兩種不同的仿真算法精度進(jìn)行更全面評(píng)估。此外，將上述蒙特卡洛仿真算法擴(kuò)展到多維、非均勻的漂移-擴(kuò)散模型或更復(fù)雜的反應(yīng)-漂移-擴(kuò)散模型及利用并行仿真方法［9］和計(jì)算設(shè)備如圖形處理單元以提高仿真的性能，都是值得深入研究的方向。

［1］M.E.J.紐曼,G.T.巴克馬.統(tǒng)計(jì)物理學(xué)中的蒙特卡羅方法（第1版）［M］.北京：世界圖書(shū)出版公司,2013.

［2］M.G.Gauthier,G.W.Slater.Building reliable lattice Monte Carlo models for real drift and diffusion problems ［J］.Physical Review E,2004,70（015763）.

［3］M.V.Chubynsky,G.W.Slater.Optimizing the accuracy of lattice Monte Carlo algorithms for simulating diffusion［J］.Physical Review E,2012,85（016709）.

［4］謝定裕.流體力學(xué)［M］.天津：南開(kāi)大學(xué)出版社,1987.

［5］E.W.Montroll,H.Scher.Random walks on lattices：IV. Continuous-time walks and influence of absorbing boundaries［J］.Journal of Statistical Physics,1973,9（2）, 101-135.

［6］花嶸,傅游,康繼昌.直接蒙特卡洛問(wèn)題的并行化方法［J］.計(jì)算機(jī)工程,2004,30（5）,40-42.

［責(zé)任編輯：胡偉］

Monte Carlo Stochastic Simulation of Drift-Diffusion Models

WANG Yu-hui

（Department of Information Science and Engineering,Hunan First Normal University,Changsha,Hunan 410205）

To investigate the accuracy of Direct Monte Carlo simulation and Random-Walk-based Monte Carlo simulation,numerical simulations are carried out to measure the relative errors introduced in the first two moments and the second central moment of the displacement of particles during the drift-diffusion process.It turns out that both the accuracy and the performance of the two approaches are almost identical.Therefore,it is concluded that for such Monte Carlo stochastic simulation of drift-diffusion models,fixed or variable time steps can be used in a flexible manner with either discrete lattice or continuous space without significant degrade of accuracy.

Monte Carlo stochastic simulation；drift-diffusion；random walk；simulation accuracy

TP391.9

A

1674-831X（2016）02-0105-04

2014-11-12

王玉輝（1979-），女，湖南長(zhǎng)沙人，碩士，湖南第一師范學(xué)校信息科學(xué)與工程學(xué)院講師，主要從事計(jì)算機(jī)仿真研究。