高效蒙特卡罗方法在偏微分方程初边值问题中的应用*

游 皎，李万爱

(中山大学中法核工程与技术学院，广东 珠海 519082)



高效蒙特卡罗方法在偏微分方程初边值问题中的应用*

游 皎，李万爱

(中山大学中法核工程与技术学院，广东 珠海 519082)

介绍了蒙特卡罗方法求解偏微分方程初边值问题的基本原理，并针对该方法在多节点、多游动次数计算中速度慢的问题，提出了一种改进方案。通过理论分析和算例验证，在相同计算条件下与传统方法相比，该改进方案不仅大大减少了计算时间，而且降低了误差，这将使蒙特卡罗方法得到更广泛的应用。

初边值问题；蒙特卡罗方法；随机游动

蒙特卡罗方法是一种利用计算机进行统计模拟的数值计算方法，已经被应用到包括物理、生物、化学和金融等众多领域[1-2]。而在科学计算中常常需要求解偏微分方程初边值问题，而该类问题传统上均使用确定性方法如有限元方法。蒙特卡罗与有限元方法的区别在于：①后者需要使用方程的弱解形式进行离散，而前者可以使用各种有效的游走规则；②后者求解大型线性方程组获得求解点值，而前者利用概率统计理论建立求解值与已知值的关系；③后者在非确定性问题上的求解需要建立相关随机模型，而前者自然地具有概率性求解的特点。故一些文献也开始探索如何有效地利用蒙特卡罗方法来求解微分方程[3-8]，为了使结果更加精确，利用蒙特卡罗法时往往需要对求解域进行细密的网格划分。而节点数增大，游走次数增多后，计算量变得很大，需要较长计算时间，这严重限制了蒙特卡罗方法的实际应用。文献[4]提出了一种改进方法，该方法已经比已有方法(文献[5-6])在计算效率上有较大改进。本文继续对文献[4]方法加以改进，方程边界条件处理采用文献[7]方法，并通过理论分析和算例结果说明，本文的改进方法在相同的条件下，不仅可以大大减少计算时间，而且降低了误差。这些改进将使蒙特卡罗方法能应用到更为占用计算资源的初边值偏微分方程计算中。

1 基本原理及改进方案

1.1 传统方案

用蒙特卡罗方法求解微分方程，首先要构造方程对应的随机概率模型。为具体起见，本文首先以二维拉普拉斯方程为例，对利用蒙特卡罗方法解决边值问题进行描述。设问题的控制方程为

(1)

其中u为要求解物理量，Ω为求解域。它满足Dirichlet边界条件

其中∂Ω为Ω的边界，up(x,y)为边界值函数。首先对求解域Ω进行网格划分，然后相应的离散格式作为蒙特卡罗方法的游动规则。这里以正规结构化网格上的差分格式为例，采用五点差分格式离散上面的偏微分方程，可得

(2)

其中(i,j)为网格点编号，在均匀网格中有

可以将这种思路推广到微分方程初边值问题中，以扩散方程为例，

(3)

其中方程的边值和初值如下

对该方程采用时间后差、空间中心差的差分格式(BTCS)并整理可得

(4)

其中

p0=λ,p1=p0σx,p2=p0σy

1.2 改进方案

上面使用蒙特卡罗方法求解偏微分方程时需要使用差分格式来对方程进行离散。在工程计算里常常需要采用较密网格来得到可靠的结果，同时由于随机游动过程具有随机性，通常m要取得很大时才能得到较为准确的结果。因此这种方法的计算量比较大，需要一些提高计算速度的方法来加速这个计算过程[1,4]。而文献[4]的方法在加快计算速度及提高精度上都更有优势，它的改进步骤为下面的随机游动3：

在这个误差减少的分析基础上，本文提出一种先计算“大网格”的方法，通过改变网格点的计算顺序以达到更好的优化。具体来说，若求解域划分为50×50的网格，可以首先计算i=10, 20, 30, 40和j=10, 20, 30, 40这4行4列的函数值，相当于把求解域分成了25个小求解域，其边界由外边界和已知近似值的点构成，这样大大缩短了其它点的马尔可夫链长度，因而缩短了计算时间。此外，在计算“大网格”时可以适当增大游动次数以获得更准确的值。本文中采用的改进版均是首先计算了间隔为5的网格点的值，再进行其他点的计算。下面将这种方法称为reorder方法。

1.3 误差分析

随机游动4其实是随机游动3在初边值问题的推广。虽然文献[4]已经验证过随机游动3的精度比随机游动1更为提高，但是没有进行数学原理上的分析证明。这里将给出随机游动3比1、随机游动4比2能够计算更准确的分析。

这里有两个误差的问题。首先是游动规则带来的误差，本文的游动规则是通过有限差分方法给出的，其空间误差为截断误差

时间误差为O(Δt)。

第二个误差为游动过程的概率统计误差。首先需要研究对于传统方法，游动次数与误差之间的关系。利用文献[6]的类似分析，以拉普拉斯方程为例，对于在求解域Ω内任意一个网格点A, 由公式可以推得满足方程的近似解为

(5)

式中s为边界点个数，up(Qj)为边界点Qj对应的函数值up(x,y)；P(A,Qj)表示粒子从A点出发后游动到边界Qj的概率。另外，定义随机变量

ξi,j的期望E(ξi,j)为

方差V(ξi,j)为

用nj表示m个游动粒子中到达边界点Qj的粒子数目，有

由契比雪夫不等式，结合式和式，对任意一给定的正数ε1>0，有

结合式和式，对任意小的正数ε2，有

由于ξi,j(i=1,2,…,m)具有相同的概率分布，并且相互独立，因此有

所以

这个结论说明一点的游动次数m越大，则该点误差越小。在离散公式(2)在均匀网格为例

一点的函数值取决于相邻四点的值。如果以该点(i,j)为起点需要游动m次，则对于传统方案，可以假定每个相邻点的值都经由m/4次游动后得到。若有一个邻点的值是已知的，即它已经经历过至少m次随机游动，直接使用邻点的已知函数值相当于过该邻点的游动次数至少为m2/4(乘法原理)。平均来说有两个邻点函数值已知则总次数为m/2+m2/2，这将大大增加随机游动次数从而减少误差。

2 数值例子

2.1 已知精确解的二维泊松方程

考虑一个温度分布T(x,y)在域Ω=[0,1]×[0,1]满足方程如下的泊松方程，

且边界条件为T(x,y)=1+x。则此方程有解析解

分别采用传统的(随机游动1，对应表 1方案0)和改进后(随机游动3+reorder方法，对应表 1方案1)的方案对上述问题进行求解，并改变参数条件，进行横向和纵向的对比。由于蒙特卡罗方法具有随机性，表 1中每个方案都进行了3次实验，结果取平均值。在表 1中，N为网格节点数，m为每个节点的游动次数，时间单位为秒，时间比为方案1的时间与方案0的时间之比，误差采用L1范数，相对误差为误差与精确解之比。可以看出，在相同的条件下，改进版不仅在时间上较传统版有巨大优势，而且在降低误差方面也表现得更好。

表1 泊松方程的蒙特卡罗求解方法的误差对比

2.2 二维非定常热传导问题

考虑如下二维非定常热传导方程

其边界条件及初始条件分别为

u(t,0,y)=u(t,1,y)=u(t,x,0)=u(t,x,1)=0,

u(x,y,0)=exp(-20(x-1/2)2-20(y-1/2)2)

由于非定常问题传统方法(即随机游动2)在二维中计算量过大，这里的传统方法采用随机游动2加上随机游动4中的时间改进方案。使用传统方法(对应表 2方案0)和改进方案(随机游动4+reorder方法，对应表 2方案1)求解t=0.05时的分布。用二阶Crank-Nicolson差分法取精细条件(N=1 000×1 000，dt=0.001)作为参考解，以此计算两种方案的误差。其计算结果及时间、误差的比较如表 2所示。可以看出改进方案不仅大程度地压缩了运行时间，而且在一定程度上减少了误差。

为了进一步分析计算结果，本文还作了温度等值线进行对比，如图 1所示。改进方案的等高线更加圆滑，说明该方案的随机误差更小，结果更稳定。

表2 二维非定常热传导问题的蒙特卡罗法误差比较(N=512,m=1 000)

图1 常热传导方程蒙特卡罗法计算的温度等值线对比Fig.1 Thermal isolines comparison using MC methods on unsteady heat conduction equation

3 总 结

本文研究了蒙特卡罗方法在微分方程初边值问题中的应用，针对传统方法计算速度较慢的问题，提出了一种适合多节点、多游动次数的改进方法并给出精度提高的分析证明过程。与较传统方法相比，改进方案无论在计算时间和结果精度上都有明显优势，适合蒙特卡罗方法在更加复杂工程问题中的应用。下一步工作将结合曲面处理技术[9],将该方法推广到一般的几何计算区域中。

[1] 刘军. 科学计算中的蒙特卡罗策略[M]. 北京：高等教育出版社，2009.

[2] 洪志敏. 基于Monte-Carlo技术的积分(微分)方程数值求解方法研究[D]. 内蒙古：内蒙古工业大学博士学位论文.

[3] 左应红，王建国. 蒙特卡罗方法在解微分方程边值问题中的应用[J]. 强激光与粒子束，2012, 24(12): 3023-3026.

[4] 周泉，陈植武，詹杰民. 蒙特卡罗法的改进：一种新型的快速算法[C]∥第二十一届全国水动力学研讨会暨第八届全国水动力学学术会议暨两岸船舶与海洋工程水动力学研讨会文集, 2008.

[5] 林建国，周俊陶. 蒙特卡罗法的一种快速计算：在势流问题中的应用[J]. 水动力学研究与进展A辑，2005, 20(3): 405-410.

[6] 史慧会，曲小钢. 蒙特卡罗方法在热传导方程中的应用[J]. 重庆工学院学报：自然科学版，2008, 22(11): 101-103.

[7]VAJARGAHBF,VAJARGAHKF.MonteCarlomethodforfindingthesolutionofDirichletpartialdifferentialequations[J].AppliedMathematicalSciences, 2007, 1(10): 453-462.

[8]LAUDUDP,BINDERK.AguildtoMonteCarlosimulationsinstatisticalphysics[M].CambridgeUniversityPress, 2009.

[9]JUNS,SUNGHWANY,NAMZC.AMonteCarlomethodforsolvingheatconductionproblemswithcomplicatedgeometry[J].NuclearEngineeringandTechnology, 2007, 29(3): 207-214.

The Efficient Monte-Carlo Method in Solving the Initial-Boundary Value Problem of Partial Differential Equations

YOUJiao,LIWanai

(Sino-French Institute of Nuclear Engineering & Technology, Sun Yat-sen university, Zhuhai 519082, China)

The primary principle of Monte Carlo method in solving the initial-boundary value problem for partial differential equations is introduced. And an improved strategy to solve the slow speed problem in computations with large grid points and random walks is proposed. Through the theoretical analysis and numerical validation, the improved strategy not only takes much less computational time, but also gives more accurate results comparing to the traditional methods under the same conditions.

initial-boundary value problem; Monte Carlo method; random walk

10.13471/j.cnki.acta.snus.2015.06.007

2015-01-17 基金项目：国家自然科学青年基金资助项目(11402313)；超级计算应用领域资助项目(141gjc10)

游皎(1992年生)，男；研究方向：核工程与技术；通讯作者：李万爱；E-mail:liwai@mail.sysu.edu.cn

O21

A

0529-6579(2015)06-0037-05