蒙特卡洛模拟

第八章 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 方法的应用有两种途径：仿真和取样。仿真是指提供实际随机现象的数学上的模仿的方法。一个典型的例子就是对中子进入反应堆屏障的运动进行仿真，用随机游动来模仿中子的锯齿形路径。取样是指通过研究少量的随机的子集来演绎大量元素的特性的方法。例如，f(x)在a?x?b上的平均值可以通过间歇性随机选取的有限个数的点的平均值来进行估计。这就是数值积分的Monte Carlo 方法。MCM已被成功地用于求解微分方程和积分方程，求解本征值，矩阵转置，以及尤其用于计算多重积分。

任何本质上属随机组员的过程或系统的仿真都需要一种产生或获得随机数的方法。这种仿真的例子在中子随机碰撞，数值统计，队列模型，战略游戏，以及其它竞赛活动中都会出现。Monte Carlo 计算方法需要有可得的、服从特定概率分布的、随机选取的数值序列。

§8.2 随机数和随机变量的产生

[5]－[10]全面的论述了产生随机数的各类方法。其中较为普遍应用的产生随机数的方法是选取一个函数g(x)，使其将整数变换为随机数。以某种方法选取x0，并按照xk?1?g(xk)产生下一个随机数。最一般的方程g(x)具有如下形式：

其中

g(x)?(ax?c)modm （8.1）

x0?初始值或种子（x0?0）

a?乘法器（a?0）

c? 增值（c?0）

m?模数

对于t数位的二进制整数，其模数通常为2。例如，对于31位的计算机m即可取2t31?1。这

里x0,a和c都是整数，且具有相同的取值范围m?a,m?c,m?x0。所需的随机数序?xn?便可由下式得

xn?1?(axn?c)modm （8.2） 该序列称为线性同余序列。例如，若x0?a?c?7且m?10，则该序列为

7，6，9，0，7，6，9，0…… （8.3） 可以证明，同余序列总会进入一个循环套；也就是说，最终总会出现一个无休止重复的数字的循环。（8.3）式中序列周期长度为4。当然，一个有用的序列必是具有相对较长周期的序列。许多作者都用术语乘同余法和混合同余法分别指代c?0和c?0时的线性同余法。选取x0,a,c和m的法则可参见[6,10]。

这里我们只关心在区间(0,1)内服从均匀分布的随机数的产生。用字符U来表示这些数字，则由式（8.2）可得

U?xn?1 （8.4） m这样U仅在数组?0,1/m,2/m,......,（对于区间（0,1）内的随机数，一(m?1)/m?中取值。种快速检测其随机性的方法是看其均值是否为0.5。其它检测方法可参见[3,6]。）产生区间

(a,b)内均匀分布的随机数X，可用下式

X?a?(b?a)U （8.5） 用计算机编码产生的随机数（利用式（8.2）和（8.4））并不是完全随机的；事实上，给定序列种子，序列的所有数字U都是完全可预测的。一些作者为强调这一点，将这种计算机产生的序列称为伪随机数。但如果适当选取a,c和m，序列U的随机性便足以通过一系列的统计检测。它们相对于真随机数具有可快速产生、需要时可再生的优点，尤其对于程序调试。

Monte Carlo 程序中通常需要产生服从给定概率分布F(x)的随机变量X。该步可用[6]，[13]－[15]中的几种方法加以实现，其中包括直接法和舍去法。

直接法（也称反演法或变换法），需要转换与随机变量X相关的累积概率函数

F(x)?prob(X?x)（即：F(x)为X?x的概率）。0?F(x)?1显然表明，通过产生(0,1)

内均匀分布随机数U，经转换我们可得服从F(x)分布的随机样本X。为了得到这样的具有概率分布F(x)的随机数X，不妨设U?F(x)，即可得

X?F(U) （8.6） 其中X具有分布函数F(x)。例如，若X是均值为?呈指数分布的随机变量，且 F(x)?1?e?x/??1 ,0?x?? （8.7）

在U?F(x)中解出X可得

1?U) （8.8） X???ln(

由于(1?U)本身就是区间(0,1)内的随机数，故可简写为

X???lnU （8.9）

有时(8.6)式所需的反函数F?1(x)不存在或很难获得。这种情况可用舍去法来处理。令

f(x)?dF(x)为随机变量X的概率密度函数。令a?x?b时的f(x)?0，且f(x)上界dx为M（即：f(x)?M），如图8.1所示。我们产生区间(0,1)内的两个随机数(U1,U2)，则 X1?a?(b?a)U1 f1?U2M （8.10） 分别为在(a,b)和(0,M)内均匀分布的随机数。若

f1?f(X1) （8.11） 则X1为X的可选值，否则被舍去，然后再试新的一组(U1,U2)。如此运用舍去法，所有位于f(x)以上的点都被舍去，而位于f(x)上或以下的点都由X1?a?(b?a)U1来产生X1。

图8.1 舍去法产生概率密度函数为

f(x)的随机变量

例8.1 设计一子程序使之产生0,1之间呈均匀分布的随机数U。用该程序产生随机变?，其概率分布由下式给定

T(?)?1(1?cos?),0???? 221解：生成U的子程序如图8.2所示。该子程序中，m?2?1?2147483647,c?0,且

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

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

0001 C********************************************************** 0002 C PROGRAM FOR GENERATING RANDOM VARIABLES 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 INTERVAL (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的子程序只是为了说明本章所介绍的一些概念。大多数计算机都有生成随机数的子程序。

为了生成随机变量?，令 U?T(?)??11(1?cos?) 2?1则有 ??T(U)?cos(1?2U)

据此，一系列具有给定分布的随机变量?便可由图8.2所示主程序中生成。

§8.3 误差计算

Monte Carlo 程序给出的解按大量的检测统计都达到了平均值。因此，该解中包含了平均值附近的浮动量，而且不可能达到100％的置信度。要计算Monte Carlo 算法的统计偏差，就必须采用与统计变量相关的各种统计方法。我们只简要介绍期望值和方差的概念，并利用中心极限定理来获得误差估计[13,16]。

设X是随机变量。则X的期望值或均值x定义为 x????? xf(x)dx （8.12）

这里f(x)是X的概率密度分布函数。如果从f(x)中取些独立的随机样本x1,x2,...,xN，那么的x估计值就表现为N个样本值的均值。

1?? xN?xn?1Nn （8.13）

?只是x的有着准确期望值的无偏估计。虽然x?的期望值x是X的真正的平均值，而x??x。因此，我们还需要x?的值在x附近的分布测度。 等于x，但x?在x附近的的值的分布，为了估计X以及x我们需要引入X的方差，其定义为X与x差的平方的期望值，即 Var(x)??2 ?(x?x)2??(x?x)2f(x)dx （8.14）

???由(x?x)2?x2?2xx?x2，故有

?2(x)??x2f(x)dx?2x?xf(x)dx?x2?f(x)dx （8.15）

?????????或者

?2(x)?x2?x2 （8.16）

?(x)?(x2?x2)1/2 （8.17）

方差的平方根称为标准差，即

?的标准差与x标准差给出了x在均值x附近的分布测度，并由此给出了误差幅度的阶数。x的标准差的关系表示为

?)? ?(x?(x)N （8.18）

?来估计x，那么结果中x?在x附近的扩该式表明，如果用根据(8.13)式由N个xn值构成的x散范围便与?(x)成比例，且随着样本数N的增加而降低。

?的扩散范围，我们定义样本方差为 为了估计x1N?)2 （8.19） S? (xn?x?N?1n?12由此式还可看出S的期望值等于?(x)。因此样本方差是?(x)的无偏估计。将（8.19）式中平方项乘出来，便可得样本标准差为

222?N? S????N?1?1/2?1N22??x?x?N?n??n?1?1/2 （8.20）

当N较大时，系数N/(N?1)可设为1。

作为获得中心极限定理的一种方法，概率论的一个基本解可考虑二次项函数 B(M)?N! pMqN?M （8.21）

M!(N?M)!该式表明N次独立随机试验中有M次成功的概率。(8.21)中，p是一次试验中成功的概率，且q?1?p。当M和N?M都较大时，便可用Stirling 公式 n!~nen?n 2?n （8.22）

因此(8.21)式可近似为正态分布[17]

?(x??x)2??)? B(M)?f(xexp??? (8.23) 2??)2??(x?2?(x)?1?)?其中x?Np且?(xNpq。因此当N??时，中心极限定理表明，描述由N点Monte

?的分布的概率密度函数是（8.23）式所示的正态分布函数f(x)。也就是Carlo算法获得的x说，大量随机变量的集合趋于呈正态分布。将（8.18）式代入（8.23）式可得

?N(x??x)2?N1?)? f(x exp??? （8.24）22??(x)2?(x)??正态（高斯）分布在工程，物理以及统计学的各类问题中都非常有用。高斯模型的显著特性

源于中心极限定理。因此，高斯模型经常用于其影响程度取决于由许多不规则的、浮动的元素构成的集合的情况。例8.2中我们给出了根据中心极限定理产生高斯随机变量的算法。

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

?落入该范围内。假设要得到x?位于x??和近估计出某一范围或区间以确保我们估计的xx??之间的概率。由定义

??x????? Prob?x???x令??x??x???)dx? （8.25） f(x??x??x可得

2/N??x?2??x???? Prob?x???x?????N/2??/???0e??d?

2 ?erf??N/2或

Prob?x?z?/2??? （8.26a） ???x?????N??x?z?/2?x?? ??1?? （8.26b）

N?其中erf(x)是误差函数，且z?/2是标准正态差的上?/2分位点。对（8.26）式可做如下解

X释：如果用来呈现独立随机观测值并构建相关随机区间x??的Monte Carlo 程序以较大的

?N????随机区间x??称N值反复运行，则这些随机区间的erf??2??x??分位点将近似等于x。

???N???为置信区间，erf??2??x??称为置信度。大多数的Monte Carlo 算法都使用误差

???是在实际均值x的标准差范围内的。由（8.26）式可得，样本均????x?/N，这表明x?位于区间x????x?/值x或三个标准差的值。如

N内的概率是0.6826或68.3％。若要求更高的置信度，可用两个

,M?1?0.6826??(x)?(x)????x?M?x Prob?x?M ???0.954, M?2 （8.27）NN???M?3?0.997,其中M是标准差的个数。

在（8.26）和（8.27）式中均假设总体标准差?为已知。由于这种情况很少出现，故必须用由（8.20）式算得的样本标准差S来估计?的值，从而使学生氏t－分布取代正态分布。我们知道当N很大，比如N?30时，t－分布近似趋于正态分布。（8.26）式等价于 Prob?x???St?/2;N?1N??x??xSt?/2;N?1? ??1?? （8.28）

N?其中t?/2;N?1为自由度是N?1的学生氏t－分布的上?/2分位点；其值在任何标准统计学课本中都有表列可查。这样置信区间的上、下限便可由下式给出 上限＝x?St?/2;N?1NSt?/2;N?1N （8.29）

下限＝x? （8.30）

Monte Carlo 算法中关于误差估计的进一步讨论，可参见[18,19]。

例8.2 利用中心极限定理产生一高斯（正态）分布的随机变量。根据中心极限定理，大量均值附近的独立随机变量的总和，无论其个体变量的分布如何，总趋向于一种高斯分布。也就是说，对于任意随机数Yi,i?1,2,....N，均值为Y，方差为Var(Y)，

Z??Yi?1Ni?NY （8.31）

NVar?Y?渐渐与N合并为零均值、单位标准差的正态分布。若Yi是均匀分布的随机变量（即Yi?Ui），

X则Y?1/2,Var(Y)?1/12，故而

Z?且变量

?Ui?1Ni?N/2 （8.32）

N/12 X??Z?? （8.33） 近似等于均值为?、方差为?的正态变量。N小于3时近似为我们所熟知的钟形高斯分布。为简化计算，通常实际中设N?12，因为这样可消除（8.32）式中的平方根项。然而N的这种取值截掉了?6?边界处的分布，且无法产生超过3?的值。对于分布曲线的两端比较重要的仿真，就必须用其它方案来产生高斯分布（参见[20]－[22]）。 这样，要产生一个均值为?、标准差为?的高斯随机变量 （1） 生成12个均匀分布的随机数U1,U2,....,U12 （2） 求得Z?，就要遵循以下步骤：

2?Ui?112i?6

（3） 令X??Z??

§8.4 数值积分

对于一维积分，现已有一些求积公式，如3.10节中所述。而对多重积分的公式则相对较少。Monte Carlo 技术对此类多重积分的重要性至少存在两方面的原因。多重积分的求积公式非常复杂，而多重积分的MCM几乎保持不变。Monte Carlo 积分的收敛性与维数无关，但对求积公式并非如此。人们已经发现，积分的统计方法是计算天线问题尤其是那些与较大结构相关的问题中的二维或三维积分的很有效的方法[23]。这里将讨论两类Monte Carlo 积分方法，即简易MCM 和具有对立变量的MCM。其它类型，如成功－失败法和控制变量法，可参见[24]－[26]。本节还将简要介绍MCM在不规范积分中的应用。

§8.4.1 简易Monte Carlo 积分

设要计算积分 I?其中R是n维空间。令X?X,X,....X随机变量，其均值由下式给出[27,28] f(X)?方差为

?Rf （8.34）

?12n?是在R内均匀分布的随机变量。则f(X)也是

I (8.35) R1R?Rf?X Var?f?X???1R??R?12f???R??R? （8.36） f???2其中 R?dX （8.37）

R如果取

的N个独立样本，即X1,X2,....,XN，它们与X具有相同的分布，且构成平均项

f?X1??f?X2??....?f?XN?1N ??f?Xi? （8.38）

NNi?1我们便会期望该平均项接近于f(X)的均值。这样，由（8.35）和（8.38）式，即可得

f?X? （8.39）?Nii?1 I?RNMonte Carlo 公式可用于有限区域R上的任何积分。为了举例说明，现将（8.39）式应用于

一维和二维积分中。 对于一维积分，设

I?由（8.39）式可得

?ba f(x)dx （8.40）

b?aN I? f?Xi? （8.41）?Ni?1其中Xi是区间?a,b?内随机变量，即

Xi?a??b?a?U,0?U?1 （8.42） 对于二维积分，有 I?则相应的Monte Carlo 公式为

???bdac fX1,X2dX1dX2 （8.43）

?(b?a)?d?c?N I? fXi1,Xi2 （8.44）?Ni?1??其中

Xi1?a?(b?a)U1,0?U1?1X?c?(d?c)U,0?U?12i22 （8.45）

（8.39）式中无偏估计值I收敛的很慢，这是由于估计值方差的级数是1/N。一种改

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

§8.4.2 具有对立变量的Monte Carlo 积分方法

术语对立变量[29,30]用来描述任一系列估计值，它们能互相抵消掉彼此的方差。方便起

X见，我们假设积分范围为区间(0,1)。设要得到下面单重积分的估计值 I?我们期望表达式

?g(U)dU （8.46）

011?g(U)?g(1?U)?与g(U)相比具有更小的方差。如果g(U)太小，那反2过来g(1?U)就很可能太大。因此，我们定义估计值

1 I?N ?2?g?U??g?1?U?? （8.47）

iii?1N1其中Ui是0和1之间的随机数。此估计值方差的级数是一个极大的进步。对于两维积分 I?则相应的估计值为

1，这是在（8.39）式基础上的4N???1100 gU1,U2dU1dU2 （8.48）

?1N1 I? gUi1,Ui2?gUi1,1?Ui2?g1?Ui1,Ui2?g1?Ui1,1?Ui2 （8.49）?Ni?14??????????根据相似性，可将该方法延拓至更高阶的积分。对于不是(0,1)的区间，可应用诸如（8.41）式到（8.45）式的转换。例如，

?f?x?dx??b?a??g?U?dU

a0b1b?aN1 ? ??g?Ui??g?1?Ui?? （8.50）

Ni?12其中g(U)?f(X)，且X?a??b?a?U。由（8.47）和（8.49）式可得，当积分维数增加时，每一维用来在简易Monte Carlo 方法上提高效率的对立变量的最小值也会增加。这样使得简易Monte Carlo 方法在多维运算中更具优越性。

§8.4.3 不规范积分 积分式

I???0 g(x)dx （8.51）

，其中f(x)在

可用Monte Carlo 仿真法进行估计[31]。对于概率密度为f(x)的随机变量区间(0,?)上的积分等于1，则有

??0?g?x?dx??g(x)dx （8.52）

0f?x?因此，为计算（8.51）式中的I值，首先要得到N个服从同一在区间(0,?)上的积分等于1的概率密度函数f(x)的独立随机变量。其样本均值

1Ng?xi? g?x?? （8.53） ?Ni?1f?xi?便是I的估计值。

例8.3 用Monte Carlo 方法计算积分 I???12?00ej??cos??d?d?

解：该积分式表示振幅和相位分别服从某一分布的圆形孔径天线的辐射状况。之所以选择该式，是因为它至少是辐射积分的一部分。且其解的闭合形式亦为可得，由此便可得Monte Carlo 解的精确性。其闭合解为 I????其中J1???是一阶Bessel 函数。

图8.3给出由（8.44）和（8.45）式计算该积分的一个简单的程序，其中a?0,b?1,c?012以及d?2?。该程序在Vax 11/780 中调用子程序RANDU来产生随机数U和U。对于不

2?J1????

同的N值，简易和对立变量Monte Carlo 方法都可用于计算辐射积分。表8.1中对??5时的结果进行了精确的比较。在应用（8.49）式时，用到了下面的对应式： U1?X1,U2?X2,1?U1?b?X1??b?a?1?U1 1?U?d?X22????d?c?1?U2

??

表8.1 例8.3辐射积分的MC方法积分结果

N 简易MCM 对立变量MCM

500 -0.2892-j0.0742 1000 -0.5737+j0.0808 2000 -0.4922-j0.0040 4000 -0.3999-j0.0345 6000 -0.3608-j0.0270 8000 -0.4327-j0.0378 10,000 -0.4229-j0.0237

-0.2887-j0.0585

-0.4982-j0.0080 -0.4682-j0.0082

-0.4216-j0.0323 -0.3787-j0.0440 -0.4139-j0.0241 -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解拉普拉斯方程

?V?0 （在区域R内） （8.54a） 满足Dirichlet边界条件

V?VP （在边界B上） （8.54b） 首先将R分成网状结构，并用其有限差当量替代?。（8.54a）在二维R内的有限差表达式可由（3.26）式给出，即

V?x,y??px?V?x??,y??px?V?x??,y??py?V?x,y????py?V?x,y??? （8.55a） 其中

px??px??py??py??221 （8.55b） 4（8.55）式中，假设网络的一个方格尺寸是?，如图8.4所示。下面给出该方程的概率解释。若随机游动粒子在某一瞬时处于?x,y?位置处，它从此点向(x??,y),(x??,y)(x,y??)及?x,y???移动的概率分别是px?,px?,py?和py?。决定粒子移动方式的方法是产生一个随机数U,0?U?1，并按下面的方式指导粒子的随机游动：

?x,y???x??,y?0?U?0.25?x,y???x??,y?0.25?U?0.5 （8.56）

?x,y???x,y???0.5?U?0.75?x,y???x,y???0.75?U?1如果不用方格而用矩形格，则有px??px?且py??py?，但px?py。在用立方体晶格表示的三维问题中，px??px??py??py??pz??pz??对区间0?U?1进行了分组。

1。两种情况中都依据概率6

图8.4 固定随机游动示意图

为了计算?x,y?处的位势，让一随机游动粒子在该点出发开始游动。粒子便开始从一点到另一点在网格中蜿蜒前行直到到达边界。此时，粒子终止游动，记下边界点的指定位势VP。

令第一个粒子游动结束时VP的值记为VP?1?，如图8.4所示。然后将第二个粒子从?x,y?释放，使其游动直到到达边界点，游动终止且该点VP的值相应的记为VP?2?。依次第三个，

??，第四个，第N个粒子从?x,y?释放重复此过程，并记下相应的指定位势VP?3?,VP?4?，

??VP?N?。根据Kakutani [32]，VP?1?,VP?2?，??VP?N?的期望值是Dirichlet问题

在?x,y?的解，即

1N V?x,y?? ?VP?i? （8.57）

Ni?1其中N较大，为游动总次数。其收敛速度随N而改变，所以需要保证许多随机游动都有精确的结果。

若要解Poisson方程

?2V??g?x,y? 在R内 （8.58a） 且V满足条件

V?VP 在B上 （8.58b） 则式（8.55）中的有限差分表示为

V?x,y??px?V?x??,y??px?V?x??,y?

?2g ?py?V?x,y????py?V?x,y???? （8.59）

4其中概率仍为式（8.55b）所示。（8.59）式的概率解释也近似于（8.55）式。但随机游动的每一步都要记录下（8.59）式中的?g/4项。如果第i次随机游动从?x,y?出发至到达边界

2需要mi步，则记录下

?2 VP?i??4由此可得，V?x,y?的Monte Carlo解为

?g?x,y? （8.60）

jjj?1mi?1N?2N?mi?1 V?x,y?? ??g?xj,yj?? （8.61）?VP?i??4N?Ni?1i?1?j?1? 对刚才所描述MCM的一个有趣的分析就是游走醉鬼问题[15,35]。我们将随机游动粒子当作一个“醉鬼”，将网格方块当作“城市街区”，结点当作“十字路口”，边界B当作“城市边界”，而且把边界B上的结点当作“警察”。尽管醉鬼想步行走回家，但他却醉的一塌

XXXX糊涂以至于在整个城市中随机游走。警察的工作是在城市边界上一看见醉鬼便将其抓获，并勒令醉鬼交付罚金VP。那么醉鬼要支付的预期罚金将会是多少呢？其答案便是（8.57）式。 在电介质边界上，指定边界条件为D1n?D2n。

?x,y?

?x?x?x?x?x?x?x?x? xxxxxxxxxxxxxxx

XXXXXXXXXXXXXXXX

[1].R.Hersch and R.J.Griego,“Brownian motion and potential theory,”Sci.Amer.,Mar.1969,pp.67-74.

本文来源：https://www.bwwdw.com/article/ocxr.html