ABSTRACT

For the Riemann-Liouville (R-L) type fractional differential equations describing sub-diffusion phenomena, a kind of explicit-implicit and implicit-explicit difference schemes for numerically solving this problem is constructed. It takes advantage of the fast computation of explicit formats and the unconditional stability of implicit formats, alternatingly applies the classical explicit format and the implicit format by time layer. Then, using the Fourier method analysis, the format is unconditionally stable and convergent. The numerical test results are consistent with the theoretical analysis results, and show that the computational efficiency of the explicit-implicit and implicit-explicit formats is better than the classical implicit format. The explicit-implicit and implicit-explicit methods are feasible to solve the R-L fractional diffusion equation.

Keywords:R-L Fractional Diffusion Equation, Explicit-Implicit and Implicit-Explicit Scheme, Fourier Method, Stability, Convergence

Riemann-Liouville型分数阶扩散方程的显–隐和隐–显差分方法及数值分析

吴立飞，孙嘉科，杨晓忠^*

华北电力大学数理学院，北京

收稿日期：2019年11月1日；录用日期：2019年11月18日；发布日期：2019年11月25日

摘 要

针对描述慢扩散现象的Riemann-Liouville (R-L)型分数阶扩散方程，构造了求解该问题的一类显–隐和隐–显差分格式。它是利用显式格式快速计算和隐式格式无条件稳定的优点，按时间层交替使用古典显式格式和隐式格式而得。使用傅里叶方法分析可知该格式为无条件稳定且收敛的。数值试验结果与理论分析结果一致，表明显–隐和隐–显格式的计算精度和计算效率均优于经典隐式格式，证实本文显–隐和隐–显格式对求解R-L型分数阶慢扩散方程是有效的。

关键词 :R-L型分数阶扩散方程，显–隐和隐–显格式，傅里叶方法，稳定性，收敛性

Copyright © 2019 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY).

http://creativecommons.org/licenses/by/4.0/

1. 引言

分数阶微分方程产生于一些反常扩散模型，有着深入的物理背景和丰富的理论意义。其特点在于它反映的不仅是局部或某个点的性质或数量，还考虑了历史以及非局部分布式的影响 [1] [2] [3]。在分数阶反常扩散模型中，时间分数阶导数主要刻画反常扩散过程具有的历史依赖性和路径依赖性，空间分数阶导数主要描述扩散过程的非局域特征，它们被用于许多物理、化学过程的数学建模中 [4]。

分数阶微分算子已成为一个描述各种复杂力学行为的重要工具。由于分数阶微分方程的解析解大多很难显式给出，因此研究解决分数阶微分方程的数值解法是必要且具有实际工程意义的 [5] [6]。近些年对于各类分数阶扩散方程的数值解法，已有不少的研究成果 [7] [8] [9] [10]。Wang Hong等对于空间分数阶扩散方程，根据所建立差分格式的特殊结构，提出了快速求解算法 [11]。崔明荣对具有非局部边界条件的时间分数阶扩散方程，考虑了两种不同的时间导数，即Caputo导数和Riemann-Liouville (R-L)导数，构造了其四阶紧致差分格式 [12]。张亚楠等提出了一种求解分数阶慢扩散方程的Crank-Nicolson型差分格式，并对其截断误差进行了详细分析 [13]。刘发旺等对一类多项分数阶波–扩散方程，提出了一些计算有效的隐式数值方法；数值结果表明了理论分析的有效性，这些方法和技术也可以推广到其他具有分数拉普拉斯的多项分数阶时空模型 [14]。吴立飞等给出时间Caputo型分数阶慢扩散方程的交替分段C-N差分格式，将经典C-N格式与四类非对称格式交替使用，提高算法的计算效率 [15]。陈昌明等人根据对分数阶导数的不同定义，利用傅里叶方法证明了显式及隐式格式的稳定性，同时结合外推技术构造了一类高精度算法 [16]。陈昌明等对具有非齐次强迫项的R-L型分数阶导数广义二阶流体的斯托克斯第一问题，构造其隐式差分格式，并分析格式的稳定性、收敛性 [17]；最后数值试验验证了算法的有效性。高广花和孙志忠对分数阶扩散方程构造了一类紧致差分格式，并用能量法证明了格式无条件稳定以及精度为时间(2-a)阶，空间四阶收敛 [18]。

分数动力学方程已被证实在描述反常扩散(慢扩散)的情况下特别适用。近几年，分数扩散方程的理论证明以及丰富的物理和生物实验验证了反常扩散现象的普遍性；对于反常慢扩散随机游走模型，以分数阶慢扩散方程代替了原来的常规扩散方程，提出描述反常扩散粒子遵循分数阶扩散方程 [1] [2] [3]：

$\frac{\partial u\left(x,t\right)}{\partial t}={}_{0}D{}_t^{1-\gamma }\left[\frac{{\partial }^{2}u\left(x,t\right)}{\partial x^2}\right]+f\left(x,t\right),t\ge 0$ . (1)

其中 $u\left(x,t\right)$ 是概率密度函数， ${}_{0}D{}_t^{1-\gamma }u$ 表示 $u\left(x,t\right)$ 的 $1-\gamma$ 阶R-L型分数阶导数，

${}_{0}D{}_t^{1-\gamma }u\left(x,t\right)=\frac{1}{\Gamma \left(\gamma \right)}\frac{\partial }{\partial t}{\displaystyle {\int }_0^t\frac{u\left(x,\tau \right)}{{\left(t-\tau \right)}^{1-\gamma }}\text{d}\tau }$ .

针对此问题，Langlands和Henry构造了隐式格式(L1近似)，并讨论了格式的稳定性和收敛性 [19]；此后，刘发旺等使用傅里叶分析法，分析了对 $0<\gamma <1$ 情形下隐式格式的稳定性和收敛性，以及全局精度 [16] [17]。但是由于分数阶导数的非局部性，分数阶扩散方程数值模拟的计算量和存储量极大，即便使用高性能计算机，也很难进行长时间历程或大计算域的模拟。对于这个问题，本文尝试使用经典显式格式和隐式格式交替使用的方法来提高数值模拟的计算效率 [20] [21]，即对分数阶次扩散模型构造显–隐(E-I)和隐–显(I-E)格式；理论分析和数值试验表明，显–隐和隐–显格式的计算效率和计算精度均优于隐式差分格式。

2. 分数阶扩散方程显–隐和隐–显方法

2.1. R-L型时间分数阶扩散方程

本文主要考虑描述慢扩散现象的R-L型分数阶扩散方程的初边值问题 [2] [3] [4]

$\frac{\partial u\left(x,t\right)}{\partial t}={}_{0}D{}_t^{1-\gamma }\left[\frac{{\partial }^{2}u\left(x,t\right)}{\partial x^2}\right]+f\left(x,t\right),0<t\le T,0<x<L$ . (2)

初始条件为： $u\left(0,t\right)=\phi \left(t\right),0\le t\le T$，

边值条件为： $u\left(1,t\right)=\psi \left(t\right),0\le t\le T,u\left(x,0\right)=w\left(x\right),0\le x\le L$。

其中 $0<\gamma <1$，$f\left(x,t\right)$ 、 $\phi \left(t\right)$ 、 $\psi \left(t\right)$ 、 $w\left(x\right)$ 都是足够光滑的函数。

给出方程(2)式两边作积分 ${}_{0}D{}_t^{\gamma -1}$ 运算，可得

${}_{0}D{}_t^{\gamma -1}\left(\frac{\partial u\left(x,t\right)}{\partial t}\right)={}_{0}D{}_t^{\gamma -1}\left\{{}_{0}D{}_t^{1-\gamma }\left[K_{\gamma }\frac{{\partial }^{2}u\left(x,t\right)}{\partial x^2}\right]+f\left(x,t\right)\right\}$ ,(3)

上式等价为

${}_0{}^{C}D{}_t^{\gamma }\frac{\partial u\left(x,t\right)}{\partial t}=K_{\gamma }\frac{{\partial }^{2}u\left(x,t\right)}{\partial x^2}-{\left\{{}_{0}D{}_t^{-\gamma }\left[K_{\gamma }\frac{{\partial }^{2}u\left(x,t\right)}{\partial x^2}\right]\right\}|}_{t=0}+{}_{0}D{}_t^{\gamma -1}f\left(x,t\right)$ , (4)

由文献 [3] 可知，如果 $u\left(x,t\right)\in C_{x,t}^{2,1}\left(\left[0,L\right]\times \left[0,T\right]\right)$，则有

${\left\{{}_{0}D{}_t^{-\gamma }\left[K_{\gamma }\frac{{\partial }^{2}u\left(x,t\right)}{\partial x^2}\right]\right\}|}_{t=0}=0$ .

即(4)式为

${}_0^C{D}_t^{\gamma }\frac{\partial u\left(x,t\right)}{\partial t}=K_{\gamma }\frac{{\partial }^{2}u\left(x,t\right)}{\partial x^2}+g\left(x,t\right),g\left(x,t\right)={}_{0}D{}_t^{\gamma -1}f\left(x,t\right)$ (5)

2.2. 显–隐和隐–显格式的构造

令 $t_k=k\tau ,k=0,1,\cdots ,N$，$x_j=jh,j=0,1,\cdots ,M$，$\tau =T/N$，$h=L/M$。定义网格函数

$u_i^n=u\left(x_i,t_n\right)$ , $f_i^n=f\left(x_i,t_n\right)$ , $g_i^n=g\left(x_i,t_n\right)$ , $\left(x_i,t_n\right)\in \left[0,L\right]\times \left[0,T\right]$ .

定义如下记号： ${\delta }_x^2{u}_i=\frac{1}{h^2}\left(u_{i+1}-2u_i+u_{i-1}\right)$。

定义离散分数阶导数算子

$D_{\tau }^{\gamma }{w}^n=\frac{{\tau }^{-\gamma }}{\Gamma \left(2-\gamma \right)}\left[w^n-{\displaystyle \underset{k=1}{\overset{n-1}{\sum }}\left(a_{n-k-1}-a_{n-k}\right)w^k}-a_{n-1}{w}^0\right]$ , $a_k={\left(k+1\right)}^{1-\alpha }-k^{1-\alpha }$ .

引理1 [6] 假设 $0<\alpha <1$，$y\in C^2\left[0,t_n\right]$，有

$\begin{array}{l}\frac{1}{\Gamma \left(1-\alpha \right)}{\displaystyle {\int }_0^{t_n}\frac{y^{\prime }\left(\xi \right)\text{d}\xi }{{\left(t-\xi \right)}^{\alpha }}}-\frac{{\tau }^{-\alpha }}{\Gamma \left(2-\alpha \right)}\left[y^n-{\displaystyle \underset{k=1}{\overset{n-1}{\sum }}\left(a_{n-k-1}-a_{n-k}\right)y^k}-a_{n-1}{y}^0\right]\\ \le \frac{1}{\Gamma \left(2-\alpha \right)}\left[\frac{1-\alpha }{12}+\frac{2^{2-\alpha }}{2-\alpha }-\left(1+2^{\alpha }\right)\right]\underset{0\le t\le t_n}{\max }\left|y^{″}\left(t\right)\right|{\tau }^{2-\alpha }\end{array}$ .

时间分数阶慢扩散方程的隐式格式为：

$\frac{{\tau }^{-\gamma }}{\Gamma \left(2-\gamma \right)}\left[u_j^k-{\displaystyle \underset{i=1}{\overset{k-1}{\sum }}\left(a_{k-i-1}-a_{k-i}\right)u_j^i}-a_{k-1}{u}_j^0\right]={\delta }_x^2{u}_j^k+g_j^k$ ,(6)

时间分数阶慢扩散方程的显式格式为：

$\frac{{\tau }^{-\gamma }}{\Gamma \left(2-\gamma \right)}\left[u_j^k-{\displaystyle \underset{i=1}{\overset{k-1}{\sum }}\left(a_{k-i-1}-a_{k-i}\right)u_j^i}-a_{k-1}{u}_j^0\right]={\delta }_x^2{u}_j^{k-1}+g_j^k,j=1,2,\cdots ,M-1$ . (7)

初始、边界条件相应离散为：

$\begin{array}{l}{u}_0^k=\phi \left(t_k\right),u_m^k=\psi \left(t_k\right),k=0,1,\cdots ,N\\ u_j^0=w\left(x_j\right),j=0,1,\cdots ,M\end{array}$

经典显式格式为条件稳定的，其优点是可以显式快速计算；而隐式格式为无条件稳定，但需要求解代数方程组，计算效率较低。结合两种差分格式的优点，在时间层上交替使用，可得显–隐和隐–显差分格式。

分数阶慢扩散方程的显–隐式格式为：

$\{\begin{array}{l}\frac{{\tau }^{-\gamma }}{\Gamma \left(2-\gamma \right)}\left[u_j^{2k+1}-{\displaystyle \underset{i=1}{\overset{2k}{\sum }}\left(a_{2k-i}-a_{2k-i+1}\right)u_j^i}-a_{2k}{u}_j^0\right]={\delta }_x^2{u}_j^{2k}+g_j^{2k+1},\\ \frac{{\tau }^{-\gamma }}{\Gamma \left(2-\gamma \right)}\left[u_j^{2k+2}-{\displaystyle \underset{i=1}{\overset{2k+1}{\sum }}\left(a_{2k-i+1}-a_{2k-i}\right)u_j^i}-a_{2k+1}{u}_j^0\right]={\delta }_x^2{u}_j^{2k+2}+g_j^{2k+2},\end{array}$ (8)

同理，可构造分数阶慢扩散方程的隐–显格式为

$\{\begin{array}{l}\frac{{\tau }^{-\gamma }}{\Gamma \left(2-\gamma \right)}\left[u_j^{2k+1}-{\displaystyle \underset{i=1}{\overset{2k}{\sum }}\left(a_{2k-i}-a_{2k-i+1}\right)u_j^i}-a_{2k}{u}_j^0\right]={\delta }_x^2{u}_j^{2k+1}+g_j^{2k+1},\\ \frac{{\tau }^{-\gamma }}{\Gamma \left(2-\gamma \right)}\left[u_j^{2k+2}-{\displaystyle \underset{i=1}{\overset{2k+1}{\sum }}\left(a_{2k-i+1}-a_{2k-i}\right)u_j^i}-a_{2k+1}{u}_j^0\right]={\delta }_x^2{u}_j^{2k+1}+g_j^{2k+2}.\end{array}$ (9)

其中 $k=0,1,2,3,\cdots ,j=1,2,\cdots ,M-1$。

3. 显–隐和隐–显方法的稳定性分析

首先分析时间分数阶慢扩散方程的显–隐格式的稳定性，其显–隐格式可以写成如下形式：

$\{\begin{array}{l}{u}_j^{2k+1}=\left(a_0-a_1-2\mu \right)u_j^{2k}+\mu u_{j+1}^{2k}+\mu u_{j-1}^{2k}+{\displaystyle \underset{i=1}{\overset{2k-1}{\sum }}\left(a_{2k-i}-a_{2k-i+1}\right)u_j^i}-a_{2k}{u}_j^0+\mu g_j^{2k+1},\\ \left(1+2\mu \right)u_j^{2k+2}-\mu u_{j+1}^{2k+2}-\mu u_{j-1}^{2k+2}=\left(a_0-a_1\right)u_j^{2k+1}+{\displaystyle \underset{i=1}{\overset{2k}{\sum }}\left(a_{2k-i+1}-a_{2k-i+2}\right)u_j^i}-a_{2k+1}{u}_j^0+\mu g_j^{2k+2},\end{array}$ (10)

其中 $\mu =\frac{{\tau }^{\gamma }}{h^2}\Gamma \left(2-\gamma \right)$。

设 $U_j^k$ 是显–隐格式的近似解，定义 ${\rho }_j^k=u_j^k-U_j^k$ $\left(k=1,2,\cdots ,N,j=1,2,\cdots ,M-1\right)$ 和 ${\rho }^k={\left[{\rho }_1^k,{\rho }_2^k,\cdots ,{\rho }_{M-1}^k\right]}^{\text{T}}$，代入(10)式可得如下误差扰动方程：

$\{\begin{array}{l}{\rho }_j^{2k+1}=\left(a_0-a_1-2\mu \right){\rho }_j^{2k}+\mu {\rho }_{j+1}^{2k}+\mu {\rho }_{j-1}^{2k}+{\displaystyle \underset{i=1}{\overset{2k-1}{\sum }}\left(a_{2k-i}-a_{2k-i+1}\right){\rho }_j^i}-a_{2k}{\rho }_j^0,\\ \left(1+2\mu \right){\rho }_j^{2k+2}-\mu {\rho }_{j+1}^{2k+2}-\mu {\rho }_{j-1}^{2k+2}=\left(a_0-a_1\right){\rho }_j^{2k+1}+{\displaystyle \underset{i=1}{\overset{2k}{\sum }}\left(a_{2k-i+1}-a_{2k-i+2}\right){\rho }_j^i}-a_{2k+1}{\rho }_j^0,\end{array}$ (11)

下面使用Fourier方法分析显–隐格式的稳定性。在此，定义网格函数

${\rho }^k\left(x\right)=\{\begin{array}{l}0,0\le x\le \frac{h}{2}\text{,}\\ {\rho }_j^k,x_j-\frac{h}{2}<x\le x_j+\frac{h}{2},j=1,2,\cdots ,M-1,\\ 0,L-\frac{h}{2}<x\le L.\end{array}$ $\left(k=1,2,\cdots ,N\right)$

${\rho }^k\left(x\right)$ 的Fourier展开形式为：

${\rho }^k\left(x\right)={\displaystyle \underset{l=-\infty }{\overset{\infty }{\sum }}{d}_k\left(l\right){\text{e}}^{i2\pi lx/L}}$ ,

其中系数 $d_k\left(l\right)=\frac{1}{L}{\displaystyle {\int }_0^L{\rho }^k\left(x\right){\text{e}}^{-i2\pi lx/L}\text{d}x}$。

另外，离散2模的定义

$\begin{array}{c}{\Vert {\rho }^k\Vert }_2^2={\displaystyle \underset{j=1}{\overset{M-1}{\sum }}h{\left|{\rho }_j^k\right|}^2}\\ ={\displaystyle {\int }_0^{\frac{h}{2}}{\left|{\rho }^k\left(x\right)\right|}^2\text{d}x}+{\displaystyle \underset{j=1}{\overset{M-1}{\sum }}{\displaystyle {\int }_{x_j-\frac{h}{2}}^{x_j+\frac{h}{2}}{\left|{\rho }^k\left(x\right)\right|}^2\text{d}x}}+{\displaystyle {\int }_{L-\frac{h}{2}}^L{\left|{\rho }^k\left(x\right)\right|}^2\text{d}x}\\ ={\displaystyle {\int }_0^L{\left|{\rho }^k\left(x\right)\right|}^2\text{d}x}\end{array}$ (12)

应用Parseval等式，有

${\displaystyle {\int }_0^L{\left|{\rho }^k(x)\right|}^2\text{d}x}=L{\displaystyle \underset{m=-\infty }{\overset{\infty }{\sum }}{\left|d_k\left(m\right)\right|}^2}$，即 ${\Vert {\rho }^k\Vert }_2^2=L{\displaystyle \underset{m=-\infty }{\overset{\infty }{\sum }}{\left|d_k\left(m\right)\right|}^2}$ .

基于上面的分析，假设(11)式有如下形式的波动解 ${\rho }_j^k=d_k{\text{e}}^{i\sigma jh}$，其中 $\sigma =2\pi m/L$。

当 $m=1$ 时(m是时间层数)，我们使用显式差分格式，有

${\rho }_j^1=\left(1-2\mu \right){\rho }_j^0+\mu {\rho }_{j+1}^0+\mu {\rho }_{j-1}^0$ ,

$d_1=\left(\mu {\text{e}}^{i\sigma h}+1-2\mu +\mu {\text{e}}^{-i\sigma h}\right)d_0$ ,

由 ${\text{e}}^{\pm i\theta }=\cos \theta \pm i\sin \theta$，可得

$\left|d_1\right|=\left|1-4\mu {\sin }^2\frac{\sigma h}{2}\right|\left|d_0\right|\le C\left|d_0\right|$ , (13)

当 $m=2$ 时，第一层使用显式格式，第二层使用隐式格式，

${\rho }_j^1=\left(1-2\mu \right){\rho }_j^0+\mu {\rho }_{j+1}^0+\mu {\rho }_{j-1}^0$ ,

$\left(1+2\mu \right){\rho }_j^2-\mu {\rho }_{j+1}^2-\mu {\rho }_{j-1}^2=\left(a_0-a_1\right){\rho }_j^1+a_1{\rho }_j^0$ ,

联立两式，可得

$\left(1+4\mu {\sin }^2\frac{\sigma h}{2}\right)d_2=\left(a_0-a_1\right)\left(1-4\mu {\sin }^2\frac{\sigma h}{2}\right)d_0+a_1{d}_0$ ,

$\left|d_2\right|=\left|\frac{\left(1-a_1\right)\left(1-4\mu {\sin }^2\frac{\sigma h}{2}\right)+a_1}{1+4\mu {\sin }^2\frac{\sigma h}{2}}\right|\left|d_0\right|\le \left|\frac{1-\left(1-a_1\right)4\mu {\sin }^2\frac{\sigma h}{2}}{1+4\mu {\sin }^2\frac{\sigma h}{2}}\right|\left|d_0\right|\le C\left|d_0\right|$ .(14)

当 $m=2k$ 时，假设 $\left|d_n\right|\le C\left|d_0\right|,n\le 2k$ 成立；当 $m=2k+2$，有

$\{\begin{array}{l}{d}_{2k+1}=\left(a_0-a_1-4\mu {\sin }^2\frac{\sigma h}{2}\right)d_{2k}+{\displaystyle \underset{i=1}{\overset{2k-1}{\sum }}\left(a_{2k-i}-a_{2k-i+1}\right)d_i}-a_{2k}{d}_0,\\ \left(1+4\mu {\sin }^2\frac{\sigma h}{2}\right)d_{2k+2}=\left(a_0-a_1\right)d_{2k+1}+{\displaystyle \underset{i=1}{\overset{2k}{\sum }}\left(a_{2k-i+1}-a_{2k-i+2}\right)d_j}-a_{2k+1}{d}_0,\end{array}$

两式结合，消去 $d_{2k+1}$ 项，可得

$\begin{array}{l}\left(1+4\mu {\sin }^2\frac{\sigma h}{2}\right)d_{2k+2}\\ =\left(a_0-a_1\right)\left[\left(a_0-a_1-4\mu {\sin }^2\frac{\sigma h}{2}\right)d_{2k}+{\displaystyle \underset{i=1}{\overset{2k-1}{\sum }}\left(a_{2k-i}-a_{2k-i+1}\right)d_i}-a_{2k}{d}_0\right]\\ +{\displaystyle \underset{i=1}{\overset{2k}{\sum }}\left(a_{2k-i+1}-a_{2k-i+2}\right)d_j}-a_{2k+1}{d}_0\end{array}$

令 $\tilde{\mu }=4\mu {\sin }^2\frac{\sigma h}{2}$，由上式可得

$\begin{array}{l}\left(1+4\mu {\sin }^2\frac{\sigma h}{2}\right)\left|d_{2k+2}\right|\\ \le \left\{\left(a_0-a_1\right)\left[\left|\left(a_0-a_1\right)+4\mu {\sin }^2\frac{\sigma h}{2}\right|+{\displaystyle \underset{i=1}{\overset{2k-1}{\sum }}\left|a_{2k-i}-a_{2k-i+1}\right|}+\left|a_{2k}\right|\right]+{\displaystyle \underset{i=1}{\overset{2k}{\sum }}\left|a_{2k-i+1}-a_{2k-i+2}\right|}+\left|a_{2k+1}\right|\right\}C\left|d_0\right|\\ \le \left\{\left(1-a_1\right)\left(1+4\mu {\sin }^2\frac{\sigma h}{2}\right)+a_1\right\}C\left|d_0\right|\\ \le \left[1+\left(1-a_1\right)4\mu {\sin }^2\frac{\sigma h}{2}\right]C\left|d_0\right|\end{array}$

$\left|d_{2k+2}\right|\le \frac{1+4\mu {\sin }^2\frac{\sigma h}{2}}{1+\left(1-a_1\right)4\mu {\sin }^2\frac{\sigma h}{2}}C\left|d_0\right|\le C\left|d_0\right|$ . (15)

综上分析，有如下定理

定理1 R-L型分数阶扩散方程的显–隐格式是无条件稳定的。

证明：由(13，14，15)式，可知

${\Vert {\rho }^k\Vert }^2=L{\displaystyle \underset{m=-\infty }{\overset{\infty }{\sum }}{\left|d_k\left(m\right)\right|}^2}\le CL{\displaystyle \underset{m=-\infty }{\overset{\infty }{\sum }}{\left|d_0\left(m\right)\right|}^2}=C{\Vert {\rho }^0\Vert }^2,k=1,2,\cdots ,N$ .

因此，分数阶扩散方程的显–隐格式是稳定的，定理1证毕。

同理，可有

定理2 R-L型分数阶扩散方程的隐–显格式是无条件稳定的。

4. 显–隐和隐–显方法的收敛性分析

引理2 [13] 假设 $0<\alpha <1$，$y\in c^2\left[0,t_n\right]$，有

$\begin{array}{l}\frac{1}{\Gamma \left(1-\alpha \right)}{\displaystyle {\int }_0^{t_n}\frac{y^{\prime }\left(\xi \right)\text{d}\xi }{{\left(t-\xi \right)}^{\alpha }}}-\frac{{\tau }^{-\alpha }}{\Gamma \left(2-\alpha \right)}\left[y^n-{\displaystyle \underset{k=1}{\overset{n-1}{\sum }}\left(a_{n-k-1}-a_{n-k}\right)y^k}-a_{n-1}{y}^0\right]\\ \le \frac{1}{\Gamma \left(2-\alpha \right)}\left[\frac{1-\alpha }{12}+\frac{2^{2-\alpha }}{2-\alpha }-\left(1+2^{\alpha }\right)\right]\underset{0\le t\le t_n}{\max }\left|y^{″}\left(t\right)\right|{\tau }^{2-\alpha }\end{array}$ (16)

引理3 [15] $D_{\tau }^{\alpha }{u}_j^n={}_0{}^{C}D{}_{\tau }^{\alpha }{u}_j^n-\theta {}_0{}^{C}D{}_{\tau }^{\alpha +1}{u}_j^n+R_1,R_1\le C{\tau }^{2-\alpha },0\le \theta \le 1$，C是常数。

首先分别分析显–隐格式(8)中，显式格式和隐式格式的截断误差。

$\frac{{\tau }^{-\gamma }}{\Gamma \left(2-\gamma \right)}\left[u_j^{2k+1}-{\displaystyle \underset{i=1}{\overset{2k}{\sum }}\left(a_{2k-i}-a_{2k-i+1}\right)u_j^i}-a_{2k}{u}_j^0\right]=\left(1-\theta \right){\delta }_x^2{u}_j^{2k}+\theta {\delta }_x^2{u}_j^{2k}+g_j^{2k+1}$ , (17)

可知，当q = 0时，为显式格式，当q = 1时，为隐式格式。定义截断误差

$R_j^{2k+1}=\frac{{\tau }^{-\gamma }}{\Gamma \left(2-\gamma \right)}\left[u_j^{2k+1}-{\displaystyle \underset{i=1}{\overset{2k}{\sum }}\left(a_{2k-i}-a_{2k-i+1}\right)u_j^i}-a_{2k}{u}_j^0\right]-\left(1-\theta \right){\delta }_x^2{u}_j^{2k}-\theta {\delta }_x^2{u}_j^{2k}-g_j^{2k+1}$ (18)

$\theta K_{\alpha }{\delta }_x^2{u}_j^k+\left(1-\theta \right)K_{\alpha }{\delta }_x^2{u}_j^{k-1}=u_{xx}+\frac{h_{}^2}{12}{u}_{xxxx}+\left(1-\theta \right)\tau u_{xxt}+O\left({\tau }^{2-\alpha }+h^2\right)$ (19)

由引理3可得

$\begin{array}{c}{R}_j^{2k+1}=\frac{{\tau }^{-\gamma }}{\Gamma \left(2-\gamma \right)}\left[u_j^{2k+1}-{\displaystyle \underset{i=1}{\overset{2k}{\sum }}\left(a_{2k-i}-a_{2k-i+1}\right)u_j^i}-a_{2k}{u}_j^0\right]-\left(1-\theta \right){\delta }_x^2{u}_j^{2k}-\theta {\delta }_x^2{u}_j^{2k}-g_j^{2k+1}\\ ={}_0{}^{C}D{}_{\tau }^{\alpha }u\left(x_i,t_{k+1}\right)-\left(1-\theta \right)\tau {}_0{}^{C}D{}_{\tau }^{\alpha +1}u\left(x_i,t_{k+1}\right)-u_{xx}+\frac{h^2}{12}{u}_{xxxx}+\left(1-\theta \right)\tau u_{xxt}+O\left({\tau }^{2-\alpha }+h^2\right)\\ =O\left({\tau }^{2-\alpha }+h^2\right)\end{array}$

因此，存在正数 $c_1$ 使得

$\left|R_j^k\right|\le c_1\left({\tau }^{2-\gamma }+h^2\right)$ , $k=1,2,\cdots ,N,j=1,2,\cdots ,M-1$ . (20)

定义

$e_j^k=u\left(x_j,t_k\right)-u_j^k$ $\left(k=1,2,\cdots ,N,j=1,2,\cdots ,M-1\right)$ .

和

$e^k={\left[e_1^k,e_2^k,\cdots ,e_{M-1}^k\right]}^{\text{T}}$ , $R^k={\left[R_1^k,R_2^k,\cdots ,R_{M-1}^k\right]}^{\text{T}}$ .

下面使用Fourier方法分析显–隐格式的收敛性。在此，定义网格函数

$e^k\left(x\right)=\{\begin{array}{l}0,0\le x\le \frac{h}{2},\\ e_j^k,x_j-\frac{h}{2}<x\le x_j+\frac{h}{2},j=1,2,\cdots ,M-1,\\ 0,L-\frac{h}{2}<x\le L.\end{array}$ $\left(k=1,2,\cdots ,N\right)$

和

$R^k\left(x\right)=\{\begin{array}{l}0,0\le x\le \frac{h}{2},\\ R_j^k,x_j-\frac{h}{2}<x\le x_j+\frac{h}{2},j=1,2,\cdots ,M-1,\\ 0,L-\frac{h}{2}<x\le L.\end{array}$ $\left(k=1,2,\cdots ,N\right)$

$e^k\left(x\right),R^k\left(x\right)$ 的Fourier展开形式为：

$e^k\left(x\right)={\displaystyle \underset{l=-\infty }{\overset{\infty }{\sum }}{\eta }_k\left(l\right){\text{e}}^{i2\pi lx/L}}$，其中系数 ${\eta }_k\left(l\right)=\frac{1}{L}{\displaystyle {\int }_0^L{e}^k\left(x\right){\text{e}}^{-i2\pi lx/L}\text{d}x}$。

$R^k\left(x\right)={\displaystyle \underset{l=-\infty }{\overset{\infty }{\sum }}{\xi }_k\left(l\right){\text{e}}^{i2\pi lx/L}}$，其中系数 ${\xi }_k\left(l\right)=\frac{1}{L}{\displaystyle {\int }_0^L{R}^k\left(x\right){\text{e}}^{-i2\pi lx/L}\text{d}x}$。

同样，关于离散2模，有

,(21)

. (22)

基于上面的分析，假设，，其中，可得

,

,

(23)

引理4 设是方程(23)的解，存在一个正数，使得

.

证明：由定义，可知，有。

另外，由(20)式和(22)式可知

(24)

(22)式等式右边的系数数列是收敛的，因此存在一个正数使得

. (25)

因此有

.

参考稳定性证明，使用数学归纳法可得

当时，

,

,

当时，

,

,

,

假设

.

由引理1，有

即

.

由于

.

则存在正常数，使得

.

综上，可有

定理3 R-L型分数阶扩散方程的显–隐格式是无条件收敛的，收敛阶是.

证明：由引理4及，可得

.

定理3证毕。

同理，可有

定理4 R-L型分数阶扩散方程的隐–显格式是无条件收敛的，收敛阶是.

5. 数值试验

数值试验基于Intel Core i5-2400 CPU@3.10GHz，在MatlabR2014a环境下进行；我们将通过数值试验验证前面的理论分析。

例1. 非齐次R-L型分数阶扩散方程 [4] [16]

.

其边界条件为：

初始条件为：.

此时解析解为：.

图1. 解析解和隐–显、显–隐格式数值解的解曲面

取时，给出解析解、显–隐和隐–显格式的解曲面图1；可以看出解析解、显–隐和

隐–显格式解曲面形状一致且光滑。解析解与差分格式解的最大模误差定义为。

对分别取0.1~0.9等不同值，取，隐式格式、显–隐格式和隐–显格式三种数值格式的最大模误差情况见表1。由表1可知，隐式格式与显–隐和隐–显格式的精度相当；对不同值隐式格式和显–隐格式的误差比较；可以看出当值大于0.5时，隐式格式的精度高于显–隐格式，当值小于0.5时，显–隐格式的精度优于隐式格式。随着值得减小，显–隐格式和隐–显格式的误差越来越小；整数阶显–隐差分格式为二阶差分格式，其精度优于隐式格式。当值接近于0时，显–隐格式和隐–显格式的特性接近整数阶的情形。

表1. 三种数值格式的精度比较分析

下面，由最大模误差定义时间收敛阶以及空间收敛阶，验证理论分析。固定空间步长h，时间收敛阶如下定义：。固定时间步长，定义空间收敛阶如下：。

表2. 三种差分格式的时间方向收敛阶(h = 64)

由表2，对不同值，隐式格式与隐–显、显–隐格式在时间方向上的收敛阶为阶；隐式格式与显–隐格式、隐–显格式的精度也相当。由表3可看出，对不同值，隐式格式和隐–显、显–隐格式在空间上的收敛阶均为2阶，与理论分析一致。

例2. 齐次R-L型分数阶扩散方程 [4] [16]

.

其边界条件为：，，。

初始条件为：

其中表示由点处热源产生的棒的温度分布。

表3. 三种差分格式的空间方向收敛阶

图2. 三种差分格式的解曲面

表4. 三种差分格式计算时间的比较(单位：s)

当初始条件为分段函数，取，，，取不同值时，隐式格式、显–隐和隐–显格式的解曲面图2；由图2可知隐式格式和显–隐(隐–显)格式解形状一致，且光滑，说明显–隐和隐–显差分格式对求解例2问题的有效性。由表4可知，隐–显格式和显–隐格式的计算时间相当，但是随着网格数的增加，隐式格式与显–隐和隐–显格式相比，计算时间增加的速度越来越快；其趋势在图3中可见。可见对于求解分数阶慢扩散模型，显–隐和隐–显格式与隐式格式相比，在计算效率方面有明显的改进。

图3. 三种差分格式的计算时间比较

6. 结论

本文对R-L型分数阶扩散方程，构造了其显–隐和隐–显差分格式；通过傅里叶方法对构造格式进行了理论分析，得出显–隐和隐–显格式均为无条件稳定和收敛的，通过数值试验证实了理论分析。交替使用显式、隐式两个古典格式，可以使算法计算复杂度得到下降。由数值试验可知，显–隐和隐–显格式的计算精度与隐式格式相当，但计算效率与隐式格式相比有着明显的提高；表明显–隐和隐–显格式对求解R-L型分数阶扩散方程是有效的。

基金项目

国家科技重大专项子课题(2017ZX07101001-01)，中央高校基金科研业务专项资金资助(2018MS168)。

文章引用

吴立飞,孙嘉科,杨晓忠. Riemann-Liouville型分数阶扩散方程的显–隐和隐–显差分方法及数值分析
Explicit-Implicit and Implicit-Explicit Difference Methods and Numerical Analysis for Riemann-Liouville Type Fractional Diffusion Equation[J]. 理论数学, 2019, 09(09): 1060-1074. https://doi.org/10.12677/PM.2019.99131

参考文献