多项分数阶非线性波动方程的数值方法及其快速实现

2022-10-20邵林馨沈卓旸马葛沁舟闵婕黄健飞

河北大学学报（自然科学版） 2022年5期

邵林馨，沈卓旸，马葛沁舟，闵婕，黄健飞

(扬州大学数学科学学院，江苏扬州 225002)

分数阶偏微分方程在反常现象和复杂系统的建模中发挥着越来越重要的作用[1-3],本文将要研究的带有空间四阶导数的多项时间分数阶非线性波动方程可以较准确地建模氧气通过毛细血管的动力学特征[3-4].大量文献表明：研究人员很难获得非线性分数阶偏微分方程的解析解[5-6]，有时候即使能获得，但由于这类解析解中往往含有无穷级数或无穷限积分，也仍然无法精确使用.因此，研究分数阶偏微分方程的数值方法越来越受到研究人员的关注[7-12]，以期为分数阶模型的实际应用和模拟提供算法支撑.

到目前为止，对于空间上具有四阶导数的时间分数阶波动方程的数值计算方法已有较多的研究.例如：Jafari等[13]利用Adomian分解方法得到了定义在有界空间域中的四阶分数阶扩散-波动方程的半解析解，并通过实例讨论了该方法的收敛性；Hu等[14]提出了一种四阶分数阶扩散波系统的有限差分方法，并证明了该方法的唯一性、稳定性和收敛性；Li等[15]对四阶分数阶扩散波动方程提出了一种带参数的5次样条方法，并用能量法严格地建立了数值方法的可解性、稳定性和收敛性；Gao等[16]给出了一类四阶时间二项分数阶波动方程的空间紧致差分格式，并证明了该差分格式的无条件稳定性和收敛性；刘新龙等[17]对时间分数阶四阶扩散方程构造了显-隐和隐-显2个差分格式，并证明了2个格式的稳定性和超收敛性，在精度一致的情况下，显-隐格式的计算效率要高于隐-显格式；Elmahdi等[18]对带空间四阶导数的时间分数阶非线性波动方程，基于等价变形的方法，构造了一个线性化差分格式和一个紧差分格式，并进行了相关的数值分析.另外，需要指出的是：对于空间上具有四阶导数的时间分数阶偏微分方程的数值方法也有很多研究，可以参考文献[19-22].

另一方面，关于多项时间分数阶偏微分方程数值方法的研究也受到了广泛关注.例如：Ren等[23]利用空间离散的紧致差分法和多项时间分数阶Caputo导数的L1 逼近，提出了一维和二维多项时间分数阶扩散波方程的有效数值格式，并加以证明无条件稳定性和全局收敛性；Zheng等[24]提出了一种基于时空谱方法求解多项时间分数扩散方程的高阶数值方法,该时空谱方法在时间和空间方向上都具有指数衰减；Abdel-Rehim等[25]通过数值计算得到了时间分数多项波动方程的近似解，包括时间分数波、强迫波和阻尼波动方程，并讨论了冯-诺依曼稳定性条件；Chen等[26]提出了一种具有可变系数的多项时间分数阶扩散和扩散波动方程的统一格式，该格式在时间方向采用有限差分近似和在空间方向采用 Legendre 谱近似；王芬玲等[27]对具有 Caputo 导数的二维多项时间分数阶扩散方程，采用线性三角元和改进 L1 格式，建立了一个全离散格式，并进行了高精度分析；Liu等[28]提出了一种二维多项时间分数混合扩散和扩散波动方程初边值问题的ADI谱方法，该方法可以处理非光滑问题，并能模拟粘弹性非牛顿流体的扩散和输运；Huang等[29]对一类二维时空多项分数阶偏微分方程设计了一个快速的ADI格式，该格式在时间方向只需要相对较低的光滑性.

从现有的文献资料来看，对带有空间四阶导数的多项时间分数阶非线性波动方程数值方法的研究还非常少.因此，本文将对上述多项分数阶非线性波动方程构造一个线性化数值方法，并进行误差分析.该线性化格式是基于等价的积分-偏微分方程来构造的，可以避免求解非线性方程组，能计算初始奇异性问题，并且还能进行快速计算.

1 问题的提出及预备知识

本文将考虑为以下形式的多项时间分数阶非线性波动方程[3-4]构造一个线性化的数值格式：

(1)

u(x,0)=φ(x),ut(x,0)=φ(x),0≤x≤L,

u(0,t)=u(L,t)=ux(0,t)=ux(L,t)=0,0

给出一个四阶导数的紧差分逼近，用来离散方程(1)中的空间方向四阶导数.

引理1[21]设u(x,·)∈C8[0,L],u(0,·)=u(L,·)=ux(0,·)=ux(L,·)=0,并令

(2)

为了便于计算，式(2)写成如下的矩阵形式

容易证明矩阵S是1个正定矩阵.

引理2[12]设u(·,t)∈C1[0,T]∩C2(0,T],且1<αK<αK-1<…<α2<α1<2，则

利用引理2，方程(1)可等价地转化为如下形式的偏微分-积分方程

(3)

引理3[29]设0<α<1,且u(·,t)∈C1[0,T],则有以下2个对Riemann-Liouville积分的逼近成立,

引理4[29]设u(·,t)∈C1[0,T]∩C2(0,T],0

在理论分析中将利用引理5，为此需要给出如下引理来说明引理3中定义的系数满足引理5的条件和结论:

2 格式的构造及其数值分析

2.1 格式的构造

假设u(·,t)∈C1[0,T]∩C2(0,T],u(x,·)∈C8[0,L],u(0,·)=u(L,·)=ux(0,·)=ux(L,·)=0.下面考虑在点(xi,tn)处来离散方程(3).具体地，用引理1来离散方程(3)中的空间四阶导数，用引理4来离散方程(3)中的时间一阶导数和时间Caputo导数，分别用引理3中的第1个逼近和第2个逼近来离散方程(3)中等号左边和右边的Riemann-Liouville积分，可得

其中，

上式两边同乘以τ，可得

(4)

(5)

从式(5)可以发现：该式是一个线性式，可以避免求解非线性方程.另外，可以采用文献[30-31]中的指数和技术来对式(5)进行快速实现，以提高计算的效率.值得一提的是：从解函数u(x,t)在时间方向的光滑性假设来看，即u(·,t)∈C1[0,T]∩C2(0,T]，格式(5)同样适用于如下初始奇异性的情况：u(·,t)=1+O(tαK),即当t→0时，u(x,t)在时间方向的二阶导数将趋于无穷.

2.2 数值分析

定理1式(5)存在唯一解.

证明：由于

因此，式(5)可表示为

(6)

对式(6)进行整理，得

下面，对式(5)进行误差估计.首先，定义1个网格函数空间

对于任意νn,wn∈V,引入下列内积和范数

在进行误差估计时，需要借助于以下结论.

(7)

其中，Bn-m是式(6)中定义的系数.

证明：用数学归纳法来证明.当J=2时，左边=B1(V0-V1)=右边，此时式(7)成立.现假设当J=Q时，式(7)成立，即

则当J=Q+1时，结合上面的假设，有

从而可得式(7)成立.

‖un-Un‖≤C(τ+h4),

其中，C是与离散参数无关的正的常数，且在不同的地方可能具有不同的值.

证明：式(4)和式(6)相减，得

{Bk}是正的单调递减序列，且函数g满足利普希茨条件，则对等式右边使用柯西-施瓦茨不等式，可得

(1+2B0) ‖en‖2-‖en-1‖2≤(B0-B1)(‖en-1‖2+‖en‖2)+

其中，L是利普希茨常数.因为

恒成立，则有

经过简单的代数运算，上式可变为

对上面不等式中的指标n从1～J进行求和，J为正整数，可得

根据式(7)和‖e0‖=0，可知

是负的，另外，注意到S是正定矩阵，则存在可逆矩阵H,使得S=HTH,从而上式可变为

由引理5可知，上式右端第1项是负的，从而上式变为

由引理3中{b(α1-1)n-m}的定义可知，

因而，上式可转化为

根据Gronwall不等式，可推得

‖eQ‖≤C(τ+h4).

证毕.

3 数值实验

该部分将给出2个数值算例来展现式(5)及其快速算法的计算效率，并验证上述误差估计的正确性.使用MATLAB 2012a来编程计算，电脑配置为：Intel(R) Core(TM) i7-8750 H CPU @ 2.20 GHz，内存为16 GB.在具体数值实验中，将使用max‖eJ‖作为计算的误差.

算例1考虑以下问题

u(x,0)=sin2(πx),ut(x,0)=0,u(0,t)=u(1,t)=ux(0,t)=ux(1,t)=0,0

其中,T=1.令上述问题的精确解为u(x,t)=(tα1+1+1)sin2(πx).显然，该解在时间方向的二阶导数连续，且满足定理2的光滑性假设.此外，取非线性项为g(u)=u2，则

首先，计算式(5)及其快速算法在时间方向的误差和数值收敛阶、平均CPU时间结果见表1.具体地，对α1,α2,α3取2组不同的组合，取一个足够小的空间步长h=0.001，这样可以忽略空间方向的计算误差.把时间步长τ从0.02开始逐步减半，计算时间方向的误差.从表1可以发现，式(5)及其快速算法在时间方向的收敛阶为1，符合定理2中的结论.此外，式(5)的快速算法比式(5)本身具有更高的计算效率，特别是当时间步长较小的时候，快速算法优势更大.

表1 在算例1中，式(5)及其快速计算在时间方向的误差、数值收敛阶和平均CPU时间Tab.1 Errors,temporal numerical convergence orders and average CPU times of scheme (5) for Example 1

其次，来计算空间方向的误差和数值收敛阶.由于空间方向具有高阶精度，因此时间方向的步长需要非常小，故取τ=1/2 000 000.显然，这个步长太小了，用快速算法来完成表2耗时约4.5 h.如果直接使用式(5)来计算，大概需要几星期，再次说明了快速算法在时间步长较小时的计算优势.另一方面，从表2可知：构造的算法在空间方向具有四阶精度，支持了定理2中的理论分析结果.

表2 算例1中，式(5)的快速计算在空间方向的误差和数值收敛阶Tab.2 Errors and spatial numerical convergence orders of the fast implementation of Scheme (5) for Example 1

算例2考虑以下问题

u(x,0)=0,ut(x,0)=0,u(0,t)=u(1,t)=ux(0,t)=ux(1,t)=0,0

其中,T=1.令该问题的精确解为u(x,t)=tα3sin2(πx),g(u)=u3，则

8π4tα3cos(2πx)-t3α3sin6(πx).

不难发现，该精确解在时间方向具有初始奇异性，但仍满足定理2中的光滑性假设，因此式(5)同样适用该问题.从表3可知：式(5)及其快速算法在时间方向的收敛阶为1，且快速算法的计算效率要明显高于式(5)本身.由于空间收敛阶的数值结果与表2较类似，这里不再列出.

表3 算例2中，格式(5)及其快速计算在时间方向的误差、数值收敛阶和平均CPU时间Tab.3 Errors,temporal numerical convergence orders and average CPU times of Scheme (5) and its fast implementation for Example 2