不规则区域上poisson方程的蒙特卡洛马尔科夫链方法求解-凯发娱乐官网

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所示的情况。

图1中， ω 表示整个求解域，并且 γ 1 , γ 2 和 γ 3 表示求解域的三个边界。下面对整个不规则求解域进行均匀网格剖分来获取网格节点。然而由于不规则边界的存在，网格剖分获得的边界点无法均匀的分布在边界上，如图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所示。

图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。

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

5. 结论

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

参考文献

notes

*通讯作者

参考文献