一类Caputo-Fabrizio型分数阶微分方程的三次B样条方法

胡行华, 蔡俊迎

(辽宁工程技术大学 理学院,辽宁 阜新 123000)

0 引 言

由于分数阶微分方程的非局部性可用来描述不同物质的“记忆性”和“遗传性”等物理特征,因此,其在众多科学和工程领域中具有重要的应用价值．分数阶导数有多种定义,较为常用的是Caputo型定义,作为弱奇异的分数阶导数定义,Caputo型定义比其他经典定义更适合于分数阶微分方程的描述,在进行Laplace变换时,其物理意义更加明确,因此,众多学者对Caputo型分数阶微分方程进行了广泛研究[1-2],但缺点是其定义中存在一个奇异核．2015年,Caputo和Fabrizio[3-4]提出了一个新的分数阶导数定义,即Caputo-Fabrizio型定义,其为指数函数与一阶导数的卷积,不存在奇异项,修正后的定义形式[5]如下:

(1)

式中,α∈(0,1),f(t)∈C[0,T]为线性光滑函数,且问题(1)存在唯一解u(t)．

然而,对于一般的右端项,通常很难获得精确解,因此在该定义下,方程数值方法的研究显得很有必要．许多学者对此展开了研究,2017年,Owolabi等[6]提出了线性和非线性具有Caputo-Fabrizio导数的分数阶微分方程的三步Adams-Bashforth格式,该格式具有条件稳定性．2020年,Cao等[14]针对Caputo-Fabrizio型分数阶常微分方程,构造了一种改进block-by-block算法．Guo等[15]基于Lagrange多项式,提出了求解Caputo-Fabrizio型分数阶微分方程的有限差分方法．2021年,Al-Smadi等[16]提出了一个Caputo-Fabrizio型非线性分数阶Abel微分方程的再生核方法,该方法具有较好的稳定性．Douaifia等[17]基于Newton插值提出了一种适用于Caputo-Fabrizio型分数阶微分方程的预测-校正数值格式．何广婷[18]提出了一种基于L2方案和递推关系的快速高阶数值方法,求解Caputo-Fabrizio型的分数阶常微分方程,该算法大大降低了存储容量和总计算成本．

通过现有文献发现,这些数值方法还存在一些不足之处:求解效率较低,求解精度有待进一步提高等．众所周知,样条插值函数非常接近被插值函数[19],而三次B样条函数具有计算简单、稳定性、收敛性有保证且易于在计算机上实现等优点,可以克服现有方法的缺点．本文基于分数阶微积分基本定理和三次B样条函数构造Caputo-Fabrizio型分数阶微分方程的数值格式．并对所构造的三次B样条方法的误差进行估计,对收敛性和稳定性进行理论证明．数值实验表明,三次B样条方法具有一定的可行性和有效性,在计算精度和计算效率上具有明显优势．

1 基 础 知 识

1.1 Caputo-Fabrizio分数阶导数的基本定义

下面,给出Caputo分数阶导数具体的定义形式．

定义1[20]函数u(t)的α阶Caputo分数阶导数定义为

(2)

其中,u(t)∈H1(0,b),b>0,α∈(0,1),Γ(·)为Gamma函数．

2015年,Caputo和Fabrizio[3]提出了一个具有光滑核的分数阶导数新定义,其定义如下．

定义2[3]令u(t)∈H1(0,b),b>0,并且α∈(0,1),通过用函数exp(-α(t-τ)/(1-α))替换核(t-τ)-α,用M(α)/(1-α)替换1/Γ(1-α),函数u(t)的α阶Caputo-Fabrizio分数阶导数定义为

(3)

其中,M(α)是一个标准化函数．与式(2)相比,新定义在t=τ时无奇异核．

Losada和Nieto[5]对Caputo-Fabrizio分数阶导数进行了修正,首先有以下定义．

定义3[5]假设u(t)∈H1(0,b),α∈(0,1),将Caputo-Fabrizio分数阶导数定义为

(4)

Losada和Nieto[5]提供了M(α)的一个显式公式:M(α)=2/(2-α)．将M(α)的显式公式代入式(4),即得式(3)修正后的公式[5]:

(5)

接下来,对于一个一般的α,α∈(0,1)阶Caputo-Fabrizio型分数阶微分方程(1),利用分数阶微积分基本定理[5]可将其转换为

(6)

由上式易得,当且仅当f(0)=0时,式(1)中的u(t)满足初始条件u(0)=u0．鉴于此,求解式(1)的数值解即转化为逼近式(6)右端项中积分的问题．

1.2 m次B样条插值函数

将区间[0,T]剖分成等距的N个小区间,对于给定的任一整数N>0,其小区间的步长为h=T/N．对于j=0,1,…,N,tj=jh,并且0=t0

(7)

(8)

(9)

在m次B样条插值函数中,相比于低次和高次的B样条函数,三次B样条函数具有需要较少的信息、精度高、计算简单和易于实现编程的优势[23]．因此,本文使用三次B样条函数求解一类Caputo-Fabrizio型分数阶微分方程．

2 三次B样条方法数值格式的构造

首先,考虑右端项与u(t)无关的线性初值问题(1),即

(10)

使用三次B样条插值函数S3(τ)来近似代替式(6)中积分里的函数f(τ),可得

[βj+1(h3+3h2(τ-tj-1)+3h(τ-tj-1)2-3(τ-tj-1)3)+βj+2(τ-tj-1)3]dτ+

(1-α)f(tk)=

其中,1≤k≤N．由此得出S3公式,并使用符号uS3(tk)表示,即

(11)

其中,β0,β1,…,βN+2为常系数,求出β的N+3个系数则得到uS3(tk)．

根据三次B样条插值函数理论,uS3(tk)满足N+1个插值条件:在节点tj上的函数值yj=f(tj)(j=0,1,…,N),且S3(tj)=yj(j=0,1,…,N)成立．此外,还需要2个边界条件才能求出β的N+3个系数,边界条件有多种类型[22]．为便于求解,本文选择普适性较强的固定边界条件S′3(t0)=f′(t0),S′3(tN)=f′(tN)．由插值条件和边界条件可得系数β满足线性方程组Mβ=f,其中f=[f′(t0),f(t0),f(t1),…,f(tN),f′(tN)]T,β=[β0,β1,…,βN+2]T,且矩阵M为

矩阵M是对角占优矩阵,因此系数β是唯一确定的,使用追赶法[22]即可获得．其他边界类型的三次B样条函数也可以类似应用,且β同样易获得．

接下来,考虑右端项与u(t)有关的线性初值问题[6]:

(12)

其中,g(t)是一个已知函数．根据式(12),有

(1-α)[u(tk)+g(tk)]=

(13)

(14)

根据上式,需要u′(0)和u′(tN)的值,它们分别由以下四阶正向差分公式和四阶反向差分公式[24]来近似,即

(15)

由式(13)—(15)即可求得u(tk)的近似数值uk．

3 三次B样条方法的误差估计、收敛性和稳定性分析

3.1 三次B样条方法的误差估计

首先,在对三次B样条方法进行误差估计之前,给出以下定理．

定理1[25]设函数f∈C4[0,T],则函数f与三次样条插值函数S3之间的距离为

定理2 若f∈C4[0,tk],且R(tk)=u(tk)-uS3(tk),则有

(16)

证明对于任意的1≤k≤N,根据定理1可得

|R(tk)|=|u(tk)-uS3(tk)|=

定理2证毕．

3.2 三次B样条方法的收敛性分析

下面给出三次B样条方法的收敛性分析．

定理3 三次B样条方法在区间[0,T]上是一致收敛的,即当h→0时,‖R3(tk)‖∞→0．

证明根据定理2中三次B样条方法的误差估计,可得

其中,f∈C4[0,T],当h→0时,有

因此,此数值格式在区间[0,T]上是一致收敛的．定理3证毕．

3.3 三次B样条方法的稳定性分析

接下来对三次B样条方法进行稳定性分析．

定理4 三次B样条方法在区间[0,T]上是无条件稳定的．

因此,三次B样条方法是无条件稳定的．定理4证毕．

4 数 值 算 例

其中,N为各数值方法在[0,1]区间分割的份数,也为待求未知量的个数,代表以上各数值方法在数值求解的计算量大小,且3种数值方法的带状矩阵皆采用Gauss消元法来计算．所有程序均在CPU为Inter Core i5 Processor 2.30 GHz、内存为8 GB的笔记本电脑环境下运行,利用MATLAB 2018a平台实现．

例1 本文考虑具有两个不同右端项的初值问题[14,18]:

情况1

(17)

情况2

(18)

两种情况的初值均为u0=0,精确解均为u(t)=t3,并且右端项中G1(t)和G2(t)分别为

G2(t)=G1(t)-t3．

可见,情况1与情况2均属于线性初值问题,情况2的右端项与u(t)有关,且G1(0)=0,G2(0)+u(0)=0．

下面,使用本文提出的三次B样条方法分别求解情况1和情况2的初值问题．

求解情况1 已知插值条件S3(tj)=G1(tj)(j=0,1,…,N)和固定边界条件S′3(0)=G′1(0),S′3(1)=

G′1(1)成立,使用三次B样条方法求解不同的分割数N(N=4,8,16,32,64,128)与不同α(α=0.3,0.7)阶的Caputo-Fabrizio型分数阶微分方程数值解,并与文献[14]中改进的block-by-block算法和文献[18]中的快速高阶数值方法进行对比．当α=0.3时,3种数值方法的最大误差和收敛阶对比如表1所示,3种数值方法的最大误差比较如图1所示．当α=0.7时,3种数值方法的最大误差和收敛阶对比如表2所示,3种数值方法的最大误差比较如图2所示．

表1 当α=0.3时,3种数值方法的最大误差和收敛阶(情况1)Table 1 Maximum errors and convergence orders of 3 numerical methods for α=0.3 (case 1)

表2 当α=0.7时,3种数值方法的最大误差和收敛阶(情况1)Table 2 Maximum errors and convergence orders of 3 numerical methods for α=0.7 (case 1)

图1 当α=0.3时,3种数值方法的最大误差(情况1) 图2 当α=0.7时,3种数值方法的最大误差(情况1)Fig.1 Maximum errors of the 3 numerical methods for α=0.3 (case 1) Fig.2 Maximum errors of the 3 numerical methods for α=0.7 (case 1)

由表1和表2可知,在不同的分数阶次下,与改进的block-by-block算法相比,本文方法的误差明显更小,数值逼近效果更佳,且收敛阶较高．与快速高阶数值方法相比,本文方法的误差明显更小,数值逼近效果更佳,收敛阶相当．此外,本文数值方法的CPU时间较短,具有可观的计算效率．显然,采用三次B样条方法求解Caputo-Fabrizio型分数阶微分方程比其他两种方法更加有效．

由图1和图2可知,改进的block-by-block算法和快速高阶数值方法在各分割数下的最大误差相差不多,与以上两种数值方法相比,本文方法的最大误差明显更小,数值逼近效果更佳．

下面,固定分割数N=100,使用三次B样条方法分别求解情况1中不同阶次α(α=0.2,0.4,0.6,0.8) 的初值问题,其各节点处的绝对误差结果如图3所示．固定分割数N=10,本文方法在阶次α取不同值时所得各节点的数值解如图4所示．

由图3可知,不同的阶次α导致本文方法的绝对误差有所不同,阶次α越小,各节点处的绝对误差越小,数值逼近的效果越佳．由图4可知,当阶次α取不同值时,各节点处的数值解均保持平稳状态,本文方法在t∈[0,1]时具有良好的稳定性．

表3 当α=0.3时,3种数值方法的最大误差和收敛阶(情况2)Table 3 Maximum errors and convergence orders of 3 numerical methods for α=0.3 (case 2)

表4 当α=0.7时,3种数值方法的最大误差和收敛阶(情况2)Table 4 Maximum errors and convergence orders of 3 numerical methods for α=0.7 (case 2)

图5 当α=0.3时,3种数值方法的最大误差(情况2) 图6 当α=0.7时,3种数值方法的最大误差(情况2)Fig.5 Maximum errors of the 3 numerical methods for α=0.3 (case 2) Fig.6 Maximum errors of the 3 numerical methods for α=0.7 (case 2)

由表3和表4可知,在不同的分数阶次下,与改进的block-by-block算法相比,本文方法的误差明显更小,数值逼近效果更佳,且收敛阶较高．与快速高阶数值方法相比,本文方法的误差明显更小,数值逼近效果更佳,收敛阶相当．此外,本文数值方法的CPU时间较短,具有可观的计算效率．显然,三次B样条方法在求解Caputo-Fabrizio型分数阶微分方程时比其他两种方法更加有效．

由图5和图6可知,改进的block-by-block算法和快速高阶数值方法在各分割数下的最大误差相差不多,与以上两种数值方法相比,本文方法的最大误差明显更小,数值逼近效果更佳．

接下来,固定分割数N=100,使用三次B样条方法分别求解情况2中不同阶次α(α=0.2,0.4,0.6,0.8) 的初值问题,其各节点处的绝对误差结果如图7所示．固定分割数N=10,本文方法在阶次α取不同值时所得各节点的数值解如图8所示．

图7 不同α值时各节点的绝对误差(情况2) 图8 不同α值时各节点的数值解(情况2)Fig.7 Absolute errors of each node with different α values (case 2) Fig.8 Numerical solutions of each node with different α values (case 2)

由图7可知,不同的阶次α导致本文方法的绝对误差有所不同,阶次α越小,各节点处的绝对误差越小,数值逼近的效果越佳．由图8可知,当阶次α取不同值时,各节点处的数值解均保持平稳状态,本文方法在t∈[0,1]时具有良好的稳定性．

例2 下面将验证三次B样条方法的稳定性,考虑如下的初值问题[14,18]:

(19)

其中,精确解为u(t)=sint,右端项满足f(0)=0．

已知插值条件S3(tj)=f(tj)(j=0,1,…,N)和固定边界条件S′3(0)=f′(0),S′3(1)=f′(1)成立,固定分割数N=10 000,使用3种数值方法分别求解不同的α(α=0.3,0.5,0.7)阶的Caputo-Fabrizio型分数阶微分方程数值解,一直计算到T=1 000．当α=0.3,0.5,0.7时,3种数值方法在不同节点处的绝对误差(|ε|=|u(tk)-uk|)和相对误差(ε′=|ε|/u(tk))结果分别如表5、表6和表7所示,3种数值方法的相对误差对比如图9、图10和图11所示．

表5 α=0.3,N=10 000时,3种数值方法在不同节点处的绝对误差和相对误差Table 5 The absolute and relative errors of the 3 numerical methods at different nodes for α=0.3,N=10 000

表6 α=0.5,N=10 000时,3种数值方法在不同节点处的绝对误差和相对误差Table 6 The absolute and relative errors of the 3 numerical methods at different nodes for α=0.5,N=10 000

表7 α=0.7,N=10 000时,3种数值方法在不同节点处的绝对误差和相对误差Table 7 The absolute and relative errors of the 3 numerical methods at different nodes for α=0.7,N=10 000

图9 当α=0.3时,3种数值方法的相对误差 图10 当α=0.5时,3种数值方法的相对误差Fig.9 Relative errors of 3 numerical methods for α=0.3 Fig.10 Relative errors of 3 numerical methods for α=0.5

图11 当α=0.7时,3种数值方法的相对误差Fig.11 Relative errors of 3 numerical methods for α=0.7

由表5所示,与改进的block-by-block算法和快速高阶数值方法相比,三次B样条方法在不同节点处的绝对误差和相对误差较小,经过长时间的计算,不同节点上的数值解均可以很好地逼近精确解,并且本文方法的CPU时间较短,计算效率可观．由图9所示,3种数值方法的相对误差在个别节点处均有起伏,但与其他两种方法相比,本文方法的平稳状态更佳．此外,3种数值方法相对误差的方差分别为1.074 4×10-10,7.603 0×10-12和5.157 2×10-13．显然,三次B样条方法在求解Caputo-Fabrizio型分数阶微分方程时比其他两种方法更加稳定．

由表6所示,与改进的block-by-block算法和快速高阶数值方法相比,三次B样条方法在不同节点处的绝对误差和相对误差较小,经过长时间的计算,不同节点上的数值解均可以很好地逼近精确解,并且本文方法的CPU时间较短,计算效率可观．由图10所示,3种数值方法的相对误差在个别节点处均有起伏,但与其他两种方法相比,本文方法的平稳状态更佳．此外,3种数值方法相对误差的方差分别为6.410 8×10-10,3.654 9×10-11和2.853 2×10-13．显然,三次B样条方法在求解Caputo-Fabrizio型分数阶微分方程时比其他两种方法更加稳定．

由表7所示,与改进的block-by-block算法和快速高阶数值方法相比,三次B样条方法在不同节点处的绝对误差和相对误差较小,经过长时间的计算,不同节点上的数值解均可以很好地逼近精确解,并且本文方法的CPU时间较短,计算效率可观．

图11(a)为当α=0.7时3种数值方法的相对误差对比,图11(b)为快速高阶数值方法和三次B样条方法的相对误差对比．由图11所示,3种数值方法在个别节点处均有起伏,但与其他两种方法相比,本文方法的平稳状态更佳．此外,3种数值方法相对误差的方差分别为4.188 4×10-9,1.514 5×10-10和1.051 7×10-10．因此,三次B样条方法在求解Caputo-Fabrizio型分数阶微分方程时比其他两种方法更加稳定．

5 结 论

本文将三次B样条方法应用于Caputo-Fabrizio定义下的分数阶微分方程的数值求解．基于分数阶微积分基本定理和三次B样条函数,构造了求解α阶线性Caputo-Fabrizio型分数阶微分方程数值解的三次B样条方法．对所构造的数值方法进行了误差估计,并对其收敛性和稳定性进行了理论证明．实验结果表明:三次B样条方法具有一定的有效性,且具有4阶收敛阶和良好的稳定性;在求解线性初值问题时,与改进的block-by-block算法和快速高阶数值方法相比,三次B样条方法明显具有较高的精度和较快的收敛速度,且计算复杂度低,计算成本小．综上,三次B样条方法在求解Caputo-Fabrizio型分数阶微分方程时具有明显优势,为一类α阶Caputo-Fabrizio型分数阶微分方程的数值求解提供了新的选择．