基于Monte Carlo方法的泊松方程求解

Yan Jia-yi;Chen Yu-mei;Zheng Yu-xia

【摘 要】基于Monte Carlo方法的主要原理,求解泊松方程第一边值问题.通过构建随机游动模型,确定统计量,抽样产生随机样本,得到泊松方程解的估计值.并给出了详细的推导步骤和算法流程,证明了统计量均值是泊松方程第一边值问题的解,为该方法在复杂问题中提供一个简单的思路. 【期刊名称】《洛阳师范学院学报》 【年(卷),期】2019(038)005 【总页数】4页(P15-18)

【关键词】MonteCarlo方法;随机游动模型;泊松方程;第一边值问题 【作 者】Yan Jia-yi;Chen Yu-mei;Zheng Yu-xia 【作者单位】;; 【正文语种】中 文 【中图分类】O242.2 0 引言

Monte Carlo方法可以用于求解确定性问题和随机问题[1-2]. 偏微分方程在自然科学、技术科学和工程科学中的应用十分广泛，比如物理学、气象学、核技术等，但很难有解析的精确解，因此研究偏微分方程数值解十分必要. 常用的数值解法包括有限差分法、有限元法等[4]，但是直接用其求解高维方程，会出现运算量大、

计算时间长的问题. 冉营丽介绍了Monte Carlo方法的优势和解题步骤，将其应用于微分方程求解[3]； 左应红和王建国将Monte Carlo方法与有限差分法结合，建立随机游动模型，并应用于拉普拉斯方程和一维热传导方程[4]. Monte Carlo方法的稳定性和收敛速度与维度无关，适用于求解高维问题，其适用范围很广[5-6]. 本文将Monte Carlo方法与五点差分格式结合，建立随机游动模型求解泊松方程第一边值问题. 1 Monte Carlo基本原理

Monte Carlo方法的思想起源可以追溯到1777年的蒲丰投针实验[7]，随着世界上第一台电子计算机的问世，Monte Carlo方法得到了开拓性的发展. 乌拉姆(S.Ulam)、冯诺依曼(J.von Neumann)和梅特罗波利斯(N.Metropolis)三位科学家为此做出了突出的贡献，分别提出了逆变换算法、取舍算法和梅特罗波利斯算法，丰富了Monte Carlo方法的理论，并在计算机上实现了模拟过程[8-9]. 随着伪随机数的发展、随机数检验方法的创新、抽样方法的完善以及降低方差技巧的改进，Monte Carlo方法的发展十分迅速[10-11].

运用Monte Carlo方法求解问题时，可以简洁地归纳为四个步骤[1]：

1) 建立概率空间(Ω、F、P). Ω是样本空间，表示全体事件； F是事件空间； P是概率函数，可以用来描述随机问题和确定性问题. 对于确定性问题，需要根据其数学特征，构造虚拟概率空间，将其转变为随机问题，再应用Monte Carlo方法求解.

2) 确定随机变量X、概率分布f和统计量h. X是样本空间上的实值函数； h是X的函数，其期望为μ，均方差为σ.

3) 随机抽样. 从概率分布f中抽样，产生样本值X，并计算统计量h的取值h(Xi)，Xi和h(Xi)独立同分布.

4) 统计估计. 计算统计量h的算术平均值 如果统计量h的期望E(h(Xi))=μ，则其

均值是统计量的无偏估计值，作为解的估计值.

Monte Carlo方法的误差是随机性误差，需要用概率论的中心极限定理来分析. 当h(Xi)独立同分布，存在期望μ和均方差σ时，绝对误差值为

相对误差值为 其中α为显著水平，1-α为置信水平，Xα为正态差[1]. 2 泊松方程求解

泊松方程是一种特殊形式的椭圆型偏微分方程，在物理学、静电学、机械工程中普遍存在. 求解如下形式的泊松方程 (1)

其中u(x,y)为未知函数，D是R2中的一个有界区域，Г是的D边界，P=(x,y)为有界区域内的点，Q为边界上的点，S(P)和f(Q)为已知函数.

偏微分方程在实际应用中很难求解其精确值. 由于Monte Carlo方法与维数无关，高维问题不会产生本质上的困难，是一个较好的选择. 2.1 建立随机游动概率模型

运用Monte Carlo方法，建立随机游动模型求解泊松方程(1).

1) 网格剖分. 把区域D用步长为δ的网格进行覆盖，区域内部网格点集记作D′，区域边界网格点集记作Г′，因此Г′是D′的边界，P∈D′和Q∈Г′. 2) 离散方程. 运用五点差分格式将方程(1)离散化得

3) 构造一条随机游动路径. 假定随机游动点P是一个随机变量，随机游动点从初始点P=P0∈D′出发，由五点差分格式可知，P0以等概率向与之邻接的四个节点P11, P12, P13, P14,随机游动一步，到达P11； 然后，再以等概率向P11四周节点随机游动一步； 重复上述操作，直到到达边界点Q=Pm∈Г′，停止随机游动过

程. 此时随机游动路径为

{ρp:P→P1→…→Pm-1→Pm=Q∈Г′} 游动步数为m，统计量为

其中Pk∈D′, k=0,1…m-1,Pm=Q∈Г1. 如果随机游动初始点P∈Г′，则初始点P不游动，此时统计量为h(P)=f(Q), P=Q∈Г.

4) 估计统计量. 重复上一步n次得 其中hi表示第i条随机游动路径得到的统计量. 2.2 理论证明

定理1 统计量h(P)的期望E(h(P))是满足泊松方程(1)的解，即E(h(P))=u(P). 证明 假设泊松方程(1)有唯一解u(P). 需要分两种情形证明. 当P=P0∈D′时，E(h(P))=u(P)； 当P=P0∈Γ′时，E(h(P))=f(P).

情形一： 当P=P0∈D′时，假设沿某一条随机游动路径ρP的概率为Pr(P)，按照数学期望的定义有 从初始点P=P0出发，向其邻接的节点P11, P12, P13, P14游动一步，一共有四条可能游动路径{P=P0→P1i，而后沿路径ρP1i游动，i=1, 2, 3, 4}，且互不相容. 对于第一步随机游动{P=P0→P1i}，统计量h(P)的增加量为 其概率为 对于后续随机游动路径ρP1i，统计量h(P)的增加量为h(P1i)，其概率为Pr(P1i)，而 因而有

由方程的意义可知E(h(P))=u(P).

情形二： 当P∈Γ′时，从初始点P出发，只有一条随机游动路径，即停在初始点，其概率为Pr(P)=1，按照数学期望的定义有

h(P)=f(P), P∈Γ′. 2.3 算法流程

随机游动模拟的核心是以等概率向其四周的点游动，每一次游动有四种可能移动的方式，构成集合L(M)={(-δ, 0), (δ, 0), (0, -δ), (0, δ)}. 随机游动模型算法如下: 1) 产生随机数U，计算M=[4U]+1，其中U为服从(0, 1)上均匀分布的随机数； 2) 计算游动点坐标Pi+1=Pi+L(M)； 3) 到达边界停止随机游动；

4) 计算统计量 一次随机游动模拟结束； 5) 重复上述过程，做n次随机游动模拟； 6) 计算解的估计值 绝对误差值 相对误差值 3 数值实验 求解二维泊松方程

其第一边值条件为

该泊松方程的精确解为u(x, y)=exsin(πy)，用Monte Carlo方法随机游动模型求该问题的近似值，所有计算结果都在MATLAB R2014a的环境下实现. 对于第i次随机游动过程，统计量的取值为

exsin(πy)，

其中P(k)=(x(k), y(k))∈D′, k=1, 2, …, m-1，Q=(x, y)∈Γ′. 设置信水平1-α=0.99，查表得正态差X0.99=2.5758. 该泊松方程的解函数在点(1, 处的精确值为1.922116，当步长时，不同模拟次数的计算结果见表1. 当步长δ分别取 模拟次数为100000，该泊松方程的解函数在点的计算结果见表2～表4.

表1 Monte Carlo方法在点处的计算结果模拟次数求解值绝对误差相对误差10-1.5155814.320894-2.8509811002.0362581.5760540.77399510002.0678550.4547370.219908100001.9393470.1405480.0724721000001.9443420.0451590.02322610000001.9410450.0142900.007362

表2 Monte Carlo方法在不同结点处的计算结果(δ, δ)/(x, y)(1/2, 1/4)(1, 1/4)(3/2, 1/4)(1/2, 1/2)(1, 1/2)(1/8,

1/8)1.1718391.9443423.1668001.6680992.748583(1/16, 1/16)1.1651101.9330573.1681311.6400922.712226(1/32, 1/32)1.1484001.9131953.1463911.6452682.708981精确解1.1658221.92211553.1690331.6487212.718282

表3 Monte Carlo方法在不同结点处的绝对误差值结果(δ, δ)/(x, y)(1/2, 1/4)(1, 1/4)(3/2, 1/4)(1/2, 1/2)(1, 1/2)(1/8,

1/8)0.0358140.0451590.0494640.0305840.038630(1/16, 1/16)0.0357920.0453950.0499010.0310620.039189(1/32, 1/32)0.0360220.0455380.0498120.0309450.039178

表4 Monte Carlo方法在不同结点处的相对误差值结果(δ, δ)/(x, y)(1/2, 1/4)(1, 1/4)(3/2, 1/4)(1/2, 1/2)(1, 1/2)(1/8, 1/8)0.0305620.0232260.015620 0.0183350.014055(1/16,

1/16)0.0307200.0234840.0157510.0189390.014449(1/32, 1/32)0.0313660.0238020.0158320.0188090.014462 表2到表4表明数值结果与理论结果一致. 4 总结

本文先简单介绍Monte Calro方法的基本原理以及误差，然后基于此，建立求解

二维泊松方程第一边值问题的随机游动模型： (1)将区域进行网格剖分； (2)用五点差分格式离散泊松方程； (3)构造随机游动模型，并确定统计量； (4)估计统计量, 并证明统计量的均值为泊松方程的解，然后给出相应的算法步骤，最后求解算例. 参考文献

【相关文献】

[1] 康崇禄. 蒙特卡罗方法理论和应用[M]. 北京:科学出版社, 2015. [2] 刘军. 科学计算中的蒙特卡罗策略[M]. 北京:高等教育出版社, 2009.

[3] 冉营丽. 探讨蒙特卡罗方法在解微分方程边值问题中的应用[J]. 梧州学院学报, 2015(3):42-46. [4] 左应红, 王建国. 蒙特卡罗方法在解微分方程边值问题中的应用[J]. 强激光与粒子束, 2012, 24(12):3023-3027.

[5] 刘芳芳, 刘播, 刘春光. 一种求解抛物型方程的Monte Carlo并行算法[J]. 高等学校计算数学学报, 2005(S1):208-212.

[6] 游皎, 李万爱. 高效蒙特卡罗方法在偏微分方程初边值问题中的应用[J]. 中山大学学报:自然科学版, 2015, 54(6):37-40.

[7] Buffon C D. Essai d’Arithmétique morale[J]. Supplément à l’ Historie Naturelle. Paris, Impimerie Royale, 1777, 4:321.

[8] Ulam S, Richtmyer R D, von Neumann J. Statistical methods in neutron diffusion[J]. Los Alamos Scientific Laboratory report LAMS-551, 1947.

[9] Metropolis N, Ulam S. The Monte Carlo method[J]. Journal of the American Statistical Association, 1949, 44:335-341.

[10] Metropolis N. The beginning of the Monte Carlo method[J]. Los Almos Science, 1987, 12:125-130.

[11] Eckhardt R. Stan Ulam, John von Neumann, and the Monte Carlo method[J]. Los Almos Science Special Issue, 1983(15):131-137.