泊松方程第一边值问题的谱配置方法
Spectral Collocation Method for Poisson Equation First Side Value Problem
DOI: 10.12677/aam.2024.135194, PDF, HTML, XML, 下载: 28  浏览: 58 
作者: 张雅津, 魏金涛:河南科技大学数学与统计学院,河南 洛阳
关键词: 泊松方程谱配置法Legendre多项式Poisson Equations Spectral Collocation Method Legendre Polynomial
摘要: 以Legendre-Gauss-Lobatto节点为配置点,利用Legendre多项式建立谱配置格式求解具有第一边值问题的泊松方程,给出算法格式,通过数值运算表明算法格式的有效性和高精度。
Abstract: Using Legendre-Gauss-Lobatto nodes as configuration points, a spectral collocation scheme is established using Legendre polynomials to solve Poisson’s equation with first boundary value problems. The algorithm format is given, and the effectiveness and high accuracy of the algorithm format are demonstrated through numerical operations.
文章引用:张雅津, 魏金涛. 泊松方程第一边值问题的谱配置方法[J]. 应用数学进展, 2024, 13(5): 2069-2077. https://doi.org/10.12677/aam.2024.135194

1. 引言

泊松方程是一类非常经典的线性椭圆偏微分方程,在电磁学、理论物理和工程计算等领域是很常见的数学问题,但是由于其计算的复杂性,泊松方程的精确解并不容易得到,随着计算机软件的蓬勃发展,求解偏微分方程的数值解变得越来越方便和重要,因此已经有不少学者投入到泊松方程第一边值问题数值解的研究中,其中有学者分别用有限元法 [1] 、有限差分法等 [2] 进行求解。本文将采用的谱配置 [3] 方法来逼近泊松方程第一边值问题,主要利用Lagrange插值 [4] 方法来逼近,并根据线性方程条件来构造数值格式,为研究泊松方程的数值解提供了一种新的方法。

2. 基于Legendre-Gauss-Lobatto节点Lagrange插值多项式及其微分矩阵

L N ( x ) ( x I = ( 1 , 1 ) ) ,表示N阶Legendre多项式, x 0 = 1 x N = 1 以及 x j = ( 1 j N 1 ) ,是 ( 1 x 2 ) x L N ( x ) = 0 的根。则以 x j 为插值节点的Lagrange插值多项式为 [5] :

ϕ j ( x ) = ( x x 0 ) ( x x 1 ) ( x x j 1 ) ( x x j + 1 ) ( x x N ) ( x j x 0 ) ( x j x 1 ) ( x j x j 1 ) ( x j x j + 1 ) ( x j x N ) , j = 0 , 1 , , N

上式等价于:

ϕ j ( x ) = ω ( x ) ( x x j ) x ω ( x j ) (2-1)

其中 ω ( x ) = ( x x 0 ) ( x x 1 ) ( x x N )

P N ( I ) 为次数不超过N的多项式集合,对于任意的 u ( x ) C ( I ) ,则有Lagrange插值多项式如下:

p N ( x ) = j = 0 N u j ϕ j ( x ) , u j = u ( x j ) , x I

根据Legendre多项式的性质 [6] :

x ( ( 1 x 2 ) x L N ( x ) ) + N ( N + 1 ) L N ( x ) = 0

我们可以将(2-1)式转化为:

ϕ j ( x ) = c ( 1 x 2 ) x L N ( x ) ( x x j ) x ( c ( 1 x j 2 ) x L N ( x j ) ) = ( x 2 1 ) x L N ( x ) N ( N + 1 ) ( x x j ) L N ( x j )

其中c为常数。

x = x k 处对 p N ( x ) 求一阶导数,可得:

x p N ( x k ) = m = 0 N d k , j u j , k = 0 , 1 , , N

其中 d k , j = x ϕ j ( x k )

D = ( d k , j ) ( N + 1 ) × ( N + 1 ) 矩阵,根据文献 [7] 得到:

d k , j = { L N ( x k ) ( x k x j ) L N ( x j ) , k j , N ( N + 1 ) 4 , k = j = 0 , N ( N + 1 ) 4 , k = j = N 0 , 0 < k = j < N , (2-2)

d k , j ( 2 ) = x 2 ϕ j ( x k ) ,则在 x = x k 处对 p N ( x ) 求二阶导数有:

x 2 p N ( x k ) = j = 0 N d k , j ( 2 ) u j , k = 0 , 1 , , N

D ^ = ( d k , m ( 2 ) ) ,则可得 [7] : D ^ = D 2

类似地,若 D ( m ) 表示m阶微分矩阵,根据文献 [8] 有: D ( m ) = D m

同理,在y轴上,用 ψ m ( y ) 表示该方向上的插值多项式基函数,并记 d ^ l , j = t ψ j ( y l ) y l 为该方向的插值节点,则可以推导出与(2-2)式类似的微分矩阵表达式。

3. 泊松方程第一边值问题的单区域谱配置法

3.1. 算法格式

泊松方程第一边值问题为:

{ Δ u = f ( s , t ) , ( s , t ) D , u ( a , t ) = f 1 ( t ) , u ( b , t ) = f 2 ( t ) , t [ c , d ] , u ( s , c ) = f 3 ( s ) , u ( s , d ) = f 4 ( s ) , s [ a , b ] . (3-1)

其中 D = ( a , b ) × ( c , d ) ,这里 Δ 是一个二维拉普拉斯算子,即 Δ u = 2 u x 2 + 2 u y 2 f ( s , t ) f 1 ( t ) f 2 ( t ) f 3 ( s ) f 4 ( s ) 为已知函数, u ( s , t ) 为待求函数。

对(3-1)式作变换: x = 2 b a s b + a b a y = 2 d c t d + c d c ,并令 x = x Ω = ( 1 , 1 ) × ( 1 , 1 ) ,则上式问题转化为 [9] :

{ 4 ( b a ) 2 x x u 4 ( d c ) 2 y y u = f ( x , y ) , ( x , y ) Ω u ( 1 , y ) = f 1 ( y ) , u ( 1 , y ) = f 2 ( y ) , t [ 1 , 1 ] u ( x , 1 ) = f 3 ( x ) , u ( x , 1 ) = f 4 ( x ) , x [ 1 , 1 ] (3-2)

P M , N ( Ω ) = P N [ 1 , 1 ] × P M [ 1 , 1 ] ,则(3-2)式的配置方法就是用数值解 [10] [11] :

u M , N ( x , y ) = n = 0 N m = 0 M u ^ m , n ϕ n ( x ) ψ m ( y ) P M , N ( Ω ) (3-3)

将(3-3)式代入(3-2)式,得到:

{ 4 ( b a ) 2 n = 0 N m = 0 M u ^ m , n ψ m ( y l ) x x ϕ n ( x k ) 4 ( d c ) 2 n = 0 N m = 0 M u ^ m , n ϕ n ( x k ) y y ψ m ( y l ) = f ( x k , y l ) , k = 1 , 2 , , N 1 ; l = 1 , 2 , , M 1 , n = 0 N m = 0 M u ^ m , n ϕ n ( x 0 ) ψ m ( y l ) = f 1 ( y l ) , n = 0 N m = 0 M u ^ m , n ϕ n ( x N ) ψ m ( y l ) = f 2 ( y l ) , l = 0 , 1 , , M n = 0 N m = 0 M u ^ m , n ϕ n ( x k ) ψ m ( y 0 ) = f 3 ( x k ) , n = 0 N m = 0 M u ^ m , n ϕ n ( x k ) ψ m ( y M ) = f 4 ( x k ) , k = 0 , 1 , , N

由微分矩阵定义,以及 u M , N ( x k , y l ) = u ^ l , k ,将上式化简为:

{ 4 ( b a ) 2 n = 0 N u ^ l , n d k , n ( 2 ) + 4 ( d c ) 2 m = 0 M u ^ m , k d ^ l , m ( 2 ) = f ( x k , y l ) , k = 1 , 2 , , N 1 ; l = 1 , 2 , M 1 u ^ l , 0 = f 1 ( y l ) , u ^ l , N = f 2 ( y l ) , l = 0 , 1 , , M u ^ 0 , k = f 3 ( x k ) , u ^ M , k = f 4 ( x k ) , k = 0 , 1 , , N . (3-4)

根据已知的边界条件(3-4)式可以等价于:

{ 4 ( b a ) 2 n = 1 N 1 u ^ l , n d k , n ( 2 ) + 4 ( d c ) 2 m = 1 M 1 u ^ m , k d ^ l , m ( 2 ) = f ( x k , y l ) 4 ( b a ) 2 u ^ l , 0 d k , 0 ( 2 ) 4 ( b a ) 2 u ^ l , N d k , N ( 2 ) 4 ( d c ) 2 u ^ 0 , k d ^ l , 0 ( 2 ) 4 ( d c ) 2 u ^ M , k d ^ l , M ( 2 ) , k = 1 , 2 , , N 1 ; l = 1 , 2 , , M 1. (3-5)

将(3-5)式按照 l = 1 , 2 , , M 1 展开,可得:

l = 1 4 ( b a ) 2 ( u ^ 1 , 1 d k , 1 ( 2 ) + u ^ 1 , 2 d k , 2 ( 2 ) + + u ^ 1 , N 1 d k , N 1 ( 2 ) ) + 4 ( d c ) 2 ( u ^ 1 , k d ^ 1 , 1 ( 2 ) + u ^ 2 , k d ^ 1 , 2 ( 2 ) + + u ^ M 1 , k d ^ 1 , M 1 ( 2 ) ) = f ( x k , y 1 ) 4 ( b a ) 2 u ^ 1 , 0 d k , 0 ( 2 ) 4 ( b a ) 2 u ^ 1 , N d k , N ( 2 ) 4 ( d c ) 2 u ^ 0 , k d ^ 1 , 0 ( 2 ) 4 ( d c ) 2 u ^ M , k d ^ 1 , M ( 2 ) ,

l = 2 4 ( b a ) 2 ( u ^ 2 , 1 d k , 1 ( 2 ) + u ^ 2 , 2 d k , 2 ( 2 ) + + u ^ 2 , N 1 d k , N 1 ( 2 ) ) + 4 ( d c ) 2 ( u ^ 1 , k d ^ 2 , 1 ( 2 ) + u ^ 2 , k d ^ 2 , 2 ( 2 ) + + u ^ M 1 , k d ^ 2 , M 1 ( 2 ) ) = f ( x k , y 2 ) 4 ( b a ) 2 u ^ 2 , 0 d k , 0 ( 2 ) 4 ( b a ) 2 u ^ 2 , N d k , N ( 2 ) 4 ( d c ) 2 u ^ 0 , k d ^ 2 , 0 ( 2 ) 4 ( d c ) 2 u ^ M , k d ^ 2 , M ( 2 ) ,

l = M 1 , 4 ( b a ) 2 ( u ^ M 1 , 1 d k , 1 ( 2 ) + u ^ M 1 , 2 d k , 2 ( 2 ) u ^ M 1 , N 1 d k , N 1 ( 2 ) ) + 4 ( d c ) 2 ( u ^ 1 , k d ^ M 1 , 1 ( 2 ) + u ^ 2 , k d ^ M 1 , 2 ( 2 ) u ^ M 1 , k d ^ M 1 , M 1 ( 2 ) ) = f ( x k , y M 1 ) 4 ( b a ) 2 u ^ M 1 , 0 d k , 0 ( 2 ) 4 ( b a ) 2 u ^ M 1 , N d k , N ( 2 ) 4 ( d c ) 2 u ^ 0 , k d ^ M 1 , 0 ( 2 ) 4 ( d c ) 2 u ^ M , k d ^ M 1 , M ( 2 ) , k = 1 , 2 , , N 1.

根据上式可得如下矩阵形式:

再将上式按照 k = 1 , 2 , , N 1 展开,得到:

令:

A = ( d 1 , 1 ( 2 ) d 1 , 2 ( 2 ) d 1 , N 1 ( 2 ) d 2 , 1 ( 2 ) d 2 , 2 ( 2 ) d 2 , N 1 ( 2 ) d N 1 , 1 ( 2 ) d N 1 , 2 ( 2 ) d N 1 , N 1 ( 2 ) ) , B = ( d ^ ( 2 ) 1 , 1 d ^ ( 2 ) 1 , 2 d ^ ( 2 ) 1 , M 1 d ^ ( 2 ) 2 , 1 d ^ ( 2 ) 2 , 2 d ^ ( 2 ) 2 , M 1 d ^ ( 2 ) M 1 , 1 d ^ ( 2 ) M 1 , 2 d ^ ( 2 ) M 1 , M 1 ) ,

F 1 = ( u ^ 1 , 0 u ^ 2 , 0 u ^ M 1 , 0 ) ( d 1 , 0 ( 2 ) d 2 , 0 ( 2 ) d N 1 , 0 ( 2 ) ) + ( u ^ 1 , N u ^ 2 , N u ^ M 1 , N ) ( d 1 , N ( 2 ) d 2 , N ( 2 ) d N 1 , N ( 2 ) )

F 2 = ( d ^ ( 2 ) 1 , 0 d ^ ( 2 ) 2 , 0 d ^ ( 2 ) M 1 , 0 ) ( u ^ 0 , 1 u ^ 0 , 2 u ^ 0 , N 1 ) + ( d ^ ( 2 ) 1 , M d ^ ( 2 ) 2 , M d ^ ( 2 ) M 1 , M ) ( u ^ M , 1 u ^ M , 2 u ^ M , N 1 ) ,

C = ( f ( x 1 , y 1 ) f ( x 2 , y 1 ) f ( x N 1 , y 1 ) f ( x 1 , y 2 ) f ( x 2 , y 2 ) f ( x N 1 , y 2 ) f ( x 1 , y M 1 ) f ( x 2 , y M 1 ) f ( x N 1 , y M 1 ) ) , X = ( u ^ 1 , 1 u ^ 1 , 2 u ^ 1 , N 1 u ^ 2 , 1 u ^ 2 , 2 u ^ 2 , N 1 u ^ M 1 , 1 u ^ M 1 , 2 u ^ M 1 , N 1 )

简化后我们可以得到以下线性矩阵方程:

4 ( b a ) 2 X A T + 4 ( d c ) 2 B X = C 4 ( b a ) 2 F 1 4 ( d c ) 2 F 2 (3-6)

3.2. 数值结果

令:

其中Y、D、H1、H2分别是矩阵X、C、F1、F2按照行进行拉长后的转置向量,因此,我们可以将(3-3)式转换成如下线性方程组:

( 4 ( b a ) 2 E M 1 A + 4 ( d c ) 2 B E N 1 ) Y = D 4 ( b a ) 2 H 1 4 ( d c ) 2 H 2

其中 E n 表示n阶单位矩阵,“ ”表示Kronecker积,则上式等价于:

Y = ( 4 ( b a ) 2 E M 1 A + 4 ( d c ) 2 B E N 1 ) 1 ( D 4 ( b a ) 2 H 1 4 ( d c ) 2 H 2 ) (3-7)

本文中我们将使用 L -误差:

E M , N = max 1 l M 1 , 1 k N 1 | u ( x k , y l ) u M , N ( x k , y l ) |

来度量数值误差。

本文我们选取的泊松方程第一边值问题的例子为:

{ Δ u = ( s 2 + t 2 ) exp ( s t ) , ( s , t ) D , u ( 2 , t ) = exp ( 2 t ) , u ( 2 , t ) = exp ( 2 t ) , t [ 2 , 2 ] , u ( s , 2 ) = exp ( 2 s ) , u ( s , 1 ) = exp ( 2 s ) , s [ 2 , 2 ] . (3-8)

其精确解为:

u ( s , t ) = exp ( s t )

根据变换关系,我们可知(1-2)式的精确解为:

u ( x , y ) = exp ( ( b a 2 x + b + a 2 ) ( d c 2 y + d + c 2 ) )

图1图2是经过数值实验得到的误差结果图。图1表示当y轴方向插值多项式次数为 M = 18 时,最大值误差 E M , N 在x轴方向 随插值多项式次数N的变化情况。由图1可以看出,误差随N的增大而呈现近似线性递减的趋势,且x轴方向和y轴方向所用的节点数量平衡。

Figure 1. Change of E M , N at M = 18 , N = 3 : 1 : 18

图1. M = 18 N = 3 : 1 : 18 时误差 E M , N 的变化

图2表示当x轴方向插值多项式次数为 N = 18 时,最大值误差 E M , N 在y轴方向 随插值多项式次数M增加时的变化情况。同样,可以得出与图1类似的结论。并且我们不难从方程(3-1)式中发现函数关于x轴和y轴高度对称,因此两个方向的数值误差阶数应该是大致相同的,而这与我们通过数值实验得到的结果是吻合的。

Figure 2. Change of E M , N at N = 18 , M = 3 : 1 : 18

图2. N = 18 M = 3 : 1 : 18 时误差 E M , N 的变化

4. 总结

本文研究了在Legendre-Gauss-Lobatto节点配置下,利用Legendre多项式建立谱配置格式来求解具有第一边值问题的泊松方程的数值解。通过数值运算表明了该算法格式的有效性和高精度。且计算过程中x轴和y轴两个方向多项式次数的平衡有效地避免数值解在某一方向上的低精度问题,保证整体数值解的准确性,平衡的节点分布有助于算法的收敛速度更快。本研究为深入理解泊松方程的数值求解方法提供了一个重要的案例,同时为进一步探索和改进相关算法奠定了坚实的基础。

参考文献

[1] 卓晓华. 二维Poisson方程有限元解法[J]. 卫生职业教育, 2004(22): 76-77.
[2] 杜书德. 基于有限差分法的泊松方程第一类边值问题求解[J]. 科技通报, 2018, 34(4): 21-24.
[3] 马亚楠, 王天军, 李冰冰. Korteweg-de Vries方程的时空谱配置方法[J]. 数值计算与计算机应用, 2021, 42(4): 351-360.
[4] Orszag, S.A. (1972) Comparison of Pseudospectral and Spectral Approximation. Studies in Applied Mathematics, 51, 253-259.
https://doi.org/10.1002/sapm1972513253
[5] Kreiss, H.O. and Oliger, J. (1972) Comparison of Accurate Methods for the Integration of Hyperbolic Equations. Tellus, 24, 199-215.
https://doi.org/10.3402/tellusa.v24i3.10634
[6] 王亚洲, 秦国良, 和文强, 包振忠. 时空耦合谱元方法求解一维Burgers方程[J]. 西安交通大学学报, 2017, 51(1): 45-50.
[7] Shen, J., Tang, T. and Wang, L.-L. (2011) Spectral Methods: Algorithms, Analysis and Applications. In Springer Series in Computational Mathematics (Volume 41), Springer, Berlin, 103-105.
https://doi.org/10.1007/978-3-540-71041-7
[8] 王天军, 殷政伟. Legendre-Gauss-Lobatto节点的一个注记[J]. 河南科技大学学报(自然科学版), 2012, 33(1): 71-74.
[9] 高文, 胡晓, 吕军亮, 等. 泊松方程第一类边值问题四阶紧差分格式数值实现[J]. 高等学校计算数学学报, 2017, 39(1): 1-16.
[10] Canuto, C. and Quarteroni, A. (1982) Approximation Results for Orthogonal Polynomials in Sobolev Spaces. Mathematics of Computation, 38, 67-86.
https://doi.org/10.1090/S0025-5718-1982-0637287-3
[11] Canuto, C. and Quarteroni, A. (1982) Error Estimates for Spectral and Pseudo-Spectral Approximations of Hyperbolic Equations. SIAM Journal on Numerical Analysis, 19, 629-642.
https://doi.org/10.1137/0719044