分数阶电报方程一类有效的普遍性差分方法

A Kind of Efficient Universal Differential Methods for Fractional Telegraph Equations

摘要: 对时间分数阶电报方程构造了一类普遍性差分方法，采用傅里叶方法分析该类差分方法的稳定性和收敛性；最后，通过数值试验验证本文方法求解分数阶电报方程的有效性。选取不同θ值进行比较分析，数值结果表明当θ取0.5附近时，数值解的精度较好；表明普遍性差分方法求解时间分数阶电报方程是有效的。

Abstract: A kind of universal difference method is constructed for the time fractional telegraph equation. The stability and convergence of the difference method are analyzed by Fourier method. Finally, the effectiveness of the method for solving the fractional telegraph equation is verified by numerical experiments. The comparison of different θ values is carried out. The numerical results show that the numerical solution is better when θ is around 0.5. Therefore, it is effective to solve the time-fractional telegraph equation by the universal difference method.

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格式无条件稳定，并且普遍性差分格式离散解较好地逼近了解析解。

Figure 1. Exact surface solution

图1. 精确解曲面

取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附近值时普遍性差分方法求解时间分数阶电报方程是可行的。

Figure 2. Universal difference scheme for surface solution when θ = 0.2

图2. 当θ = 0.2时普遍性差分格式解曲面

Figure 3. Universal difference scheme for surface solution when θ = 0.4

图3. 当θ = 0.4时普遍性差分格式解曲面

Figure 4. Universal difference scheme for surface solution when θ = 0.5

图4. 当θ = 0.5时普遍性差分格式解曲面

Table 1. Error analysis of solutions of universal difference schemes (M = 200, N = 100)

表1. 普遍性差分格式解的误差分析(取M = 200, N = 100)

Figure 5. Solution of C-N (θ = 0.5) format at N = 200 and M = 50

图5. N = 200，M = 50时C-N (θ = 0.5)格式解

Figure 6. Solution of C-N (θ = 0.5) format at N = 400 and M = 200

图6. N = 400，M = 200时C-N (θ = 0.5)格式解

Figure 7. Implicit (θ = 1) scheme solutions for N = 200 and M = 50

图7. N = 200，M = 50时隐式(θ = 1)格式解

Figure 8. Implicit (θ = 1) scheme solutions for N = 400 and M = 200

图8. N = 400，M = 200时隐式(θ = 1)格式解

Figure 9. Error surface of universal difference schemes when theta takes different values

图9. 当θ取不同值时普遍性差分格式的误差曲面

5. 结论

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

基金项目

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