第二章数值积分和MonteCarlo方法第一节数值积分令,则零阶近似一

第二章 数值积分和Monte Carlo 方法 第一节 数值积分 ()b

a S f x dx =? 令 10,

,k k n h x x x a x b +=-==， 则

()()()1

1

0(),

''()

k k

n k k k k k k k k x x S f x dx

f x f x x f f f x f f x +-=='=+-+

==∑

?

()x f

a

k x b x

零阶近似

()()h f x f k O +=

()()

()

∑∑-=-=O +=O +=1

1

0n k k n k k h f h h f h S

一阶近似

()()

()

21h h

f f x x f x f k

k k k O +--+=+ ∵

()()?

+=---=-++1

212212

122k k

x x k k k k k k h x x x x

x dx x x

∴

()()

∑-=+O +???

?

?-+?=1

02121n k k k k h h f f h f S

()()

∑

-=+O ++=1

0212

1

n k k k h f f h 从直观看，用()11

2

k k f f ++近似()f x 比只用k f 或1k f +好。这方法也称

Trapezoid 方法。 这样的数值积分方法的

优点：

● 简单直观，误差可以控制

缺点：

● “平均主义”，

在()0≈x f 的区域，()k f x x ?对S 贡献很小，但消耗

同等的机时。在多自由度系统这弱点尤为特出。

问题： 直观地看，零级近似和一级近似的差别在哪？ 习题： 编程序数值计算高斯积分。

第二节 Monte Carlo 方法 如何用随机方法求积分？

例如，可用‘抛石子’方法。但这方法不比简单的数值积分有效。 1．简单抽样的Monte Carlo 方法

均匀地随机地选取[b a ,]中{}k M x 个点，显然，

(1

1()1/

M

k

k S f x M

==

+O ∑

当M 足够大，当然可以得到足够好的积分值。

问题：为什么误差是(1/O ？

答 ：不妨把这看成一个M 次测量的实验，假设每次测量都是独立

的，由涨落理论，误差应为(O 。 比较误差：

Monte Carlo 方法 1/S ?∝ 数值积分Trapezoid 方法 2

1/(1/)S M h M ?∝=

对单自由度而言，数值积分方法要有效得多。

对多自由度，例如d 自由度，

Monte Carlo 方法 1/S ?∝ ！！ 数值积分Trapezoid 方法 2/1/(1/)d

d S M h M ?∝= 当d 非常大，数值积分方法根本没法和Mont

e Carlo 方法比较。 我们当然可以再改进数值积分方法的精度，但这种改进的量级没法和d 的大小比拟。在多体系统的数值模拟中，d 通常至少是410！

Monte Carlo 方法真的是完美了吗？ 当然不是。

● ‘平均主义’的弱点其实还没改进

下面我们引入所谓的重要抽样Monte Carlo 方法 ● 当引入重要抽样方法后，每次抽样的样品可能不独立

如何取得独立的抽样，是Monte Carlo 方法的重点所在！

2．重要抽样的Monte Carlo 方法

如果被积函数f(x)是均匀的函数，则简单抽样方法已经可以得到 相当准确的积分值。

如果被积函数f(x)不是均匀的函数----这在高维积分十分常见，则必须引入重要抽样。我们希望在对积分贡献大的区域多取样品，在贡献小的区域少取样品。

设积分为

()()()b

a S f x f x W x dx =<>=?

其中f(x)是均匀函数，W(x)是不均匀函数，而且()0W x ≥，可以归一化，给予概率分布的含义。

假设我们可以按照分布W(x)得到{}l M x 个点，则

(1

1()1/

M

l

l S f x M

==

+O ∑

注意：现在W(x)不出现在求和式子中，而是体现在{}l x 的分布里。 证明：如第一节方法，把[b a ,]分为n 个小区间， [],l k k x x x h +落在区间的数目约为

()k W x h M ?

对[],l k k x x x h ∈+

()()

()()()1

11

l k M

n l k k l

k f x f x f x f x W x h M

M

M

-=≈∴

≈

∑

∑

显然

()()()(

)1lim lim M n

l

k

k

M n l

k

f x f x W x h f x S M

→∞→∞

==

=∑∑

()W x ，产生M 个{}k x （即上面的{}l x ）

问题：可以用产生随机数的方法产生()W x 吗？ 答 ：多自由度，有相互作用时不可能 3． Markov 过程 产生{}k x 的方法：

构造一个Markov 过程，即给出一个动力学规则，由t x 随机地产生1t x +。那么，给一个0x ，可产生

0001,,

,

M M M x x x x +

如果随机过程满足一定条件，则当0M 足够大时，即达到“平衡

态”时，00

1,M M M x x ++按()W x 分布。

? 各态历经

从概率上说，在有限时间内t x 可走遍[b a ,]

? 细致平衡

Markov 过程由从t x 到1t x +的转移概率定义。

设()T x x '→为过程从x 转移（或跃迁）到x '的概率，()T x x '→为从x '到x 的概率。则细致平衡条件为

()()()

()

T x x W x T x x W x ''→='→

记忆：从x 转移到x '的概率正比与()W x '

证明：设(),W x t 为x 的“非平衡态”分布

()

()()()(),,,x x dW x t T x x W x t T x x W x t dt ''

'''=-→+→∑∑ ()()()

,,,0t W x t W x dW x t dt

→∞

→∞→

()()()(),,x x T x x W x T x x W x '

'

'''∴

→∞=→∞∑∑

细致平衡条件显然保证()W x 满足上述方程，所以，

()(),W x W x ∞=

不过，细致平衡条件是充分条件。 第三节 Metropolis 算法和Heat －bath 算法

Markov 过程的全部信息包含在转移概率()T x x '→中

● 细致平衡条件是平衡态()W x 的要求 ● 是否各态历经常常可以直观判断

()T x x '→必须给出从任意 x 到任意 x ’ 的概率，

但这并不一定要求从x 到任意x ’ 的概率都是零。各态历经只要求在有限时间内，即()T x x '→的有限次作用，能到达x ’。

1． Metropolis 算法

Metropolis (1915-1999) The Paper was cited 7500 times from 1988 to 2003 令

()

()()

()()(')/()1

(1,(')/())

m T x x W x W x W x W x W x W x Min W x W x '→'?

设 ()()W x W x '<

()()(')/()(')

1()

m m T x x W x W x W x T x x W x '→=='→

设 ()()W x W x '≥

()()1(')

()/(')()

m m T x x W x T x x W x W x W x '→=='→

即()m T x x '→满足细致平衡条件。

但是，注意到()T x x '→是转移概率，所以应有归一化条件

()'

1x T x x '→=∑ 上面的()m T x x '→还不满足归一化条件。 所以，完整的Metropolis 算法的转移概率是

()()()m T x x S x x T x x '''→=→→

其中()S x x '→是从x 选中x ’ 的概率，而且()(')S x x S x x '→=→。这样的转移概率仍然满足细致平衡条件。当然，归一化条件也能保证。

例如，可选

1． ()S x x '→=常数，即均匀选取x ’

2． ()'[,]

'[,]

x x x S x x x x x εεεε∈-+?'→=?

?-+?常数

归纳起来，Metropolis 算法包括如下步骤： 1） 设t x x =，按()S x x '→选取尝试x ’，

2） 以概率()m T x x '→取1't x x +=

以概率 ()1m T x x '-→取1t x x += （！！）

第二步要记住，初学者容易误解，以为一旦x ’ 不被采纳，重新选取尝试x ’。

不难理解，无论是1或2方案的()S x x '→，Metropolis 算法都满足各态历经条件，因为只要(')0W x ≠，在有限时间取到x ’ 的概率不为零。

Metropolis 算法非常普遍。

一个重要的例子，在物理学中，常取()W x 为正则分布，

()/()H x kT W x e -∝

其中H(x) 代表能量，T 是温度，k 是Boltzmann 常数。

所以，跃迁概率为

()()()()()()()()

()()()//1

(1,)

H x H x H x H x kT

m kT H x H x H x H x e

T x x Min e '-'

---'>'≤??'→=?

??=

如何在计算机上实现？

① 设已到达t x x =点，按均匀分布随机地取

[],,

x x x εεε'∈-+为一‘恰当’的数

例如，取[]1,0∈r 为计算机上常用的“随机”数，则

()21x x r ε'=+-

② 以概率()m T x x '→取1't x x +=

以概率 ()1m T x x '-→取1t x x +=

例如，计算 ()()x H x H -'，取随机数[]1,0∈r

如果 (

)(

)()/H x H x kT e

r '-

->，取 1't x x +=，否则1t x x +=

()()()/H x H x kT

e '

--

③

如何取0M M 和？

这与具体系统有关。通常通过尝试改变0M 以判断0M 是否够大，即随机过程是否已到达“平衡”态。

一般 010M M >。

不过，对单自由度问题，0M 不重要。 练习： ① ()22

1

x x H σ=

计算2x <>，和解析结果比较

② ()424

1

21x x x H λσ+=

， 计算2x <> ● 讨论ε的作用

● 对没有解析解的系统，如何判断数值结果的可靠性和正确性？

● 有兴趣的同学还可以和数值积分比较，会发现Monte Carlo 方法在单自由度情形没有优势

2． Heat-bath 算法

Heat-bath 算法没有Metropolis 算法那么普遍 但也是多体系统研究中的典型方法之一 Heat-bath 算法直接取转移概率为

()(')h T x x W

x '→= 这算法自然是各态历经并且满足细致平衡条件的。

● Heat-bath 算法的特点是()h T x x '→与x 无关 在某些情形会显示出优越性

● 当(')W x 较复杂，在计算机上不一定能够找到简单的实现 事实上，对单自由度情形，已经是直接产生平衡态分布()W x

例如，如果x 只取 x1和x2, 可取

()(1)

'1

(2)

'2

h W x x x T x x W x x x =?'→=?

=?

计算机上的实现：取随机数[]1,0∈r

如果 (1)W x r >，取 11t x x +=，否则12t x x += 这一算法可简单推广到x 取多个分立值情形。

当x 是连续变量时，不一定能实现Heat-bath 算法。但我们可以结合Metropolis 算法和Heat-bath 算法。

例如，取()()1()h T x x S x x T x x ''→=→→

()1111()/(()())

'()/(()())'0h

W x W x W x x x T x x W x W x W x x x otherwise

+=??'→=+=???

这算法Metropolis 的特征多些。

归纳起来，Metrpolis 算法或Heat-bath 算法可以相当普遍地实现重要抽样的Monte Carlo 模拟方法。但是，对多体系统，抽样得到的样品可能不独立，即有关联。设关联时间为 τ, 即每 τ 个样品才

有一个是独立的。那么，误差为

S O ?∝

当 τ 很大时，Monte Carlo 方法遇到极大困难。对物理学家而言，这是Monte Carlo 算法需要解决的重要问题。

第四节 Langevin 方程

问题：系统的Hamiltonian 量为S ，

平衡态分布为 )exp(S -，（这里温度已吸收到S ）。

假设系统0=t 时处于一初始状态，系统如何演化至平衡态？ 单自由度

()

实数:x x S S =

Langevin 方程

()()()t x

x S dt t dx η+??-=

()

()()

():0

2t t t t t η

η

ηηηηδ''==-高斯随机数

对固定t

()ηP ~ 2

ησ-e

()??

-∞

+∞

--==??=2

2

1ησηση

ηηηηe

d Z e

d Z

t

如果没有随机力，平衡态为

()()

0dx t S x dt x

?=-=?，即能量取极小值。如果存在随机力，体系会被推离能量极小，处于某种能量较高的平衡态。

由于随机力的存在，Langevin 方程存在复杂性，因为我们必须考虑对随机力带来的奇异性。为了简单起见，我们对 时间分立化 在数值模拟中应用

较直观

()()t t t t t '?=

'δηηη2，???'

≠'

=='t t t t t t 0

1

δ Z ~

2

1

ησ

=

σ

σ21

ln =??-

Z ∴ 4

t

?=σ

Langevin 方程

()()()()

()()S x t x t x t t x t t t t x

η??=+?-=-

?+?? 令

ηη→?2

t

∴ ()(')

t t t t ηηδ'=

()()()()

()t t t t x t x S t x η?+???-

=?2 方程的解 ()t x η是随机变量，在数值模拟中给定初始值()t x x η,0还不确定，与随机力有关。也就是说，在t 时刻，x 遵从一个分布

();P x t 。

物理量()x O 的平均值

()()()();x t dx P x t x ηη=O ?

问题：()()x t ηη的含义？ 答：必须对t 之前的所有随机力做平均。 Langevin 方程其实也是一种Monte Carlo 方法 设 ')(,

)(x t t x x t x =?+=

()]

4/))('(exp[]

2/)(exp['22t t x

x S x x t x x P ????+--∝-∝→η

()]

4/)')'('(ex p[]

2/)(ex p['22t t x x S x x t x x P ????+--∝-∝→η

注意到x 和x ’之差是无穷小，

)]

()'((exp[]

)

()'(exp[)'()'(x S x S x x S x x x x P x x P --=??--=→→

Fokker-Planck 方程的推导和平衡态

()()()()

()()

()

()()()()

()() +??O ?+

??O ?=

O ?=O ?2

22

21t x t x t x t x t x t x t x t x ηηηηηηη

ηη

η

()()()()()()()()()()()()()()η

ηηηηη

ηηηηηη2

22

22

12???

? ???+???-?O ?+???

? ???+???-?O ?=t t t t x t x S t x t x t t t t x t x S t x t x

∵()()()以及更早的随机力有关只与无关，

与t t t t x ?-ηηη

∴

()()

()

()()()()()()()↑

?????O ?-

=??O ?t t x t x S t x t x t x t x t x ηηηηηηη

∵ ()0=t η

又 ∵ ()()()

22

2t t t x ?O +?=? η

η

∴

()()η

η

η2

2x x S x t

t x ?O

?+

???O ?-=?O ?

()()()()t x P x

t x P x S x x x dx x x S x t x P dx ;;;2222还作用于分步积分

??

↑

↑????

??????+??O =?

???

???O ?+???O ?-=??

这里做分步积分时，假设0);(=±∞

t P 另一方面

()()

()

()

???O =?O ?t

t x P x dx t

t x ,η

η Fokker-Planck 方程

()()()()0;,

;;→??∞→???? ????+????-=-=??t

t x P t x x S x x H t x P H t

t x P FP FP 当 ∴

()0;=∞x P H FP

显然 ()∞;x P ~ ()

x S e -