两项时间分数阶扩散方程的一种数值解法

余跃玉

(四川文理学院 数学学院，四川 达州 635000)

0 引言

大量研究发现，分数阶微分方程可以更好地描述物质的记忆和遗传性，反映复杂的物理过程及动力学行为.比如，它们能够很好地模拟地下水运输、热传导、空气污染、地震等复杂动力系统[1-4]，在物理、金融、工程、医学、社会科学等领域也有相当广泛的运用[5-6].虽然分数阶微分方程的物理意义比整数阶方程更直观，模型更简单，但分数阶导数本身的非局部性，给理论分析和数值求解都带来了很多的困难.通常，其解析解带有极其复杂的特殊函数，比如Wright函数和Mittag-Leffler函数等，或者根本就难以获得，所以寻求精度较高的数值解就显得尤为重要[7-9].近年来，受到广大学者关注的分数阶模型主要有分数阶水波模型、分数阶反应扩散模型、分数阶Cable模型、分数阶扩散波动模型等[10-14].

分数阶扩散模型在物理领域常用来模拟异常扩散问题，其粒子的扩散不同于布朗运动的无规则性，它是不同物质间的一种相互扩散，其分子因热运动从高浓度区域向低浓度区域迁移.有限差分法是分数阶扩散方程的一种重要的数值解法，目前已有大量的研究成果，如Yao等[14]考虑了Neumann边界条件下空间四阶时间分数阶慢扩散问题，Cui[15]构造了第一类Dirichlet边界条件下的分数阶慢扩散方程.为了更精准地描述复杂环境下的扩散问题，经典扩散方程或单项分数阶扩散方程往往是不够的，我们需要引入多项时间分数阶导数[16]刻画多种复杂材料中发生的介质不规律扩散现象.

考虑两项时间分数阶扩散方程的初边值问题

(1)

这里的Г(·)为Gamma函数.

1 差分格式的建立及其可解性

对γ=α和γ=β阶Caputo型时间分数阶导数做如下近似：

对空间二阶导数采用中心差分，可得一种隐式格式：

令

显然，pk>0,k=0,1,2,…,j.

将pk代入上述差分格式，并化简可得：

当j=0时，

当j>0时，

则(6)式的矩阵表示形式为

AU1=p0U0+F1,

(8)

(7)式的矩阵表示形式为

由于矩阵A为严格对角占优矩阵，且aii>0,aij≤0(i≠j)，所以矩阵A的行列式|A|≠0，因此线性方程组(8)和(9)都有唯一解，由此可得：

定理1两项时间分数阶扩散方程(1)的差分格式(6)和(7)是唯一可解的.

2 差分格式的稳定性

当j=0时，

当j>0时，

其中，i=1,2,…,m-1;j=0,1,2,…,n-1.

利用引理1、引理2及数学归纳法可得：

定理2||Ej||∞≤||E0||∞,j=1,2,3,….

所以，||E1||∞≤||E0||∞,j=1,2,3,….

所以，||Ej+1||∞≤||E0||∞,j=1,2,3,….】

定理3多项时间分数阶扩散波动方程(1)的差分格式(6)和(7)是无条件稳定的.

3 数值格式的收敛性

当j=0时，

(14)

当j>0时，

利用泰勒定理及积分中值定理，易得此格式的截断误差为

由于矩阵A为严格对角占优矩阵，|A|≠0，即A可逆，所以(16)和(17)式又可化为

所以

假设||ek||2≤C(τμ+h2),k

其中C为正常数.

4 数值实验

令两项时间分数阶扩散方程的初边值问题(1)中φ(x)=x2-x3,ψ(x)=0,λ=0.8,L=T=1，源项

α=0.7,β=0.6，则该方程的精确解为

u(x,t)=(t2+1)(x2-x3).

取时间和空间步长分别为τ=0.1,h=0.1;τ=0.05,h=0.05;τ=0.025,h=0.025;τ=0.01,h=0.01，分别计算这4种情况下两项时间分数阶扩散方程在t=1时的数值解，以及该方程在t=1时的精确解，其结果见表1.从表1可见，网格剖分越细，数值解越接近精确解，这说明本文数值格式是有效的.

表1 t=1时不同步长下数值解与精确解的比较

表2给出了不同时间步长和不同空间步长在时间t分别取0.1,0.2,…,1时的L∞-误差，结果表明，数值格式的收敛性为o(τmin{2-α,2-β}+h2)，实验结果与理论分析高度一致.

表2 t时刻每个时间层上的L∞-误差