漂移—擴散模型的蒙特卡洛隨機仿真

王玉輝

（湖南第一師范學院 信息科學與工程學院，湖南 長沙 410205）

漂移—擴散模型的蒙特卡洛隨機仿真

王玉輝

（湖南第一師范學院 信息科學與工程學院，湖南 長沙 410205）

為比較直接蒙特卡洛仿真與基于隨機游走的蒙特卡洛仿真的精度，使用數值仿真對粒子漂移和擴散過程中的一階、二階和二階中心矩的相對誤差進行了測試。實驗結果表明，兩種方法的精度和性能大致相當。因此，在漂移—擴散模型的蒙特卡洛隨機仿真中，固定或變化的時間步長可以根據需要和連續空間或離散網格靈活組合使用，并不會嚴重影響仿真的精度。

蒙特卡洛隨機仿真；漂移—擴散；隨機游走；仿真精度

漂移—擴散是各種物理、化學、生物系統中一種常見的現象。對于這類漂移—擴散或更復雜的反應—漂移—擴散模型，當問題的解析解難以求得或明確需要采用單個粒子追蹤等方法模擬粒子的微觀、隨機運動過程時，蒙特卡洛隨機仿真就成了求解這類模型的一種重要方法［1］。

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

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

1蒙特卡洛隨機仿真算法

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

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

任意時刻，相應的均值和方差為：

原則上，蒙特卡洛隨機仿真的目標是精確模擬和復現上述公式（3）和（2）描述的粒子宏觀運動特征和微觀分布狀態。

2.1直接蒙特卡洛仿真算法

直接蒙特卡洛仿真的思想比較簡單：根據公式（2）描述的粒子概率密度分布函數及給定的時間步長，可推斷下一時刻粒子的位置分布服從正態分布。因此，可以直接使用正態分布進行隨機采樣，得到一個隨機的跳躍步長，再使用該步長更新粒子的位置。對所有粒子重復上述過程，直至仿真結束，即可實現蒙特卡洛隨機仿真。

顯然，上述隨機跳躍步長采樣將導致粒子在連續空間的位置更新。若設定一個網格尺寸，可以將粒子在連續空間的位置映射到網格格點中，但同時仍需保留在連續空間中的位置狀態變量。否則，可能影響和降低仿真的精度。

原則上，直接蒙特卡洛仿真可以使用固定或者變化的時間步長。步長的大小只影響仿真的時間分辨率，而不會對仿真的精度產生太大的影響。因而，該方法具有簡單、直觀、靈活的優點。

2.2基于隨機游走的網格蒙特卡洛仿真算法

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

對于公式（1）中的漂移—擴散模型，假設網格步長為λ，則一個連續時間隨機游走粒子的首次跳出當前網格的時間概率密度分布可以用以下解析形式給出：

向左右兩側轉移的時間積分概率為：

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

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

使用Taylor展開公式，可以得到首次跳出時間的一階和二階矩：

假設使用一個固定時間步長隨機游走模型，其首次跳出時間的概率分布密度函數為：

其中，τ為時間步長，p0為停留概率，δ為Dirac函數。易知，公式（10）對應的Laplace變換為：

同樣使用公式Taylor展開，可得：

比較公式（12）和公式（7）中的系數，可得以下關系：

根據公式（13）和（14）即可構造固定步長的網格蒙特卡洛仿真算法：首先根據給定的網格步長λ、飄移速度q和擴散系數D，計算停留概率p0和時間步長τ，在每個時間步長上，對每個粒子隨機產生一個在［0，1］區間均勻分布的隨機數，根據公式（13）計算和確定該粒子是否跳躍；若是，再產生一個在［0，1］區間均勻分布的隨機數，根據公式（5）計算和確定轉移到的網格格點。重復上述過程，直至仿真結束。

實際上，為了提高效率，可以將公式（13）和公式（5）結合，每次只采樣產生一個隨機數，即可確定其下一時刻的位置。

對于固定步長的網格蒙特卡洛仿真算法，當網格步長和飄移—擴散參數給定時，即可唯一地確定時間步長、停留概率和轉移概率。使用其他時間步長或者變化時間步長都將損害仿真的精度，甚至導致不正確的結果。

2.數值結果

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

根據第2節中的蒙特卡洛仿真算法，使用標準C++語言編程，利用boost軟件庫提供的隨機數生成器，不難實現上述的直接蒙特卡洛仿真算法和固定時間步長的網格蒙特卡洛仿真算法。

為了進行數值仿真和便于比較分析，實驗中選擇使用一個固定的擴散系數和相對變化的其他參數，具體參數設置如表1中所示。

表1　飄移—擴散模型仿真參數設置

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

圖1　平均相對誤差（粒子數目N=104）

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

圖2　平均相對誤差（粒子數目N=10^5）

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

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

圖3　不同時刻的粒子分布概率密度

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

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

值得注意的是，對于蒙特卡洛仿真的精度評估還有其他度量指標和方法。但是對于分析比較兩種蒙特卡洛仿真的精度而言，上述基于一階、二階和二階中心矩的平均相對誤差以及不同時刻的粒子分布概率密度的仿真結果已經能夠基本上證明兩者具有大致相當的精度。

此外，在仿真實驗中測試得到的誤差不僅包含算法自身帶來的誤差，還包括偽隨機數產生器的誤差。進一步的研究還應考慮不同的偽隨機數產生器對仿真結果和精度的影響，從而更準確地評估不同仿真算法的真實性能。

結論

蒙特卡洛隨機仿真是求解漂移—擴散模型的一類重要方法，其精度主要取決于算法自身帶來的誤差和偽隨機數產生器等因素，而算法的精度直接決定著仿真結果的可信度。

針對用于漂移—擴散模型的直接蒙特卡洛仿真方法和基于隨機游走的蒙特卡洛仿真方法，使用數值仿真對粒子漂移和擴散過程中的一階、二階和二階中心矩的平均相對誤差進行了測試，分析和比較了兩種不同方法的精度。實驗結果證實兩種方法具有大致相當的仿真精度和時間性能。因此，在漂移-擴散模型的蒙特卡洛隨機仿真中，固定或變化的時間步長可以根據實際需要和連續空間或離散網格靈活組合使用，而不會嚴重影響仿真的精度。

進一步的研究應充分考慮偽隨機數產生器在蒙特卡洛仿真中帶來的誤差以及影響，因為一旦正確設計和實現了蒙特卡洛仿真算法，其精度將主要取決于偽隨機數產生器的好壞。而在蒙特卡洛仿真中，通過嚴格隨機性指標測試的偽隨機數產生器并不一定意味著帶來較高的仿真精度。除了采用粒子漂移和擴散過程中的一階、二階和二階中心矩的平均相對誤差進行仿真精度分析評估外，還可以考慮其他可能的度量指標和評估方法，對兩種不同的仿真算法精度進行更全面評估。此外，將上述蒙特卡洛仿真算法擴展到多維、非均勻的漂移-擴散模型或更復雜的反應-漂移-擴散模型及利用并行仿真方法［9］和計算設備如圖形處理單元以提高仿真的性能，都是值得深入研究的方向。

［1］M.E.J.紐曼,G.T.巴克馬.統計物理學中的蒙特卡羅方法（第1版）［M］.北京：世界圖書出版公司,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］謝定裕.流體力學［M］.天津：南開大學出版社,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］花嶸,傅游,康繼昌.直接蒙特卡洛問題的并行化方法［J］.計算機工程,2004,30（5）,40-42.

［責任編輯：胡偉］

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-），女，湖南長沙人，碩士，湖南第一師范學校信息科學與工程學院講師，主要從事計算機仿真研究。