当前位置:首页 >> 数学 >>

第五章 蒙特卡罗方法(三)


蒙特卡罗方法(三)

5.3 Metropolis 算法
目的:在积分变量 X 的空间(可能是多维的)内产生一组按

概率密度 ω (X ) 分布的点

Metropolis 算法:假想有一个随机行走者在 X 的空间中运动。
该随机行走过程相继各步经过的点产生出一个序列:X0,

X1, …;随着行走的路程越长,这些点的分布会逼近所要求
的分布 ω (X ) 。

注意在这里,行走不是完全无规的。每一步可能达到的位置 依赖且仅依赖于上一步的位置。这种随机行走可以用下列跃

迁概率描述

状态
A

A
1/2

B
1/4

C
0

D
1/4

B
C

0
0

1/3
1

2/3
0

0
0

D

1/2

0

1/2

0

考虑大量行走者从不同的初始点出发在 X 空间独立的随机行
走。若 Nn(X ) 是 n 步后这些行走者在 X 点的密度, Nn(Y) 是

n 步后这些行走者在 Y 点的密度,那么在下一步从 X 点走到 Y 点的行走者净数目为
X 点上的一个行走者转移到Y 的概率



时,行走者的布局不会有

净变化,系统将会达到平衡。这个条件称为细致平衡条件。

为了使行走者的分布达到给定的概率密度分布ω (X ) ,需要规定
适当的跃迁概率

细致平衡条件
一的定出 跃迁概率

并不能唯

有多种具体方案给出跃迁概率,其中最重要的是 Metropolis方案

Metropolis 的方案
Metropolis提出的随机行走的规则 1. 设行走者处于序列中的 Xn 点上,为了产生 Xn+1,行走者迈 出试探性的一步到一新点 Xt 。 该新点可以用任何方便的方 法选取,如可以在点 Xn 周围的一个边长为 d 的多维立方 体中均匀地随机选取。

2. 然后按比值

来决定是“接受”还是“拒绝”该试探步。

?如果 r ≥1,接受该步(即取 X n+1 = X t );

?如果 r <1,则以概率 r 接受这一步:把 r 和一个在[0, 1]区
间上均匀分布的随机数 h 比较,若 h <r 就接受这一步(即取 X
n+1 =

X t ) ,否则就舍弃该步(即取 X n+1 = X n )

这样产生出 X n+1 之后,可以再从X n+1 出发迈出一个试验步,

按照同样的过程产生 X n+2 ,不断重复以上步骤,得到整个随
机行走序列。

证明 Metropopis方案能导致平衡分布 ω (X ) 根据Metropopis方案,从 X 到 Y 的转移概率为

其中 T 是从 X 试探到 Y 的概率,A 是接受这一试探步的概率。 首先试探概率满足

因此,Metropopis随机行走者的平衡分布满足

若 ω (X ) > ω (Y ),则 A (Y→X )=1 , 若 ω (X ) < ω (Y ), 则 A (X→Y )=1 , 在两种情况下,Metropolis 随机行走者的平衡分布都满足

所以行走者最终确实会达到分布 ω 。

其它的选择 细致平衡条件并没有定出跃迁概率的具体形式,因此除了 Metropolis 方案外,还存在其它一些选择,例如 Barker 方法

容易证明 Barker 方法确实满足细致平衡条件。

Metropolis 方法要注意的地方 1. 随机行走从何出发,即选择何处为 X0 。原则上,任何位置

都合适,但实际中,合适的起始点是 ω 值大的地方。
2. 步长 d 如何选取? 太小了关联太强,需要很多步才能达到一个统计独立的构型。 太大了会使大多数尝试都失败,导致更新缓慢。 经验的做法是使得大约的一半试探步入选。

3. 构成随机行走的点 X0 , X1 , X2 , …,由于产生的方法,彼此 不是独立的。

关联函数 关联函数定义为

实际做蒙特卡罗模拟时,需要在生成的点列中进行采样。相 邻两个样本点之间应有足够大的的间隔,使得 C(l) 足够小。

自相关函数图例

应用实例:单原子气体 Metropolis 算法最早被用在经典流体的的模拟中 ( Metropolis, 1953)。这里讨论最简单的经典流体—— 单原子气体。

考虑一定体积,温度为T 的单原子气体,这是一个正则系综,
其物理量 A 的期望值为

其中 R=(r1, r2, …, rN ) 为3N 维坐标矢量,假设粒子间为对势

这里没有考虑 U 对速度的依赖,因为速度的分布由麦克斯韦
分布给出。 任何依赖于速度的物理量可以直接用该分布解析 的计算。下面只考虑依赖于位形的物理量, 如总势能等。

用蒙特卡洛方法计算<A> 的方法为

其中Ri 为相空间中依据下面的分布函数抽样得到的样本点

构型的更新 很多时候,更新一次构型的只改变一个粒子的坐标效率较高, 特别是当系统非常接近平衡态或接近相变时。

随机的选择第 i 个粒子,其坐标被更新为

其中 hk 为第 k 个方向的步长,α 、β、 ? 为 [0, 1] 区间的随机数

这步更新被接受的概率为

新构型中,只有第 i 个粒子的坐标变化了,因此没有必要为了 得到概率 p 而计算整个 ω (Rn+1)

我们将 Δ U 表示为新旧构型的能量差

单原子气体的具体算法

(1). 构造系统的一个初始态 (2). 随机选择一个粒子 i,产生一个试探步 r’i=ri+δ (3). 计算这一试探位移引起的能量变化 ? E (4). 如果 ? E <0, 接受试探步;执行第 (2) 步 (5).如果 ? E >0,生成一个随机数r, 满足0<r<1

(6). 如果 r< exp(- ? E/KT),接受这一试探步;否则拒绝
(7). 执行第 (2) 步。

势函数的截断 为了有效的降低计算量,需要对位势引入截断, rij ≤rc 。 例如简单流体中典型的相互作用为 Lennard-Jones 势

其截断长度可取为 rc =3σ,其中 σ 为势函数为 0 时的距离

长程相互作用, 例如库仑势,不能被截断,需要额外的方法来处理。

步长的选择

?步长 h 如果太小: 粒子只能移动很小的距离, 所以需要
很多步才能达到一个统计独立的构型。 ?步长 h 如果太大:粒子虽然平均来说可以移动更大的距离, 但是绝大多数移动对系统构型的改变过大,以至于能量会显 著的增加。 这导致拒绝率过高,从而构型的更新缓慢。

一个经验法则是接受率平均在 0.4 和 0.6 之间。 对于硬球, 接受率应该低些,约为 0.1。 用来取平均的样本点之间典型 的应有 10-15 个间隔样本点。

应用实例——Ising 模型

经典的二维 Ising 模型:在一个二维正方格子上,每个格点 i 上
有一个自旋,可以取值 +1 或 -1。相邻自旋通过一个交换耦合 能 J 相互作用,此外还存在一个外磁场 B。

系统的哈密顿量为

其中 <i j> 代表对最近邻求和。 伊辛模型最早被用来研究磁相变,另一个有趣的应用是二元合金。

交换耦合能 J

零温下,当没有外磁场时 ? J > 0 : 如果所有自旋都朝同一方向,系统能量会最低, 对应铁磁态 ? J<0:如果相邻两个自旋朝向相反,系统能量会最低。对 应反铁磁态

二维 Ising 模型示意图

周期性边条件

二维 Ising 模型不同温度下典型的自旋构型

二维 Ising 模型的解析结果 Lars Onsager 于1944 年得到二维 Ising 模型的精确解,并

证明存在相变点,这是统计物理学发展过程中的里程碑。其
临界指数为

三维伊辛模型仍然缺乏精确解

需要计算的物理量 磁化强度

内能

比热

磁化率

正则系综下物理量的计算 正则系综下磁化强度的平均值被定义为

其中

为构型 σ 的平均自旋。

磁化强度可以通过蒙特卡罗方法计算

其中 σ=1,2,…,M 根据下面的分布函数取样

随机的在所有格点生成 1 或 -1, 然后随机的选择一个格点
来更新它, 接受率为

其中

j 对 i 的所有邻居求和。
注意与每一个格点相联系的 可以储存并在模拟中更新。 另外,接受概率只有 5 种可能的值,可以事先计算并储存起来,

以避免重复进行指数运算。

二维 Ising 模型的具体算法

(1). 随机生成一个初始构型 (2). 随机选择一个格点 i (3). 计算如果将点 i 的自旋翻转引起的能量变化? E (4). 如果 ? E <0, 接受试探步;执行第 (2) 步 (5).如果 ? E >0,生成一个随机数r, 满足0<r<1 (6). 如果 r< exp(- ? E/KT),则翻转格点 i 的自旋 (7). 执行第二步。

比热随温度变化曲线

20×20 格子

磁化强度随温度变化曲线


赞助商链接
相关文章:
更多相关标签: