数值求解一维时间—空间分数阶扩散方程

李玉山（甘肃政法学院 网络空间安全学院，甘肃 兰州，730070）

1 引 言

进入21世纪以来,随着分数阶微积分在反常扩散、非牛顿流体力学、多孔介质力学、粘弹性材料力学、地球物理、生物医学工程、经济学等诸多领域的建模及广泛应用[1-6],分数阶微积分的理论及应用开始受到越来越多的关注,相关的结果也陆续出现,参见Podlubny, Kilbas等人的专著[7、8]以及中文专著[9、10]。分数阶反常扩散方程是由经典扩散方程推广而得到的。相比整数阶扩散方程，分数阶扩散方程更适合描述反常扩散现象，因为分数阶导数能够描述具有记忆和遗传性质的非均匀物质[12、13]。

时间-空间分数阶扩散方程

02α＜＜被运用在反常扩散模型中，通常,若考虑时间相关性或扩散的记忆性，就得到时间分数阶扩散方程；若考虑空间相关性或非局域性，则得到空间分数阶扩散方程；如果既考虑时间相关性，又考虑空间整体性，则得到所谓的时间-空间分数阶扩散方程。对于时间分数阶扩散方程的数值求解，已经有了很多的研究成果[15-20]。对于时间-空间扩散方程的数值求解，相对来说，研究的比较少，刘发旺等[21,22]给出了空间分数阶拉普拉斯算子的处理办法，但是时间项依然是整数阶的。

本文将探讨一维时间-空间扩散方程的数值求解，利用矩阵转换技术处理空间分数阶拉普拉斯算子和有限差分法处理时间分数阶项，得到代数方程组，同时利用分离变量法得到解析解表达式，并且给出数值算例作比较。

2 问题的提出

考虑如下的齐次Dirichlet边界条件的时间-空间分数阶扩散方程的初边值问题：

其中Γ为伽马函数。

为阶数为α(0 ＜α≤2)的空间分数阶拉普拉斯算子，由拉普拉斯算子-（）Δ的谱分解来定义，定义如下：

设

那么对于任意的fγ∈F，定义分数阶拉普拉斯算子为：

本文考虑给定f(x),p(t),g(x)和方程(1)-(3)，数值求解u(x,t)。

3 数值方法

利用矩阵转换技巧求解齐次Dirichlet边界条件的时间-空间分数阶扩散方程的初边值问题。先考虑标准的一维扩散方程：

取 正 整 数 1N＞ ,空 间 步 长 1/hN= ,ixih=(0iN≤≤)，引进标准的中心差分离散二阶空间导数：

利用文献[21,22]的结果，方程（14）-（16）可以离散为如下矩阵形式

取正整数1M＞,时间步长(0jM≤≤)，用隐式差分格式离散时间分数阶项[20]：

将（19）写为矩阵形式，代入（17）得：

利用Matlab软件求解代数方程组（20），即可得到问题（1）-（3）的近似解。

为了验证本文的方法，利用分离变量法，给出问题（1）-（3）的解析解如下：

其中,Eβ β为广义的Mittag-Leffler函数，定义如下：

在计算（21）的过程中，只计算无穷级数的前50项来近似得到u(x,t)，利用Matlab程序[23]计算Mittag-Leffler函数时取精度为10-6。

4 数值例子

例:取T=1,初始值u(x,0)=g(x)=x(1-x),f(x)=sin(x),p(t)=t,u(0,t)=u(1,t)=0,γ=1,M=N=100。

图1 α =1.2, β =0.3时的数值解、解析解以及误差图

图2 α =1.8, β =0.3时的数值解、解析解以及误差图

图3 α =1.2, β =0.7时的数值解、解析解以及误差图

图4 α =1.8, β =0.7时的数值解、解析解以及误差图

表1 数值解和精确解CPU占用时间对比 单位（s）

图1-图4给出了例1中α=1.2,1.8和β=0.3,0.7的数值结果，并与解析表达式（3.15）做了比较，可以看出，本文中的方法很有效，并且具有较高的精度。由表1可以看出，随着时间和空间离散点的增加，解析解的CPU占用时间远高于本文提出的数值方法的CPU占用时间，这是由于利用解析表达式（3.15）计算时，牵涉到计算两个无穷级数、数值积分以及Mittag-Leffler函数，在时效上有很大的缺陷，而本文提出的方法避免了这一缺陷，只需要在每一时间层解一个线性方程组，有较高的精度，且占用极少的CPU时间。

5 结 语

本文研究了一维具有齐次Dirichlet边界条件的时间-空间分数阶扩散方程的数值求解，采用矩阵转换技巧求解且给出了数值例子，并和解析解做了比较。数值例子表明，1）本文的方法是很有效的，具有较高的计算精度，可以处理实际问题；2）矩阵转换技巧求解只需要求解线性方程组，而利用解析表达式求解计算比较复杂，且占用CPU时间较长；3）利用矩阵转换技巧可以推广到第二类、第三类边界条件以及非齐次边界条件的问题。