Riesz空间分布阶的分数阶扩散方程的数值模拟

陈景华，陈雪娟

(1.集美大学理学院，福建 厦门 361021；2.集美大学理学院数字福建大数据建模与智能计算研究所，福建 厦门 361021)

0 引言

分数阶微分方程是广义的非整数阶的微分方程。由于分数阶算子的非局部性，分数阶微分方程非常适合用来描述现实世界中具有记忆及遗传性质的材料[1-3]。 然而，还有很多问题不能采用分数阶系统来解决，如非均质多孔和裂隙介质中的溶质迁移试验[4-5]。由于复杂流体流动的力学性质会随时空尺度而发生一定的变化，因此污染物等溶质扩散将产生所谓的多尺度效应。 近年来，研究者们[6-9]采用Riesz分布阶的分数阶扩散方程(其中分布阶导数在数学上可以看成是对分数阶导数在给定范围内关于其阶数的一个积分，即分数阶导数是在一个单位区间上的积分)来描述这种反常扩散。Diethelm等[10]研究了Caputo分布阶常微分方程的数值解并进行了收敛性分析；Atanackovic等[11]证明了时间分布阶柯西问题解存在性，并求出解析解；Luchko等[12]考虑了一个开区间上的分布阶时间分数阶扩散方程，并证明了这个问题的解的唯一性和存在性；Jiao等[9]给出了分布阶导数在信号处理中的分布阶的滤波器和控制系统中的最优分布阶阻尼两个应用。 Podlubny等[13]用矩阵方法离散分布阶导数和积分,给出数值例子验证算法的精确性，但是没有具体的数值理论分析；Gorenflo等[14]研究了一种分布阶时间分数阶扩散波动方程，并利用傅里叶变换和拉普拉斯变换技术得到了此问题的基本解。现有文献关于分布阶微分方程数值方法的研究还比较少。2014年，Katsikadelis[15]研究了Caputo分布阶线性与非线性常微分方程的数值解，但没有进行理论分析。2018年，Semarya等[16]运用移位的分数阶积分切比雪夫算子矩阵与修正的Picard数值方法求解分布阶线性分数阶常微分方程，但没有进行理论分析。分布阶微分算子作为分数阶微分算子的发展形式，已经越来越受到研究人员的重视。相对于固定阶导数，这类算子的应用领域更为广阔，但由于算子定义的复杂性，相应的数值算法更为复杂，并不能简单照搬分数阶的数值算法。本文考虑用有效的数值方法求解空间分布阶的分数阶扩散方程。

考虑以下Riesz空间分布阶的分数阶扩散方程[17-18]：

(1)

边界条件：

u(xL,t)=0,u(xR,t)=0,t∈(0,T],

(2)

初始条件：

u(x,0)=φ(x),x∈[xL,xR],

(3)

这里P(α)是非负加权函数[10]，满足：

(4)

(5)

这里cα=-1/(2cos(πα/2)),且

(7)

1 积分项的离散

首先，把积分区间划分成若干个子区间：1=α0<α1<…<αS=2,S∈N，并记

Δαi=αi-αi-1=1/S=σ,αi-1/2=(αi-1+αi)/2,i=1,2,…,S。

(8)

(9)

(10)

则Riesz空间分布阶的分数阶扩散方程可以离散成具有多项分数阶导数的微分方程：

∂u(x,t)/∂t=Dxu+f(x,t)+O(σ4)

(11)

具有边界和初始条件(2)～(3)。

2 数值离散格式

首先对空间和时间进行离散。取正整数N,设时间步长τ=T/N,tn=nτ,0≤n≤N。记时间区域Ωτ={tn|0≤n≤N}。取正整数M，并设空间步长h=(xR-xL)/M,xi=xL+ih,0≤i≤M。空间区域Ωh={xi|0≤i≤M}。空间平均差分算子Aα定义为：

平均算子A定义为：

(12)

(13)

(14)

(15)

(16)

其中,

(17)

(18)

且

(19)

(20)

及差分算子δxu：

(21)

(22)

其中,

(23)

由方程(9)～(16)得到：

A(Dxu(xi,t))=δxu(xi,t)+O(h4)。

(24)

考虑方程(1)在点(xi,t)处的方程为：

ut(xi,t)=Dxu(xi,t)+f(xi,t)+O(σ4),

(25)

在方程(25)两端取平均算子A，则有：

Aut(xi,t)=ADxu(xi,t)+Af(xi,t)+O(σ4),xi∈Ωh。

(26)

由式(24)可得：

Aut(xi,t)=δxu(xi,t)+Af(xi,t)+O(σ4+h4),xi∈Ωh。

(27)

下面对时间离散。对于方程(27)，在t=tn和t=tn+1处取平均，且利用泰勒展开，可得到：

(28)

(29)

(30)

得到以下四阶拟紧差分格式：

(31)

(32)

为了逼近Riesz空间分布阶的分数阶扩散方程，将方程(31)改成：

(33)

线性系统(33)的系数矩阵A=(as,t)写成下式：

(34)

则有以下的定理1。

定理1 差分格式(31)～(32)是唯一可解的。

3 稳定性和收敛性的证明

为了证明差分格式的稳定性和收敛性，有以下引理。

引理1A是自相似的，即对任意的u,v∈Vh,成立(Au,v)=(u,Av)。

证明

(35)

即证得。

证明

(36)

且

(37)

可得

(38)

即证得。

引理3 对任意的v∈Vh, 成立(δxv,v)≤0。

证明

(39)

使用引理3，容易得到以下定理2，其证明类似文献[18]的定理2。

(40)

(41)

则有

(42)

由定理2可以直接推出如下定理3。

定理3 差分格式(31)～(32)按初值和右端项f是无条件稳定的。

现在考虑差分格式 (31)～(32)的收敛性。

定理4 差分格式(31)～(32)是收敛的，且收敛阶为O(τ2+σ4+h4)。

(43)

将式(28)和式(29)分别减去式(31)和式(32)，可得到误差方程

(44)

(45)

根据定理2得到：

(46)

证毕。

4 数值例子

例1 考虑以下Riesz空间分布阶扩散方程：

(47)

边界条件：

u(0,t)=0,u(1,t)=0,t∈(0,1],

(48)

初始条件：

u(x,0)=x4(1-x)4,x∈[0,1]。

(49)

这里，

P(α)=-2Γ(9-α)cos(πα/2)，

(50)

且

(51)

方程(47)～(49)的精确解为u(x,t)=e-tx4(1-x)4。

图1显示了T=1.0时刻数值解与精确解的比较,取σ=1/1 000、τ=h2=1/400二者非常吻合。表 1 显示σ=1/1 000、τ=h2、T=1.0时数值解与精确解的最大误差及收敛。误差率接近:Rε=lg(ε(h1)/ε(h2))/lg 2≈4，这与定理4中差分格式的收敛阶是O(τ2+σ4+h4)是一致的。从图1和表1可以看出，数值结果与理论分析结果是一致的。

表1 T=1.0时刻数值方法的最大误差和收敛阶(σ=1/1 000，τ=h2)

5 结论

本文发展了一个Riesz空间分布阶的分数阶扩散方程的数值方法。将分布阶方程离散为含有多项分数阶导数的微分方程，利用有限差分法对得到的微分方程进行数值求解。证明差分格式是稳定的和无条件收敛的。此外，本文的数值方法的收敛阶为O(τ2+σ4+h4)。最后，给出了数值例子来证明本文数值方法与理论结果是一致的。这种方法和分析技术可用于求解和分析其他类型的分数阶偏微分方程。