利用Laplace变换求解时间分数阶对流扩散方程

Solving Time Fractional Convection-Diffusion Equation Using Laplace Transform

摘要: 本文提出了求解时间分数阶对流扩散方程两种高效数值算法。首先基于Laplace变换及指数变换将原问题转化为整数阶扩散问题；然后采用Crank-Nicolson格式并分别结合二阶中心差分和四阶紧致差分方法，设计出两种求解时间分数阶对流扩散方程的高精度差分格式，并利用Fourier方法证明两种差分格式都是稳定的。数值实验验证了两种格式的有效性。

Abstract: This paper proposes two efficient numerical algorithms for solving the time fractional convection-diffusion equation. First, the original problem is transformed into an integer-order diffusion problem based on Laplace transform and exponential transform. Then, using Crank-Nicolson format and combining second order central difference and fourth order compact difference respectively, two kinds of high precision difference formats for solving time fractional order convection-diffusion equations are designed, and both schemes are proved to be stable by using Fourier method. Numerical experiments verify the effectiveness of the two formats.

1. 引言

在研究分子扩散、粘性流体运动及物质运输过程等问题时，经常要根据实际情况建立起相应的对流扩散方程来求解。特别是分数阶对流扩散方程，有时更能准确刻画物质的运动规律、描述一些物理现象。因此，设计出比较稳定的数值格式来求解时间分数阶对流扩散方程一直备受人们的关注，具有极为重要的理论意义和实际意义。

本文研究如下的时间分数阶对流扩散方程：

$\{\begin{array}{l}\frac{{\partial }^{\alpha }u}{\partial t^{\alpha }}=\frac{{\partial }^{2}u}{\partial x^2}-\beta \frac{\partial u}{\partial x}+f\left(x,t\right),0\le x\le L,0\le t\le T,0<\alpha <1,\beta \in Z^{+},\\ \begin{array}{l}u\left(x,0\right)=g\left(x\right),0\le x\le L,\hfill \\ u\left(0,t\right)=\varphi \left(t\right),u\left(L,t\right)=\psi \left(t\right),0\le t\le T,\hfill \end{array}\end{array}$ (1)

国内外学者们对分数阶导数、整数阶对流扩散方程及分数阶对流扩散方程进行了许多重要研究。如陈焕贞等 [1] 详细梳理了分数阶扩散模型方面的理论和数值模拟方面的进展和动态；Ren等 [2] 提出了一种用Laplace变换近似分数阶项的方法；李娴娟等 [3] 把谱方法用于求解分数阶偏微分方程在空间和时间方向上的离散，还用最优误差估计证明了该方法的收敛性，并用数值结果验证了理论分析；王玉娇等 [4] 讨论了分数阶导数的数值解法以及其在生物力学方面的一些应用；叶俊杰等 [5] 研究了一种求解Riemann-Liouville型分数阶微分方程的微分变换方法；董建强等 [6] 使用样条插值和帕德逼近设计了一种求解一维对流扩散方程具有高精度的紧致差分格式；郭非凡等 [7] 研究了变系数分数阶对流扩散方程的一种有限差分方法；刘扬等 [8] 构造了一种严格对角占优的Crank-Nicolson有限差分格式来求解一维非定常对流扩散方程；黄素珍等 [9] 通过指数变换、反指数变换和待定系数法，建立了一种三层差分格式来求解一维对流扩散方程；宋光珍 [10] 针对时间分数阶扩散方程初边值问题给出了两种高精度数值格式；余跃玉 [11] 针对变系数时间分数阶对流扩散方程设计出一种隐式有限差分格式，并给出格式差分解存在性与唯一性的证明；王学彬等 [12] 给出了两种常见分数阶导数的Laplace变换公式，同时给出了具体实例来说明如何使用Laplace变换；覃平阳等 [13] 提出一种隐式差分格式；张艳等 [14] 通过Laplace变换得出了分数阶微分方程Mittag-Leffler函数形式解析解；沈淑君等 [15] 讨论了时间、空间及时空方向的分数阶对流扩散方程的解法；常娟等 [16] 设计出一种高阶紧致指数型差分格式来求解对流扩散方程，并采用具有并行性的AGE迭代法对其求解。

本文利用Laplace变换逼近时间分数阶对流扩散方程的分数阶项，将其转化为整数阶问题，然后分别通过二阶中心差分格式和四阶紧致差分格式求方程的数值解。

2. 时间分数阶对流扩散方程的数值解

2.1. 利用Laplace变换将分数阶问题转化为整数阶问题

针对方程(1)，首先使用Laplace变换逼近方程的分数阶导数：

$L\left\{\frac{{\partial }^{\alpha }u\left(x,t\right)}{\partial t^{\alpha }}\right\}=s^{\alpha }\hat{u}\left(x,s\right)-s^{\alpha -1}\left(x,0\right)=s^{\alpha }\left[\hat{u}\left(x,s\right)-s^{-1}u\left(x,0\right)\right]$ (2)

其中， $\hat{u}\left(x,s\right)$ 是 $u\left(x,t\right)$ 的Laplace变换。由于 $\alpha \in \left(0,1\right)$，故可以将 $s^{\alpha }$ 线性化：

$s^{\alpha }\approx \alpha s^1+\left(1-\alpha \right)s^0=\alpha s+\left(1-\alpha \right)$ (3)

将(3)式代入到(2)式中可得(4)式：

$\begin{array}{c}L\left\{\frac{{\partial }^{\alpha }u\left(x,t\right)}{\partial t^{\alpha }}\right\}\approx \left[\alpha s+\left(1-\alpha \right)\right]\left[\hat{u}\left(x,s\right)-s^{-1}u\left(x,0\right)\right]\\ =\alpha s\left[\hat{u}\left(x,s\right)-s^{-1}u\left(x,0\right)\right]+\left(1-\alpha \right)\left[\hat{u}\left(x,s\right)-s^{-1}u\left(x,0\right)\right]\end{array}$ (4)

再由Laplace逆变换，可得：

$\frac{{\partial }^{\alpha }u\left(x,t\right)}{\partial t^{\alpha }}\approx \alpha \frac{\partial u\left(x,t\right)}{\partial t}+\left(1-\alpha \right)\left[u\left(x,t\right)-u\left(x,0\right)\right]$ (5)

于是，原时间分数阶对流扩散方程可转化为如下的整数阶方程(6)：

$\alpha \frac{\partial u\left(x,t\right)}{\partial t}+\left(1-\alpha \right)\left[u\left(x,t\right)-g\left(x\right)\right]=\frac{{\partial }^{2}u}{\partial x^2}-\beta \frac{\partial u}{\partial x}+f\left(x,t\right)$ (6)

引入指数变换：

$u\left(x,t\right)={\text{e}}^{\frac{\beta }{2}x}v\left(x,t\right)$ (7)

将(7)带入方程(6)整理可得：

$\alpha \frac{\partial v}{\partial t}=\frac{{\partial }^{2}v}{\partial x^2}-\left(1-\alpha +\frac{{\beta }^2}{4}\right)v\left(x,t\right)+{\text{e}}^{-\frac{\beta }{2}x}\left[f\left(x,t\right)+\left(1-\alpha \right)g\left(x\right)\right]$ (8)

令 $\lambda =1-\alpha +\frac{{\beta }^2}{4}$， $G\left(x,t\right)={\text{e}}^{-\frac{\beta }{2}x}\left[f\left(x,t\right)+\left(1-\alpha \right)g\left(x\right)\right]$，有：

$\alpha \frac{\partial v}{\partial t}=\frac{{\partial }^{2}v}{\partial x^2}-\lambda v\left(x,t\right)+G\left(x,t\right)$ (9)

2.2. 二阶中心差分格式

以 $\tau$ 表示时间步长，以h表示空间步长。网格点为 $\left(x_i,t_k\right)$，其中 $x_i=ih\left(0\le i\le N\right)$， $t_k=k\tau \left(0\le k\le M\right)$。

考虑方程(9)在 $k+1/2$ 时刻的值并对于时间和空间导数均利用中心差分，可得：

$\alpha {\delta }_t{v}_i^{k\text{+}\frac{1}{2}}={\delta }_x^2{v}_i^{k\text{+}\frac{1}{2}}-\lambda v_i^{k\text{+}\frac{1}{2}}+G_i^{k\text{+}\frac{1}{2}}+O\left({\tau }^2+h^2\right)$ (10)

式中：

${\delta }_t{v}_i^{k\text{+}\frac{1}{2}}=\frac{1}{\tau }\left(v_i^{k+1}-v_i^k\right)$ (11)

${\delta }_x^2{v}_i=\frac{1}{h^2}\left(v_{i+1}-2v_i+v_{i-1}\right)$ (12)

项中 $v_i^{k\text{+}\frac{1}{2}}$ 以第k + 1层和第k层的算数平均值代替，代入(10)式并略去高阶无穷小项，得到二阶中心差分格式：

$\begin{array}{l}-\frac{r}{2}{v}_{i-1}^{k+1}+\left(\alpha +r+\frac{\lambda \tau }{2}\right)v_i^{k+1}-\frac{r}{2}{v}_{i+1}^{k+1}\\ =\frac{r}{2}{v}_{i-1}^k+\left(\alpha -r-\frac{\lambda \tau }{2}\right)v_i^k+\frac{r}{2}{v}_{i+1}^k+\tau G_i^{k+\frac{1}{2}},1\le i\le N-1,0\le k\le M-1.\end{array}$ (13)

2.3. 四阶紧致差分格式

为了进一步改进上述算法的空间精度，我们引入如下的padé逼近：

$v^{″}\left(x_i\right)=\frac{{\delta }_x^2}{1+\frac{h^2}{12}{\delta }_x^2}v\left(x_i\right)+O\left({\tau }^2+h^4\right)$ (14)

考虑方程(9)在 $k+1/2$ 时刻的值并对于时间和空间导数项分别利用中心差分和上述padé逼近,可得：

$\alpha {\delta }_t{v}_i^{k\text{+}\frac{1}{2}}=\frac{{\delta }_x^2}{1+\frac{h^2}{12}{\delta }_x^2}{v}_i^{k\text{+}\frac{1}{2}}-\lambda v_i^{k\text{+}\frac{1}{2}}+G_i^{k\text{+}\frac{1}{2}}+O\left({\tau }^2+h^4\right)$ (15)

并略去高阶无穷小项，得到如下四阶紧致差分格式：

$\begin{array}{l}\left(\frac{\alpha }{12}+\frac{\lambda \tau }{24}-\frac{r}{2}\right)v_{i-1}^{k+1}+\left(\frac{5}{6}\alpha +\frac{5\lambda \tau }{12}+r\right)v_i^{k+1}+\left(\frac{\alpha }{12}+\frac{\lambda \tau }{24}-\frac{r}{2}\right)v_{i+1}^{k+1}\\ =\left(\frac{\alpha }{12}-\frac{\lambda \tau }{24}+\frac{r}{2}\right)v_{i-1}^k+\left(\frac{5}{6}\alpha -\frac{5\lambda \tau }{12}-r\right)v_i^k+\left(\frac{\alpha }{12}-\frac{\lambda \tau }{24}+\frac{r}{2}\right)v_{i+1}^k\\ +\frac{\tau }{12}\left(G_{i-1}^{k\text{+}\frac{1}{2}}+10G_i^{k\text{+}\frac{1}{2}}+G_{i+1}^{k\text{+}\frac{1}{2}}\right),1\le i\le N-1,0\le k\le M-1.\end{array}$ (16)

3. 稳定性分析

下面我们将利用Fourier分析法，给出上述两种格式的稳定性分析。

定理1. 对于任意 $\tau \ge 0$，差分格式(13)是无条件稳定的。

证明：利用Fourier分析法，令 $v_j^n=w^n{\text{e}}^{ikjh}$，得：

$w^{n+1}=\frac{\frac{r}{2}{\text{e}}^{ikh}+\left(\alpha -\frac{\lambda \tau }{2}-r\right)+\frac{r}{2}{\text{e}}^{-ikh}}{-\frac{r}{2}{\text{e}}^{ikh}+\left(\alpha +\frac{\lambda \tau }{2}+r\right)-\frac{r}{2}{\text{e}}^{-ikh}}$

由此得到增长因子：

$Gr\left(\tau ,k\right)=\frac{\frac{r}{2}{\text{e}}^{ikh}+\left(\alpha -\frac{\lambda \tau }{2}-r\right)+\frac{r}{2}{\text{e}}^{-ikh}}{-\frac{r}{2}{\text{e}}^{ikh}+\left(\alpha +\frac{\lambda \tau }{2}+r\right)-\frac{r}{2}{\text{e}}^{-ikh}}=\frac{\left(\alpha -\frac{\lambda \tau }{2}-r\right)+r\cos \left(kh\right)}{\left(\alpha +\frac{\lambda \tau }{2}+r\right)-r\cos \left(kh\right)}=\frac{\alpha -\frac{\lambda \tau }{2}-2r{\sin }^2\frac{kh}{2}}{\alpha +\frac{\lambda \tau }{2}+2r{\sin }^2\frac{kh}{2}},$

于是：

$\left|Gr\left(\tau ,k\right)\right|\le 1\iff 4\alpha \left(\frac{\lambda \tau }{2}+2r{\sin }^2\frac{kh}{2}\right)\ge 0.$

定理2. 对于任意 $\tau \ge 0$，差分格式(22)是无条件稳定的。

证明：利用Fourier分析法，令 $v_j^n=w^n{\text{e}}^{ikjh}$，得：

$w^{n+1}=\frac{\left(\frac{\alpha }{12}-\frac{\lambda \tau }{24}+\frac{r}{2}\right){\text{e}}^{ikh}+\left(\frac{5}{6}\alpha -\frac{5\lambda \tau }{12}-r\right)+\left(\frac{\alpha }{12}-\frac{\lambda \tau }{24}+\frac{r}{2}\right){\text{e}}^{-ikh}}{\left(\frac{\alpha }{12}\text{+}\frac{\lambda \tau }{24}-\frac{r}{2}\right){\text{e}}^{ikh}+\left(\frac{5}{6}\alpha +\frac{5\lambda \tau }{12}+r\right)+\left(\frac{\alpha }{12}\text{+}\frac{\lambda \tau }{24}-\frac{r}{2}\right){\text{e}}^{-ikh}},$

由此得到增长因子：

$\begin{array}{c}Gr\left(\tau ,k\right)=\frac{\left(\frac{\alpha }{12}-\frac{\lambda \tau }{24}+\frac{r}{2}\right){\text{e}}^{ikh}+\left(\frac{5}{6}\alpha -\frac{5\lambda \tau }{12}-r\right)+\left(\frac{\alpha }{12}-\frac{\lambda \tau }{24}+\frac{r}{2}\right){\text{e}}^{-ikh}}{\left(\frac{\alpha }{12}\text{+}\frac{\lambda \tau }{24}-\frac{r}{2}\right){\text{e}}^{ikh}+\left(\frac{5}{6}\alpha +\frac{5\lambda \tau }{12}+r\right)+\left(\frac{\alpha }{12}\text{+}\frac{\lambda \tau }{24}-\frac{r}{2}\right){\text{e}}^{-ikh}}\\ =\frac{\left(\frac{5}{6}\alpha -\frac{5\lambda \tau }{12}-r\right)+\left(\frac{\alpha }{6}-\frac{\lambda \tau }{12}+r\right)\cos \left(kh\right)}{\left(\frac{5}{6}\alpha +\frac{5\lambda \tau }{12}+r\right)+\left(\frac{\alpha }{6}+\frac{\lambda \tau }{12}-r\right)\cos \left(kh\right)}\\ =\frac{\frac{\alpha }{6}\cos \left(kh\right)+\frac{5}{6}\alpha -\left[\frac{\lambda \tau }{2}+\left(2r-\frac{\lambda \tau }{6}\right){\sin }^2\frac{kh}{2}\right]}{\frac{\alpha }{6}\cos \left(kh\right)+\frac{5}{6}\alpha +\left[\frac{\lambda \tau }{2}+\left(2r-\frac{\lambda \tau }{6}\right){\sin }^2\frac{kh}{2}\right]},\end{array}$

于是：

$\begin{array}{l}\left|Gr\left(\tau ,k\right)\right|\le 1\\ \iff \left|\frac{\alpha }{6}\cos \left(kh\right)+\frac{5}{6}\alpha -\left[\frac{\lambda \tau }{2}+\left(2r-\frac{\lambda \tau }{6}\right){\sin }^2\frac{kh}{2}\right]\right|\\ \le \left|\frac{\alpha }{6}\cos \left(kh\right)+\frac{5}{6}\alpha +\left[\frac{\lambda \tau }{2}+\left(2r-\frac{\lambda \tau }{6}\right){\sin }^2\frac{kh}{2}\right]\right|\\ \iff 4\left[\frac{\alpha }{6}\cos \left(kh\right)+\frac{5}{6}\alpha \right]\left[\frac{\lambda \tau }{2}+\left(2r-\frac{\lambda \tau }{6}\right){\sin }^2\frac{kh}{2}\right]\ge 0,\end{array}$

易证 $\frac{\lambda \tau }{2}+\left(2r-\frac{\lambda \tau }{6}\right){\sin }^2\frac{kh}{2}\ge 0$，那么恒有 $\left|Gr\left(\tau ,k\right)\right|\le 1$。从而此差分格式无条件稳定。

4. 数值算例

在本章，我们通过具体的数值算例来验证两种差分格式的有效性和准确性。首先，定义如下的最大相对误差：

$\text{MRE}\left(N,M\right)=\underset{1\le i\le N-1,1\le k\le M}{\max }\left|\frac{u\left(x_i,t_k\right)-u_i^k}{u\left(x_i,t_k\right)}\right|$

其中， $u\left(x_i,t_k\right)$ 为方程的精确解， $u_i^k$ 为方程的数值解( $1\le i\le N-1,1\le k\le M$ )。

4.1. 算例一

考虑如下时间分数阶对流扩散方程：

$\{\begin{array}{l}\frac{{\partial }^{\alpha }u\left(x,t\right)}{\partial t^{\alpha }}=\frac{{\partial }^{2}u\left(x,t\right)}{\partial x^2}-2\frac{\partial u\left(x,t\right)}{\partial x}+f\left(x,t\right),\\ \begin{array}{l}u\left(x,0\right)=\sin \left(\pi x\right),0\le x\le 1,\hfill \\ u\left(0,t\right)=0,u\left(1,t\right)=0,0\le t\le 1,\hfill \end{array}\end{array}$

其中，

$\{\begin{array}{l}u\left(x,t\right)=\left(t^3+1\right)\sin \left(\pi x\right),\\ f\left(x,t\right)=\frac{\sin \left(\pi x\right)}{\Gamma \left(1-\alpha \right)}\left(\frac{3}{3-\alpha }-\frac{6}{2-\alpha }+\frac{3}{1-\alpha }\right)t^{3-\alpha }+2\pi \left(1+t^3\right)\cos \left(\pi x\right)+{\pi }^2\left(1+t^3\right)\sin \left(\pi x\right),\end{array}$

表1，表2给出了不同时空步长、不同 $\alpha$ 下的两种格式的最大相对误差。对比数据可以发现当 $\alpha$ 接近1或0时，数值解精度越高，而且在剖分数较少的情况下，格式2比格式1具有更高的精度。

Table 1. Example 1: Maximum relative error of the second-order central difference scheme

表1. 算例一：二阶中心差分格式最大相对误差

Table 2. Example 1: Maximum relative error of the fourth-order central difference scheme

表2. 算例一：四阶中心差分格式最大相对误差

图1给出了两种格式在 $\alpha$ 分别取0.3和0.8时的真实解图像和误差图像。

Figure 1. When alpha is 0.3 and 0.8, the numerical solution image of the first equation and the error image of the two formats

图1. $\alpha$ 分别取0.3和0.8时算例一方程的数值解图像和两种格式的误差图像

4.2. 算例二

考虑如下时间分数阶对流扩散方程：

$\{\begin{array}{l}\frac{{\partial }^{\alpha }u\left(x,t\right)}{\partial t^{\alpha }}=\frac{{\partial }^{2}u\left(x,t\right)}{\partial x^2}-2\frac{\partial u\left(x,t\right)}{\partial x}+f\left(x,t\right),\\ \begin{array}{l}u\left(x,0\right)=0.5x^4+1,0\le x\le 1,\hfill \\ u\left(0,t\right)=t^3+1,u\left(1,t\right)=1.5\left(t^3+1\right),0\le t\le 1,\hfill \end{array}\end{array}$

其中，

$\{\begin{array}{l}u\left(x,t\right)=\left(t^3+1\right)\left(0.5x^4+1\right),\\ f\left(x,t\right)=\frac{0.5x^4+1}{\Gamma \left(1-\alpha \right)}\left(\frac{3}{3-\alpha }-\frac{6}{2-\alpha }+\frac{3}{1-\alpha }\right)t^{3-\alpha }-\left(-\text{4}{x}^3+6x^2\right)\left(t^3+1\right),\end{array}$

表3，表4给出了不同时空步长、不同α下的两种格式的最大相对误差。通过数据对比可以得到类似例一的结论。

Table 3. Example 2: Maximum relative error of the second-order central difference scheme

表3. 算例二：二阶中心差分格式最大相对误差

Table 4. Example 2: Maximum relative error of the fourth-order central difference scheme

表4. 算例二：四阶中心差分格式最大相对误差

图2给出了两种格式在α分别取0.3和0.8时的真实解图像和误差图像。

Figure 2. When alpha is 0.3 and 0.8, the numerical solution image of the second equation and the error image of the two formats

图2. $\alpha$ 分别取0.3和0.8时算例二方程的数值解图像和两种格式的误差图像

5. 结论

基于Laplace变换和有限差分方法，本文提出了两种求解时间分数阶对流扩散方程的高效数值算法，证明了格式的无条件稳定性。数值实验表明，四阶紧致差分格式比二阶中心差分格式具有更高的计算效率和精度优势。与此同时，Laplace变换的使用限制了算法的精度，我们下一步将对此变换做进一步的改进。

基金项目

国家自然科学基金青年项目(11701196，11701197)资助。