蒙特卡洛模拟

《蒙特卡洛模拟》由会员分享，可在线阅读，更多相关《蒙特卡洛模拟（21页珍藏版）》请在装配图网上搜索。

1、第八章 Monte Carlo 法8.1 概述Monte Carlo 法不同于前面几章所介绍的确定性数值方法，它是用来解决数学和物理问 题的非确定性的（概率统计的或随机的）数值方法。Monte Carlo方法（MCM）,也称为统 计试验方法，是理论物理学两大主要学科的合并：即随机过程的概率统计理论（用于处理布 朗运动或随机游动实验）和位势理论，主要是研究均匀介质的稳定状态1。它是用一系列随 机数来近似解决问题的一种方法，是通过寻找一个概率统计的相似体并用实验取样过程来获 得该相似体的近似解的处理数学问题的一种手段。运用该近似方法所获得的问题的解 in spirit 更接近于物理实验结果，而不是

2、经典数值计算结果。普遍认为我们当前所应用的MC技术，其发展约可追溯至1944年，尽管在早些时候仍 有许多未解决的实例。MCM的发展归功于核武器早期工作期间Los Alamos （美国国家实验 室中子散射研究中心）的一批科学家。 Los Alamos 小组的基础工作刺激了一次巨大的学科 文化的迸发，并鼓励了 MCM在各种问题中的应用2-4“Monte Carlo”的名称取自于Monaco （摩纳哥）内以赌博娱乐而闻名的一座城市。Monte Carlo 方法的应用有两种途径：仿真和取样。仿真是指提供实际随机现象的数学 上的模仿的方法。一个典型的例子就是对中子进入反应堆屏障的运动进行仿真，用随机游动

3、 来模仿中子的锯齿形路径。取样是指通过研究少量的随机的子集来演绎大量元素的特性的方法。例如，f （x）在a x 0）00a=乘法器（a 0） c=增值（c 0） m =模数对于t数位的二进制整数，其模数通常为2t。例如，对于31位的计算机m即可取231-1。这 里x ,a和c都是整数，且具有相同的取值范围m a,m c,m x。所需的随机数序x 0 0 n 便可由下式得x 二(ax + c)modm(8.2)n+1n该序列称为线性同余序列 例如，若x = a = c = 7且m io，则该序列为07, 6, 9, 0, 7, 6, 9, 0 (8.3)可以证明，同余序列总会进入一个循环套；也就

4、是说，最终总会出现一个无休止重复的数字 的循环。(8.3)式中序列周期长度为 4。当然,一个有用的序列必是具有相对较长周期的序 列。许多作者都用术语乘同余法和混合同余法分别指代c二0和c工0时的线性同余法。选 取x0, a, c和m的法则可参见6,10。这里我们只关心在区间(0,1)内服从均匀分布的随机数的产生。用字符U来表示这些 数字,则由式(8.2)可得xU =(8.4)m这样U仅在数组0,1/m,2/m,(m-1)/m中取值。(对于区间(0,1)内的随机数，一 种快速检测其随机性的方法是看其均值是否为 0.5。其它检测方法可参见3,6。)产生区间 (a,b)内均匀分布的随机数X，可用下式

5、X = a + (b a)U(8.5)用计算机编码产生的随机数(利用式(8.2)和(8.4)并不是完全随机的；事实上, 给定序列种子，序列的所有数字U都是完全可预测的。一些作者为强调这一点，将这种计 算机产生的序列称为伪随机数但如果适当选取a,c和m，序列U的随机性便足以通过一 系列的统计检测。它们相对于真随机数具有可快速产生、需要时可再生的优点,尤其对于程 序调试。Monte Carlo程序中通常需要产生服从给定概率分布F(x)的随机变量X。该步可用 6， 1315中的几种方法加以实现，其中包括直接法和舍去法。直接法(也称反演法或变换法) ，需要转换与随机变量 X 相关的累积概率函数F(x)

6、二prob(X x)(即：F(x)为X x的概率)。0 F(x) x b时的f (x)二0，且f (x)上界dx为M (即：f (x) M)，如图8.1所示。我们产生区间(0,1)内的两个随机数(U,U2)，则X = a + (b a)Uf = U M(8 10)1 1 1 2分别为在(a,b)和(0,M)内均匀分布的随机数。若f1 f(X1)(8.11)则X1为X的可选值，否则被舍去，然后再试新的一组严2)。如此运用舍去法，所有位于f (x)以上的点都被舍去，而位于f (x)上或以下的点都由X = a + (b a)U来产生X1。图 8.1舍去法产生概率密度函数为 f ( x ) 的随机变量

7、例8.1设计一子程序使之产生0,1之间呈均匀分布的随机数U。用该程序产生随机变0 , 其概率分布由下式给定1T(6)二(1 cos6),0 6 兀解：生成U的子程序如图8.2所示。该子程序中，m = 221 1 = 2147483647,c = 0,且 a二75二16807。应用种子数(如1234)，主程序中每调用一次子程序，就会生成一个随机 数U。种子数可取1到m间的任一整数。000200030004000500060007000800090010001100120013001400150001000200030004x-、*1* *1* *1* *1* *1* *1* *1* *1* *1

8、* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* vt*、t*、t*、t*、t* *1* *1* *1* *1* *1* rT kt*rT* rT* rTxC PROGRAM FOR GENERATING RANDOM VARIABLES WITHC A GIVEN PROBABILITY DISTRIBUTIONx-、*1* *1

9、* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* vt*、t*、t*、t*、t* *1* *1* *1* *1* *1* rT kt*rT* rT* rTxDOUBLE PRECISION ISEEDISEED=1234.D010DO 10 I=1,100CALL RANSOM

10、(ISEED,R) THETA=ACOSD(1.0-2.0*R) WRITE(6,*)I,THETA CONTINUESTOPENDx-、*1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* vt*、t*、t*、t*、t* *1* *1* *1* *1* *1* rT

11、kt*rT* rT* rTxC SUBROUTINE FOR GENERATING RANDOM NUMBERS INC THE INTERVAL (0,1)x-、*1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* vt*、t*、t*、t*、t* *1* *1* *1

12、* *1* *1* rT kt*rT* rT* rTx00050006SUBROUTINE RANDOM (ISEED,R)0007DOUBLE PRECISION ISDDE,DEL,A0008DATA DEL,A/2147483647.D0,16807.D0/00090010ISDDE=DMOD(A*ISDDE,DEL)0011R=ISDDE/DEL 0012 RETURN 0013 END图 8.2 例 8.1 的随机数生成器图 8.2 的子程序只是为了说明本章所介绍的一些概念。大多数计算机都有生成随机数 的子程序。为了生成随机变量0，令1U = T (0) = (1 - cos 0)2

13、则有0 二 T-i (U)二 cos-i (1 U)据此，一系列具有给定分布的随机变量0便可由图8.所示主程序中生成。8.3 误差计算Monte Carlo 程序给出的解按大量的检测统计都达到了平均值。因此，该解中包含了平 均值附近的浮动量，而且不可能达到100的置信度。要计算 Monte Carlo 算法的统计偏差， 就必须采用与统计变量相关的各种统计方法。我们只简要介绍期望值和方差的概念，并利用 中心极限定理来获得误差估计13,16。设X是随机变量。则X的期望值或均值x定义为xf ( x)dxg这里f (x)是X的概率密度分布函数。如果从f (x)中取些独立的随机样本x , xx ,1 2

14、 N 那么的x估计值就表现为N个样本值的均值。入1 yx = x(8.13)Nnn=1x是X的真正的平均值，而x只是x的有着准确期望值的无偏估计。虽然x的期望值 等于x，但x丰x。因此，我们还需要x的值在x附近的分布测度。为了估计X以及x在x附近的的值的分布，我们需要引入X的方差，其定义为X与x 差的平方的期望值，即Var(x) =G2 = (x x)2 =(x x)2f (x)dx(8.14)g由(x x)2 = x2 2xx + 无2，故有8.15)g2(x) =x2f (x)dx 2xjg xf(x)dx + x2卜 f (x)dxggg或者g 2(x) = x2 x28.16)方差的平

15、方根称为标准差即 G (x) = (x2 x2)1/28.17)标准差给出了x在均值x附近的分布测度，并由此给出了误差幅度的阶数。x的标准差与x 的标准差的关系表示为g (x)8.18)该式表明，如果用根据(8.13)式由N个x值构成的x来估计x，那么结果中x在x附近的扩n散范围便与g (x)成比例，且随着样本数N的增加而降低。为了估计x的扩散范围，我们定义样本方差为S 2 = 艺(x x )2N 1 nn=18.19)由此式还可看出S 2的期望值等于g 2(x)。因此样本方差是G 2(x)的无偏估计。将(8.19)式中平方项乘出来，便可得样本标准差为S=(N1/21 yN入乙 x 2 x 2

16、 Nnn=1N 1J1/28.20)当N较大时，系数N /(N 一 1)可设为1。作为获得中心极限定理的一种方法，概率论的一个基本解可考虑二次项函数B (M) =pMqN-MM!(N -M)!该式表明N次独立随机试验中有M次成功的概率。（8.21）中，p是一次试验中成功的概率,且q - 1 - P。当M和N 一 M都较大时，便可用Stirling公式8.22)因此（8.21）式可近似为正态分布17(8.23)(X X )22g 2 (x)其中x二化且.（X）丽。因此当N *时，中心极限定理表明，描述由N点MonteCarlo算法获得的X的分布的概率密度函数是（8.23）式所示的正态分布函数f

17、（x）。也就是说，大量随机变量的集合趋于呈正态分布。将（8.18）式代入（8.23）式可得exp2G 2(X)8.24)正态（高斯）分布在工程，物理以及统计学的各类问题中都非常有用。高斯模型的显著特性 源于中心极限定理。因此，高斯模型经常用于其影响程度取决于由许多不规则的、浮动的元 素构成的集合的情况。例8.2中我们给出了根据中心极限定理产生高斯随机变量的算法。由于样本数N是有限的，所以Monte Carlo算法不可能达到绝对的确定性。故而在X附 近估计出某一范围或区间以确保我们估计的X落入该范围内。假设要得到X位于X8和 X + 8之间的概率。由定义8.25)Pr obX 8 X X +x-

18、8(X - X)8 /g)、:兀 08.26a)Pr ob X za/2- X X + z -Na/28.26b)其中erf （X）是误差函数，且Za/2是标准正态差的上a/2分位点。对（8.26）式可做如下解释：如果用来呈现独立随机观测值并构建相关随机区间x 5的Monte Carlo程序以较大的 N值反复运行，则这些随机区间的erf /xj分位点将近似等于X。随机区间X 5称 为置信区间，erf彳茨X）称为置信度。大多数的 Monte Carlo算法都使用误差 5=b（x）/N，这表明X是在实际均值x的标准差范围内的。由（8.26）式可得，样本均值X位于区间Xb（x）/jN内的概率是0.6

19、826或68.3%。若要求更高的置信度，可用两个或三个标准差的值。如Pr ob X 一 M 印 X X + M 1vN0.6826,二 30时，t 分布近似趋于正态分布。（8.26）式等价于f St_ StPr obi X /2N-1 X = 1 a（8.28）IvNI其中t 为自由度是N -1的学生氏t分布的上a/2分位点；其值在任何标准统计学/2; N1课本中都有表列可查。这样置信区间的上、下限便可由下式给出_ St上限=X + a d j（8.29）JNSt下限=x 1（8.30）PNMonte Carlo算法中关于误差估计的进一步讨论，可参见18,19。例8.2 利用中心极限定理产生一

20、高斯（正态）分布的随机变量 。根据中心极限定理， 大量均值附近的独立随机变量的总和，无论其个体变量的分布如何，总趋向于一种高斯分布。也就是说，对于任意随机数Y, i = 1,2,.N，均值为Y ,方差为Var（Y）,i Y NY8.31)ivWVar(Y )渐渐与N合并为零均值、单位标准差的正态分布。若Y是均匀分布的随机变量（即Y = U ）， iii则 Y 二 1/2,Var（Y）二（卫，故而兰 U -N/2iZ = 4（8.32）JN/12且变量X +卩（8.33）近似等于均值为卩、方差为b 2的正态变量。N小于3时近似为我们所熟知的钟形高斯分布。为简化计算，通常实际中设N二12，因为这样

21、可消除（8.32）式中的平方根项。然而N的这种取值截掉了 士 6b边界处的分布，且无法产生超过3b的值。对于分布曲线的两端比较重要的仿真，就必须用其它方案来产生高斯分布（参见2022）。这样，要产生一个均值为卩、标准差为b的高斯随机变量，就要遵循以下步骤：（1）生成12个均匀分布的随机数U ,U U1 2 12（2）求得Z =兰U - 6ii=1（3）令 X =bZ +卩8.4 数值积分对于一维积分，现已有一些求积公式，如3.10 节中所述。而对多重积分的公式则相对 较少。Monte Carlo技术对此类多重积分的重要性至少存在两方面的原因。多重积分的求积 公式非常复杂，而多重积分的MCM几乎

22、保持不变。Monte Carlo积分的收敛性与维数无关， 但对求积公式并非如此。人们已经发现，积分的统计方法是计算天线问题尤其是那些与较大 结构相关的问题中的二维或三维积分的很有效的方法23。这里将讨论两类Monte Carlo积 分方法，即简易MCM和具有对立变量的MCM。其它类型，如成功一失败法和控制变量法， 可参见24 26。本节还将简要介绍MCM在不规范积分中的应用。8.4.1 简易 Monte Carlo 积分设要计算积分I = J f（8.34）R其中R是n维空间。令X =（X1, X 2,.Xn）是在R内均匀分布的随机变量。则f （ X ）也是随机变量，其均值由下式给出27,28

23、IR(8.35)方差为其中如果取Var 0(X 讥肯 jRf 2 -声 “ 1|R 二 j dXR8.36)8.37)的N个独立样本，即X1, X2，., XN，它们与X具有相同的分布，且构成平均项f (X 1)+ f (X 2)+ . + f &丄1 艺 f(X)NNii=18.38)我们便会期望该平均项接近于f （X）的均值。这样，由（8.35 ）和（8.38）式，即可得8.39)Monte Carlo公式可用于有限区域R上的任何积分。为了举例说明，现将（8.39）式应用于 一维和二维积分中。对于一维积分，设I = jb f (x)dxa由( 8.39)式可得8.40)b-aii=18.4

24、1)其中X是区间（a,b）内随机变量，即i8.42)X = a +(b - ab,0 U 1i对于二维积分，有I = jb jd f (Xi, X2XidX28.43)则相应的 Monte Carlo 公式为aci=18.44)其中X i = a + （b a）U i,0 U i 1 i（8.45）X 2 = c + （d c）U 2,0 U 2 1 i（8.39）式中无偏估计值I收敛的很慢，这是由于估计值方差的级数是1/N。一种改 良的方法，即控制变量法，可通过减小估计值方差来提高其准确性和收敛性。 842具有对立变量的Monte Carlo积分方法术语对立变量 29,30用来描述任一系列估

25、计值，它们能互相抵消掉彼此的方差。方便起见，我们假设积分范围为区间(0,1)。设要得到下面单重积分的估计值I = f1 g (U )dU(8.46)0我们期望表达式2 g (U) + g (1 - U)与g (U)相比具有更小的方差。如果g (U)太小，那反过来g (1 - U)就很可能太大。因此，我们定义估计值1=；迟 2 也)+ giJ(847)i=11其中U .是0和1之间的随机数。此估计值方差的级数是_，这是在(8.39)式基础上的iN 4一个极大的进步。对于两维积分I = J*1 J*1 g U 1,U 2 u 1 dU 2(8.48)00则相应的估计值为8.49)I =丄迓-Ig

26、U 1,U 2)+ g U 1,1 - U 2)+ g ( U 1, U 2)+ g ( U 1,1 - U 2N 4i iiii iiii=1根据相似性，可将该方法延拓至更高阶的积分。对于不是(0,1)的区间，可应用诸如(8.41) 式到(8.45)式的转换。例如，Jbf (x )7x = (b - a )J1 g(U )dUa08.50)畔为 2 ig 0)+ g n- ui)i=1其中g(U) = f (X)，且X = a +(b a技。由(8.47)和(8.49)式可得，当积分维数增加时，每一维用来在简易Monte Carlo方法上提高效率的对立变量的最小值也会增加。这样 使得简易 M

27、onte Carlo 方法在多维运算中更具优越性。8.4.3不规范积分积分式I =代 g (x)dx(8.51)0可用Monte Carlo仿真法进行估计31。对于概率密度为f (x)的随机变量，其中f (x)在= J g ( x)dx08.52)区间(0,)上的积分等于1,则有J0因此，为计算(8.51)式中的I值，首先要得到N个服从同一在区间(0严)上的积分等于1 的概率密度函数 f(x) 的独立随机变量。其样本均值1 yN g (x ) 而=n yf(nii=1便是 I 的估计值。例8.3 用Monte Carlo方法计算积分00解：该积分式表示振幅和相位分别服从某一分布的圆形孔径天线的

28、辐射状况。之所以选择该 式，是因为它至少是辐射积分的一部分。且其解的闭合形式亦为可得，由此便可得 Monte Carlo 解的精确性。其闭合解为I （x）= SJ1G）a其中JC）是一阶Bessel函数。图8.3给出由（8.44）和（8.45）式计算该积分的一个简单的程序，其中a = 0,b = 1,c = 0以及d = 2兀。该程序在Vax 11/780中调用子程序RANDU来产生随机数U1和U2。对于不同的N值，简易和对立变量Monte Carlo方法都可用于计算辐射积分。表8.1中对a = 5时 的结果进行了精确的比较。在应用（8.49）式时，用到了下面的对应式：U1 三 X 1,U 2

29、 三 X 2,1 - U1 三 b - X1 =（ - a X - U1）(d-U2表8.1 例8.3辐射积分的MC方法积分结果N简易MCM对立变量MCM500-0.2892-j0.0742-0.2887-j0.05851000-0.5737+j0.0808-0.4982-j0.00802000-0.4922-j0.0040-0.4682-j0.00824000-0.3999-j0.0345-0.4216-j0.03236000-0.3608-j0.0270-0.3787-j0.04408000-0.4327-j0.0378-0.4139-j0.024110,000-0.4229-j0.023

30、7-0.4121-j0.0240精确解：-0.4116+j08.5 位势问题的解位势理论与布朗运动（或随机游动）的关系是由Kakutani在1944年首次提出的32。 自此，所谓的概率位势理论便在诸如热传导3338、静电学3946以及电力工程等许 多学科领域得到广泛应用。不同方程式的概率解或 MC 解的一个基本概念就是随机游动。 不同类型的随机游动应用不同的Monte Carlo方法。最常见的类型是固定随机游动和自动定00010002000300040005000600070008000900100011001200130014001500160017001800190020002100220

31、02300240025002600270028002900300031003200330034003500360037003800390040位随机游动。其它不常见类型有迁离法、缩减边界法、内切形法以及表面密度法。C INTEGRATION USING CRUDE MONTE CARLOC AND ANTITHETIC METHODSC ONLY FEW LINES NEED BE CHANGED TO USE THIS PROGRAMC FOR ANY MULTI DIMENSIONAL INTEGRATION*1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *

32、1* *1* vt* 、t* 、t* 、t* 1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* 、t* 、t* 、t* 、t* 、t* 、t* 、t* 、t* 、t* 、t* 、t* 、t* 、t* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* *1* ZT* rT*DATA IS1,IS2,IS3,IS4/1234,5678,9012,3456/DATA A,B,C/0.0,1.0,0.0/C SPECIFY THE INTEG

33、RAND F(RHO,PHI)=RHO*CEXP(J*ALPHA*RHO*COS(PHI)J=(0.0,1.0)ALPHA=5.0PIE=3.1415927 D=2.0*PIE DO 30 NRUN=500,10000,500SUM1=(0.0,0.0)SUM2=(0.0,0.0)DO 10 I=1,NRUNCALL RANDU(IS1,IS2,U1)CALL RANDU(IS3,IS4,U2) X1=A+(B-A)*U1X2=C+(D-C)*U2X3=(B-A)*(1.0-U1)X4=(D-C)*(1.0-U2) SUM1=SUM1+F(X1,X2)SUM2=SUM2+F(X1,X2)+F(

34、X1,X4)+F(X3,X2)+F(X3,X4)10 CONTINUE AREA1=(B-A)*(D-C)*SUM1/FLOAT(NRUN) AREA2=(B-A)*(D-C)*SUM2/(4.0*FLOAT(NRUN)PRINT *,NRUN,AREA1,AREA2 WRITE(6,*) NRUN,AREA1,AREA2WRITE(6,20) NRUN,AREA1,AREA220FORMAT(2X,NRUN=,I5,3X,AREA1=,F12.6,3X,F12.6,AREA2=,1 F12.6,3X,F12.6,/)30 CONTINUESTOPEND图 8.3 例 8.3 中用 Monte

35、 Carlo 方法计算二维积分的程序8.5.1 固定随机游动为具体起见，假设用固定随机游动的MCM解拉普拉斯方程V 2V二0(在区域R内)满足 Dirichlet 边界条件V = V(在边界B 上)P8.54b)首先将R分成网状结构，并用其有限差当量替代V 2。(8.54a)在二维R内的有限差表达式可由(3.26)式给出，即V(x,y)= p V(x + A,y)+ p V(x-A,y)+ p V(x,y + A)+ p V(x,y-A) (8.55a)y-x+其中1p 二p 二p 二p 二丁x+x-y +y -48.55b)(8.55)式中，假设网络的一个方格尺寸是A，如图8.4所示。下面给

36、出该方程的概率解释。 若随机游动粒子在某一瞬时处于(x, y)位置处，它从此点向(x + A, y), (x - A, y) (x, y + A)及(x, y -A)移动的概率分别是p , p , p和p 。决定粒子移动方式的方法是产生一个 x+ x-y+y -随机数U,0 U 1，并按下面的方式指导粒子的随机游动:(x, y)t (x + A, y)0 U 0.25(x,y)6-A,y)0.25 U 0.5(x, y)t (x, y + A)0.5 U 0.75(x,y)(x,y-A)0.75U 18.56)如果不用方格而用矩形格，则有px+二p且p 二p ，但p鼻p。在用立方体晶x-y +

37、y -xy1格表示的三维问题中，p二p二p二p二p二p二。两种情况中都依据概率x+x -y +y -z +z - 6对区间0 U 1进行了分组。图 8.4 固定随机游动示意图为了计算(x, y)处的位势，让一随机游动粒子在该点出发开始游动。粒子便开始从一点到另一点在网格中蜿蜒前行直到到达边界。此时，粒子终止游动，记下边界点的指定位势Vp。令第一个粒子游动结束时V的值记为V 6),如图8.4所示。然后将第二个粒子从(x, y)释 PP放，使其游动直到到达边界点，游动终止且该点V的值相应的记为V (2)。依次第三个，PP第四个，第N个粒子从(x, y)释放重复此过程，并记下相应的指定位势V (3)

38、,V (4)，PPV (N )。 根据 Kakutani 32, V (L), V (2)， V (N)的期望值是 Dirichlet 问题P P P P 在(x, y)的解，即V (x, y )=丄兰V (i)N pi=18.57)其中N较大，为游动总次数。其收敛速度随而改变，所以需要保证许多随机游动都有 精确的结果。若要解 Poisson 方程V 2V = g (x, y )在 R 内8.58a)且 V 满足条件8.58b)则式(8.55)中的有限差分表示为V(x, y)= p V(x + A, y)+ p V(x A, y)x+ p V (x ,y+y+A)+ p V(x, y A)+A

39、2gy 48.59)其中概率仍为式(8.55b)所示。(8.59)式的概率解释也近似于(8.55)式。但随机游动的 每一步都要记录下(8.59)式中的Ag/4项。如果第i次随机游动从(x,y)出发至到达边界 需要m步，则记录下i8.60)V 0+孚芦 g (x , y )P 4j jj=1由此可得，V(x, y )的Monte Carlo解为V , y )=N N c)+4N 迓艺g(x ,y )(8.61)j jj=1对刚才所描述MCM的一个有趣的分析就是游走醉鬼问题15,35。我们将随机游动粒子 当作一个“醉鬼”，将网格方块当作“城市街区”结点当作“十字路口”，边界B当作“城 市边界”，而且把边界 B 上的结点当作“警察”。尽管醉鬼想步行走回家，但他却醉的一塌 糊涂以至于在整个城市中随机游走。警察的工作是在城市边界上一看见醉鬼便将其抓获，并勒令醉鬼交付罚金Vp。那么醉鬼要支付的预期罚金将会是多少呢？其答案便是（8.57）式。 在电介质边界上，指定边界条件为D二D 。1n2 n(x,y )x x x x x x x x x x x x x xXXXXXXXXXXXXXXXX1.R.Hersch and R.J.Griego,“Brownian motion and potential theory,” Sci.Amer.,Mar.1969,pp.67-74.