DOI: 10.12677/AAM.2019.89181, PDF, HTML, XML,  被引量    国家科技经费支持
作者: 吴立飞, 杨晓忠：华北电力大学数理学院，北京
关键词: 时间分数阶电报方程；普遍性差分方法；傅里叶方法；稳定性；数值试验；Time Fractional Telegraph Equation； Universal Difference Method； Fourier Method； Stability； Numerical Experiment

文章引用：吴立飞, 杨晓忠. 分数阶电报方程一类有效的普遍性差分方法[J]. 应用数学进展, 2019, 8(9): 1544-1555. https://doi.org/10.12677/AAM.2019.89181

1. 引言

随着分数阶微积分在工程和其他方面的重要作用，最近几年，更多的学者研究分数阶微积分的求解方法，其中就有很多关于分数阶电报方程的研究 [1] [2] [3] [4] 。目前最为常用的当属分离变量法，通过该方法或者与其他方法相结合可以一定条件下给出分数阶电报方程的解析解，章红梅等 [5] 、J. Chen等 [6] 、王学彬和刘发旺 [7] [8] 等人就是借助了该方法求解了不同边界条件的包括一维、二维、三维情况在内的各类时间分数阶电报方程的解析解或近似解，部分解析解还可以表示成级数形式，极大地推进了分数阶电报方程的工程应用。此外R. F. Camargo等人 [9] 借助Laplace和Fourier变换讨论了广义上的空间–时间分数阶电报方程的解的性质，但未能给出具体的解析解。

随着分数阶电报方程在各个领域中的应用不断增强，并且一般情况下解析解很难给出，因此分数阶电报方程的数值算法的研究也随之不断发展，当前已有许多学者和专家给出了不同的求解分数阶电报方程数值解的方法。综合来看，目前为止最有效、最常用的方法有差分方法、有限元方法、谱方法、同伦分析法等。Ford [10] 等采用有限差分法求解分数阶电报方程的近似数值解；文献 [11] 通过组合正交勒让德多项式与Tau方法方式来求解分数阶电报方程的数值解；文献 [12] 用同伦分析法求解空间分数电报方程的数值解。文献 [13] 则利用切比雪夫多项式求解电报方程的数值解。陈善镇 [14] 采用有限差分法构造了Riesz空间分数阶电报方程的隐式差分格式，通过数值算例验证了差分解法的可行性，同时差分格式稳定性和收敛性的证明则通过矩阵分析的方法给出。

针对时间分数阶电报方程，本文构造了一类普遍性差分方法，采用傅里叶方法分析差分格式的稳定性和收敛性；通过数值算例验证理论分析，说明普遍性差分方法求解分数阶电报方程的有效性。

2. 分数阶电报方程的普遍性差分格式

2.1. 分数阶电报方程

考虑如下时间分数阶电报方程：

$\frac{{\partial }^{2\alpha }u\left(x,t\right)}{\partial t^{2\alpha }}+\frac{{\partial }^{\alpha }u\left(x,t\right)}{\partial t^{\alpha }}=k\frac{{\partial }^{2}u\left(x,t\right)}{\partial x^2}+f\left(x,t\right),\text{0}<x<L,t>0,\frac{\text{1}}{\text{2}}<\alpha \le 1$ (1)

其初始条件为： $u\left(x,0\right)={\varphi }_0\left(x\right)$ ， $\frac{\partial u\left(x,0\right)}{\partial t}={\varphi }_1\left(x\right)$ ，边界条件为： $u\left(0,t\right)=0$ ， $u\left(L,t\right)=0.$

2.2. 普遍性差分格式的构造

为了得到分数阶电报方程的显–隐差分格式，首先对求解区域作一网格剖分：取空间步长 $h=L/M$ 和时间步长 $\tau =T/N$ ，其中M，N都是自然数；有 $x_i=ih\left(i=1,2,\cdots ,M\right),Mh=L$ ； $t_k=k\tau \left(k=1,2,\cdots ,N\right),N\tau =T$ 网格节点为 $\left(x_i,t_k\right).$

时间分数阶导数可采用L1公式 [3] ：

$\frac{{\partial }^{\alpha }u\left(x_i,t_{k+1}\right)}{\partial t^{\alpha }}=\frac{{\tau }^{-\alpha }}{\Gamma \left(2-\alpha \right)}{\displaystyle \underset{j=0}{\overset{k}{\sum }}{c}_j{\nabla }_{t}u\left(x_i,t_{k-j}\right)}+o\left({\tau }^{2-\alpha }\right)$ (2)

和

$\frac{{\partial }^{2\alpha }u\left(x_i,t_{k+1}\right)}{\partial t^{2\alpha }}=\frac{{\tau }^{-2\alpha }}{\Gamma \left(3-2\alpha \right)}{\displaystyle \underset{j=0}{\overset{k\text{-}1}{\sum }}{d}_j{\delta }_t^{2}u\left(x_i,t_{k-j}\right)}+2d_k\left(u\left(x_i,t_1\right)-u\left(x_i,t_0\right)-\tau u_t\left(x_i\right)\right)+o\left({\tau }^{3-2\alpha }\right)$ (3)

当 $j=0$ 或者 $j=k$ 时上式会出现项 $u_i^{-1}$ ，用初始边界条件去逼近 $u_i^{-1}$ ， $u_i^{-1}=u_i^1-2k{u}_1\left(x_i\right)+O\left({\tau }^3\right)$ ，误差忽略不计。

空间二阶导数采用中心差分离散：

$\frac{{\partial }^{2}u\left(x_i,t_{k+1}\right)}{\partial x^2}={\delta }_x^2\left[\theta u\left(x_i,t_{k+1}\right)+\left(1-\theta \right)u\left(x_i,t_k\right)\right]+o\left(h^2\right)$ (4)

其中 ${\delta }_x^{2}u\left(x_i,t_k\right)=u\left(x_{i-1},t_k\right)-2u\left(x_i,t_k\right)+u\left(x_{i+1},t_k\right)$ ， ${\nabla }_{t}u\left(x_i,t_k\right)=u\left(x_i,t_{k+1}\right)-u\left(x_i,t_k\right)$ ， ${\delta }_t^{2}u\left(x_i,t_k\right)=u\left(x_i,t_{k-1}\right)-2u\left(x_i,t_k\right)+u\left(x_i,t_{k+1}\right)$ ， $c_j={\left(j+1\right)}^{1-\alpha }-j^{1-\alpha }$ ， $d_j={\left(j+1\right)}^{2-2\alpha }-j^{2-2\alpha }$ ， $j=0,1,2,\cdots ,k-1.$

把式(2, 3, 4)代入分数阶电报方程(1)得到其普遍性差分格式为：

$\frac{{\tau }^{-2\alpha }}{\Gamma \left(3-2\alpha \right)}{\displaystyle \underset{j=0}{\overset{k\text{-}1}{\sum }}{c}_j{\delta }_t^{2}u\left(x_i,t_{k-j}\right)}+\frac{{\tau }^{-\alpha }}{\Gamma \left(2-\alpha \right)}{\displaystyle \underset{j=0}{\overset{k}{\sum }}{d}_j{\nabla }_{t}u\left(x_i,t_{k-j}\right)}=\frac{k}{h^2}{\delta }_x^2\left[\theta u\left(x_i,t_{k+1}\right)+\left(1-\theta \right)u\left(x_i,t_k\right)\right]+f\left(x_i,t_{k+1}\right)$ (5)

其中， $0\le \theta \le 1$ 。令 $u_i^k$ 表示函数 $u\left(x,t\right)$ 在点 $\left(x_i,t_k\right)$ 的值； $f_i^k$ 表示函数 $f\left(x,t\right)$ 在点 $\left(x_i,t_k\right)$ 的值；设 $a=\frac{{\tau }^{-2\alpha }}{\Gamma \left(3-2\alpha \right)}$ ， $b=\frac{{\tau }^{-\alpha }}{\Gamma \left(2-\alpha \right)}$ ， $r=\frac{k}{h^2}$ ，则普遍性差分格式可写为如下形式：

$a{\displaystyle \underset{j=0}{\overset{k-1}{\sum }}{d}_j{\delta }_t^2{u}_i^{k-j}}+2d_k\left[u\left(x_i,t_1\right)-u\left(x_i,t_0\right)-\tau u_t\left(x_i\right)\right]+b{\displaystyle \underset{j=0}{\overset{k-1}{\sum }}{c}_j{\nabla }_t{u}_i^{k-j}}=r{\delta }_x^2\left(\theta u_i^{k+1}+\left(1-\theta \right)u_i^k\right)+f_i^{k+1}.$

普遍性差分格式改成矩阵形式：

$\{\begin{array}{l}\left[\left(2a+b\right)E+r\theta G\right]U^1=\left[\left(2a+b\right)E-r\left(1-\theta \right)G\right]U^0+2a\tau H+F^1,\\ A{U}^{k+1}=B{U}^k+q_1{U}^{k-1}+\cdots +q_{k-2}{U}^2+q_{k-1}{U}^1+q_k{U}^0+2a{d}_k\tau H+F^{k+1}.\end{array}$ (6)

其中， $q_0=a\left(2d_0-d_1\right)+b\left(c_0-c_1\right)$ ， $q_j=a\left(-d_{j-1}+2d_j-d_{j+1}\right)+b\left(c_j-c_{j+1}\right),j=1,2,3,\cdots ,k-2$ ， $q_{k-1}=a\left(-d_{k-2}+2d_{k-1}-2d_k\right)+b\left(c_{k-1}-c_k\right)$ ， $q_k=a\left(2d_k-d_{k-1}\right)+b{c}_k$ ， $U^k={\left[u_1^k,u_2^k,\cdots ,u_{m-1}^k\right]}^{\prime }$ ， $U^0={\left[{\phi }_0\left(x_1\right),{\phi }_0\left(x_2\right),\cdots ,{\phi }_0\left(x_{m-1}\right)\right]}^{\prime }$ ， $H={\left[{\phi }_1\left(x_1\right),{\phi }_1\left(x_2\right),\cdots ,{\phi }_1\left(x_{m-1}\right)\right]}^{\prime }$ ， $F={\left[f_1,f_2,\cdots ,f_{m-1}\right]}^{\prime }$ ， $A=\left(a{d}_0+b{c}_0\right)E+r\theta G$ ， $B=\left(a\left(2d_0-d_1\right)+b\left(c_0-c_1\right)\right)E-r\left(1-\theta \right)G$ ，

$G=\left[\begin{array}{ccccc}2& -1& & & \\ -1& 2& -1& & \\ & & \ddots & \ddots & \ddots \\ & & \ddots & \ddots & -1\\ & & & -1& 2\end{array}\right].$

由矩阵A定义可知，A是严格对角占优矩阵，所以分数阶电报方程的普遍性差分格式是唯一可解的。

3. 普遍性差分格式解的稳定性和收敛性分析

3.1. 差分格式稳定性分析

假设 ${\tilde{u}}_i^n$ 是普遍性差分格式的近似解，定义误差 ${\epsilon }_i^n=u_i^n-{\tilde{u}}_i^n$ ， ${\epsilon }_0^n={\epsilon }_M^n=0$ ， ${\epsilon }^n=\left({\epsilon }_1^n,{\epsilon }_2^n,\cdots ,{\epsilon }_{M-1}^n\right)$ 。满足如下误差扰动方程：

$\{\begin{array}{l}\left(2a+b+2r\theta \right){\epsilon }_i^1-r\theta \left({\epsilon }_{i+1}^1+{\epsilon }_{i-1}^1\right)=\left(2a+b-2r\left(1-\theta \right)\right){\epsilon }_i^0+\left(1-\theta \right)r\left({\epsilon }_{i+1}^0+{\epsilon }_{i-1}^0\right),\\ \left(a+b+2r\theta \right){\epsilon }_i^{k+1}-r\theta \left({\epsilon }_{i+1}^{k+1}+{\epsilon }_{i-1}^{k+1}\right)=\left[a\left(2d_0-d_1\right)+b\left(c_0-c_1\right)-2r\left(1-\theta \right)\right]{\epsilon }_i^k\\ +r\left(1-\theta \right)\left({\epsilon }_{i+1}^k+{\epsilon }_{i-1}^k\right)+{\displaystyle \underset{j=1}{\overset{k-1}{\sum }}{q}_j{\epsilon }_i^{k-j}+q_k{\epsilon }_i^0.}\end{array}$ (7)

下面使用Fourier方法分析普遍性差分格式的稳定性。在此，定义网格函数

${\epsilon }^n\left(x\right)=\{\begin{array}{l}{\epsilon }_l^n,x\in \left(x_{l-1/2},x_{l+1/2}\right],l=1,2,3,\cdots ,M-1,\\ 0,x\in \left[L_1,h/2\right]\cup \left(L_2-h/2,L_2\right].\end{array}$

对 ${\epsilon }^n\left(x\right)$ 进行傅立叶展开：

${\epsilon }^n\left(x\right)={\displaystyle \underset{j=-\infty }{\overset{+\infty }{\sum }}{v}_j^n}{\text{e}}^{i\frac{2j\pi x}{L}}\left(i=\sqrt{-1}\right),n=0,1,2,\cdots ,N$

其中 $v_j^n=\frac{1}{L}{\displaystyle {\int }_{L_1}^{L_2}{\epsilon }^n{\text{e}}^{i\frac{2j\pi x}{L}}}\text{d}x,j=\pm 0,\pm 1,\cdots .$

定义范数： ${\Vert {\epsilon }^n\left(x\right)\Vert }^2=\sqrt{{\displaystyle \underset{j=1}{\overset{M-1}{\sum }}h{\left|{\epsilon }_j^n\right|}^2}}\left({\epsilon }_1^n,{\epsilon }_2^n,\cdots ,{\epsilon }_{M-1}^n\right)$ 。由于 ${\epsilon }_0^n={\epsilon }_M^n=0$ ，有 ${\Vert {\epsilon }^n\Vert }^2={\displaystyle {\int }_{L_1}^{L_2}{\left|{\epsilon }^n\left(x\right)\right|}^2\text{d}x}={\Vert {\epsilon }^n\left(x\right)\Vert }^2.$

从而利用Parseval定理 ${\Vert {\epsilon }^n\left(x\right)\Vert }^2=L{\displaystyle \underset{j=-\infty }{\overset{+\infty }{\sum }}{\left|v_j^n\right|}^2}$ ，得出

${\Vert {\epsilon }^n\left(x\right)\Vert }^2={\displaystyle \underset{j=1}{\overset{M-1}{\sum }}h{\left|{\epsilon }_j^n\right|}^2}=L{\displaystyle \underset{j=-\infty }{\overset{+\infty }{\sum }}{\left|v_j^n\right|}^2},n=0,1,\cdots ,N$

基于以上分析，可以假设方程误差解的形式为 ${\epsilon }_l^n=v^n{\text{e}}^{i\lambda lh}\left(\lambda =\frac{2j\pi }{L}\right).$

当 $k=1$ 时，

$\left(2a+b+2r\theta -r\theta {\text{e}}^{i\lambda h}-r\theta {\text{e}}^{-i\lambda h}\right)v^1=\left(2a+b-2r\left(1-\theta \right)+\left(1-\theta \right)r\left({\text{e}}^{i\lambda h}+{\text{e}}^{-i\lambda h}\right)\right)v^0$ (8)

因 ${\text{e}}^{\pm i\lambda h}=\cos \lambda h\pm i\sin \lambda h$ ， ${\sin }^2\frac{\lambda h}{2}=-\frac{1}{4}\left({\text{e}}^{i\lambda h}-2+{\text{e}}^{-i\lambda h}\right)$ ，可得

$\left|v^1\right|=\frac{\left|\text{2}a+b-4r\left(1-\theta \right){\sin }^2\frac{\lambda h}{2}\right|}{\left|2a+b+4r\theta {\sin }^2\lambda h\right|}\left|v^0\right|.$

由 $\frac{\left|\text{2}a+b-4r\left(1-\theta \right){\sin }^2\frac{\lambda h}{2}\right|}{\left|2a+b+4r\theta {\sin }^2\lambda h\right|}\le 1$ ，可得 $2r\left(1-2\theta \right)\le 2a+b$ ；因此，普遍性差分格式的稳定条件为：当 $0\le \theta <\frac{1}{2}$ 时， $r\le \frac{2a+b}{2\left(1-2\theta \right)}$ ；当 $\frac{1}{2}\le \theta \le 1$ 时，无条件稳定。

当 $k>1$ 时，

$\begin{array}{l}\left(a+b+2r\theta -r\theta \left({\text{e}}^{i\lambda h}+{\text{e}}^{-i\lambda h}\right)\right)v^{k+1}\\ =\left[a\left(2d_0-d_1\right)+b\left(c_0-c_1\right)-2r\left(1-\theta \right)+r\left(1-\theta \right)\left({\text{e}}^{i\lambda h}+{\text{e}}^{-i\lambda h}\right)\right]v^k+{\displaystyle \underset{j=1}{\overset{k-1}{\sum }}{q}_j{v}^{k-j}+q_k{v}^0}\end{array}$ (9)

$\begin{array}{l}\left(a+b+4r\theta {\sin }^2\lambda h\right)\left|v^{k+1}\right|\\ =\left[a\left(2d_0-d_1\right)+b\left(c_0-c_1\right)-4r\left(1-\theta \right){\sin }^2\lambda h\right]\left|v^k\right|+{\displaystyle \underset{j=1}{\overset{k-1}{\sum }}{q}_j\left|v^{k-j}\right|+q_k\left|v^0\right|}\\ \le \left\{\left|\left[a\left(2d_0-d_1\right)+b\left(c_0-c_1\right)-4r\left(1-\theta \right){\sin }^2\lambda h\right]\right|+{\displaystyle \underset{j=1}{\overset{k-1}{\sum }}\left|q_j\right|}+\left|q_k\right|\right\}\left|v^0\right|\\ \le \left\{\left|\left[a\left(2d_0-d_1\right)+b\left(c_0-c_1\right)-4r\left(1-\theta \right){\sin }^2\lambda h\right]\right|+b{c}_1+a\left(d_0-d_1\right)\right\}\left|v^0\right|\\ \le \left|\left[a+b-4r\left(1-\theta \right){\sin }^2\lambda h\right]\right|\left|v^0\right|\end{array}$

由 $\frac{\left|a+b-4r\left(1-\theta \right){\sin }^2\frac{\lambda h}{2}\right|}{\left|a+b+4r\theta {\sin }^2\lambda h\right|}\le 1$ ，可得 $2r\left(1-2\theta \right)\le 2a+b$ ；此时普遍性差分格式的稳定条件为：当 $0\le \theta <\frac{1}{2}$ 时， $r\le \frac{a+b}{2\left(1-2\theta \right)}$ ；当 $\frac{1}{2}\le \theta \le 1$ 时，无条件稳定。

所以，有如下定理：

定理2：当 $\frac{1}{2}\le \theta \le 1$ 时，分数阶电报方程的普遍性差分格式为无条件稳定的；当 $0\le \theta <\frac{1}{2}$ 且 $r\le \frac{a+b}{2\left(1-2\theta \right)}$ 时，分数阶电报方程的普遍性差分格式是稳定的。

3.2. 差分格式收敛性分析

定义精确解和数值解之间的误差为

$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}}.$

则初始条件和边界条件为

$\begin{array}{l}{e}_i^0=0,i=0,1,\cdots ,M,\\ e_0^k=e_M^k=0,k=0,1,\cdots ,N\end{array}$

下面使用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模，有

${\Vert e^k\Vert }_2^2={\displaystyle \underset{j=1}{\overset{M-1}{\sum }}h{\left|e_j^k\right|}^2}={\displaystyle \underset{l=-\infty }{\overset{\infty }{\sum }}{\left|{\eta }_k\left(l\right)\right|}^2},$ (10)

${\Vert R^k\Vert }_2^2={\displaystyle \underset{j=1}{\overset{M-1}{\sum }}h{\left|R_j^k\right|}^2}={\displaystyle \underset{l=-\infty }{\overset{\infty }{\sum }}{\left|{\xi }_k\left(l\right)\right|}^2}.$ (11)

基于上面的分析，假设 $e_j^k={\eta }_k{\text{e}}^{i\sigma jh}$ ， $R_j^k={\xi }_k{\text{e}}^{i\sigma jh}$ ，其中 $\sigma =2\pi m/L.$

$\{\begin{array}{l}\left(2a+b+2r\theta \right)e_i^1-r\theta \left(e_{i+1}^1+e_{i-1}^1\right)=\left(2a+b-2r\left(1-\theta \right)\right)e_i^0+\left(1-\theta \right)r\left(e_{i+1}^0+e_{i-1}^0\right)+R_i^1,\\ \left(a+b+2r\theta \right)e_i^{k+1}-r\theta \left(e_{i+1}^{k+1}+e_{i-1}^{k+1}\right)=\left[a\left(2d_0-d_1\right)+b\left(c_0-c_1\right)-2r\left(1-\theta \right)\right]e_i^k\\ +r\left(1-\theta \right)\left(e_{i+1}^k+e_{i-1}^k\right)+{\displaystyle \underset{j=1}{\overset{k-1}{\sum }}{q}_j{e}_i^{k-j}}+q_k{e}_i^0+R_i^{k+1}.\end{array}$ (12)

其中 $\left|R_i^k\right|\le K\left({\tau }^{2-\alpha }+h^2\right).$

引理1：设 ${\eta }_k\left(k=1,2,\cdots ,N\right)$ 是方程(10)的解，存在一个正数 $c_2$ ，使得

$\left|{\eta }_k\right|\le c_2{a}_{k-1}^{-1}{\xi }_1,k=1,2,\cdots ,N.$

证明由定义，可知 $e^0=0$ ，有 ${\eta }_0\equiv {\eta }_0\left(m\right)=0.$

另外，有

${\Vert R^k\Vert }_2\le c_1\sqrt{L}\left({\tau }^{2-\alpha }+h^2\right),$ (13)

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

$\left|{\xi }_k\right|\equiv \left|{\xi }_k\left(m\right)\right|\le c_2{\xi }_1\equiv c_2\left|{\xi }_1\left(m\right)\right|,k=1,2,\cdots ,N.$ (14)

因此有

$\left|{\xi }_k\right|\le c_2{\xi }_1,k=1,2,\cdots ,N.$

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

当 $k=1$ 时，有

$\left(2a+b+2r\theta \right)e_i^1-r\theta \left(e_{i+1}^1+e_{i-1}^1\right)=\left(2a+b-2r\left(1-\theta \right)\right)e_i^0+\left(1-\theta \right)r\left(e_{i+1}^0+e_{i-1}^0\right)+R_i^1,$

$\left|{\eta }_1\right|=\frac{1}{\left|2a+b+4r\theta {\sin }^2\lambda h\right|}{\xi }_1\le q_0^{-1}{\xi }_1.$

当 $k>1$ 时，

$\begin{array}{l}\left(a+b+2r\theta \right)e_i^{k+1}-r\theta \left(e_{i+1}^{k+1}+e_{i-1}^{k+1}\right)\\ =\left[a\left(2d_0-d_1\right)+b\left(c_0-c_1\right)-2r\left(1-\theta \right)\right]e_i^k+r\left(1-\theta \right)\left(e_{i+1}^k+e_{i-1}^k\right)+{\displaystyle \underset{j=1}{\overset{k-1}{\sum }}{q}_j{e}_i^{k-j}}+q_k{e}_i^0+R_i^{k+1}\end{array}$

$\begin{array}{l}\left(a+b+4r\theta {\sin }^2\lambda h\right)\left|{\eta }^{k+1}\right|=\left[a\left(2d_0-d_1\right)+b\left(c_0-c_1\right)-4r\left(1-\theta \right){\sin }^2\lambda h\right]\left|{\eta }^k\right|+{\displaystyle \underset{j=1}{\overset{k-1}{\sum }}{q}_j\left|{\eta }^{k-j}\right|+\left|R^k\right|}\\ \le \left\{\left|\left[a\left(2d_0-d_1\right)+b\left(c_0-c_1\right)-4r\left(1-\theta \right){\sin }^2\lambda h\right]\right|q_k^{-1}+{\displaystyle \underset{j=1}{\overset{k-1}{\sum }}\left|q_j\right|q_{k-j}^{-1}}+1\right\}{c}_2\left|{\xi }_1\right|\\ \le \left\{\left|\left[a\left(2d_0-d_1\right)+b\left(c_0-c_1\right)-4r\left(1-\theta \right){\sin }^2\lambda h\right]\right|+{\displaystyle \underset{j=1}{\overset{k-1}{\sum }}\left|q_j\right|}+q_k\right\}{c}_2{q}_k^{-1}\left|{\xi }_1\right|\\ \le \left\{\left|\left[a\left(2d_0-d_1\right)+b\left(c_0-c_1\right)-4r\left(1-\theta \right){\sin }^2\lambda h\right]\right|+b{c}_1+a\left(d_0-d_1\right)\right\}{c}_2{q}_k^{-1}\left|{\xi }_1\right|\\ \le \left|\left[a+b-4r\left(1-\theta \right){\sin }^2\lambda h\right]\right|c_2{q}_k^{-1}\left|{\xi }_1\right|\end{array}$

即

$\left|{\eta }^{k+1}\right|\le \frac{\left|a+b-4r\left(1-\theta \right){\sin }^2\lambda h\right|}{a+b+4r\theta {\sin }^2\lambda h}{c}_2{q}_k^{-1}\left|{\xi }_1\right|\le c_2{q}_k^{-1}\left|{\xi }_1\right|,$

$\underset{n\to \infty }{\lim }\frac{c_n^{-1}}{n^{\alpha }}=\underset{n\to \infty }{\lim }\frac{n^{-\alpha }}{n^{1-\alpha }-{\left(n-1\right)}^{1-\alpha }}=\underset{n\to \infty }{\lim }\frac{n^{-1}}{1-{\left(1-\frac{1}{n}\right)}^{1-\alpha }}=\frac{1}{1-\alpha },$

由 $q_k=a\left(2d_k-d_{k-1}\right)+c_k=\frac{{\tau }^{-2\alpha }}{\Gamma \left(3-2\alpha \right)}\left(2d_k-d_{k-1}\right)+\frac{{\tau }^{-\alpha }}{\Gamma \left(2-\alpha \right)}{c}_k$ ，可知 $\frac{q_k^{-1}}{n^{\alpha }}\le \frac{c_k^{-1}}{n^{\alpha }}$ ， $\underset{n\to \infty }{\lim }\frac{q_k^{-1}}{n^{\alpha }}\le \underset{n\to \infty }{\lim }\frac{c_k^{-1}}{n^{\alpha }}=\frac{1}{1-\alpha }.$

则存在正常数 $C_1$ ，使得

$\left|{\eta }^{k+1}\right|\le c_2{q}_k^{-1}\left|{\xi }_1\right|\le \frac{{\left(k\tau \right)}^{\alpha }}{\Gamma \left(2-\alpha \right)}\frac{c_k}{k^{\alpha }}{c}_2\left|{\xi }_1\right|\le \frac{T^{\alpha }}{\left(1-\alpha \right)\Gamma \left(2-\alpha \right)}{c}_2\left|{\xi }_1\right|=C_1\left|{\xi }_1\right|,k=1,2,\cdots ,N.$

其中 $C_1=\frac{T^{\alpha }}{\left(1-\alpha \right)\Gamma \left(2-\alpha \right)}{c}_2.$

同稳定性分析，由 $\frac{\left|a+b-4r\left(1-\theta \right){\sin }^2\lambda h\right|}{a+b+4r\theta {\sin }^2\lambda h}\le 1$ 可知，普遍性差分格式的收敛条件为：

当 $0\le \theta <\frac{1}{2}$ 时， $r\le \frac{a+b}{2\left(1-2\theta \right)}$ ；当 $\frac{1}{2}\le \theta \le 1$ 时，无条件收敛。

综上，可有

定理3：当 $\frac{1}{2}\le \theta \le 1$ 时，分数阶电报方程的普遍性差分格式是无条件收敛的，收敛阶是 $O\left({\tau }^{2-\alpha }+h^2\right)$ ；当 $0\le \theta <\frac{1}{2}$ 时，分数阶电报方程的普遍性差分格式收敛的条件是 $r\le \frac{a+b}{2\left(1-2\theta \right)}.$

4. 数值试验

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

例1考虑如下时间分数阶电报方程 [6]

$\frac{{\partial }^{1.6}u\left(x,t\right)}{\partial t^{1.6}}+\frac{{\partial }^{0.8}u\left(x,t\right)}{\partial t^{0.8}}=\frac{{\partial }^{2}u\left(x,t\right)}{\partial x^2}+f\left(x,t\right),\text{0}<x<1,t>0$

其中 $f\left(x,t\right)=x\left(1-x\right)\left(\frac{2}{\Gamma \left(1.4\right)}{t}^{0.4}+\frac{2}{\Gamma \left(2.2\right)}{t}^{1.2}\right)+2\left(1+t^2\right)$ ，

初边值条件： $u\left(x,0\right)=x\left(1-x\right),u_t\left(x,0\right)=0,0<x<1$ ； $u\left(0,t\right)=u_1\left(t\right)=0,u\left(L,t\right)=u_2\left(t\right)=0.$

则该方程有精确解： $u\left(x,t\right)=\left(1+t^2\right)x\left(1-x\right).$

取 $M=200$ ， $N=100$ 时，其解析解曲面见图1；此时 $r=1.0e+04$ ，令 $r_c=\frac{a+b}{2\left(1-2\theta \right)}$ ，当 $\theta =0.\text{2}$ 时， $r_c=4.57537e+03<r$ ，不满足普遍性差分格式稳定条件，其数值解呈指数增长趋势，其数值解曲面见图2；当 $\theta =0.\text{4}$ 时， $r_c=1.37261e+04>r$ ，满足普遍性差分格式稳定条件，其数值解曲面见图3~图4可以很好逼近解析解曲面、且光滑；当 $\theta \ge 0.\text{5}$ 时，普遍性差分格式无条件稳定，其解析解曲面与解析解曲面相近，且光滑，可以很好逼近解析解；与前面普遍性差分格式的稳定性分析一致。

当 $\theta \ge 0.\text{5}$ 时，普遍性差分格式无条件稳定。考察 $\theta =0.\text{5}$ 和1.0时普遍性差分格式的稳定性，分别取 $M=200$ ， $N=50$ 和 $M=400$ ， $N=200$ 两种情况，由图5~图8可知，隐式格式和C-N格式无条件稳定，并且普遍性差分格式离散解较好地逼近了解析解。

取M = 200，N = 100，差分格式解与解析解的误差情况见表1。由表1可知，当θ = 0.5时(C-N格式)的精度比θ = 1 (隐式格式)的精度要好，且数值解均小于解析解；当θ = 0.4时数值解的绝对误差与θ = 0.5时基本相同，但此时数值解大于解析解的值。因此当θ取0.4至0.5中间某值得时候精度可取的最好。如图9所以，当θ = 0.45时的误差要小于其他三种情况。由表1和图9可以看出，当θ取值在0.5附近时精度较高，对于例1而言，当θ = 0.45时普遍性差分格式的精度最好。对方程进行离散时网格划分越密集，所得的差分格式解越平滑越接近于解析解，采用普遍性差分格式，尤其是θ = 0.5附近值时普遍性差分方法求解时间分数阶电报方程是可行的。

5. 结论

本文对时间分数阶电报方程构造了一类普遍性差分方法，分析了该种差分方法解的存在唯一性，稳定性和收敛性，最后进行了数值实验；对比解析解与数值解，结果表明该类方法尤其是当θ取0.5附近时，对数值求解时间分数阶电报方程是有效的、可行的。

基金项目

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

参考文献

[1]	郭柏林, 蒲学科, 黄凤辉. 分数阶偏微分方程及其数值解[M]. 北京: 科学出版社, 2001.
[2]	周激流. 分数阶微积分原理及其在现代信号分析与处理中的应用[M]. 北京: 科学出版社, 2010.
[3]	孙志忠, 高广花. 分数阶微分方程的有限差分方法[M]. 北京: 科学出版社, 2015.
[4]	刘发旺, 庄平辉, 刘青霞. 分数阶偏微分方程数值方法及其应用[M]. 北京:科学出版社, 2015.
[5]	章红梅, 刘发旺. 时间分数阶电报方程的一种解技巧[J]. 厦门大学学报: 自然科学版, 2006, 46(1): 10-13.
[6]	Chen, J., Liu, F. and Anh, V. (2008) Analytical Solution for the Time-Fractional Telegraph Equation by the Method of Separating Variables. Journal of Mathematical Analysis and Applications, 338, 1364-1377.
[7]	王学彬, 刘发旺. 二维和三维的时间分数阶电报方程的解析解[J]. 山东大学学报(理学版), 2012, 47(8): 114-121.
[8]	王学彬, 刘发旺. Dirich1et边界条件下时间分数阶电报方程的解析解[J]. 山东大学学报: 自然科学版, 2012, 47(8): 1-9.
[9]	Camargo, R.F., Chiacchio, A.O. and Oliveira, E.C. (2008) Dif-ferentiation to Fractional Orders and the Fractional Telegraph Equation. Journal of Mathematical Physics, 49, Article ID: 033505. https://doi.org/10.1063/1.2890375
[10]	Ford, N.J., Rodrigues, M.M., Xiao, J. and Yan, Y. (2013) Numerical Analysis of a Two-Parameter Fractional Telegraph Equation. Journal of Computational and Applied Mathematics, 249, 95-106. https://doi.org/10.1016/j.cam.2013.02.009
[11]	Saadatmandi, A. and Mohabbati, M. (2015) Numerical Solution of Fractional Telegraph Equation via the Tau Method. Mathematical Reports, 17, 155-166.
[12]	Kumar, S. (2014) A New Analytical Modelling for Fractional Telegraph Equation via Laplace Transform. Applied Mathematical Modelling, 38, 3154-3163. https://doi.org/10.1016/j.apm.2013.11.035
[13]	牛变玲, 李灯熬, 赵富强, 解加全. 分数阶电报方程的Chebyshev多项式数值解法研究[J]. 工程数学学报, 2018, 35(1): 79-87.
[14]	陈善镇. 两类空间分数阶偏微分方程模型有限差分逼近的若干研究[D]: [博士学位论文]. 济南: 山东大学数学学院, 2015.