带有弱奇异核的四阶偏积分微分方程的Sinc-Galerkin方法
The Sinc-Galerkin Method for the Fourth-Order Partial Integro-Differential Equation with Weakly Singular Kernels
摘要: 本文我们利用Sinc-Galerkin方法求解带有弱奇异核的四阶偏积分微分方程。首先在时间上借助L1格式和梯形卷积求积公式分别离散分数阶导数和积分,其次在空间上利用Sinc-Galerkin方法近似四阶偏导项,得到方程的全离散格式。最后推导出数值格式的收敛阶并通过数值算例来验证该方法的准确性和有效性。
Abstract: The Sinc-Galerkin method is considered and analyzed for solving the fourth-order partial integro-differential equation with weakly singular kernels. At first, for the temporal direction, we use L1 scheme to approximate Caputo derivative and the trapezoidal convolution quadrature rule to discretize the Riemann-Liouville fractional integral term. Then for space, we utilize Sinc-Galerkin method to deal with fourth-order partial derivative and obtain the fully discrete scheme. Finally, we deduce the convergence order of the numerical scheme and verify the accuracy and effectiveness of the proposed method through a numerical example.
文章引用:贾毅. 带有弱奇异核的四阶偏积分微分方程的Sinc-Galerkin方法[J]. 应用数学进展, 2022, 11(2): 651-663. https://doi.org/10.12677/AAM.2022.112072

1. 引言

近年来,带有弱奇异核的积分微分方程在连续介质力学、热力学、电动力学和种群生物学等领域得到了广泛的应用 [1],因此受到广大学者的关注和研究。针对这类方程,Hu等 [2] 提出向后有限差分方法,并证明了数值格式的稳定性和收敛性,Mohebbi等 [3] 提出高阶紧致差分格式,Chen等 [4] 提出了积分微分方程的二阶BDF交替方向隐式差分格式,Slodicka等 [5] 利用向后差分和线性有限元方法研究带有弱奇异型记忆项的抛物型方程。

本文考虑如下带有弱奇异核的分数阶积分微分方程:

{ D 0 C t β u ( x , t ) + I ( α ) u x x x x ( x , t ) = f ( x , t ) x ( a , b ) , t ( 0 , T ] , u ( x , 0 ) = u 0 ( x ) , x ( a , b ) , u ( a , t ) = u ( b , t ) = u x ( a , t ) = u x ( b , t ) = 0 , t ( 0 , T ] . (1)

其中 0 < β , α < 1 u x x x x = 4 u x 4 D 0 C t β u ( x , t ) 是关于t的分数阶Caputo导数:

D 0 C t β u ( x , t ) = 1 Γ ( 1 β ) 0 t u ( x , s ) s ( t s ) β d s . (2)

I ( α ) u ( x , t ) 是关于t的分数阶 Riemann-Liouville积分:

I ( α ) u ( x , t ) = 1 Γ ( α ) 0 t u ( x , s ) ( t s ) α 1 d s . (3)

在本文中,我们利用Sinc-Galerkin方法离散方程(1)中的四阶偏导项。Sinc-Galerkin方法最初是Stenger [6] 在研究Whittaker基数的基础上提出的,并应用该方法求解二阶微分方程。Smith等 [7] 深入研究了该方法,证明了该方法的收敛速率为 O ( e k M ) ,给出相关矩阵的性质,并利用该方法求解了四阶微分方程。此后有关Sinc-Galerkin方法的研究与应用大量涌现出来。Zarebnia等 [8] 利用该方法研究Troesch问题,并验证了Sinc-Galerkin的指数收敛速度。El-Gamel等 [9] 利用Sinc-Galerkin方法求解六阶两点边值问题,并验证了算法的有效性。Rashidinia等 [10] 在研究非线性两点边值问题中比较Sinc-Galerkin方法和Sinc配置方法,根据数值算例发现两种方法精度基本一致。Qiu等 [11] 利用该方法求解四阶偏积分微分方程,并将相关性质进行了推广。

2. 预备知识

将区域[0,T]均匀剖分,设 τ = T / N 为时间步长,令 t n = n τ n = 0 , 1 , 2 , , N 。为了方便,假设本文中的 C , C 1 , C 2 , 均为不依赖变量的正常数,符号 C ( a , b ) 表示与参数 a , b 有关的正常数。记 exp ( y ) = e y 。接下来给出分数阶Caputo导数的L1格式。

引理2.1 [12] 假设 0 < β < 1 v ( t ) [ 0 , T ] ,则有

D 0 C t β v ( t n ) = τ β Γ ( 2 β ) [ a 0 ( β ) v ( t n ) i = 1 n 1 ( a n i 1 ( β ) a n i ( β ) ) v ( t i ) a n 1 ( β ) v ( t 0 ) ] + R 1 n , (4)

其中 a n i ( β ) = ( n i + 1 ) 1 β ( n i ) 1 β | R 1 n | C τ 2 β

为了离散 t n 处的分数阶Riemann-Liouville积分 I ( α ) φ ,下面给出梯形卷积求积公式 [13]

Q ( α ) ( φ ) = τ α p = 0 n w p ( α ) φ n p + τ α w ˜ n ( α ) φ 0 , (5)

其中权 w p ( α ) 可以通过 ( 2 ( 1 z ) 1 + z ) α = p = 0 w p ( α ) z p 得到,修正权 w ˜ n ( α ) = n α Γ ( α + 1 ) p = 0 n w p ( α )

引理2.2 [4] 令 0 ϑ < 1 n 1 φ ( t ) 在区间 0 < t T 连续可导, B ϑ ( 0 , T ] = { z | | z | ϑ = sup 0 < t T t ϑ | z ( t ) | < } 。当 φ t t B ϑ ( 0 , T ] ,则存在常数C,使得

| I ( α ) φ ( t n ) Q ( α ) ( φ ) | = | R 2 n | C τ 2 ( t n α 1 | φ t ( 0 ) | + t n α ϑ | φ t t | ϑ ) . (6)

引理2.3 [11] 假设 0 < α < 1 n 0 ,则对于 w p ( α ) w ˜ n ( α )

( 1 ) w n ( α ) = 2 α j = 0 n κ j ( α ) σ n j ( α ) , κ j ( α ) = Γ ( α + j ) Γ ( α ) Γ ( j + 1 ) , σ n j ( α ) = Γ ( α + 1 ) Γ ( α n + j + 1 ) Γ ( n j + 1 ) , j 0 , ( 2 ) | w n ( α ) | = O ( τ t n α 1 ) , | w ˜ n ( α ) | C ( α ) τ t n α 1 , n 1 , ( 3 ) p = 0 n | w p ( α ) | C ( α , T ) , 0 n N . (7)

接下来,介绍Sinc函数及其相关性质。在实轴 上,Sinc函数的定义为

Sinc ( x ) = { sin ( π x ) π x , x 0 , 1 , x = 0. (8)

给定步长h,第k次平移的Sinc函数为

S ( k , h ) ( x ) = Sinc ( x k h h ) , k = 0 , ± 1 , ± 2 , . (9)

令v为定义在实轴上的可解析的有界函数,若级数 r ( v , h ) ( x ) = k = + v ( k h ) S ( k , h ) ( x ) 收敛,则 r ( v , h ) ( x ) 为函数v的Whittaker基数展开式,且有 r ( v , h ) ( k h ) = v ( k h ) 。为了近似有限区间 Ω = ( a , b ) 上的函数 v ( x ) ,引入映射 ϕ

ϕ ( x ) = log ( x a b x ) . (10)

且有 ϕ ( a ) = ϕ ( b ) = + 。映射 ϕ 将区域

D = { x = m + i n : | arg ( x a b x ) | < d π 2 } (11)

映射到能够推导出Whittaker基数展开式性质的区域

D d = { w = u + i v : | v | < d π 2 } . (12)

因此有限区间上Sinc方法的基函数可以表示为

S k ( x ) : = S ( k , h ) ° ( ϕ ( x ) ) = Sinc ( ϕ ( x ) k h h ) , k = 0 , ± 1 , ± 2 , . (13)

区间 ( a , b ) 上的Sinc方法网格节点为

x j = ϕ 1 ( j h ) = a + b exp ( j h ) 1 + exp ( j h ) . (14)

应用Sinc-Galerkin方法求解问题会得出加权内积项,接下来我们介绍Sinc梯形求积公式的定义和相关性质。

定义2.1 [7] 令函数 ψ ϕ 的逆映射, B ( D ) 为在区域D上可解析的函数F的集合,且满足

ψ ( x + L ) | F ( z ) d z | 0 , x ± , (15)

其中 L = i v : | v | < d π 2 ,且在边界上(记为 D )满足

H ( F ) = D | F ( z ) d z | < . (16)

定理2.1 [7] 令 F B ( D ) i = 1 ,给定充分小的h,则有

Ω F ( z ) d z h j = + F ( z j ) ϕ ( z j ) = i 2 D F ( z ) k ( ϕ , h ) ( z ) sin ( π ϕ ( z ) / h ) I F , (17)

其中

k ( ϕ , h ) ( z ) = exp ( i π h ϕ ( z ) sgn ( Im ( ϕ ( z ) ) ) ) . (18)

利用 | k ( ϕ , h ) ( z ) | | z D = exp ( π d h ) ,可以得到 | I F | exp ( π d / h ) 2 sinh ( π d / h ) H ( F )

定理2.2 [7] 假设 η γ L ˜ 为正常数,函数F和 ϕ 的定义与定理2.1一致,则有

| F ( z ) ϕ ( z ) | L ˜ { exp ( η | ϕ ( z ) | ) , z ψ ( ( , 0 ) ) , exp ( γ | ϕ ( z ) | ) , z ψ ( ( 0 , + ) ) . (19)

因此有限项Sinc梯形求积公式可改写为

| Ω F ( z ) d z h j = M N F ( z j ) ϕ ( z j ) | L ˜ ( 1 η exp ( η M h ) + 1 γ exp ( γ N h ) ) + | I F | , (20)

其中 | I F | C exp ( π d / h ) 。令 h = ( π d ) / ( η M ) N = η M / γ + 1 (这里 表示下取整),则有

Ω F ( z ) d z = h j = M N F ( z j ) ϕ ( z j ) + O ( exp ( π d η M ) ) . (21)

引理2.4 [7] 令 S k ( x ) = S ( k , h ) ° ( ϕ ( x ) ) δ k j ( l ) = h l [ d l d ϕ l [ S ( j , h ) ° ϕ ( x ) ] ] | x = x j l = 0 , 1 , 2 , 3 , 4 ,则Sinc导函数在节点处的值为

δ k j ( 0 ) = { 0 , j k , 1 , j = k , δ k j ( 1 ) = { ( 1 ) j k j k , j k , 0 , j = k , δ k j ( 2 ) = { 2 ( 1 ) j k ( j k ) 2 , j k , π 2 3 , j = k , δ k j ( 3 ) = { ( 1 ) j k ( j k ) 3 [ 6 π 2 ( j k ) 2 ] , j k , 0 , j = k , δ k j ( 4 ) = { 4 ( 1 ) j k ( j k ) 4 [ 6 π 2 ( j k ) 2 ] , j k , π 4 5 , j = k . (22)

引理2.5 [7] 令 d > 0 k ,对于 h > 0 l = 0 , 1 , 2 , 3 , 4 ,有

| S k ( l ) ( z ) sin ( π ϕ ( z ) / h ) | z D C l ( h , d ) , (23)

其中

C 0 ( h , d ) = h π d , C 1 ( h , d ) = 1 d C 0 ( h , d ) + 1 d tanh ( π d h ) , C 2 ( h , d ) = 2 d C 1 ( h , d ) + π h d , C 3 ( h , d ) = 3 d C 2 ( h , d ) + π 2 h 2 d tanh ( π d h ) , C 4 ( h , d ) = 4 d C 3 ( h , d ) + π 3 h 3 d . (24)

3. Sinc方法的离散格式构造及收敛阶分析

3.1. 半离散格式

u p ( x ) = u ( x , t p ) p = 0 , 1 , 2 , , n f n ( x ) = f ( x , t n ) x Ω 0 n N 。利用引理2.1和引理2.2则方程(1)在 t n 的半离散格式为

L u n ( x ) = τ β Γ ( 2 β ) [ a 0 ( β ) u n ( x ) i = 1 n 1 ( a n i 1 ( β ) a n i ( β ) ) u n ( x ) a n 1 ( β ) u 0 ( x ) ] + τ α p = 0 n w p ( α ) u x x x x n p ( x ) + τ α w ˜ n ( α ) u x x x x 0 ( x ) + R t n = f n ( x ) . (25)

且很容易得到 R t n = R 1 n + R 2 n = C τ 2 β + C τ 1 + α

3.2. 全离散格式

Sinc方法的逼近格式为

(26)

其中未知系数 { u j } 通过对残差正交化得出,即

L u M n ( x ) f n , S k = 0 , M k N . (27)

这里 , 为加权内积,其中权函数 w ( x ) 的选取与偏微分方程的边界条件,定义的区域有关 [7]。为了利用Sinc求积公式逼近积分项,这里采取另一种分析方法:

0 L u n f n , S k = τ β a 0 ( β ) Γ ( 2 β ) u n , S k τ β Γ ( 2 β ) p = 1 n 1 ( a n i 1 ( β ) a n i ( β ) ) u p , S k τ β a n 1 ( β ) Γ ( 2 β ) u 0 , S k + τ α p = 0 n w p ( α ) u x x x x n p , S k + w ˜ n ( α ) u x x x x 0 , S k + R t n , S k f n , S k . (28)

因此上式中的加权内积项可以被定理2.1和定理2.2精确逼近。

我们首先讨论 u x x x x n , S k ,利用四次分部积分得到

(29)

其中 ( , ) 表示 L 2 内积,且边界项 B T

B T = { u x x x n ( w S k ) u x x n ( w S k ) x + u x n ( w S k ) x x u n ( w S k ) x x x } ( x ) | a b (30)

应用定理2.1中的求积公式可得

u x x x x n , S k = h j = u n ( x j ) ϕ ( x j ) ( w S k ) x x x x ( x j ) + B T + I F ( 1 ) , (31)

其中积分误差项 I F ( 1 ) 是将式(17)中的F替换为 u n ( w S k ) x x x x 得到的。这里为了计算方便,将无限项求和截取为从-M到N的有限项求和。假设 1 n N ,选取合适的权函数,使得 B T = 0 ,利用 δ k j ( l ) 的定义,可以得到

u x x x x n , S k = h j = M N u n ( x j ) ϕ ( x j ) [ l = 0 4 δ k j ( l ) h l g ( l ) ( x j ) ] + I F ( 1 ) . (32)

其中

{ g ( 0 ) ( x ) = w ( 4 ) ( x ) , g ( 1 ) ( x ) = { 4 ϕ w + 6 ϕ w + 4 ϕ w + ϕ w } ( x ) , g ( 2 ) ( x ) = { 6 [ ϕ ] 2 w + 12 ϕ ϕ w + 4 ϕ ϕ w + 3 [ ϕ ] 2 w } ( x ) , g ( 3 ) ( x ) = { 4 [ ϕ ] 3 w + 6 [ ϕ ] 2 ϕ w } ( x ) , g ( 4 ) ( x ) = { [ ϕ ] 4 } ( x ) . (33)

同理得到剩下的积分的逼近格式

G , S k = h w ( x k ) G ( x k ) ϕ ( x k ) + I F ( 2 ) , (34)

其中G为(28)中的 u n u p u 0 f n R t n ,积分项误差 I F ( 2 ) 是分别将式(17)中F替换为 u n w S k u p w S k u 0 w S k f n w S k R t n w S k 得到的。

假设 u M n ( x ) 为方程(1)的近似解,记 u j n u ( x j , t n ) ,令 u ( x j , t n ) k = M , M + 1 , , N 1 n N 。令 μ = τ β / Γ ( 2 β ) ,则确定未知系数 { u j n } 的Sinc-Galerkin全离散格式为

μ a 0 ( β ) w ( x k ) ϕ ( x k ) u k n μ p = 1 n 1 ( a n p 1 ( β ) a n p ( β ) ) w ( x k ) ϕ ( x k ) u p n μ a n 1 ( β ) w ( x k ) ϕ ( x k ) u 0 n = τ α p = 0 n w p ( α ) j = M N u j n p ϕ ( x j ) [ l = 0 4 δ k j ( l ) h l g ( l ) ( x j ) ] w ˜ n ( α ) j = M N u j 0 ϕ ( x j ) [ l = 0 4 δ k j ( l ) h l g ( l ) ( x j ) ] + w ( x k ) ϕ ( x k ) f n ( x k ) . (35)

D ( ) 为对角线上元素为 s ( x j ) M + N + 1 阶对角矩阵, j = M , , N U n = [ u M , , u N ] T F n = [ f M , , f N ] T I ( l ) = [ δ k j ( l ) ] M × M 。令 B = D ( w ϕ ) P = l = 0 4 I ( l ) h l D ( g ( l ) ( x j ) ϕ ( x j ) ) A = μ a 0 ( β ) B + τ α w 0 ( α ) P 。则可以得到式(35)的矩阵形式:

A U n = μ p = 1 n 1 ( a n p 1 ( β ) a n p ( β ) ) B U p + μ a n 1 ( β ) B U 0 τ α p = 0 n 1 w n p ( α ) P U p w ˜ n ( α ) P U 0 + B F n . (36)

3.3. 收敛阶分析

给定函数 H ( ) 满足定义2.1,函数 C l ( h , d ) 满足引理2.5, G ( ) 满足(34),对于(28),利用定理2.1和引理2.5可以得到

| u x x x x n , S k h j = u n ( x j ) ϕ ( x j ) ( w S k ) x x x x ( x j ) | = | I F ( 1 ) | exp ( π d / h ) 2 D l = 0 4 | u n ( x ) g ( l ) ( x ) | | S k ( l ) ( x ) sin ( π ϕ ( x ) / h ) | | d x | exp ( π d / h ) 2 [ l = 0 4 C l ( h , d ) H ( u n g ( l ) ) ] C 1 exp ( π d / h ) , (37)

| G , S k h w ( x k ) G ( x k ) ϕ ( x k ) | = | I F ( 0 ) | 1 2 C 0 ( h , d ) H ( w G ) exp ( π d / h ) . (38)

在实际计算中,我们需要将式(37)中的无限项求和截断为有限项,即

| u x x x x n , S k h j = M N l = 0 4 u n ( x j ) ϕ ( x j ) h l δ k j ( l ) g ( l ) | | I F ( 1 ) | + h | j = M 1 l = 0 4 u n ( x j ) ϕ ( x j ) h l δ k j ( l ) g ( l ) | + h | j = N + 1 l = 0 4 u n ( x j ) ϕ ( x j ) h l δ k j ( l ) g ( l ) | . (39)

接着继续分析(39)中级数的截断项,令 0 n N 0 l 4 ,通过定理2.2可以得到

| u n ( x ) g ( l ) ( x ) ϕ ( x ) | L ˜ { exp ( η | ϕ ( x ) | ) x ( a , a + b 2 ) , exp ( γ | ϕ ( x ) | ) x [ a + b 2 , b ) . (40)

通过引理2.4 和 δ k j ( l ) 的定义,得到如下不等式

| δ k j ( 0 ) | 1 , | δ k j ( 1 ) | 1 , | δ k j ( 2 ) | π 2 3 , | δ k j ( 3 ) | 2 π 2 3 4 , | δ k j ( 4 ) | π 4 5 . (41)

结合(39),(40)和(41),可以得到

h | j = M 1 l = 0 4 u n ( x j ) ϕ ( x j ) h l δ k j ( l ) g ( l ) | h j = M + 1 | l = 0 4 u n ( x j ) ϕ ( x j ) g ( l ) ( x j ) | | δ k , j ( l ) h l | L ˜ h ( 1 + 1 h + π 2 3 h 2 + 2 π 2 3 4 h 3 + π 4 5 h 4 ) j = M + 1 e η j h L ˜ η C 2 e η M h , (42)

其中 C 2 = ( 1 + 1 h + π 2 3 h 2 + 2 π 2 3 4 h 3 + π 4 5 h 4 ) ,同理对另一侧截断级数分析可以得到

h | j = N + 1 l = 0 4 u n ( x j ) ϕ ( x j ) h l δ k j ( l ) g ( l ) | L ˜ γ C 2 e η N h . (43)

因此,将(42)和(43)代入(39),得到

| u x x x x n , S k h j = M M u n ( x j ) ϕ ( x j ) ( l = 0 4 δ k j ( l ) h l g ( l ) ( x j ) ) | L ˜ C 2 ( exp ( η M h ) η + exp ( γ N h ) γ ) + C 1 exp ( π d / h ) χ ( h , M ) . (44)

这里选取 h = ( π d ) / ( η M ) N = η M / γ + 1 ,注意到 C 1 = 1 2 [ l = 0 4 C l ( h , d ) H ( u n g ( l ) ) ] ,结合引理2.5有

C l ( h , d ) C h 1 l = C ( η π d ) l 1 2 M l 1 2 C M l 1 2 . (45)

因此 C 1 C M 3 2 ,且有

χ ( h , M ) C 3 M 2 exp ( π d η M ) + C M 3 2 exp ( π d η M ) C M 2 exp ( π d η M ) , (46)

其中 C 3 = L ˜ [ π 4 5 ( η π d ) 2 + 2 π 2 3 4 ( η π d ) 3 / 2 + π 2 3 ( η π d ) + ( η π d ) 1 / 2 + 1 ] ( 1 γ + 1 η ) 。最后讨论带有权和修正权的误差估计,结合引理2.3,(44)和(46),可以得到

| τ α p = 0 n w p ( α ) [ u x x x x n , S k h j = M N u n ( x j ) ϕ ( x j ) ( l = 0 4 δ k j ( l ) h l g ( l ) ( x j ) ) ] | τ α p = 0 n | w p ( α ) | C M 2 exp ( π d η M ) C ( α , T ) τ α M 2 exp ( π d η M ) , (47)

| w ˜ n ( α ) [ u x x x x 0 , S k h j = M N u 0 ( x j ) ϕ ( x j ) ( l = 0 4 δ k j ( l ) h l g ( l ) ( x j ) ) ] | | w ˜ n ( α ) | C M 2 exp ( π d η M ) C τ t n α 1 M 2 exp ( π d η M ) C τ α n 1 α M 2 exp ( π d η M ) C τ α M 2 exp ( π d η M ) . (48)

注意到 a 0 β a n β p = 1 n 1 ( a n p 1 ( β ) a n p ( β ) ) 有界,结合上述证明,可以得到(28)中其余积分项的估计。为了表述方便,其结果由下述定理给出。

定理3.1令映射 ϕ ,Sinc网格节点 x j ,函数 H ( ) ,函数 C l ( h , d ) ,集合 B ( D ) 的定义与前文一致,给定l, l = 0 , 1 , 2 , 3 , 4 1 n N ,选择 h = ( π d ) / ( η M ) N = η M / γ + 1 ,对于式(28)中的每一个积分项有如下估计:

1) 令 w u p B ( D ) p = 1 , 2 , , n 1 ,则有

τ β Γ ( 2 β ) | p = 1 n 1 ( a n p 1 ( β ) a n p ( β ) ) [ u p , S k h w ( x k ) u p ( x k ) ϕ ( x k ) ] | C N β ( a 0 ( β ) a n 1 ( β ) ) H ( w u p ) T β 2 Γ ( 2 β ) M 1 2 e π d η M . (49)

2) 令 w u 0 B ( D ) ,则有

τ β a 0 ( β ) Γ ( 2 β ) | u 0 , S k h w ( x k ) u 0 ( x k ) ϕ ( x k ) | C N β a 0 ( β ) H ( w u 0 ) T β 2 Γ ( 2 β ) M 1 2 e π d η M . (50)

3) 令 w u n 1 B ( D ) ,则有

τ β a n 1 ( β ) Γ ( 2 β ) | u n 1 , S k h w ( x k ) u n 1 ( x k ) ϕ ( x k ) | C N β a n 1 ( β ) H ( w u n 1 ) T β 2 Γ ( 2 β ) M 1 2 e π d η M . (51)

4) 令 w f n B ( D ) ,则有

| f n , S k h w ( x k ) f n ( x k ) ϕ ( x k ) | C H ( w f n ) 2 M 1 2 e π d η M . (52)

5) 令 w R t n B ( D ) δ = min { 1 + α , 2 β } ,则有

| R t n , S k | | h w ( x k ) R t n ( x k ) ϕ ( x k ) | + C 0 ( h , d ) 2 H ( w R t n ) e π d / h C τ δ + C H ( w R t n ) 2 M 1 2 e π d η M . (53)

6) 令 u n g ( l ) B ( D ) ,则有

| τ α p = 0 n w p ( α ) [ u x x x x n , S k h j = M N u n ( x j ) ϕ ( x j ) ( l = 0 4 δ k j ( l ) h l g ( l ) ( x j ) ) ] | C ( α , T ) τ α M 2 e π d η M , (54)

| w ˜ n ( α ) [ u x x x x 0 , S k h j = M N u 0 ( x j ) ϕ ( x j ) ( l = 0 4 δ k j ( l ) h l g ( l ) ( x j ) ) ] | C τ α M 2 e π d η M . (55)

定理3.1 给出方程(1)的Sinc-Galerkin全离散格式的精确估计。令 δ = min { 1 + α , 2 β } ,则全离散格式的收敛阶为 O ( τ δ + ( N β M 1 2 + M 2 ) exp ( π d η M ) )

4. 数值实验

定义误差的无穷范数为

E r r o r ( M , N ) = max | u ( x j , t n ) u j n | .

时间和空间上的收敛阶分别为

R a t e t = log N 2 N 1 ( E r r o r ( M , N 1 ) E r r o r ( M , N 2 ) ) , R a t e x = log M 2 M 1 ( E r r o r ( M 1 , N ) E r r o r ( M 2 , N ) ) . (56)

下面通过一个数值算例验证收敛阶。给定区域 Ω = ( 0 , 1 ) ,令 d = π / 2 T = 1 ,取 η = γ = 1 h = π / 2 M ,则空间上的误差为 O ( exp ( π M / 2 ) ) 。Sinc-Galerkin方法的权函数为 w ( x ) = 1 / ( ϕ ( x ) ) 3 / 2 ,记 B ( p , q ) 为Beta函数,令精确解为 u ( x , t ) = ( 1 + t 1 + α ) x 5 2 ( log x ) 3 ,对应的源项为

f ( x , t ) = ( 1 + α ) B ( 1 + α , 1 β ) Γ ( 1 β ) t 1 + α β x 5 2 ( log x ) 3 + ( t α Γ ( α + 1 ) + Γ ( α + 2 ) t 2 α + 1 Γ ( 2 α + 2 ) ) 3 [ 112 log x 16 ( log x ) 2 5 ( log x ) 3 + 128 ] 16 x 3 / 2 . (57)

首先计算 L 范数下的时间误差和收敛阶。在空间上令 M = 32 ,选择不同的 N 进行验证,结果如图1表1所示。在表1中分别选择 α = 0.25 β = 0.6 α = 0.7 β = 0.35 得到 L 范数下的时间误差和收敛阶。在图1中,分别画出了数值模拟和理论估计在时间上的收敛阶图像。结合表1的数据和图1的结果来看,在时间上收敛阶与理论分析一致。接着令 N = 10000 ,选择不同的M得到相应的空间Sinc点,完成数值模拟。接着选择四组不同的 α , β 验证,在表2中得到四组 L 范数下的空间误差和收敛阶。选择 α = 0.45 β = 0.6 ,在图2中得到空间上的收敛阶曲线并与理论估计值比较。综合来看,在空间上收敛阶与理论分析也一致,均为指数收敛。

Table 1. The maximum-norm errors, temporal convergence rates with different α , β when M = 32

表1. M = 32 L 范数下的时间误差和收敛阶

Table 2. The maximum-norm errors, spatial convergence rates with different α , β when N = 10000

表2. N = 10000 L 范数下的空间误差和收敛阶

Figure 1. The temporal convergence rates when M = 32

图1. M = 32 时得到的时间收敛阶

Figure 2. The spatial convergence rate when N = 10000

图2. N = 10000 时得到的空间收敛阶

5. 结论

本文中,我们利用Sinc-Galerkin方法求解带有弱奇异核的四阶偏积分微分方程。在时间上利用L1格式和梯形卷积求积公式分别离散分数阶导数和积分。在空间上针对四阶导数的复杂性,利用Sinc-Galerkin方法去近似它。相较有限差分方法,Sinc-Galerkin方法通过较少的点近似高阶导数并可以达到指数收敛。最后我们通过数值算例证明方法的准确性和有效性。结果表明,在空间上利用较少的点即可达到较高的精度,且收敛阶与理论估计一致,再一次证明方法的准确性和有效性。相比之下时间上的精度还不够理想,在今后的研究中我们将考虑分数阶算子的其他差分格式,使时间上的精度得到提高。

参考文献

[1] 刘发旺, 庄平辉, 刘青霞. 分数阶偏微分方程数值方法及其应用[M]. 北京: 科学出版社, 2015.
[2] Hu, S.F., Qiu, W.L. and Chen, H.B. (2019) A Backward Euler Difference Scheme for the Integro-Differential Equations with the Multi-Term Kernels. International Journal of Computer Mathematics, 6, 1254-1267.
https://doi.org/10.1080/00207160.2019.1613529
[3] Mohebbi, A. (2017) Compact Finite Difference Scheme for the Solution of a Time Fractional Partial Integro-Differential Equation with a Weakly Singular Kernel. Mathematical Methods in the Applied Sciences, 40, 7627-7639.
https://doi.org/10.1002/mma.4549
[4] Chen, H.B., Xu, D. and Peng, Y.L. (2017) A Second Order BDF Alternating Direction Implicit Difference Scheme for the Two-Dimensional Fractional Evolution Equation. Applied Mathematical Modelling, 41, 54-67.
https://doi.org/10.1016/j.apm.2016.05.047
[5] Slodicka, M. (1997) Numerical Solution of a Parabolic Equation with a Weakly Singular Positive-Type Memory Term. Electronic Journal of Differential Equations, 1997, 1-12.
[6] Stenger, F. (1979) A Sinc-Galerkin Method of Solution of Boundary Value Problems. Mathematics of Computation, 33, 85-109.
https://doi.org/10.1090/S0025-5718-1979-0514812-4
[7] Smith, R.C., Bogar, G.A., Bowers, K.L. and Lund, J. (1991) The Sinc-Galerkin Method for Fourth-Order Differential Equations. SIAM Journal on Numerical Analysis, 28, 760-788.
https://doi.org/10.1137/0728041
[8] Zarebnia, M. and Sajjadian, M. (2012) The Sinc-Galerkin Method for Solving Troesch’s Problem. Mathematical and Computer Modelling, 56, 218-228.
https://doi.org/10.1016/j.mcm.2011.11.071
[9] El-Gamel, M., Cannon, J.R. and Zayed, A.I. (2004) Sinc-Galerkin Method for Solving Linear Sixth-Order Boundary-Value Problems. Mathematics of Computation, 73, 1325-1343.
https://doi.org/10.1090/S0025-5718-03-01587-4
[10] Rashidinia, J. and Nabati, M. (2013) Sinc-Galerkin and Sinc-Collocation Methods in the Solution of Nonlinear Two-Point Boundary Value Problems. Computational & Applied Mathematics, 32, 315-330.
https://doi.org/10.1007/s40314-013-0021-y
[11] Qiu, W.L., Xu, D. and Guo, J. (2020) The Crank-Nicolson-Type Sinc-Galerkin Method for the Fourth-Order Partial Integro-Differential Equation with a Weakly Singular Kernel. Applied Numerical Mathematics, 159, 239-258.
https://doi.org/10.1016/j.apnum.2020.09.011
[12] Sun, Z. and Wu, X. (2006) A Fully Discrete Difference Scheme for a Diffusion-Wave System. Applied Numerical Mathematics, 56, 193-209.
https://doi.org/10.1016/j.apnum.2005.03.003
[13] Cuesta, E. and Palencia, C. (2003) A Fractional Trapezoidal Rule for Integro-Differential Equations of Fractional Order in Banach Spaces. Applied Numerical Mathematics, 45, 139-159.
https://doi.org/10.1016/S0168-9274(02)00186-1