(完整版)蒙特卡洛算法详讲

Monte Carlo 法
§8.1 概述
Monte Carlo 法不同于前面几章所介绍的确定性数值方法，它是用来解决数
学和物理问题的非确定性的（概率统计的或随机的）数值方法。

Monte Carlo 方
法（MCM ），也称为统计试验方法，是理论物理学两大主要学科的合并：即随机
过程的概率统计理论（用于处理布朗运动或随机游动实验）和位势理论，主要是
研究均匀介质的稳定状态[1]。

它是用一系列随机数来近似解决问题的一种方法，
是通过寻找一个概率统计的相似体并用实验取样过程来获得该相似体的近似解
的处理数学问题的一种手段。

运用该近似方法所获得的问题的解in spirit 更接近
于物理实验结果，而不是经典数值计算结果。

普遍认为我们当前所应用的MC 技术，其发展约可追溯至1944年，尽管在
早些时候仍有许多未解决的实例。

MCM 的发展归功于核武器早期工作期间Los
Alamos （美国国家实验室中子散射研究中心）的一批科学家。

Los Alamos 小组
的基础工作刺激了一次巨大的学科文化的迸发，并鼓励了MCM 在各种问题中的
应用[2]-[4]。

“Monte Carlo ”的名称取自于Monaco （摩纳哥）内以赌博娱乐而闻名
的一座城市。

Monte Carlo 方法的应用有两种途径：仿真和取样。

仿真是指提供实际随机
现象的数学上的模仿的方法。

一个典型的例子就是对中子进入反应堆屏障的运动
进行仿真，用随机游动来模仿中子的锯齿形路径。

取样是指通过研究少量的随机
的子集来演绎大量元素的特性的方法。

例如，)(x f 在b x a <<上的平均值可以
通过间歇性随机选取的有限个数的点的平均值来进行估计。

这就是数值积分的
Monte Carlo 方法。

MCM 已被成功地用于求解微分方程和积分方程，求解本征
值，矩阵转置，以及尤其用于计算多重积分。

任何本质上属随机组员的过程或系统的仿真都需要一种产生或获得随机数
的方法。

这种仿真的例子在中子随机碰撞，数值统计，队列模型，战略游戏，以
及其它竞赛活动中都会出现。

Monte Carlo 计算方法需要有可得的、服从特定概
率分布的、随机选取的数值序列。

§8.2 随机数和随机变量的产生
[5]－[10]全面的论述了产生随机数的各类方法。

其中较为普遍应用的产生随
机数的方法是选取一个函数)(x g ，使其将整数变换为随机数。

以某种方法选取
0x ，并按照)(1k k x g x =+产生下一个随机数。

最一般的方程)(x g 具有如下形式：
m c ax x g mod )()(+=
（8.1）
其中
=0x 初始值或种子（00>x ） =a 乘法器（0≥a ）
=c 增值（0≥c ）
=m 模数
对于t 数位的二进制整数，其模数通常为t 2。

例如，对于31位的计算机m 即可
取1312-。

这里a x ,0和c 都是整数，且具有相同的取值范围0,,x m c m a m >>>。

所需的随机数序{}n x 便可由下式得
m c ax x n n m od )(1+=+
（8.2）
该序列称为线性同余序列。

例如，若70===c a x 且10=m ，则该序列为
7，6，9，0，7，6，9，0…… （8.3）
可以证明，同余序列总会进入一个循环套；也就是说，最终总会出现一个无休止
重复的数字的循环。

（8.3）式中序列周期长度为4。

当然，一个有用的序列必是
具有相对较长周期的序列。

许多作者都用术语乘同余法和混合同余法分别指代
0=c 和0≠c 时的线性同余法。

选取c a x ,,0和m 的法则可参见[6,10]。

这里我们只关心在区间)1,0(内服从均匀分布的随机数的产生。

用字符U 来
表示这些数字，则由式（8.2）可得
m
x U n 1-=
（8.4）
这样U 仅在数组{}m m m m /)1(,......,/2,/1,0-中取值。

（对于区间（0,1）内的随机
数，一种快速检测其随机性的方法是看其均值是否为0.5。

其它检测方法可参见
[3,6]。

）产生区间),(b a 内均匀分布的随机数X ，可用下式
U a b a X )(-+=
（8.5）
用计算机编码产生的随机数（利用式（8.2）和（8.4））并不是完全随机的；
事实上，给定序列种子，序列的所有数字U 都是完全可预测的。

一些作者为强调
这一点，将这种计算机产生的序列称为伪随机数。

但如果适当选取c a ,和m ，序
列U 的随机性便足以通过一系列的统计检测。

它们相对于真随机数具有可快速产
生、需要时可再生的优点，尤其对于程序调试。

Monte Carlo 程序中通常需要产生服从给定概率分布)(x F 的随机变量X 。

该步可用[6]，[13]－[15]中的几种方法加以实现，其中包括直接法和舍去法。

直接法（也称反演法或变换法），需要转换与随机变量X 相关的累积概率函
数)()(x X prob x F ≤=（即：)(x F 为x X ≤的概率）。

1)(0≤≤x F 显然表明，通
过产生(0,1)内均匀分布随机数U ，经转换我们可得服从)(x F 分布的随机样本X 。

为了得到这样的具有概率分布)(x F 的随机数X ，不妨设)(x F U =，即可得
)(1U F X -= （8.6）
其中X 具有分布函数)(x F 。

例如，若X 是均值为μ呈指数分布的随机变量，且 ∞<<-=-x e x F x 0,1)(/μ （8.7）
在)(x F U =中解出X 可得
)1ln(U X --=μ （8.8）
由于)1(U -本身就是区间(0,1)内的随机数，故可简写为
U X ln μ-= （8.9）
有时(8.6)式所需的反函数)(1x F -不存在或很难获得。

这种情况可用舍去法来处理。

令dx
x dF x f )()(=为随机变量X 的概率密度函数。

令b x a >>时的0)(=x f ，且)(x f 上界为M （即：M x f ≤)(），如图8.1所示。

我们产生区间(0,1)内的两个随机数),(21U U ，则
11)(U a b a X -+= M U f 21= （8.10）
分别为在(a,b)和(0,M)内均匀分布的随机数。

若
)(11X f f ≤ （8.11）
则1X 为X 的可选值，否则被舍去，然后再试新的一组),(21U U 。

如此运用舍去法，所有位于)(x f 以上的点都被舍去，而位于)(x f 上或以下的点都由11)(U a b a X -+=来产生1X 。

图8.1 舍去法产生概率密度函数为)(x f 的随机变量
例8.1 设计一子程序使之产生0,1之间呈均匀分布的随机数U 。

用该程序产生随机变Θ，其概率分布由下式给定
πθθθ<<-=0),cos 1(2
1)(T 解：生成U 的子程序如图8.2所示。

该子程序中，,0,21474836471221==-=c m 且1680775==a 。

应用种子数（如1234），主程序中每调用一次子程序，就会生成一个随机数U 。

种子数可取1到m 间的任一整数。

0001
C**********************************************************
0002 C PROGRAM FOR GENERATING RANDOM V ARIABLES WITH
0003 C A GIVEN PROBABILITY DISTRIBUTION
0004
C**********************************************************
0005
0006 DOUBLE PRECISION ISEED
0007
0008 ISEED=1234.D0
0009 DO 10 I=1,100
0010 CALL RANSOM(ISEED,R)
0011 THETA=ACOSD(1.0-2.0*R)
0012 WRITE(6,*)I,THETA
0013 10 CONTINUE
0014 STOP
0015 END
0001
C**********************************************************
0002 C SUBROUTINE FOR GENERATING RANDOM NUMBERS IN
0003 C THE INTERV AL (0,1)
0004
C**********************************************************
0005
0006 SUBROUTINE RANDOM (ISEED,R)
0007 DOUBLE PRECISION ISDDE,DEL,A
0008 DATA DEL,A/2147483647.D0,16807.D0/
0009
0010 ISDDE=DMOD(A*ISDDE,DEL)
0011 R=ISDDE/DEL
0012 RETURN
0013 END
图8.2 例8.1的随机数生成器
图8.2的子程序只是为了说明本章所介绍的一些概念。

大多数计算机都有生成随机数的子程序。

为了生成随机变量Θ，令
)cos 1(2
1)(Θ-=Θ=T U 则有 )21(cos )(11U U T -==Θ--
据此，一系列具有给定分布的随机变量Θ便可由图8.2所示主程序中生成。

§8.3 误差计算
Monte Carlo 程序给出的解按大量的检测统计都达到了平均值。

因此，该解中包含了平均值附近的浮动量，而且不可能达到100％的置信度。

要计算Monte Carlo 算法的统计偏差，就必须采用与统计变量相关的各种统计方法。

我们只简要介绍期望值和方差的概念，并利用中心极限定理来获得误差估计[13,16]。

设X 是随机变量。

则X 的期望值或均值x 定义为
⎰∞
∞-=dx x xf x )( （8.12）
这里)(x f 是X 的概率密度分布函数。

如果从)(x f 中取些独立的随机样本N x x x ,...,,21，那么的x 估计值就表现为N 个样本值的均值。

∑==N n n x
N x 11ˆ
（8.13）
x 是X 的真正的平均值，而x
ˆ只是x 的有着准确期望值的无偏估计。

虽然x ˆ的期望值等于x ，但x x
≠ˆ。

因此，我们还需要x ˆ的值在x 附近的分布测度。

为了估计X 以及x
ˆ在x 附近的的值的分布，我们需要引入X 的方差，其定义为X 与x 差的平方的期望值，即
⎰∞
∞--=-==dx x f x x x x x Var )()()()(222σ （8.14） 由2222)(x x x x x x +-=-，故有
⎰⎰⎰∞
∞-∞∞-∞∞-+-=dx x f x dx x xf x dx x f x x )()(2)()(222σ （8.15）
或者 222)(x x x -=σ （8.16）
方差的平方根称为标准差，即 2/122)()(x x x -=σ （8.17）
标准差给出了x 在均值x 附近的分布测度，并由此给出了误差幅度的阶数。

x
ˆ的标准差与x 的标准差的关系表示为
N x x )
()ˆ(σσ=
（8.18）
该式表明，如果用根据(8.13)式由N 个n x 值构成的x
ˆ来估计x ，那么结果中x ˆ在x 附近的扩散范围便与)(x σ成比例，且随着样本数N 的增加而降低。

为了估计x
ˆ的扩散范围，我们定义样本方差为 ∑=--=N n n x x N S 122
)ˆ(11 （8.19）
由此式还可看出2S 的期望值等于)(2x σ。

因此样本方差是)(2x σ的无偏估计。

将（8.19）式中平方项乘出来，便可得样本标准差为
2/11222/1ˆ11⎥⎦⎤⎢⎣⎡-⎪⎭
⎫ ⎝⎛-=∑=N n n x x N N N S
（8.20） 当N 较大时，系数)1/(-N N 可设为1。

作为获得中心极限定理的一种方法，概率论的一个基本解可考虑二次项函数 M N M q p M N M N M B --=
)!
(!!)( （8.21）
该式表明N 次独立随机试验中有M 次成功的概率。

(8.21)中，p 是一次试验中成
功的概率，且p q -=1。

当M 和M N -都较大时，便可用Stirling 公式
n e n n n n π2~!- （8.22）
因此(8.21)式可近似为正态分布[17]
⎥⎦⎤⎢⎣⎡--=≈)ˆ(2)ˆ(exp 2)ˆ(1
)ˆ()(22x x x x x f M B σπσ (8.23) 其中p N x =且Npq x =)ˆ(σ。

因此当∞→N 时，中心极限定理表明，描述由N 点
Monte Carlo 算法获得的x
ˆ的分布的概率密度函数是（8.23）式所示的正态分布函数)(x f 。

也就是说，大量随机变量的集合趋于呈正态分布。

将（8.18）式代入（8.23）式可得
⎥⎦
⎤⎢⎣⎡--=)(2)ˆ(exp )(12)ˆ(22x x x N x N x f σσπ （8.24）
正态（高斯）分布在工程，物理以及统计学的各类问题中都非常有用。

高斯模型的显著特性源于中心极限定理。

因此，高斯模型经常用于其影响程度取决于由许多不规则的、浮动的元素构成的集合的情况。

例8.2中我们给出了根据中心极限定理产生高斯随机变量的算法。

由于样本数N 是有限的，所以Monte Carlo 算法不可能达到绝对的确定性。

故而在x 附近估计出某一范围或区间以确保我们估计的x
ˆ落入该范围内。

假设要得到x
ˆ位于δ-x 和δ+x 之间的概率。

由定义 {}⎰+-=+<<-δ
δδδx x x d x
f x x x ob ˆ)ˆ(ˆPr （8.25） 令()()
x N x x σλ/2ˆ-=可得 {}()
()⎰-=+<<-σδλλπδδ/2/02
2
ˆPr N d e x x x ob ()⎪⎪⎭⎫ ⎝
⎛=x N erf σδ2/ （8.26a ）
或
ασσαα-=⎭⎬⎫⎩
⎨⎧+≤≤-1ˆPr 2/2/N z x x N z x ob
（8.26b ）
其中)(x erf 是误差函数，且2/αz 是标准正态差的上2/α分位点。

对（8.26）式可做如下解释：如果用来呈现独立随机观测值并构建相关随机区间δ±x 的Monte
Carlo 程序以较大的N 值反复运行，则这些随机区间的()⎪⎪⎭
⎫ ⎝⎛x N erf σδ2分位点将
近似等于x ˆ。

随机区间δ±x 称为置信区间，()⎪⎪⎭
⎫ ⎝⎛x N erf σδ2称为置信度。

大多数
的Monte Carlo 算法都使用误差()N x /σδ=，这表明x ˆ是在实际均值x 的标准差范围内的。

由（8.26）式可得，样本均值x
ˆ位于区间()N x x /ˆσ±内的概率是0.6826或68.3％。

若要求更高的置信度，可用两个或三个标准差的值。

如
⎪⎩⎪⎨⎧=⎭⎬⎫⎩⎨⎧+<<-,997.0,954.0,6826.0)(ˆ)(Pr N x M x x N x M x ob σσ 3
21===M M M （8.27）
其中M 是标准差的个数。

在（8.26）和（8.27）式中均假设总体标准差σ为已知。

由于这种情况很少出现，故必须用由（8.20）式算得的样本标准差S 来估计σ的值，从而使学生氏t －分布取代正态分布。

我们知道当N 很大，比如30>N 时，t －分布近似趋于正态分布。

（8.26）式等价于
ααα-=⎭⎬⎫⎩
⎨⎧+≤≤---1ˆPr 1;2/1;2/N St x x N St x ob N N （8.28）
其中1;2/-N t α为自由度是1-N 的学生氏t －分布的上2/α分位点；其值在任何标准统计学课本中都有表列可查。

这样置信区间的上、下限便可由下式给出 上限＝N St x N 1;2/-+α （8.29）
下限＝N St x N 1
;2/--α
（8.30）
Monte Carlo 算法中关于误差估计的进一步讨论，可参见[18,19]。

例8.2 利用中心极限定理产生一高斯（正态）分布的随机变量X 。

根据中心极限定理，大量均值附近的独立随机变量的总和，无论其个体变量的分布如何，
总趋向于一种高斯分布。

也就是说，对于任意随机数N i Y i ,....2,1,=，均值为Y ，方差为)(Y Var ，
()
Y Var N Y N Y Z N i i -=
∑=1 （8.31） 渐渐与N 合并为零均值、单位标准差的正态分布。

若i Y 是均匀分布的随机变量（即i i U Y =），则12/1)(,2/1==Y Var Y ，故而
12/2/1N N U Z N i i -=
∑=
（8.32）
且变量 μσ+=Z X （8.33）
近似等于均值为μ、方差为2σ的正态变量。

N 小于3时近似为我们所熟知的钟形高斯分布。

为简化计算，通常实际中设12=N ，因为这样可消除（8.32）式中的平方根项。

然而N 的这种取值截掉了σ6±边界处的分布，且无法产生超过σ3的值。

对于分布曲线的两端比较重要的仿真，就必须用其它方案来产生高斯分布（参见[20]－[22]）。

这样，要产生一个均值为μ、标准差为σ的高斯随机变量X ，就要遵循以下步骤：
（1） 生成12个均匀分布的随机数1221,....,,U U U
（2） 求得612
1-=∑=i i U Z
（3） 令μσ+=Z X
§8.4 数值积分
对于一维积分，现已有一些求积公式，如3.10节中所述。

而对多重积分的公式则相对较少。

Monte Carlo 技术对此类多重积分的重要性至少存在两方面的原因。

多重积分的求积公式非常复杂，而多重积分的MCM 几乎保持不变。

Monte Carlo 积分的收敛性与维数无关，但对求积公式并非如此。

人们已经发现，积分
的统计方法是计算天线问题尤其是那些与较大结构相关的问题中的二维或三维积分的很有效的方法[23]。

这里将讨论两类Monte Carlo 积分方法，即简易MCM 和具有对立变量的MCM 。

其它类型，如成功－失败法和控制变量法，可参见[24]－[26]。

本节还将简要介绍MCM 在不规范积分中的应用。

§8.4.1 简易Monte Carlo 积分
设要计算积分
⎰=R
f I （8.34）
其中R 是n 维空间。

令()n X X X X ,....,21=是在R 内均匀分布的随机变量。

则)(X f 也是随机变量，其均值由下式给出[27,28]
R I f R X f R ==
⎰1)( (8.35)
方差为
()()2
211⎪⎪⎭⎫ ⎝⎛-=
⎰⎰R R f R f R X f Var （8.36） 其中
⎰=R
dX R （8.37） 如果取X 的N 个独立样本，即N X X X ,....,,21，它们与X 具有相同的分布，且构成平均项
()()()()∑==+++N
i i
N X f N N X f X f X f 1211.... （8.38）
我们便会期望该平均项接近于)(X f 的均值。

这样，由（8.35）和（8.38）式，即可得
()∑==N
i i
X f N R I 1 （8.39）
Monte Carlo 公式可用于有限区域R 上的任何积分。

为了举例说明，现将（8.39）式应用于一维和二维积分中。

对于一维积分，设
⎰=b
a
dx x f I )(
（8.40）
由（8.39）式可得
()∑=-=N
i i X f N a b I 1
（8.41）
其中i X 是区间()b a ,内随机变量，即
()10,<<-+=U U a b a X i （8.42）
对于二维积分，有
()
⎰⎰=b a
d c
dX dX X X
f I 2121
,
（8.43）
则相应的Monte Carlo 公式为
()()
∑=--=N
i i i X X f N c d a b I 1
21,)( （8.44）
其中
1
0,)(10,)(2
2
2111<<-+=<<-+=U U c d c X U U a b a X i
i
（8.45） （8.39）式中无偏估计值I 收敛的很慢，这是由于估计值方差的级数是N /1。

一种改良的方法，即控制变量法，可通过减小估计值方差来提高其准确性和收敛性。

§8.4.2 具有对立变量的Monte Carlo 积分方法 术语对立变量[29,30]用来描述任一系列估计值，它们能互相抵消掉彼此的方差。

方便起见，我们假设积分范围为区间(0,1)。

设要得到下面单重积分的估计值 ⎰=1
0)(dU U g I
（8.46） 我们期望表达式
[])1()(2
1
U g U g -+与)(U g 相比具有更小的方差。

如果)(U g 太小，那反过来)1(U g -就很可能太大。

因此，我们定义估计值
()()[]∑=-+=
N
i i i U g U g N
I 112
1
1
其中i U 是0和1之间的随机数。

此估计值方差的级数是4
1
N ，这是在（8.39）式基础上的一个极大的进步。

对于两维积分 (
)
⎰⎰=101
2121,dU dU U U g I
（8.48）
则相应的估计值为
()()()()[]
∑=--+-+-+=
N
i i i i i i i i i U U g U U g U U g U U g N
I 1212121211,1,11,,4
1
1
（8.49）
根据相似性，可将该方法延拓至更高阶的积分。

对于不是(0,1)的区间，可应用诸如（8.41）式到（8.45）式的转换。

例如，
()()()⎰⎰-=1
dU U g a b dx x f b a
()()[]∑=-+-≈N i i i U g U g N a b 112
1
（8.50）
其中)()(X f U g =，且()U a b a X -+=。

由（8.47）和（8.49）式可得，当积分维数增加时，每一维用来在简易Monte Carlo 方法上提高效率的对立变量的最小值也会增加。

这样使得简易Monte Carlo 方法在多维运算中更具优越性。

§8.4.3 不规范积分 积分式
⎰∞
=0)(dx x g I
（8.51）
可用Monte Carlo 仿真法进行估计[31]。

对于概率密度为)(x f 的随机变量X ，其中)(x f 在区间(0,∞)上的积分等于1，则有
()()
⎰⎰
∞∞
=00
)(dx x g dx x f x g （8.52）
因此，为计算（8.51）式中的I 值，首先要得到N 个服从同一在区间(0,∞)上的积分等于1的概率密度函数)(x f 的独立随机变量。

其样本均值
()()
()
∑
==
N
i i i x f x g N
x g 1
1
便是I 的估计值。

例8.3 用Monte Carlo 方法计算积分 ⎰
⎰
=1020
cos πφαρφρρd d e I j
解：该积分式表示振幅和相位分别服从某一分布的圆形孔径天线的辐射状况。

之所以选择该式，是因为它至少是辐射积分的一部分。

且其解的闭合形式亦为可得，由此便可得Monte Carlo 解的精确性。

其闭合解为 ()()
α
απα12J I =
其中()α1J 是一阶Bessel 函数。

图8.3给出由（8.44）和（8.45）式计算该积分的一个简单的程序，其中
0,1,0===c b a 以及π2=d 。

该程序在Vax 11/780 中调用子程序RANDU 来产
生随机数1U 和2U 。

对于不同的N 值，简易和对立变量Monte Carlo 方法都可用于计算辐射积分。

表8.1中对5=α时的结果进行了精确的比较。

在应用（8.49）式时，用到了下面的对应式：
()()111221111,,U a b X b U X U X U --=-≡-≡≡ ()()22211U c d X d U --=-≡-
表8.1 例8.3辐射积分的MC 方法积分结果
N 简易MCM 对立变量MCM 500 -0.2892-j0.0742
-0.2887-j0.0585
1000 -0.5737+j0.0808
-0.4982-j0.0080
2000 -0.4922-j0.0040
-0.4682-j0.0082
4000 -0.3999-j0.0345
-0.4216-j0.0323
6000 -0.3608-j0.0270
-0.3787-j0.0440
8000 -0.4327-j0.0378
-0.4139-j0.0241
10,000 -0.4229-j0.0237
-0.4121-j0.0240 精确解：-0.4116+j0
§8.5 位势问题的解
位势理论与布朗运动（或随机游动）的关系是由Kakutani在1944年首次提出的[32]。

自此，所谓的概率位势理论便在诸如热传导[33]－[38]、静电学[39]－[46]以及电力工程等许多学科领域得到广泛应用。

不同方程式的概率解或MC解的一个基本概念就是随机游动。

不同类型的随机游动应用不同的Monte Carlo方法。

最常见的类型是固定随机游动和自动定位随机游动。

其它不常见类型有迁离法、缩减边界法、内切形法以及表面密度法。

0001
C***************************************************************** 0002 C INTEGRATION USING CRUDE MONTE CARLO
0003 C AND ANTITHETIC METHODS
0004 C ONLY FEW LINES NEED BE CHANGED TO USE THIS PROGRAM
0005 C FOR ANY MULTI—DIMENSIONAL INTEGRATION
0006
C**************************************************************** 0007
0008 DATA IS1,IS2,IS3,IS4/1234,5678,9012,3456/
0009 DATA A,B,C/0.0,1.0,0.0/
0010
0011 C SPECIFY THE INTEGRAND
0012 F(RHO,PHI)=RHO*CEXP(J*ALPHA*RHO*COS(PHI))
0013
0014 J=(0.0,1.0)
0015 ALPHA=5.0
0016 PIE=3.1415927
0017 D=2.0*PIE
0018 DO 30 NRUN=500,10000,500
0019 SUM1=(0.0,0.0)
0020 SUM2=(0.0,0.0)
0021 DO 10 I=1,NRUN
0022 CALL RANDU(IS1,IS2,U1)
0023 CALL RANDU(IS3,IS4,U2)
0024 X1=A+(B-A)*U1
0025 X2=C+(D-C)*U2
0026 X3=(B-A)*(1.0-U1)
0027 X4=(D-C)*(1.0-U2)
0028 SUM1=SUM1+F(X1,X2)
0029 SUM2=SUM2+F(X1,X2)+F(X1,X4)+F(X3,X2)+F(X3,X4)
0030 10 CONTINUE
0031 AREA1=(B-A)*(D-C)*SUM1/FLOAT(NRUN)
0032 AREA2=(B-A)*(D-C)*SUM2/(4.0*FLOAT(NRUN)) 0033 PRINT *,NRUN,AREA1,AREA2 0034 WRITE(6,*) NRUN,AREA1,AREA2 0035 WRITE(6,20) NRUN,AREA1,AREA2 0036
20
FORMAT(2X,’NRUN=’,I5,3X,’AREA1=’,F12.6,3X,F12.6,’AREA2=’, 0037 1 F12.6,3X,F12.6,/) 0038 30 CONTINUE 0039 STOP 0040 END
图8.3 例8.3中用Monte Carlo 方法计算二维积分的程序 §8.5.1 固定随机游动
为具体起见，假设用固定随机游动的MCM 解拉普拉斯方程
02=∇V （在区域R 内） （8.54a ）
满足Dirichlet 边界条件
P V V =
（在边界
B 上）
（8.54b ）
首先将R 分成网状结构，并用其有限差当量替代2∇。

（8.54a ）在二维R 内的有限差表达式可由（3.26）式给出，即
()()()()()∆-+∆++∆-+∆+=-+-+y x V p y x V p y x V p y x V p y x V y y x x ,,,,, （8.55a ） 其中
4
1
=
===-+-+y y x x p p p p （8.55b ）
（8.55）式中，假设网络的一个方格尺寸是∆，如图8.4所示。

下面给出该方程的概率解释。

若随机游动粒子在某一瞬时处于()y x ,位置处，它从此点向),(),,(y x y x ∆-∆+),(∆+y x 及()∆-y x ,移动的概率分别是+-+y x x p p p ,,和-y p 。

决定粒子移动方式的方法是产生一个随机数10,<<U U ，并按下面的方式指导粒子的随机游动：
()()
()()
()()()()
∆-→∆+→∆-→∆+→y x y x y x y x y x y x y x y x ,,,,,,,,
1
75.075.05.05.025.025
.00<<<<<<<<U U U U
（8.56）
如果不用方格而用矩形格，则有-+=x x p p 且-+=y y p p ，但y x p p ≠。

在用立方体晶格表示的三维问题中，6
1
======-+-+-+z z y y x x p p p p p p 。

两种情况中都依据概率对区间10<<U 进行了分组。

图8.4 固定随机游动示意图
为了计算()y x ,处的位势，让一随机游动粒子在该点出发开始游动。

粒子便开始从一点到另一点在网格中蜿蜒前行直到到达边界。

此时，粒子终止游动，记下边界点的指定位势P V 。

令第一个粒子游动结束时P V 的值记为()1P V ，如图8.4所示。

然后将第二个粒子从()y x ,释放，使其游动直到到达边界点，游动终止且该点P V 的值相应的记为()2P V 。

依次第三个，第四个， ，第N 个粒子从()y x ,释放重复此过程，并记下相应的指定位势()()4,3P P V V ， ()N V P 。

根据Kakutani [32]，()()2,1P P V V ， ()N V P 的期望值是Dirichlet 问题在()y x ,的解，即
()()∑==
N
i P i V N
y x V 1
1,
（8.57）
其中N 较大，为游动总次数。

其收敛速度随N 而改变，所以需要保证许多随机游动都有精确的结果。

若要解Poisson 方程
()y x g V ,2-=∇
在
R 内
（8.58a ） 且V 满足条件
P V V =
在
B 上
（8.58b ）
则式（8.55）中的有限差分表示为
()()()y x V p y x V p y x V x x ,,,∆-+∆+=-+
()()4
,,2g
y x V p y x V p y y ∆+∆-+∆++-+
（8.59）
其中概率仍为式（8.55b ）所示。

（8.59）式的概率解释也近似于（8.55）式。

但随机游动的每一步都要记录下（8.59）式中的4/2g ∆项。

如果第i 次随机游动从
()y x ,出发至到达边界需要i m 步，则记录下
()()∑=∆+
i
m j j
j
P y x g i V 1
2,4
（8.60）
由此可得，()y x V ,的Monte Carlo 解为
()()()∑∑∑=-==⎥⎦⎤⎢⎣⎡∆+=N i m j j j N i P i y x g N i V N y x V 11
1
21,41,
（8.61）
对刚才所描述MCM 的一个有趣的分析就是游走醉鬼问题[15,35]。

我们将随机游动粒子当作一个“醉鬼”，将网格方块当作“城市街区”，结点当作“十字路口”，边界B 当作“城市边界”，而且把边界B 上的结点当作“警察”。

尽管醉鬼想步行走回家，但他却醉的一塌糊涂以至于在整个城市中随机游走。

警察的工作是在城市边界上一看见醉鬼便将其抓获，并勒令醉鬼交付罚金P V 。

那么醉鬼要支付的预期罚金将会是多少呢？其答案便是（8.57）式。

在电介质边界上，指定边界条件为n n D D 21=。

()y x ,
x
ˆx ˆx ˆx ˆx ˆx ˆx ˆx ˆx ˆ x x x x x x x x x x x x x x
X X X X X X X X X X X X X X X X X X X X
[1].R.Hersch and R.J.Griego,“Brownian motion and potential theory,”
Sci.Amer.,Mar.1969,pp.67-74.。