危国华

（福建广播电视大学 三明分校，福建 三明 365000）

近年来研究者们发现分数阶微分方程比经典的整数阶微分方程更为准确的描述现实生活中的各种现象，尤其分数阶扩散方程已被广泛用来描述特征方差尺度为时间的幂次的反常扩散过程。随后学者们讨论分数阶微分方程的数值解法[1-9]。本文尝试采用差分/多项式基点插值配置法处理多项时间-空间分数阶微分方程，提出半离散格式以及全离散格式，并通过数值例子验证方法的有效性。

讨论如下多项的时间-空间分数阶微分方程：

带有初边值条件：

其中d1，d2，…，dk，γ1，γ2，…，γk，κβ，κ0，κ1都是非负常数，且表示函数 u（x，t）对时间变量t的γl阶Caputo分数阶导数[10]，即表示函数u（x，t）对空间变量x的β（1＜β＜2）阶Riemann-Liouville分数阶导数[10]，即

这里Γ（·）为伽马函数。

本文假设时间-空间分数阶微分方程(1-3)的解

1 离散格式

1.1 半离散格式

首先离散时间变量，令 tn=nτ，n=0，1，…，N,τ=T/N，其中τ为时间步长。

对于 0＜γ＜1，引入记号…，在t=tn+1时间点离散各项，利用[3]

得到下面格式：

其中

以及

其中C为正的常数。

经整理，（4）式可改写为

令 un=un（x）为 u（x，tn）的数值近似，fn=fn（x）表示f（x，tn），略去误差项 τRn+1（x），可得到半离散格式：

1.2 多项式点插值

接下来介绍本文用到的离散空间变量的多项式点插值。

首先，在问题域[a，b]内生成场节点，节点可以随机生成，也可以人工加入节点。设有内部节点x1,x2,…xM-1,和边界节点 x0，xM。

对每一个场节点，都做一个包含该节点的支持域，在实际计算中，一般取节点平均间距的二到三倍大小[11]。

设第i个节点xi的支持域为Ωi，支持域Ωi内包含 mi个节点 xi1，xi2，ximi。为书写方便起见，记Di={i1,i2,…，imi}，这就意味着，当 l∈ Di时，xl∈ Ωi。

接下来，在节点的支持域内对函数u（x）利用基函数进行逼近。为了描述方便，假设计算点的支持域内包含m个场节点x1,x2,…,xm，连续函数u（x）可以近似表示为：

其中pj（x）为空间坐标变量的单项式，被称为多项式基，m是多项式基的个数，a=[a1,a2,…am]T,aj为待定常数。一般地，可取一维的m-1次完备多项式为

为了确定待定系数 aj（j=1,2,… ,m ），令（i=1,2,… ,m ），即

写成矩阵形式为

其中 Us=[u（x1）u（x2）… [u（xm）]是节点上的函数值向量，以及

由于x1,x2,…,xm互不相等，则矩阵Pm是m阶可逆方阵.由(11)，有：

将（12）代入（9）,可得

其中 ΦT（x）是形函数向量，其定义为

容易得到，形函数向量ΦT（x）的k阶(k为正整数)导数为

形函数向量ΦT（x）的β阶 (β为正实数)左侧Riemann-Liouville导数为：

于是，我们有

1.3 全离散格式

接下来，我们将对方程（8）进行空间离散。

考虑内部节点 xi（i=1,2,… ,m ），每个节点 xi对应的支持域为 Ωi，支持域 Ωi内包含 mi个节点 xi1，xi2，ximi，

为方便起见，记假设方程（8）的近似解为

将该式代入方程（8），并令 x=xi，于是，可得

边界条件为

由幂函数的分数阶导数公式[10]：

得到(14)式中形函数 φj（x）（j=1,2,… ,n）在 xi的分数阶导数表达式为：

其中 Pi，m为节点 xi的支持域｛xl，l=i1，i2，… ,in｝形成的方阵，

2 数值结果

为了避免在计算中产生Runge现象，在构造场节点时选择6个临近节点作为局部支持域时k=5，此时，基函数为5次完备多项式基。

例1考虑下列方程：

带有边值条件

以及初始条件：

其中

该方程的精确解为 u（x，t）=x2t2。

首先，将问题域[0,1]用规则节点 xi=ih（i=0，1，2，…，100）划分，采用全离散格式(14)进行求解（计算中取m=5），从图1看出，精确解和数值解吻合得很好。

图1 规则点划分的数值解与精确解Figure 1 Numerical and exact solutions for regular distributed nodes

表1中列出时间步长变化时，数值解与精确解的误差以及误差阶。

表1 采用规则点划分求解的数值解与精确解的误差以及误差阶Table 1 Error of the numerical and exact solutions for regular distributed nodes

接下来将问题域用不规则节点划分，采用全离散格式(14)进行求解。图2给出此时精确解与数值解。

图2 不规则点划分的数值解与精确解Figure 2 Numerical and exact solutions for irregular distributed nodes

表2列出针对不规则点划分空间变量，时间步长变化时，数值解与精确解的误差以及误差阶。由此看出，本节提出的数值方法仍然适用于不规则点划分。

表2 采用不规则点划分求得的数值解与精确解的误差以及误差阶Table 2 Error of the numerical and exact solutions for irregular distributed nodes

例2考虑下列时间-空间分数阶微分方程

带有边值条件

以及初始条件

其中 γ1=0.9，γ2=0.6，β=1.8，且

该问题的精确解为 u（x，t）=t2x2（1-x）2。

同样，将问题域[0，1]分别用等距节点划分以及非等距节点划分，采用全离散格式(14)进行求解。从图3和4看出，无论是规则点，还是不等距节点，精确解和数值解吻合得很好。

图3 规则点划分时的数值解与精确解Figure 3 Numerical and exact solutions for regular distributed nodes

图4 非等距节点划分时的数值解与精确解Figure 4 Numerical and exact solutions for irregular distributed nodes

表3和4中给出规则点和不规则点划分空间变量，时间步长变化时的误差以及误差阶。

表3 规则点划分时求得的数值解与精确解的误差以及误差阶Table 3 Error of the numerical and exact solutions for regular distributed nodes

表4 非等距节点划分时求得的数值解与精确解的误差以及误差阶Table 4 Error of the numerical and exact solutions for irregular distributed nodes

[1]DENTZ M,BERKOWIZ B.Transport behavior of a passive solute in continuous time random walks and multirate mass transfer[J].Water Resources Research,2003,39(5)：1111-1130.

[2]HARVEY C,GORELICK S.M.Rate-limited mass transer or macrodispersion：Which dominates plume evolution at the Macrodispersion Experiment (MADE)site [J].Water Resources Research,2000,36(3)：637-650.

[3]SUN Z Z,WU X N.A fully discrete difference scheme for a diffusion-wave system[J].Apply Numerical Mathematics,2006(56)：193–209.

[4]ZHANG Y,MEERSCHAERT M M,BAEUMER B.Particle tracking for time-fractional diffusion[J].Physical Review E,2008(78)：036705.

[5]MEERSCHAERT M M,ZHANG Y,BAEUMER B.Particle tracking for fractional diffusion with two time scales[J].Computers and Mathematics with Applications,2010(59)：1078-1086.

[6]ZHUANG P,LIU F,ANH V,et al.New solution and analytical techniques of the implicit numerical method for the anomaloussubdiffusion equation[J].Siam Journal on Numerical Analysis,2008,46(2)：1079-1095.

[7]LIU Q,LIU F,TURNER I,et al.Numerical simulation for the 3D seepage flow with fractional derivatives in porous media[J].Ima Journal of Applied Mathematics,2008,74(2)：201-229(29).

[8]LI C,ZENG F.Finite difference methods for fractional differential equations[J].International Journal of Bifurcation and Chaos,2012,22(4)：1230014.

[9]LI C,DENG W H,SHEN X Q.Exact Solutions and Their Asymptotic Behaviors for the Averaged Generalized Fractional Elastic Models[J].Communications in Theoretical Physics,2014,62(4)：443-450.

[10]PODLUBNY.Fractional differential equations[M].Salt Lake：Academic Press,1999.

[11]LIU G R.Mesh free methods：moving beyond the finite element method[M].Boca Raton：CRCPress,2005.