Advances in Applied Mathematics
Vol. 08  No. 02 ( 2019 ), Article ID: 28760 , 7 pages
10.12677/AAM.2019.82020

The Method of Monte Carlo Markov Chain for Solving the Poisson Equation on Irregular Domain

Chengfeng Zheng1,2, Zhizhong Yan1,2*

1School of Mathematics and Statistics, Beijing Institute of Technology, Beijing

2Beijing Key Laboratory on MCAACI, Beijing Institute of Technology, Beijing

Received: Jan. 10th, 2019; accepted: Jan. 25th, 2019; published: Feb. 1st, 2019

ABSTRACT

In this paper, the Monte Carlo Markov chain method for solving Poisson equation in irregular domain is studied. In irregular domain, the numerical calculation of differential equations is usually difficult. Based on the idea of finite difference, an absorbable Markov chain is constructed to solve differential equations in irregular domain. The numerical experiments in irregular domain show that the method is feasible and effective. This method provides a new idea for solving Poisson’s equation in irregular region and keeps the second order convergence order.

Keywords:Monte Carlo Chain Method, Poisson Equation, Irregular Domain

不规则区域上Poisson方程 的蒙特卡洛马尔科夫链方法求解

郑成丰1,2,闫志忠1,2*

1北京理工大学数学与统计学院,北京

2北京理工大学复杂信息数学表征分析与应用北京市重点实验室,北京

收稿日期:2019年1月10日;录用日期:2019年1月25日;发布日期:2019年2月1日

摘 要

本文研究了不规则区域上泊松方程的蒙特卡洛马尔科夫链方法数值求解。在不规则区域上,微分方程的数值计算通常是困难的。基于有限差分的思想,本文构造了可吸收马尔科夫链来求解不规则区域上的微分方程。不规则区域的数值实验结果表明了该方法的可行性和有效性。该方法为求解不规则区域上的泊松方程提供了一种新的计算思路,并保持了空间方向上二阶收敛阶。

关键词 :蒙特卡洛马尔科夫链方法,Poisson方程,不规则区域

Copyright © 2019 by authors and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY).

http://creativecommons.org/licenses/by/4.0/

1. 引言

众所周知,Poisson方程是一种经典的偏微分方程,广泛地应用于静电学,机械工程和理论物理等领域 [1] [2] [3] [4] [5]。由于方程求解域逐渐复杂化,Poisson方程难以得到解析解。传统上人们采用有限差分方法和有限元方法来得到Poisson方程的数值解,但是这两种各有优缺点。有限差分法是一种早期成熟的经典数值方法,在复杂边界的情况下,该方法并不灵活,而且对网格生成的要求非常严格 [6]。而有限元法是一种常用的数值方法,对不同求解域的微分方程具有较强的适应性,但是该方法的计算过程抽象而且复杂 [7]。近年来,蒙特卡洛马尔科夫链方法也越来越引起人们的关注 [8] [9] [10]。蒙特卡洛马尔科夫链方法在求解微分方程时具有良好的边界适应性,并且易于编程实现,是一种非常具有发展前景的数值方法。鉴于以上分析,本文考虑不规则区域上Poisson方程的蒙特卡洛马尔科夫链方法。

首先通过网格剖分构建不规则区域Poisson方程的有限差分方法,进而建立蒙特卡洛随机模型获得Poisson方程的近似解,根据蒙特卡洛随机模型引入马尔科夫链从而获得Poisson方程的数值解。

2. 蒙特卡洛随机模型

本文考虑迪利克雷条件下的Poisson方程

{ Δ u = f in Ω u = g on Γ (1)

这里u 表示的是待求解的实函数,x,y为实自变量,并且 Γ 为求解域 Ω 的边界。由于本文考虑的是不规则区域上的Poisson方程,不失一般性,求解域不妨设为图1所示的情况。

Figure 1. The irregular solution domain of the Poisson equation

图1. Poisson方程的不规则的求解域

图1中, Ω 表示整个求解域,并且 Γ 1 , Γ 2 Γ 3 表示求解域的三个边界。下面对整个不规则求解域进行均匀网格剖分来获取网格节点。然而由于不规则边界的存在,网格剖分获得的边界点无法均匀的分布在边界上,如图2所示。

Figure 2. Grid points obtained by mesh generation for irregular domain

图2. 不规则区域剖分获得的网格节点

图2中,蓝色点(见联机版)表示内部网格节点 Ω k 并且k表示内部网格节点的序列编号。红色点(见联机版)表示由网格剖分获得的边界点 Γ b ,这里b表示边界点的序列编号。对于x方向剖分的次数为 M 1 ,该方向的剖分步长为 Δ x i ( i = 1 , 2 , , M ) 。并且y方向的剖分次数为 N 1 ,其对应的剖分步长为 Δ y j ( j = 1 , 2 , , N ) 。由于生成网格节点的不均匀性,方程(1)不能直接通过二阶中心差分进行离散,下面本文给出了非均匀离散形式的二阶中心差分。

2 Δ x i 1 ( U i + 1 , j U i , j ) + Δ x i ( U i 1 , j U i , j ) Δ x i 1 Δ x i ( Δ x i + Δ x i 1 ) + 2 Δ y j 1 ( U i , j + 1 U i , j ) + Δ y j ( U i , j 1 U i , j ) Δ y j 1 Δ y j ( Δ y j + Δ y j 1 ) = f i , j (2)

由方程(2)可以得到:

1 Δ x i ( Δ x i + Δ x i 1 ) U i + 1 , j ( 1 Δ x i ( Δ x i + Δ x i 1 ) + 1 Δ x i 1 ( Δ x i + Δ x i 1 ) ) U i , j + 1 Δ x i 1 ( Δ x i + Δ x i 1 ) U i 1 , j + 1 Δ y j ( Δ y j + Δ y j 1 ) U i , j + 1 ( 1 Δ y j ( Δ y j + Δ y j 1 ) + 1 Δ y j 1 ( Δ y j + Δ y j 1 ) ) U i , j + 1 Δ y j 1 ( Δ y j + Δ y j 1 ) U i , j 1 = 1 2 f i , j (3)

λ 1 = 1 Δ x i ( Δ x i + Δ x i 1 ) λ 2 = 1 Δ x i 1 ( Δ x i + Δ x i 1 ) λ 3 = 1 Δ y j ( Δ y j + Δ y j 1 ) λ 4 = 1 Δ y j 1 ( Δ y j + Δ y j 1 ) λ = λ 1 + λ 2 + λ 3 + λ 4 ,可以得到:

U i , j = λ 1 λ U i + 1 , j + λ 2 λ U i 1 , j + λ 3 λ U i , j + 1 + λ 4 λ U i , j 1 + 1 2 λ f i , j (4)

由方程(4)可知,求解域中的任意内点 Ω k 的函数值与其四个相邻点的函数值有关。由此可以对其他内点可以建立相同的近似方程,通过简化,可以得到所有内点与边界点的关系,从而可以得到任意内点 Ω k 的函数值 U ( Ω k ) 。如果有 N p 个粒子在 Ω k 点处释放,则这些粒子以方程(4)中对应的概率在网格上随机游

动。网格节点 Ω k 上粒子的移动方向可根据方程(4)中以 [ λ 1 k λ k , λ 2 k λ k , λ 3 k λ k , λ 4 k λ k ] 为概率生成的随机数R确定,R在 [ 1 , 2 , 3 , 4 ] 取值。其中0表示向左移动,1表示向右移动,2表示向上移动,3表示向下移动,粒子随机游动直至被边界点吸收,这样便构成了蒙特卡洛随机模型。那么任意内点 Ω k 的函数值可以近似的表示为:

U ( Ω k ) = b g ( Γ b ) P b N p + 1 N p i = 1 N p m = 1 n i 1 1 2 λ m f ( Ω m ) (5)

这里 P b 表示被边界点 Γ b 吸收的粒子的个数, n i 表示第i个粒子到达边界所运动的步数,即第i个粒子经过了 n i 1 个内点在第 n i 次运动后被边界点吸收。 λ m 表示在内点 Ω m 处对应的 λ 的值。

3. Poisson方程的蒙特卡洛马尔科夫链方法

为了详细介绍马尔科夫链,下面将对求解域中的内点和边界点进行编号,如图3所示。

Figure 3. Serial numbers of interior points and boundary points in irregular domain

图3. 不规则区域中的内点和边界点的编号

图3中,内点的个数为 N f ,边界点的个数为 N b 。所有的内点和边界点按照图3的顺序排列成一列,然后每个粒子按照相应的概率移动直到被吸收,显然这是一个可吸收的马尔可夫链,其中马尔可夫链具有 N b 个可吸收态和 N f 个不可吸收态。在方程(4)中,在 U i , j 处释放的粒子到 [ U i + 1 , j , U i 1 , j , U i , j + 1 , U i , j 1 ] 的概率为它们各自的系数。假设 [ U i + 1 , j , U i 1 , j , U i , j + 1 , U i , j 1 ] 的编号为 [ m , n , p , q ] 。对于编号为k的内点 U i , j 概率分布可以表示成下面的形式:

P k = [ 0 , , λ 1 k λ k , 0 , , λ 2 k λ k , 0 , λ 3 k λ k , 0 , λ 4 k λ k , 0 , , 0 ] 1 × ( N f + N b ) (6)

这里 [ λ 1 k λ k , λ 2 k λ k , λ 3 k λ k , λ 4 k λ k ] 分别处在 P k [ m , n , p , q ] 的位置。由于每个网格点都可以确定其概率分布,那么整个马尔可夫链的概率传递矩阵可以得到,记为:

P = [ I 0 R Q ] ( N f + N b ) × ( N f + N b ) (7)

其中, P i , j 表示粒子从状态i到状态j的概率,矩阵 R ( N f × N b ) 表示粒子从不可吸收状态到可吸收状态的概率,矩阵 Q ( N f × N f ) 表示粒子从不可吸收状态到不可吸收状态的概率,矩阵 I ( N b × N b ) 表示粒子从可吸收状态到可吸收状态的概率,并且它是一个单位矩阵。矩阵0的元素都是0,表示粒子从可吸收状态到不可吸收状态的概率。对于可吸收马尔可夫链,矩阵 E Q 具有可逆矩阵,其中E是与Q同阶的单位矩阵。

N = ( E Q ) 1

N i , j 具有一定的实际意义,表示从状态i到状态j的传输路径数。那么对应的可吸收概率矩阵为:

B = N R

那么对于所有内点的函数值可以表示为下面的形式:

U Ω = B g Γ + N F Ω (8)

这里 U Ω 表示所有内点的函数值, g Γ 表示所有边界点的函数值, F Ω = 1 2 λ Ω f Ω ,并且 f Ω 表示所有内点对应的f的函数值, λ Ω 为相应的内点所对应的 λ 的值。这里值得注意的是 U Ω 所对应的内点的函数值的编号与内点的编号是一致的。

4. 数值实验

为了进一步证明本文提出方法的有效性,本文将使用下面的形式来计算数值收敛阶。

E N ( h ) = | U ( x c , y j ) u ( x c , y j ) | , r a t e h = log 2 ( E N ( h ) E N ( 2 h ) )

这里h表示的是在x和y方向上的剖分次数 M = N 时, h = Δ x = Δ y E N ( h ) 表示数值解和精确解之间的最大范数误差, r a t e h 表示空间收敛阶,并且 c = r o u n d ( M / 2 )

考虑二维不规则区域上Poisson方程的迪利克雷问题,可以描述为下面的形式:

{ 2 u x 2 + 2 u y 2 = 2 e x + y , ( x , y ) Ω u = e x + y , ( x , y ) Γ

其中 Ω Γ = { ( x , y ) | x 2 4 + y 2 1 , ( x 1 ) 2 + y 2 1 4 , ( x + 1 ) 2 + y 2 1 4 } ,该问题的解析解为 u = e x + y

Figure 4. The comparison between the exact solution (right panel) and the numerical solution (left panel) with M = N = 80

图4. 在 M = N = 80 下的解析解(右边)和数值解(左边)的对比

Figure 5. The error obtained by Monte Carlo Markov chain method with M = N = 80

图5. 在 M = N = 80 下的数值解与解析解的误差

Table 1. The numerical convergence order by Monte Carlo Markov chain method

表1. 蒙特卡洛马尔科夫链方法数值收敛阶情况

图4展示了在 M = N = 80 下不规则区域上Poisson方程得到的精确解和数值解的比较。可以看出,两者的结果是一致的。图5展示了不规则区域上Poisson方程的数值解与解析解之间的误差,可以清楚的看出两者之间的误差很小。通过表1可以看出由蒙特卡洛马尔科夫链方法获得的数值收敛阶始终在2阶左右,其效果是显著的。

5. 结论

本文引入一种新的蒙特卡罗马尔可夫链方法求解不规则区域上的Poisson方程。数值结果表明了该方法的有效性和适应性。该方法为求解不规则区域上的微分方程提供了一种新的思路。此外,该方法简单,易于编程,还可推广到其它更复杂的微分方程。

文章引用

郑成丰,闫志忠. 不规则区域上Poisson方程的蒙特卡洛马尔科夫链方法求解
The Method of Monte Carlo Markov Chain for Solving the Poisson Equation on Irregular Domain[J]. 应用数学进展, 2019, 08(02): 181-187. https://doi.org/10.12677/AAM.2019.82020

参考文献

  1. 1. 王忆锋, 唐利斌. 利用有限差分和MATLAB矩阵运算直接求解二维泊松方程[J]. 红外技术, 2010, 32(4): 213-216.

  2. 2. 邵肖伟, 刘政凯, 李厚强. 一种基于Poisson方程的分离型图像修复方法[J]. 电路与系统学报, 2008, 13(6): 1-6.

  3. 3. 张建桥. 基于泊松方程的数字图像无缝拼合[J]. 现代电子技术, 2010, 33(17): 139-141.

  4. 4. 张琦, 周冉辉, 刘睿, 等. 基于泊松方程的磁罗盘磁域自差修正[J]. 舰船电子工程, 2011, 31(9): 50-53.

  5. 5. Nakamura, T. and Yabe, T. (1999) Cubic Interpolated Propagation Scheme for Solving the Hy-per-Dimensional Vlasov-Poisson Equation in Phase Space. Computer Physics Communications, 120, 122-154. https://doi.org/10.1016/S0010-4655(99)00247-7

  6. 6. Frey, W.H. (1977) Flexible Finite-Difference Stencils from Isoparametric Finite Elements. International Journal for Numerical Methods in Engineering, 11, 1653-1665. https://doi.org/10.1002/nme.1620111103

  7. 7. Frind, E.O. and Pinder, G.F. (1979) A Collocation Finite Element Method for Potential Problems in Irregular Domains. International Journal for Numerical Methods in Engineering, 14, 681-701. https://doi.org/10.1002/nme.1620140505

  8. 8. Farnoosh, R. and Ebrahimi, M. (2008) Monte Carlo Method for Solving Fredholm Integral Equations of the Second Kind. Applied Mathematics & Computation, 195, 309-315. https://doi.org/10.1016/j.amc.2007.04.097

  9. 9. Vajargah, B.F. and Vajargah, K.F. (2007) Monte Carlo Method for Finding the Solution of Dirichlet Partial Differential Equations. Applied Mathematical Sciences, 1, 453-462.

  10. 10. Gu, K. and Sadiku, M.N.O. (2000) Absorbing Markov Chain Solution for Possion’s Equation. Pro-ceedings of the IEEE Southeast Con 2000 Preparing for the New Millennium, Nashville, TN, 9-9 April 2000.

NOTES

*通讯作者

期刊菜单