动力学蒙特卡洛方法在计算物理教学中的引入

陈含爽

（安徽大学物理与材料科学学院，安徽合肥230601）

随着计算机技术的快速发展，计算物理作为一门独立的学科崭露头角，成为联系理论物理和实验物理的桥梁。计算物理也是大多数本科高等院校物理以及相近专业开设的一门本科生专业课。大多数计算物理教材都会涉及蒙特卡洛方法的介绍，如圆周率和定积分的计算、任意分布随机数的产生等[1]。蒙特卡洛方法是发展最为成熟的计算机模拟方法之一，最早是在1957年由Metropolis和Ulam等针对中子输运问题时提出的[2]。蒙特卡洛方法基本原理是基于多次随机采样来得到数值结果的计算机模型方法，已被广泛地应用于数值优化和数值积分等问题。蒙特卡洛方法在许多领域均有广泛的应用，如统计物理、计算生物学、金融学以及人工智能等领域[3]。

动力学蒙特卡洛方法是模拟一些实际过程时间演化的有效的蒙特卡罗方法，其应用范围十分广泛，包括晶体生长、表面反应和扩散等[4]。动力学蒙特卡洛方法的提出是为了解决常规蒙特卡洛方法效率不高的问题。例如，在模拟自旋模型的Metropolis算法中，自旋翻转的尝试可能不被执行，这些尝试被称为空事件，尤其在低温情况下，空事件出现的可能性很大，这时蒙特卡洛模拟的效率就不高，动力学蒙特卡洛方法可以避免空事件的发生[5]。动力学蒙特卡洛方法和Gillespie算法[6]本质上是相同的，其算法实现是要回答两个问题：一是下一个事件是什么时候发生，二是下一次发生的是哪一个事件。然而，大多数计算物理教科书都没有涉及动力学蒙特卡洛方法的介绍。为此，本文尝试通过简洁的数学语言介绍动力学蒙特卡洛方法的基本原理以及实现方法。最后，通过一个简单的例子——Susceptible-Infected-Susceptible：SIS传播模型[7]，来阐述动力学蒙特卡洛方法实现的过程，以及如何进行结果分析和讨论。

1 动力学蒙特卡洛方法的基本原理

现有n个独立无关的随机过程，假设过程本身发生的时间忽略不计，相继两次过程之间的等待时间τ满足指数分布（泊松过程），即

其中λi是第i个过程的发生速率。定义存活概率函数其含义是第i个过程的等待时间大于τ的概率，那么所有过程的存活概率为因此，整个过程的等待时间满足概率密度函数是可以看出整个过程仍是泊松过程，任意一个过程发生的速率为所有过程发生速率之和，即因为存活概率是0到1之间的随机数，所以等待时间可以写做

其中r是0到1之间均匀分布的随机数。第i个随机过程发生的概率密度为

所以第i个随机过程发生的概率为

基于上述理论，动力学蒙特卡洛方法可以总结为以下两步：

（1）产生0到1之间的均匀分布的随机数r，计算下一个过程发生的时间间隔更新时间

（2）产生另一个0到1之间的均匀分布的随机数u，若满足不等式则发生第j个随机过程，更新粒子数和反应速率。

2 动力学蒙特卡洛方法的应用

作为一个具体例子，下面将动力学蒙特卡洛方法应用到SIS模型。给定一个大小为N的网络，网络上每个节点放置一个粒子，粒子的状态要么是易感态（susceptible），用S表示，要么是感染态（infected），用I表示。一方面，一个S粒子与一个I粒子接触，即这两个粒子是邻居节点，那么S粒子被感染变成I粒子的速率为λ；另一方面，一个I粒子自发地恢复成S粒子的速率为μ。这两个过程可以表示成下列两个反应：

SIS模型的动力学蒙特卡洛方法的实施如下：

（1）初始条件设置：随机选择一部分粒子作为感染种子节点（例如10%的节点），即这些种子节点的初始状态设为I，剩下节点的初始状态设为S。设定初始时间t=0，产生初始I粒子和SI边的列表，记NI和NSI为I粒子和SI边的数目。

（2）产生均匀分布的随机数r∈[0,1]，计算下一个过程发生的时间间隔更新时间

（4）更新I粒子和SI边的列表。回到（2）。

3 结果分析

定义有效感染系数β=λ μ。图1给出了ER随机网络[8]上SIS模型的感染范围ρ(t)=NI(t)N的3条时间序列，对应于3个不同的β值。当β较小时，最终的感染范围变为零，称之为健康态，该状态是吸收态。吸收态的意思是一旦体系达到这个状态，会永远停留在此态[9]。对SIS模型来说，一旦所有节点的状态都是S态，(5)式两个反应都无法进行。当β较大时，最终的感染范围会维持在某个值附近。由于有限大小网络原因，感染范围会在稳态值附近涨落，网络规模越小，涨落越大。

图1 3个不同感染速率下的ER网络上感染范围ρ(t)的时间序列。网络大小N=1 000，平均度 k=10

图2 给出了β=0.15时两次相继反应之间的等待时间Δt的分布P( )Δt，如果将纵坐标设置成对数值，如图2中的插图所示，可以看出满足指数分布，即这与之前泊松过程假设是一致的。

图2 β=0.15时两次相继反应之间的等待时间Δt的分布P( )Δt

图3 给出了感染范围的稳态平均值 ρ随有效感染系数β的变化图，其中是采样时间长度。当β＜βc时，ρ=0；当β＞βc时，ρ＞0。体系在传播速率阈值β=βc时发生了相变。然而，当β稍稍大于βc时，感染范围并不大，由于有限尺度涨落原因，体系很容易达到吸收态。因此，很难通过图3来准确地定出传播速率阈值，为此，采取准静态方法来确定βc的值。准静态方法的基本思想是：一旦体系达到吸收态，随机地赋予体系的一个采用过的非吸收态构型。通过计算ρ的涨落有关的量，定义为极大值的位置给出传播速率阈值βc的值，如图4所示，βc=0.106。

图3 感染范围的稳态平均值 ρ随有效感染系数β的变化图

图4 χ随有效感染系数β的变化图

4 结论

本文介绍了动力学蒙特卡洛方法的基本原理以及程序实现的基本步骤，并以一个简单的例子——网络上SIS疾病传播过程来说明动力学蒙特卡洛方法的实际应用。通过计算传播范围与传播速率之间的关系来阐述该模型相变的性质，并通过计算涨落相关的物理量来确定相变的位置，即传播速率的阈值。本文对动力学蒙特卡洛方法的基本原理的阐述是考虑了本科生对数学和物理知识掌握的实际程度，相信只要本科生掌握了一些微积分和概率论的知识，本文的推导是完全可以看懂的。此外，本文讨论的是泊松过程，若随机过程的等待时间不满足指数分布，即非泊松过程，也可以使用类似的方法推导[11]。