A Class of Semi-Implicit Parallel Difference Method for Time Fractional Diffusion Equations

Lifei Wu, Jiake Sun, Xiaozhong Yang^*

Mathematics and Physics Department, North China Electric Power University, Beijing, China.

DOI: 10.4236/jamp.2020.81012                

Abstract

In this paper, we construct a class of semi-implicit difference method for time fractional diffusion equations—the group explicit (GE) difference scheme, which is a difference scheme with good parallelism constructed using Saul’yev asymmetric scheme. The stability and convergence of the GE scheme of time fractional diffusion equation are analyzed by mathematical induction. Then, the theoretical analysis is verified by numerical experiments, which shows that the GE scheme is effective for solving the time fractional diffusion equation.

Keywords

Time Fractional Diffusion Equation, Group Explicit Scheme, Stability, Parallel Computation, Numerical Experiment

1. Introduction

The fractional anomalous diffusion model has a profound physical background and rich theoretical connotation. It is widely used in the fields of fluid mechanics, signal processing and information recognition, fractal theory, etc. It has become an important tool for describing various complex mechanical behaviors [1] [2] [3]. However, the analytical solutions of fractional differential equations are mostly difficult to give explicitly. It is necessary and important to study the numerical solution of fractional differential and integral equations [4] [5] [6].

In recent years, there have been many research results on numerical solutions of fractional differential equations, such as spectral methods [7], finite element methods [8] [9]. However, difference methods are still dominated [10] [11] [12]. Zhuang and Liu [13] constructed a class of implicit difference schemes with unconditional stability and convergence for time fractional diffusion equations. Yuste [14] constructed a weighted average finite difference method for the time fractional slow diffusion problem by G-L approximation and proved its stability. Tadjeran and Meerschaert [15] constructed a numerical method with second-order precision in time and space for a class of variable value initial step fractional diffusion equations. The method is to combine the classical CN scheme with the extrapolation technique. Gao and Sun [16] gave a compact difference scheme for the slow diffusion equation, which has fourth-order accuracy in space. However, due to the memory and non-locality of fractional derivatives, the computational and storage quantities of numerical calculations of fractional differential equations are very large. When we simulate practical problems, the requirements of computational resources will be very high. However, the existing serial algorithm has a large amount of calculation and relatively low computational efficiency. We mainly study the fast numerical algorithm of fractional differential equation to improve the numerical simulation efficiency of fractional order modeling this paper.

With the rapid development of multi-core and cluster technology, parallel algorithms have become one of the mainstream technologies to improve the efficiency of numerical calculations [17] [18]. Zhang and Gu [19] proposed a piecewise implicit scheme for the integer order diffusion equations in an asymmetrical scheme, and used alternating techniques to construct multiple explicit-implicit and implicit alternating parallel methods. This kind of parallel method has been widely used in integer order evolution equations. For fractional differential equations, Gong and Bao [20] [21] performed parallel computation on the explicit difference schemes of fractional reaction-diffusion equations. The core content of their parallelization is parallel calculation of matrix and vector product, vector and vector addition. Sweilam and Moharram [22] constructed a parallel C-N scheme for time fractional parabolic equations. The core of the method is to solve the equation $Ax=b$ using the preconditioned conjugate gradient method. We do not study the parallel algorithm of equations from the perspective of numerical algebra, but based on the parallelization of traditional differential schemes for solving fractional diffusion equations numerically [23].

For the time fractional diffusion equation, the Saul’yev asymmetric scheme is given, and then the Saul’yev asymmetric scheme is used to construct the group explicit (GE) scheme of the time fractional diffusion equation with parallel nature. The mathematical induction method is used to prove that the GE of the time fractional diffusion equation is unconditionally stable and convergent. Finally, the theoretical analysis is verified by numerical experiments, which shows that the GE scheme is very effective for solving time fractional diffusion equations.

2. GE Scheme of Fractional Diffusion Equation

2.1. Fractional Diffusion Equation

We consider the initial-boundary value problem of the fractional diffusion equation

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

with the initial conditions

$u\left(x,0\right)=g\left(x\right),0\le x\le L.$ (2)

and the boundary conditions

$u\left(0,t\right)=u\left(L,t\right)=0,0\le t\le T.$ (3)

where $0<\alpha <1$.

Fractional derivative $\frac{{\partial }^{\alpha }u\left(x,t\right)}{\partial t^{\alpha }}$ is Caputo derivative

$\frac{{\partial }^{\alpha }u\left(x,t\right)}{\partial t^{\alpha }}=\frac{1}{\Gamma \left(1-\alpha \right)}{\displaystyle {\int }_0^t\frac{\partial u\left(x,\xi \right)}{\partial \xi }}\frac{\text{d}\xi }{{\left(t-\xi \right)}^{\alpha }},0<\alpha <1.$

When $\alpha =1$, Equation (1) is a well-known diffusion equation (Markovian process):

$\frac{\partial u\left(x,t\right)}{\partial t}=\frac{{\partial }^{2}u\left(x,t\right)}{\partial x^2}+f\left(x,t\right),0\le x\le L,0\le t\le T.$ (4)

2.2. Construction of GE Scheme for Fractional Diffusion Equation

Define $t_k=k\tau ,k=0,1,2,\cdots ,n$, $x_i=ih,i=0,1,2,\cdots ,m$, where $\tau =\frac{T}{n}$, $h=\frac{L}{m}$

are the grid sizes in time and space respectively. Let $u_i^k$ be the numerical approximation to $u\left(x_i,t_k\right)$. The time fractional derivative term is approximated by the following scheme:

$\begin{array}{l}\frac{{\partial }^{\alpha }u\left(x_i,t_{k+1}\right)}{\partial t^{\alpha }}\\ =\frac{1}{\Gamma \left(1-\alpha \right)}{\displaystyle \underset{j=0}{\overset{k}{\sum }}{\displaystyle {\int }_{j\tau }^{\left(j+1\right)\tau }\frac{\partial u\left(x,\xi \right)}{\partial \xi }\frac{\text{d}\xi }{{\left(t_{k+1}-\xi \right)}^{\alpha }}}}\\ \approx \frac{1}{\Gamma \left(1-\alpha \right)}{\displaystyle \underset{j=0}{\overset{k}{\sum }}\frac{u\left(x_i,t_{k+1}\right)-u\left(x_i,t_k\right)}{\tau }{\displaystyle {\int }_{j\tau }^{\left(j+1\right)\tau }\frac{\text{d}\xi }{{\left(t_{k+1}-\xi \right)}^{\alpha }}}}\\ =\frac{{\tau }^{-\alpha }}{\Gamma \left(2-\alpha \right)}\left[u\left(x_i,t_{k+1}\right)-u\left(x_i,t_k\right)\right]+\frac{{\tau }^{-\alpha }}{\Gamma \left(2-\alpha \right)}{\displaystyle \underset{j=1}{\overset{k}{\sum }}\left[b_{j}u\left(x_i,t_{k+1-j}\right)-u\left(x_i,t_{k-j}\right)\right]},\end{array}$

here $b_j={\left(j+1\right)}^{\alpha }-j^{\alpha },j=0,1,2,\cdots ,n$.

$L_{h,\tau }^{\alpha }u\left(x_i,t_{k+1}\right)=\frac{{\tau }^{-\alpha }}{\Gamma \left(2-\alpha \right)}{\displaystyle \underset{j=0}{\overset{k}{\sum }}\left[b_{j}u\left(x_i,t_{k+1-j}\right)-u\left(x_i,t_{k-j}\right)\right]},$

and $\left|\frac{{\partial }^{\alpha }u\left(x_i,t_{k+1}\right)}{\partial t^{\alpha }}-L_{h,\tau }^{\alpha }u\left(x_i,t_{k+1}\right)\right|\le C{\tau }^{2-\alpha }$, where C is a constant.

Construct the Saul’yev asymmetric scheme of Equation (1):

$L_{h,\tau }^{\alpha }u\left(x_i,t_{k+1}\right)=\frac{1}{h^2}\left(u_{i+1}^k-u_i^k-u_i^{k+1}+u_{i-1}^{k+1}\right)+f_i^{k+1}.$ (5)

$L_{h,\tau }^{\alpha }u\left(x_i,t_{k+1}\right)=\frac{1}{h^2}\left(u_{i+1}^{k+1}-u_i^{k+1}-u_i^k+u_{i-1}^k\right)+f_i^{k+1}.$ (6)

where $i=1,2,\cdots ,m-1$, $k=0,1,2,\cdots ,n-1$. Let $\mu =\frac{{\tau }^{\alpha }}{h^2}$, $r=\mu \Gamma \left(2-\alpha \right)$, the above two types of Saul’yev asymmetric schemes can be rewritten as follows:

When $k=0$,

$\left(1+r\right)u_i^1-r{u}_{i-1}^1=\left(1-r\right)u_i^0+r{u}_{i+1}^0+f_i^1,$ (7)

$-r{u}_{i+1}^1+\left(1+r\right)u_i^1=\left(1-r\right)u_i^0+r{u}_{i-1}^0+f_i^1.$ (8)

When $k>0$,

$\left(1+r\right)u_i^{k+1}-r{u}_{i-1}^{k+1}=\left(1-r-b_1\right)u_i^k+r{u}_{i+1}^k+{\displaystyle \underset{j=1}{\overset{k}{\sum }}{u}_i^{k-j}\left[b_j-b_{j+1}\right]}+b_k{u}_i^0+f_i^{k+1},$ (9)

$-r{u}_{i+1}^{k+1}+\left(1+r\right)u_i^{k+1}=\left(1-r-b_1\right)u_i^k+r{u}_{i-1}^k+{\displaystyle \underset{j=1}{\overset{k}{\sum }}{u}_i^{k-j}\left[b_j-b_{j+1}\right]}+b_k{u}_i^0+f_i^{k+1},$ (10)

The GE scheme of the time fractional diffusion equation is designed as follows: the four “○” in Figure 1 indicate the establishment of the Saul’yev asymmetric format (10) at point $\left(i,k+1\right)$, four “□” means that the Saul’yev asymmetric format (9) is established at point $\left(i+1,k+1\right)$. The grouping explicit scheme is a combination of these two asymmetric schemes, and the difference equation is the following $2\times 2$ equations.

when $k=0$,

$\{\begin{array}{l}-r{u}_{i+1}^1+\left(1+r\right)u_i^1=\left(1-r\right)u_i^0+r{u}_{i-1}^0+f_i^1,\\ \left(1+r\right)u_{i+1}^1-r{u}_i^1=\left(1-r\right)u_{i+1}^0+r{u}_{i+2}^0+f_{i+1}^1.\end{array}$ (11)

when $k>0$,

$\{\begin{array}{l}-r{u}_{i+1}^{k+1}+\left(1+r\right)u_i^{k+1}=\left(1-r\right)u_i^k+r{u}_{i-1}^k-b_1{u}_i^k+{\displaystyle \underset{j=1}{\overset{k-1}{\sum }}{u}_i^{k-j}\left[b_j-b_{j+1}\right]}+b_k{u}_i^0+f{,}_i^{k+1}\\ \left(1+r\right)u_{i+1}^{k+1}-r{u}_i^{k+1}={\left(1-r\right)}_1{u}_{i+1}^k+r{u}_{i+2}^k-b_1{u}_{i+1}^k+{\displaystyle \underset{j=1}{\overset{k-1}{\sum }}{u}_{i+1}^{k-j}\left[b_j-b_{j+1}\right]}+b_k{u}_{i+1}^0+f_{i+1}^{k+1}.\end{array}$ (12)

Simplified, when $k=0$,

$\{\begin{array}{l}{u}_i^1=\frac{1}{1+2r}\left[r\left(1+r\right)u_{i-1}^0+\left(1-r^2\right)u_i^0+r\left(1-r\right)u_{i+1}^0+r^2{u}_{i+2}^0+\left(1+r\right)f_i^1+r{f}_{i+1}^1\right],\\ u_{i+1}^1=\frac{1}{1+2r}\left[r^2{u}_{i-1}^0+r\left(1-r\right)u_i^0+\left(1-r^2\right)u_{i+1}^0+r\left(1+r\right)u_{i+2}^0+r{f}_i^1+\left(1+r\right)f_{i+1}^1\right].\end{array}$ (13)

Figure 1. Diagram of the GE scheme segmentation processing point.

when $k>0$,

$\{\begin{array}{l}{u}_i^{k+1}=\frac{1}{1+2r}\left[r\left(1+r\right)u_{i-1}^k+\left(1-r^2\right)u_i^k+r\left(1-r\right)u_{i+1}^k+r^2{u}_{i+2}^k+\left(1+r\right)h_i^k+r{h}_{i+1}^k+\left(1+r\right)f_i^1+r{f}_{i+1}^1\right],\\ u_{i+1}^{k+1}=\frac{1}{1+2r}\left[r^2{u}_{i-1}^k+r\left(1-r\right)u_i^k+\left(1-r^2\right)u_{i+1}^k+r\left(1+r\right)u_{i+2}^k+r{h}_i^k+\left(1+r\right)h_{i+1}^k+r{f}_i^1+\left(1+r\right)f_{i+1}^1\right].\end{array}$ (14)

where $h_i^k=h\left(u_i^k\right)=-b_1{u}_i^k+{\displaystyle \underset{j=1}{\overset{k-1}{\sum }}{u}_i^{k-j}\left[b_j-b_{j+1}\right]}+b_k{u}_i^0$.

From Equations (13) and (14), we can know that the value of two points $u_i^{k+1}$ and $u_{i+1}^{k+1}$ on the $k+1$ ^th time layer can be calculated explicitly by function values known in the previous k layers. Therefore, it can be seen that the method formed by Equations (9) and (10) is called a group explicit method, or GE method, and this method is easy to calculate in parallel.

According to the parity of the number of points m, the GE method has different forms.

When the number of points m is even, then $m-1$ is an odd number, there are $\left(m-2\right)/2$ GE groups, and a single point. This single point is either a right single point or a left single point.

1) GER method

It has a right single point and $\left(m-2\right)/2$ GE groups. From $\left(1,k+1\right)$ to $\left(m-2,k+1\right)$ using GE scheme $\left(m-2\right)/2$ times consecutively, the value of right single point $\left(m-1,k+1\right)$ is calculated by the asymmetric scheme (9), see figure. The GER method can be expressed as the following matrix form:

$\left(I+r{G}_1\right)U^{k+1}=\left(I-r{G}_2\right)U^k+h^k+f^{k+1}+b_1,k=1,2,3,\cdots$

$G_1={\left[\begin{array}{ccccc}G& & & & \\ & G& & & \\ & & \ddots & & \\ & & & G& \\ & & & & 1\end{array}\right]}_{\left(m-1\right)\times \left(m-1\right)}$, $G=\left[\begin{array}{cc}1& -1\\ -1& 1\end{array}\right]$, $b_1={\left[r{u}_0^k,0,\cdots ,0,r{u}_m^{k+1}\right]}^{\text{T}}$,

$h^k={\left[h_1^k,h_2^k,\cdots ,h_{m-1}^k\right]}^{\text{T}}$.

2) GEL method

It has a left single point, except that the value of the point $\left(1,k+1\right)$ is calculated by the asymmetric scheme (10), from $\left(2,k+1\right)$ to $\left(m-1,k+1\right)$ using GE scheme $\left(m-2\right)/2$ times consecutively.

$\left(I+r{G}_2\right)U^{k+1}=\left(I-r{G}_1\right)U^k+h^k+f^{k+1}+b_2,k=1,2,3,\cdots$,

where $b_2={\left[r{u}_0^{k+1},0,\cdots ,0,r{u}_m^k\right]}^{\text{T}}$.

When the number of points m is odd, and $m-1$ is even, there are $\left(m-1\right)/2$ or $\left(m-3\right)/2$ GE groups, including two single-point calculations, that is, there are both right single points and left single points.

3) GEU method

Except for $\left(m-3\right)/2$ GE groups, it has both right single points and left single points. GEU (group explicit with both ungrouped ends) scheme:

$\left(I+r{\hat{G}}_1\right)U^{k+1}=\left(I-r{\hat{G}}_2\right)U^k+h^k+f_i^{k+1}+b_3,k=1,2,3,\cdots$,

${\hat{G}}_1={\left[\begin{array}{ccccc}1& & & & \\ & G& & & \\ & & \ddots & & \\ & & & G& \\ & & & & 1\end{array}\right]}_{\left(m-1\right)\times \left(m-1\right)}$, ${\hat{G}}_2={\left[\begin{array}{ccc}G& & \\ & \ddots & \\ & & G\end{array}\right]}_{\left(m-1\right)\times \left(m-1\right)}$, $b_3={\left[r{u}_0^{k+1},0,\cdots ,0,r{u}_m^{k+1}\right]}^{\text{T}}$.

4) GEC method

It includes $\left(m-1\right)/2$ times GE scheme on pairwise points consecutively, so there is GEC (group explicit complete) method:

$\left(I+r{\hat{G}}_2\right)U^{k+1}=\left(I-r{\hat{G}}_1\right)U^k+h^k+f^{k+1}+b_4,k=1,2,3,\cdots$,

$b_4={\left[r{u}_0^k,0,\cdots ,0,r{u}_m^k\right]}^{\text{T}}$.

Using the properties of function $g\left(x\right)=x^{1-\alpha }\left(x\ge 1\right)$, we can draw the following conclusions:

$\{\begin{array}{l}1=b>b_1>b_2>\cdots \to 0,\\ {\displaystyle \underset{j=1}{\overset{k}{\sum }}{c}_j=1-b_k},\\ {\displaystyle \underset{j=1}{\overset{\infty }{\sum }}{c}_j=1,1>2-2^{1-\alpha }=c_1>c_2>c_3>\cdots \to 0.}\end{array}$

3. Stability Analysis of GE Scheme

In this section, we analyze the stability of the GE scheme. Taking GER scheme as an example, the stability of GER scheme is analyzed. Assume ${\tilde{u}}_i^k\left(i=0,1,2,\cdots ,m;k=0,1,2,\cdots ,n\right)$ is the numerical solution of GER scheme. Error ${\tilde{\epsilon }}_i^k={\tilde{u}}_i^k-u_i^k$ satisfy the following equations:

when $k=0$,

$\left(I+r{G}_2\right)E^1=\left(I-r{G}_1\right)E^0$,

when $k>0$,

$\left(I+r{G}_1\right)E^{k+1}=\left[\left(1+c_1\right)I-r{G}_2\right]E^k+c_2{E}^{k-1}+\cdots +c_k{E}^1+b_k{E}^0$,

where $E^k={\left[{\epsilon }_1^k,{\epsilon }_2^k,\cdots ,{\epsilon }_{m-1}^k\right]}^{\prime }$. Let $\left|{\epsilon }_{l_k}^k\right|=\underset{1\le i\le m-1}{\max }\left|{\epsilon }_i^k\right|$,

when $k=0$, for GE scheme (13), if $r\le 1$ then error satisfy:

$\begin{array}{c}{\Vert {\epsilon }_i^1\Vert }_{\infty }=\frac{1}{1+2r}{\Vert r\left(1+r\right){\epsilon }_{i-1}^0+\left(1-r^2\right){\epsilon }_i^0+r\left(1-r\right){\epsilon }_{i+1}^0+r^2{\epsilon }_{i+2}^0\Vert }_{\infty }\\ \le \frac{1}{1+2r}\left[r\left(1+r\right)+\left(1-r^2\right)+r\left(1-r\right)+r^2\right]{\Vert {\epsilon }_{l_0}^0\Vert }_{\infty }\\ \le {\Vert {\epsilon }_{l_0}^0\Vert }_{\infty },\end{array}$

$\begin{array}{c}{\Vert {\epsilon }_{i+1}^1\Vert }_{\infty }=\frac{1}{1+2r}{\Vert r^2{\epsilon }_{i-1}^0+r\left(1-r\right){\epsilon }_i^0+\left(1-r^2\right){\epsilon }_{i+1}^0+r\left(1+r\right){\epsilon }_{i+2}^0\Vert }_{\infty }\\ \le \frac{1}{1+2r}\left[r^2+r\left(1-r\right)+\left(1-r^2\right)+r\left(1+r\right)\right]{\Vert {\epsilon }_{l_0}^0\Vert }_{\infty }\\ \le {\Vert {\epsilon }_{l_0}^0\Vert }_{\infty }.\end{array}$

In summary, we have ${\Vert E^1\Vert }_{\infty }\le {\Vert E^0\Vert }_{\infty }$.

Assume that ${\Vert E^j\Vert }_{\infty }\le {\Vert E^0\Vert }_{\infty },j=1,2,\cdots ,k$, $\left|{\epsilon }_{l_{k+1}}^{k+1}\right|=\underset{1\le i\le m-1}{\max }\left|{\epsilon }_i^{k+1}\right|$, which error satisfy the following equations:

For GE scheme (13), if $r\le 1$ then error satisfy:

$\begin{array}{c}{\Vert {\epsilon }_i^{k+1}\Vert }_{\infty }=\frac{1}{1+2r}{\Vert r\left(1+r\right){\epsilon }_{i-1}^k+\left(1-r^2\right){\epsilon }_i^k+r\left(1-r\right){\epsilon }_{i+1}^k+r^2{\epsilon }_{i+2}^k+\left(1+r\right){\tilde{h}}_i^k+r{\tilde{h}}_{i+1}^k\Vert }_{\infty }\\ \le \frac{1}{1+2r}\left(\left(1+2r\right)+\left(1+r\right)\left(-b_1+c_2+\cdots +c_k+b_k\right)+r\left(-b_1+c_2+\cdots +c_k+b_k\right)\right){\Vert {\epsilon }_{l_0}^0\Vert }_{\infty }\\ \le \left(1-b_1+c_2+\cdots +c_k+b_k\right){\Vert {\epsilon }_{l_0}^0\Vert }_{\infty }\\ \le {\Vert {\epsilon }_{l_0}^0\Vert }_{\infty },\end{array}$

$\begin{array}{c}{\Vert {\epsilon }_{i+1}^{k+1}\Vert }_{\infty }=\frac{1}{1+2r}{\Vert r^2{\epsilon }_{i-1}^k+r\left(1-r\right){\epsilon }_i^k+\left(1-r^2\right){\epsilon }_{i+1}^k+r\left(1+r\right){\epsilon }_{i+2}^k+r{\tilde{h}}_i^k+\left(1+r\right){\tilde{h}}_{i+1}^k\Vert }_{\infty }\\ \le \frac{1}{1+2r}\left(\left(1+2r\right)+r\left(-b_1+c_2+\cdots +c_k+b_k\right)+\left(1+r\right)\left(-b_1+c_2+\cdots +c_k+b_k\right)\right){\Vert {\epsilon }_{l_0}^0\Vert }_{\infty }\\ \le \left(1-b_1+c_2+\cdots +c_k+b_k\right){\Vert {\epsilon }_{l_0}^0\Vert }_{\infty }\\ \le {\Vert {\epsilon }_{l_0}^0\Vert }_{\infty }.\end{array}$

where ${\tilde{h}}_i^k=-b_1{\epsilon }_i^k+{\displaystyle \underset{j=1}{\overset{k-1}{\sum }}{\epsilon }_i^{k-j}\left[b_j-b_{j+1}\right]}+b_k{\epsilon }_i^0$.

In summary, we have ${\Vert E^{k+1}\Vert }_{\infty }\le {\Vert E^0\Vert }_{\infty }$. Therefore, when $r\le 1$, GE scheme of the fractional diffusion equation is stable. Then, we can get Theorem 1.

Theorem 1. When $r\le 1$, GE scheme of fractional diffusion equation is stable.

4. Error Analysis of GE Scheme

Let $u\left(x_i,t_k\right)\left(i=1,2,\cdots ,m-1;k=1,2,\cdots ,n\right)$ is analytical solution of time fractional diffusion equation at the grid point $\left(x_i,t_k\right)$. Define $e_i^k=u\left(x_i,t_k\right)-u_i^k$ and $e^k={\left(e_1^k,e_2^k,\cdots ,e_{m-1}^k\right)}^{\text{T}}$, assuming no initial errors, that is $e^0=0$. The accuracy of the GE scheme is discussed below, bring $e^k$ into GE scheme, then we have:

when $k=0$,

$\left(I+r{G}_2\right)e^1=R^1.$

when $k>0$,

$\left(I+r{G}_1\right)e^{k+1}=\left(I-r{G}_2\right)e^k-b_1{e}^k+{\displaystyle \underset{j=1}{\overset{k-1}{\sum }}{e}^{k-j}\left(b_j-b_{j+1}\right)}+R^{k+1}.$

For GE scheme, firstly we analyze the accuracy of the two types of Saul’yev asymmetric schemes.

When Saul’yev asymmetric scheme $I_1$ :

when $k=0$,

$\left(1+r\right)e_i^1-r{e}_{i+1}^1=R_i^1,$

when $k>0$,

$\begin{array}{l}\left(1+r\right)e_i^{k+1}-r{e}_{i+1}^{k+1}\\ =r{e}_{i-1}^k+\left(1-r-b_1\right)e_i^k+{\displaystyle \underset{j=1}{\overset{k}{\sum }}{c}_{j+1}{e}_i^{k-j}}+R_i^{k+1}.\end{array}$

$\begin{array}{c}{R}_i^{k+1}=L_{h,\tau }^{\alpha }u\left(x_i,t_{k+1}\right)-\frac{1}{h^2}\left[u\left(x_{i+1},t_{k+1}\right)-u\left(x_i,t_{k+1}\right)-u\left(x_i,t_k\right)+u\left(x_{i-1},t_k\right)\right]\\ =u\left(x_i,t_{k+1}\right)-u\left(x_i,t_k\right)+{\displaystyle \underset{j=1}{\overset{k}{\sum }}{b}_j\left[u\left(x_i,t_{k+1-j}\right)-u\left(x_i,t_{k-j}\right)\right]}\\ -\frac{\mu \Gamma \left(2-\alpha \right)}{2}\left[u\left(x_{i+1},t_{k+1}\right)-u\left(x_i,t_{k+1}\right)-u\left(x_i,t_k\right)+u\left(x_{i-1},t_k\right)\right]\\ ={\displaystyle \underset{j=0}{\overset{k}{\sum }}{b}_j\left[u\left(x_i,t_{k+1-j}\right)-u\left(x_i,t_{k-j}\right)\right]}\\ -\frac{{\tau }^{\alpha }\Gamma \left(2-\alpha \right)}{2h^2}\left[u\left(x_{i+1},t_{k+1}\right)-u\left(x_i,t_{k+1}\right)-u\left(x_i,t_k\right)+u\left(x_{i-1},t_k\right)\right],\end{array}$

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

$\begin{array}{l}\frac{1}{h^2}\left[u\left(x_{i+1},t_{k+1}\right)-u\left(x_i,t_{k+1}\right)-u\left(x_i,t_k\right)+u\left(x_{i-1},t_k\right)\right]\\ =\frac{{\partial }^{2}u\left(x_i,t_{k+1}\right)}{\partial x^2}+\left(\frac{\tau }{h}\frac{{\partial }^{2}u\left(x_i,t_{k+1}\right)}{\partial x\partial t}-\frac{{\tau }^2}{2h}\frac{{\partial }^{3}u\left(x_i,t_{k+1}\right)}{\partial x\partial t^2}+\frac{\tau h}{6}\frac{{\partial }^{4}u\left(x_i,t_{k+1}\right)}{\partial x^3\partial t}\right)\\ -\frac{\tau }{2}\frac{{\partial }^{3}u\left(x_i,t_{k+1}\right)}{\partial x^2\partial t}+\frac{h^2}{24}\frac{{\partial }^{4}u\left(x_i,t_{k+1}\right)}{\partial x^4}+C_2\left({\tau }^{2-\alpha }+h^2\right),\end{array}$

$\begin{array}{c}{R}_i^{k+1}={\tau }^{\alpha }\Gamma \left(2-\alpha \right)\{\left[\frac{{\partial }^{\alpha }u\left(x_i,t_{k+1}\right)}{\partial t^{\alpha }}-\frac{{\partial }^{2}u\left(x_i,t_{k+1}\right)}{\partial x^2}\right]\\ -\left(\frac{\tau }{h}\frac{{\partial }^{2}u\left(x_i,t_{k+1}\right)}{\partial x\partial t}-\frac{{\tau }^2}{2h}\frac{{\partial }^{3}u\left(x_i,t_{k+1}\right)}{\partial x\partial t^2}+\frac{\tau h}{6}\frac{{\partial }^{4}u\left(x_i,t_{k+1}\right)}{\partial x^3\partial t}\right)\\ -\frac{\tau }{2}\frac{{\partial }^{3}u\left(x_i,t_{k+1}\right)}{\partial x^2\partial t}+\frac{h^2}{24}\frac{{\partial }^{4}u\left(x_i,t_{k+1}\right)}{\partial x^4}\}+C\left({\tau }^2+h^2{\tau }^{\alpha }\right).\end{array}$

When Saul’yev asymmetric scheme $I_2$ :

$\begin{array}{l}\frac{1}{h^2}\left[u\left(x_{i+1},t_k\right)-u\left(x_i,t_k\right)-u\left(x_i,t_{k+1}\right)+u\left(x_{i-1},t_{k+1}\right)\right]\\ =\frac{{\partial }^{2}u\left(x_i,t_{k+1}\right)}{\partial x^2}+\left(-\frac{\tau }{h}\frac{{\partial }^{2}u\left(x_i,t_{k+1}\right)}{\partial x\partial t}+\frac{{\tau }^2}{2h}\frac{{\partial }^{3}u\left(x_i,t_{k+1}\right)}{\partial x\partial t^2}-\frac{\tau h}{6}\frac{{\partial }^{4}u\left(x_i,t_{k+1}\right)}{\partial x^3\partial t}\right)\\ -\frac{\tau }{2}\frac{{\partial }^{3}u\left(x_i,t_{k+1}\right)}{\partial x^2\partial t}+\frac{h^2}{24}\frac{{\partial }^{4}u\left(x_i,t_{k+1}\right)}{\partial x^4}+C_2\left({\tau }^{2-\alpha }+h^2\right),\end{array}$

$\begin{array}{c}{R}_i^{k+1}={\tau }^{\alpha }\Gamma \left(2-\alpha \right)\{\left[\frac{{\partial }^{\alpha }u\left(x_i,t_{k+1}\right)}{\partial t^{\alpha }}-\frac{{\partial }^{2}u\left(x_i,t_{k+1}\right)}{\partial x^2}\right]\\ -\left(-\frac{\tau }{h}\frac{{\partial }^{2}u\left(x_i,t_{k+1}\right)}{\partial x\partial t}+\frac{{\tau }^2}{2h}\frac{{\partial }^{3}u\left(x_i,t_{k+1}\right)}{\partial x\partial t^2}-\frac{\tau h}{6}\frac{{\partial }^{4}u\left(x_i,t_{k+1}\right)}{\partial x^3\partial t}\right)\\ -\frac{\tau }{2}\frac{{\partial }^{3}u\left(x_i,t_{k+1}\right)}{\partial x^2\partial t}+\frac{h^2}{24}\frac{{\partial }^{4}u\left(x_i,t_{k+1}\right)}{\partial x^4}\}+C\left({\tau }^2+h^2{\tau }^{\alpha }\right).\end{array}$

It can be seen that there are three terms in the Saul’yev asymmetric scheme $I_1$ and $R_i^{k+1}$ in and of which have the same form but opposite signs ( $\frac{\tau }{h}\frac{{\partial }^{2}u\left(x_i,t_{k+1}\right)}{\partial x\partial t}$ etc.), so inner errorscan be eliminated by simultaneous $I_1$ and $I_2$, so there is no term of $\left(\frac{\tau }{h}\right)$, that is $R_i^{k+1}\le C\left({\tau }^2+{\tau }^{\alpha }{h}^2\right)$.

In summary, for the GE scheme we have $R_i^k\le C\left({\tau }^2+{\tau }^{\alpha }{h}^2\right),k=1,2,\cdots ,m$.

Lemma 1. ${\Vert e^k\Vert }_{\infty }\le C{b}_{k-1}^{-1}\left({\tau }^2+{\tau }^{\alpha }{h}^2\right),k=1,2,\cdots ,n$. Where ${\Vert e^k\Vert }_{\infty }=\underset{1\le i\le m-1}{\max }\left|e_i^k\right|$, C is a constant.

Proof. Using mathematical induction. When $k=1$, let ${\Vert e^1\Vert }_{\infty }=\left|e_l^1\right|=\underset{1\le i\le m-1}{\max }\left|e_i^1\right|$, We will consider $\left|e_l^1\right|$ in two cases.

When Saul’yev asymmetric scheme $I_1$ :

$\begin{array}{c}\left|e_l^1\right|\le \left(1+r\right)\left|e_l^1\right|-r\left|e_{l+1}^1\right|\\ \le \left|\left(1+r\right)e_l^1-r{e}_{l+1}^1\right|\\ =\left|R_i^1\right|\le C{b}_0^{-1}\left({\tau }^2+{\tau }^{\alpha }{h}^2\right).\end{array}$

Similarly, when Saul’yev asymmetric scheme $I_2$, we have $\left|e_l^1\right|\le C{b}_0^{-1}\left({\tau }^2+{\tau }^{\alpha }{h}^2\right)$.

Assume that ${\Vert e^j\Vert }_{\infty }\le C{b}_j^{-1}\left({\tau }^2+{\tau }^{\alpha }{h}^2\right),j=1,2,\cdots ,k$ and $\left|e_l^{k+1}\right|=\underset{1\le i\le m-1}{\max }\left|e_i^{k+1}\right|$. From the property of $b_j$ we can know that $b_j^{-1}\le b_k^{-1},j=0,1,2,\cdots ,k$.

Similarly we consider $e_l^{k+1}$ in two cases. For Saul’yev asymmetric scheme $I_1$ :

$\begin{array}{c}\left|e_l^{k+1}\right|\le \left(1+r\right)\left|e_l^{k+1}\right|-r\left|e_{l+1}^{k+1}\right|\\ \le \left|\left(1+r\right)e_l^{k+1}-r{e}_{l+1}^{k+1}\right|\\ =\left|\left(1-r-b_1\right)e_l^k+r{e}_{l-1}^k+{\displaystyle \underset{j=1}{\overset{k-1}{\sum }}{c}_{j+1}{e}_l^{k-j}}+R_l^{k+1}\right|\\ \le \left|\left(1+r-b_1\right)e_l^k+r{e}_l^k+{\displaystyle \underset{j=1}{\overset{k-1}{\sum }}{c}_{j+1}{e}_l^{k-j}}+R_l^{k+1}\right|\end{array}$

$\begin{array}{c}\le \left(1-b_1\right)\left|e_l^k\right|+{\displaystyle \underset{j=1}{\overset{k-1}{\sum }}{c}_{j+1}\left|e_l^{k-j}\right|}+\left|R_l^{k+1}\right|\\ \le \left(1-b_1\right){\Vert e^k\Vert }_{\infty }+{\displaystyle \underset{j=1}{\overset{k-1}{\sum }}{c}_{j+1}{\Vert e^{k-j}\Vert }_{\infty }}+C\left({\tau }^2+{\tau }^{\alpha }{h}^2\right)\\ \le C{b}_k^{-1}\left({\tau }^2+{\tau }^{\alpha }{h}^2\right).\end{array}$

Similarly, when Saul’yev asymmetric scheme $I_2$, we have $\left|e_l^{k+1}\right|\le C{b}_k^{-1}\left({\tau }^2+{\tau }^{\alpha }{h}^2\right)$.

So for both types of Saul’yev asymmetric schemes, we have $\left|e_l^{k+1}\right|\le C{b}_k^{-1}\left({\tau }^2+{\tau }^{\alpha }{h}^2\right)$, $k^{\alpha }{\tau }^{\alpha }\le T$ and

$\underset{k\to \infty }{\lim }\frac{b_k^{-1}}{k^{\alpha }}=\underset{k\to \infty }{\lim }\frac{k^{-\alpha }}{{\left(k+1\right)}^{1-\alpha }-k^{1-\alpha }}=\underset{k\to \infty }{\lim }\frac{k^{-1}}{{\left(1+\frac{1}{k}\right)}^{1-\alpha }-1}=\underset{k\to \infty }{\lim }\frac{k^{-1}}{\left(1-\alpha \right)k^{-1}}=\frac{1}{1-\alpha }.$

Therefore, we have

${\Vert e^k\Vert }_{\infty }\le \tilde{C}{k}^{\alpha }\left({\tau }^2+{\tau }^{\alpha }{h}^2\right)=\tilde{C}{k}^{\alpha }{\tau }^{\alpha }\left({\tau }^{2-\alpha }+h^2\right)\le \tilde{C}{T}^{\alpha }\left({\tau }^{2-\alpha }+h^2\right)$,

where C is a constant.

Theorem 3. Assume $u_i^k$ is an approximate solution calculated from the GE scheme of the fractional diffusion equation, then

$\left|u_i^k-u\left(x_i,t_k\right)\right|\le C\left({\tau }^{2-\alpha }+h^2\right)$.

5. Numerical Experiment

In this section, we will compare and analyze the analytic solutions of GE scheme and the classical implicit scheme by numerical examples. It shows that GE scheme in this paper is effective for solving time fractional diffusion equations. GE scheme can also be applied to other types of time fractional equations.

Consider the following time fractional diffusion equation [13] :

$\frac{{\partial }^{\alpha }u\left(x,t\right)}{\partial t^{\alpha }}=\frac{{\partial }^{2}u\left(x,t\right)}{\partial x^2},0\le x\le 2,0<t.$ (15)

The boundary conditions are: $u\left(0,t\right)=u\left(2,t\right)=0$, and the initial conditions are:

$u\left(x,0\right)=g\left(x\right)=\{\begin{array}{l}2x,0\le x\le 0.5,\\ \left(4-2x\right)/3,0.5\le x\le 2.\end{array}$ (16)

The function $g\left(x\right)$ indicates that the heat source is at point $x=0.5$.

Through finite sine transform and Laplace transform, we can obtain the analytical solution of the time fractional diffusion Equation (15) under the above boundary condition (16):

$u\left(x,t\right)=\frac{2}{L}{\displaystyle \underset{n=1}{\overset{\infty }{\sum }}{E}_{\alpha }\left(-a^2{n}^2{t}^{\alpha }\right)\sin \left(anx\right){\displaystyle {\int }_0^{L}g\left(r\right)\sin \left(anr\right)\text{d}r}}$, (17)

where $E_{\alpha }\left(z\right)$ is Mittag-Leffler function $E_{\alpha }\left(-z\right)={\displaystyle \underset{k=0}{\overset{\infty }{\sum }}\frac{z^k}{\Gamma \left(\alpha k+1\right)}}$, $a=\pi /L$.

Take $t=0.5$ and $\alpha =0.5$, the analytical solution is compared with the numerical solution of implicit scheme and GE scheme given in this paper. Although there is an analytical solution for this type of fractional diffusion equation, it can be seen from the form that the analytic solution is very complicated. For the convenience of calculation, only the first 20 terms are taken for the analytical solution formula (3). For the numerical solution, only 40 and 1000 of space and time are considered. Compare and analysis of analytical solution and two numerical solutions are shown in Table 1, whose curves are show in Figure 2.

It can be seen from Figure 2 that the surface of GE scheme solution is smooth, which is same as the surfaces of analytical solution and the implicit scheme solution, which shows that GE scheme is feasible to solve the fractional diffusion equation. In terms of accuracy, it can be seen from Table 1 that the difference scheme solution is very close to the analytical solution (in the case of the first 20 terms). It can be seen that difference scheme is effective for solving time fractional diffusion equation. In terms of calculation time (CPU time), it can be seen that when the first 20 terms of the analytical solution formula are taken, the calculation time is 21.5483 s, which shows that the calculation amount is very large. The computation time of the two difference schemes is less than 1 s. It can be seen that the difference method is effective numerical methods for solving the time fractional diffusion equation.

By comparing GE scheme of this paper with the analytical solution and the implicit difference scheme, it can be seen from Table 1 that the absolute error of the numerical solution of two schemes is between 10⁻³ - 10⁻⁴ , but the calculation amount (CPU time) of GE scheme given in this paper is only 57.1% of implicit difference scheme. Because GE schemein this paper has the property of parallel computing, compared with the implicit difference scheme, GE scheme of the fractional diffusion equation improves the computing efficiency by about 43% when the calculation accuracy is equivalent. When performing long term calculations, the advantages of parallel computing in GE scheme will be more obvious.

Table 1. Compare and analysis of two difference scheme solutions and analytical solutions.

Figure 2. Curves of solutions of analytic solution and two difference schemes.

Table 2. Order of time convergence of GE scheme ( $m=100$ ).

Table 3. Order of convergence of GE scheme ( $n=2000$ ).

The convergence of GE scheme and its definition is verified in following. Define $E_{\infty }\left(h,\tau \right)=\underset{0\le i\le M}{\max }\left|u\left(x_i,t_N\right)-u_i^N\right|$. A spatial division of 80 is selected to analyze the time convergence of GE scheme, a time division of 2000 is selected to verify the convergence of GE scheme in the spatial direction. The calculation results are shown in Table 2 and Table 3. From Table 2, we can see that GE scheme converges linearly in the time direction. From Table 3, we can see that GE scheme converges squarely in the spatial direction. The numerical results are consistent with the theoretical analysis.

6. Conclusion

In this paper, the group explicit (GE) scheme of the time fractional diffusion equation is constructed by applying the Saul’yev asymmetric scheme. We analyzed the stability and convergence of GE scheme. GE scheme has the property of parallel computing, and its computation efficiency is nearly 60% less than that of the classic implicit scheme. The numerical experimental results are consistent with the theoretical analysis. GE scheme has $2-\alpha$ order of convergence in time and second order of convergence in space. Theoretical analysis and numerical experiments show that GE scheme is effective for solving time fractional diffusion equations. Especially for long term history and large computational domain problems, the advantages of GE scheme for parallel computing will be more obvious.

Acknowledgements

The research was supported by National Science and Technology Major Special Subproject (2017ZX07101001-1) and the Fundamental Research Funds of the Central Universities (2018MS168).

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Uchaikin, V.V. (2013) Fractional Derivatives for Physicists and Engineers, Volume II: Applications. Springer, Berlin. https://doi.org/10.1007/978-3-642-33911-0
[2]	Diethelm, K. (2010) The Analysis of Fraction Differential Equations. Springer, Berlin. https://doi.org/10.1007/978-3-642-14574-2
[3]	Chen, W., Sun, H.G., Li, X.C., et al. (2010) Fractional Derivative Modeling of Mechanics and Engineering Problems. Science Press, Beijing. (In Chinese)
[4]	Guo, B.L., Pu, X.K., Huang, F.H., et al. (2015) Fractional Partial Differential Equations and Their Numerical Solutions. Science Press, Beijing. https://doi.org/10.1142/9543
[5]	Liu, F.W., Zhuang, P.H. and Liu, Q.X. (2015) Numerical Methods and Its Application of Fractional Partial Differential Equation. Science Press, Beijing. (In Chinese)
[6]	Sun, Z.Z. and Gao, G.H. (2015) Finite Difference Methods for Fractional Differential Equations. Science Press, Beijing. (In Chinese)
[7]	Li, X.J. and Xu, C.J. (2009) A Space-Time Spectral Method for the Time Fractional Diffusion Equation. SIAM Journal on Numerical Analysis, 47, 2108-2131. https://doi.org/10.1137/080718942
[8]	Zeng, F.H., Li, C.P., Liu, F.W., et al. (2013) The Use of Finite Difference/Element Approaches for Solving the Time-Fractional Sub-Diffusion Equation. SIAM Journal on Scientific Computing, 35, A2976-A3000. https://doi.org/10.1137/130910865
[9]	Neville, J.F., Xiao, J.Y. and Yan, Y.B. (2014) A Finite Element Method for Time Fractional Partial Differential Equations. Fractional Calculus & Applied Analysis, 14, 454-474.
[10]	Wang, H. and Treena, S.B. (2012) A Fast Finite Difference Method for Two-Dimensional Space-Fractional Diffusion Equations. SIAM Journal Scientific Computing, 34, A2444-A2458. https://doi.org/10.1137/12086491X
[11]	Tan, P.Y. and Zhang, X.D. (2008) Numerical Solution of Space-Time Fractional Convection-Diffusion Equation. Computational Mathematics, 30, 305-310. (In Chinese)
[12]	Liu, F.W., Zhuang, P.H., Anh, V., et al. (2007) Stability and Convergence of the Difference Methods for the Space-Time Fractional Advection-Diffusion Equation. Applied Mathematics and Computation, 191, 12-20. https://doi.org/10.1016/j.amc.2006.08.162
[13]	Zhuang, P.H. and Liu, F.W. (2006) Implicit Difference Approximation for the Time Fractional Diffusion Equation. Journal of Applied Mathematics and Computing, 22, 87-99. https://doi.org/10.1007/BF02832039
[14]	Yuste, S.B. (2006) Weighted Average Finite Difference Methods for Fractional Diffusion Equations. Journal of Computational Physics, 216, 264-274. https://doi.org/10.1016/j.jcp.2005.12.006
[15]	Tadjeran, C., Meerschaert, M.M. and Scheffler, H.P. (2006) A Second-Order Accurate Numerical Approximation for the Fraction Diffusion Equation. Journal of Computational Physics, 213, 205-123. https://doi.org/10.1016/j.jcp.2005.08.008
[16]	Gao, G.H. and Sun, Z.Z. (2011) A Compact Finite Difference Scheme for the Fractional Sub-Diffusion Equations. Journal of Computational Physics, 230, 586-595. https://doi.org/10.1016/j.jcp.2010.10.007
[17]	Petter, B. and Mitchell, L. (2000) Parallel Solution of Partial Differential Equations. Springer-Verlag, New York.
[18]	Diethelm, K. (2011) An Efficient Parallel Algorithm for the Numerical Solution of Fractional Differential Equations. Fractional Calculus & Applied Analysis, 14, 475-490. https://doi.org/10.2478/s13540-011-0029-1
[19]	Zhang, B.L., Gu, T.X. and Mo, Z.Y. (1999) Principles and Methods of Numerical Parallel Computation. National Defence Industry Press, Beijing. (In Chinese)
[20]	Gong, C.Y., Bao, W.M., Tang, G.J., et al. (2014) An Efficient Parallel Solution for Caputo Fractional Reaction-Diffusion Equation. The Journal of Supercomputing, 68, 1521-1537. https://doi.org/10.1007/s11227-014-1123-z
[21]	Gong, C.Y., Bao, W.M. and Tang, G.J. (2013) A Parallel Algorithm for the Riesz Fraction Reaction-Diffusion Equation with Explicit Finite Difference Method. Fractional Calculus & Applied Analysis, 16, 654-669. https://doi.org/10.2478/s13540-013-0041-8
[22]	Sweilam, N.H., Moharram, H., Abdel Moniem, N.K., et al. (2014) A Parallel Crank-Nicolson Finite Difference Method for Time-Fractional Parabolic Equation. Journal of Numerical Mathematics, 22, 363-382. https://doi.org/10.1515/jnma-2014-0016
[23]	Liu, W. (2012) The Actual Combat Matlab Parallel Programming. Beihang University Press, Beijing. (In Chinese)