Alternating Segment Explicit-Implicit and Implicit-Explicit Parallel Difference Method for Time Fractional Sub-Diffusion Equation

Lifei Wu, Yadi Zhao, Xiaozhong Yang^*

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

DOI: 10.4236/jamp.2018.65089                

Abstract

The fractional diffusion equations can accurately describe the migration process of anomalous diffusion, which are widely applied in the field of natural science and engineering calculations. This paper proposed a kind of numerical methods with parallel nature which were the alternating segment explicit-implicit (ASE-I) and implicit-explicit (ASI-E) difference method for the time fractional sub-diffusion equation. It is based on the combination of the explicit scheme, implicit scheme, improved Saul’yev asymmetric scheme and the alternating segment technique. Theoretical analyses have shown that the solution of ASE-I (ASI-E) scheme is uniquely solvable. At the same time the stability and convergence of the two schemes were proved by the mathematical induction. The theoretical analyses are verified by numerical experiments. Meanwhile the ASE-I (ASI-E) scheme has the higher computational efficiency compared with the implicit scheme. Therefore it is feasible to use the parallel difference schemes for solving the time fractional diffusion equation.

Keywords

Time Fractional Diffusion Equation, ASE-I, ASI-E, Stability, Parallel Computing

1. Introduction

Fractional differential equations arise from some anomalous diffusion models and can be very useful in describing the memory and heritability of various complex substances. Because of its deep physical background and rich theoretical significance, it has been widely used in various fields such as fluid mechanics, signal processing and information recognition [1] [2] [3] . Due to the success in the analysis of a discrete non-Markovian Random Walk Approximation for the anomalous diffusion and its close connection with fractional calculation, the anomalous diffusion has become a more interesting research direction in the field of complex systems [4] . The time fractional diffusion equation is a class of diffusion equations where the first-order time derivative is replaced by the time fractional order derivative α. Due to the numerical calculation and storage capacity of the analytic solution of such equations are very large, the studying on numerical algorithm for solving the models has become one of the main subjects in anomalous diffusion [5] [6] [7] [8] .

In recent years, there have been many research achievements on the numerical algorithm of fractional diffusion equations. Xu Chuanju et al. applied the spectral method directly to the solution of time fractional derivatives, and proved the convergence by providing a priori error estimate [9] . Liu Fawang et al. proposed the implicit RBF method for solving the time fractional diffusion equation. Since no grid meshing is consistent with the definition of fractional derivatives, it has been regarded as a promising research direction [10] . However in the existing numerical algorithms, the finite difference method is still dominant. A class of unconditionally stable and convergent implicit difference approximate method for fractional diffusion equation was constructed by Zhuang Pinghui et al., which was second-order in space and $2-\alpha$ order in time [11] . Yuste established the forward Euler difference scheme and put forward the weighted average finite difference scheme by G-L (Grunwald-Letnikov) approximation method for the time fractional sub-diffusion model, as well as proved the stability of the scheme [12] . Pu Hai et al. examined a C-N (Crank-Nicolson) difference method to solve a class of sub-diffusion equations with variable coefficients [13] . The approach was proved unconditionally stable by the energy method and convergent to temporally $\min \{2-\alpha /2,1+\alpha \}$ order and spatially second order. Gao Guanghua et al. derived a compact finite difference scheme for the sub-diffusion equation, which is fourth-order accuracy approximation for the space derivative [14] .

However the computational complexity of existing serial algorithms is relatively high and the computation efficiency is low. With the rapid development of multi core and cluster technology, the parallel algorithm of diffusion equation has also widely used in numerical calculations [15] [16] . Zhang Baolin et al. proposed the idea of using Saul’yev asymmetric scheme to construct segment implicit scheme, and using the alternative technique to establish a variety of explicit-implicit and implicit alternative parallel methods [17] . Yuan Guangwei et al. put forward an efficient parallel method which not only keeps the conservation of the implicit scheme, but also maintains the required accuracy and unconditional stability by taking the prediction correction method [18] . In this paper, we will study the application of the parallel difference method which was used in the integer-order equations for solving the time fractional order diffusion equation [19] [20] [21] .

The structure of this paper is arranged as follow. In Section 2, the alternating segment explicit-implicit (ASE-I) parallel difference method is constructed. The unconditional stability and convergence are analyzed. In Section 3, we give the alternating segment implicit-explicit (ASI-E) parallel difference method. In Section 4, numerical experiments are presented to support our theoretical analysis and indicate that the ASE-I (ASI-E) scheme is effective for solving time fractional sub-diffusion equation.

2. ASE-I Parallel Difference Method

2.1. Time Fractional Sub-Diffusion Equation

The fractional sub-diffusion equation is considered as follows

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

Initial boundary conditions:

$u\left(0,t\right)=u\left(L,t\right)=0,u\left(x,0\right)=f\left(x\right).$

$u\left(x,t\right)$ indicates the diffusion concentration in point x at time t, fractional order derivative α in Equation(1) is Caputo fractional derivative defined by

$\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$ (2)

By taking the finite sine transform and Laplace transform, the exact solution for the Equation (1) with the boundary conditions as above is obtained as Equation (3)

$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^{t}f\left(r\right)\sin \left(anr\right)\text{d}r}}$ (3)

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=\frac{\text{π}}{L}$ .

2.2. Construction of ASE-I Scheme

Define $t_k=k\tau ,k=0,1,\cdots ,n$ , $x_i=ih,i=0,1,\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 $L_{\text{1}}$ interpolation approximation of the time fractional derivative is defined as follows

$\begin{array}{c}\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 }}}}\\ =\frac{{\tau }^{1-\alpha }}{\Gamma \left(2-\alpha \right)}{\displaystyle \underset{j=0}{\overset{k}{\sum }}\frac{u\left(x_i,t_{k+1-j}\right)-u\left(x_i,t_{k-j}\right)}{\tau }\left[{\left(j+1\right)}^{1-\alpha }-j^{1-\alpha }\right]}\end{array}$

Let $b_j={\left(j+1\right)}^{1-\alpha }-j^{1-\alpha },j=0,1,2,\cdots ,n$ , the approximate form can be rewritten as

$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 }}{b}_j\left[u\left(x_i,t_{k+1-j}\right)-u\left(x_i,t_{k-j}\right)\right]}$ (4)

Moreover the spatial derivative can be discretized in the following four schemes: first, the classical implicit scheme is

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

Second, the classical explicit scheme is

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

At last, we present the two improved Saul’yev asymmetric schemes

$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)$ , (7)

$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)$ . (8)

Substituting Equation (4-8) into Equation (1), let $\mu ={\tau }^{\alpha }/h^2,r=\mu \Gamma \left(2-\alpha \right)$ , we can respectively derive the four schemes of Equation (1):

When $k=0$ ,

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

When $k>0$ ,

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

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

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

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

As the scheme we constructed above, the classical implicit scheme (5) is absolutely stable but it is inconvenience to efficiently obtain the results because of needing to solve three diagonal matrixes. The classical explicit scheme (6) has ideal parallelism but it is conditionally stable. The improved Saul’yev asymmetric schemes (7) and (8) are convenient to parallel computing, but are conditionally stable. So the ASE-I scheme which we constructed is combined with the advantages of the above schemes and the design is as follows

Let $m-1=Bl$ , here B is a positive odd number, w is a positive integer, $B\ge 3,l\ge 3$ . We divide the points on each time level into B sections, order as $S_1,S_2,\cdots ,S_B$ . And on the even level, we arrange the computation according to the rule of “the explicit segment-the implicit segment-the explicit segment”. When it turns to the odd level, the rule changes into “the implicit segment-the explicit segment-the implicit segment” that makes the implicit segment and the explicit segment doing alternatively at different time level. In this format, for $i_0\ge 0,i_0=l,2l,\cdots ,\left(B-2\right)l$ , we consider the calculation of the implicit segment point $\left(i_0+i,k+l\right),i=1,2,\cdots ,l$ . The left boundary point $\left(i_0+1,k+1\right)$ of the implicit segment is calculated with the improved Saul’yev scheme (7), the right boundary point $\left(i_0+l,k+1\right)$ is calculated with the improved Saul’yev scheme (8), and the “interior point” $\left(i_0+i,k+1\right),i=2,3,\cdots ,l-1$ are calculated with the classical implicit scheme (5). See Figure 1, indicates the classical implicit scheme, indicates the improved Saul’yev schemes, thus we get the following implicit segment Equation (13).

$\begin{array}{l}\left[\begin{array}{ccccc}1+r& -r& & & \\ -r& 1+2r& -r& & \\ & \ddots & \ddots & \ddots & \\ & & -r& 1+2r& -r\\ & & & -r& 1+r\end{array}\right]\left[\begin{array}{c}{u}_{i_0+1}^{k+1}\\ u_{i_0+2}^{k+1}\\ ⋮\\ u_{i_0+l-1}^{k+1}\\ u_{i_0+l}^{k+1}\end{array}\right]\\ =\left[\begin{array}{c}\left(1-r-b_1\right)u_{i_0+1}^k+r{u}_{i_0}^k\\ \left(1-b_1\right)u_{i_0+2}^k\\ ⋮\\ \left(1-b_1\right)u_{i_0+l-1}^k\\ \left(1-r-b_1\right)u_{i_0+l}^k+r{u}_{i_0+l+1}^k\end{array}\right]+\left[\begin{array}{c}{\displaystyle \underset{j=1}{\overset{k-1}{\sum }}\left(b_j-b_{j+1}\right)u_{i_0+1}^{k-j}}+b_k{u}_{i_0+1}^0\\ {\displaystyle \underset{j=1}{\overset{k-1}{\sum }}\left(b_j-b_{j+1}\right)u_{i_0+2}^{k-j}}+b_k{u}_{i_0+2}^0\\ ⋮\\ {\displaystyle \underset{j=1}{\overset{k-1}{\sum }}\left(b_j-b_{j+1}\right)u_{i_0+l-1}^{k-j}}+b_k{u}_{i_0+l-1}^0\\ {\displaystyle \underset{j=1}{\overset{k-1}{\sum }}\left(b_j-b_{j+1}\right)u_{i_0+l}^{k-j}}+b_k{u}_{i_0+l}^0\end{array}\right].\end{array}$ (13)

In order to improve the calculation accuracy, the left and right boundary point of implicit segment will be replaced by the implicit scheme when $i_0=0$ and $i_0=\left(B-1\right)l$ .

At the same time the explicit segment scheme is

$\begin{array}{c}\left[\begin{array}{c}{u}_{i_0+1}^{k+1}\\ u_{i_0+2}^{k+1}\\ ⋮\\ u_{i_0+l-1}^{k+1}\\ u_{i_0+l}^{k+1}\end{array}\right]=\left[\begin{array}{ccccc}1-2r-b_1& r& & & \\ r& 1-2r-b_1& r& & \\ & \ddots & \ddots & \ddots & \\ & & r& 1-2r-b_1& r\\ & & & r& 1-2r-b_1\end{array}\right]\left[\begin{array}{c}{u}_{i_0+1}^k\\ u_{i_0+2}^k\\ ⋮\\ u_{i_0+l-1}^k\\ u_{i_0+l}^k\end{array}\right]\\ +\left[\begin{array}{c}r{u}_{i_0}^k\\ 0\\ ⋮\\ 0\\ r{u}_{i_0+l+1}^k\end{array}\right]+\left[\begin{array}{c}{\displaystyle \underset{j=1}{\overset{k-1}{\sum }}\left(b_j-b_{j+1}\right)u_{i_0+1}^{k-j}}+b_k{u}_{i_0+1}^0\\ {\displaystyle \underset{j=1}{\overset{k-1}{\sum }}\left(b_j-b_{j+1}\right)u_{i_0+2}^{k-j}}+b_k{u}_{i_0+2}^0\\ ⋮\\ {\displaystyle \underset{j=1}{\overset{k-1}{\sum }}\left(b_j-b_{j+1}\right)u_{i_0+l-1}^{k-j}}+b_k{u}_{i_0+l-1}^0\\ {\displaystyle \underset{j=1}{\overset{k-1}{\sum }}\left(b_j-b_{j+1}\right)u_{i_0+l}^{k-j}}+b_k{u}_{i_0+l}^0\end{array}\right].\end{array}$ (14)

Figure 1. Schematic of segment implicit.

We use ○ to denote the classical explicit scheme, ● to denote the classical implicit scheme, the remainder are two improved asymmetric formats.

Let $m=26,l=5,B=5$ and give the schematic of the ASE-I scheme (see Figure 2).

Thus the ASE-I scheme can be written as follows

$\{\begin{array}{l}\left(I+r{G}_1\right)U^{k+1}=\left(c_{1}I-r{G}_2\right)U^k+b^k+{\displaystyle \underset{j=2}{\overset{k}{\sum }}{c}_j{U}^{k-j+1}+b_k{U}^0}\\ \left(I+r{G}_2\right)U^{k+2}=\left(c_{1}I-r{G}_1\right)U^{k+1}+b^{k+2}+{\displaystyle \underset{j=2}{\overset{k+1}{\sum }}{c}_j{U}^{k-j+2}+b_{k+1}{U}^0}\end{array},k=0,2,\cdots$ (15)

where $b^k={\left[r{u}_{\text{0}}^k,0,\cdots ,0,r{u}_m^k\right]}^{\text{T}}$ , $c_j=b_{j-1}-b_j$ , $U^k={\left[u_1^k,u_2^k,\cdots ,u_{m-1}^k\right]}^{\text{T}}$ for any non-negative k was established, $U^0=f={\left[f\left(x_{\text{1}}\right),\cdots ,f\left(x_{m-1}\right)\right]}^{\text{T}}$ ,

$G_1={\left[\begin{array}{cccccccc}{Q}_l& & & & & & & \\ & G_l^{\left(1\right)}& & & & & & \\ & & Q_l& & & & & \\ & & & G_l^{\left(2\right)}& & & & \\ & & & & \ddots & & & \\ & & & & & Q_l& & \\ & & & & & & G_l^{\left(\frac{B-1}{2}\right)}& \\ & & & & & & & Q_l\end{array}\right]}_{\left(m-1\right)\times \left(m-1\right)}$ ,

$G_2={\left[\begin{array}{cccccccc}{\tilde{G}}_{l+1}^{\left(1\right)}& & & & & & & \\ & Q_{l-2}& & & & & & \\ & & G_{l+2}^{\left(2\right)}& & & & & \\ & & & Q_{l-2}& & & & \\ & & & & \ddots & & & \\ & & & & & G_{l+2}^{\left(\frac{B-1}{2}\right)}& & \\ & & & & & & Q_{l-2}& \\ & & & & & & & {\bar{G}}_{l+1}^{\left(\frac{B+1}{2}\right)}\end{array}\right]}_{\left(m-1\right)\times \left(m-1\right)}$ ,

${\tilde{G}}_{l+1}^{\left(1\right)}={\left[\begin{array}{ccccc}2& -1& & & \\ -1& 2& -1& & \\ & \ddots & \ddots & \ddots & \\ & & -1& 2& -1\\ & & & -1& 1\end{array}\right]}_{\left(l+1\right)\times \left(l+1\right)}$ , ${\tilde{G}}_{l+1}^{\left(\frac{B+1}{2}\right)}={\left[\begin{array}{ccccc}1& -1& & & \\ -1& 2& -1& & \\ & \ddots & \ddots & \ddots & \\ & & -1& 2& -1\\ & & & -1& 2\end{array}\right]}_{\left(l+1\right)\times \left(l+1\right)}$ ,

Figure 2. Schematic diagram of the ASE-I scheme.

$G_{l^{*}}^{\left(i\right)}={\left[\begin{array}{ccccc}1& -1& & & \\ -1& 2& -1& & \\ & \ddots & \ddots & \ddots & \\ & & -1& 2& -1\\ & & & -1& 1\end{array}\right]}_{l^{*}\times l^{*}}\left(l^{*}=l\text{or}l+2,i=1,2,\cdots ,\frac{B-1}{2}\right),$

$Q_{l^{″}}$ is a zero matrix with $l^{″}\times l^{″}$ order ( $l^{″}=l$ or $l-2$ ).

Using the properties of the function $g\left(x\right)=x^{1-\alpha }\left(x\ge 1\right)$ , a set of conclusions can be obtained:

$\{\begin{array}{l}1=b_0>b_1>\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>\cdots \to 0}\end{array}$ (16)

2.3. Existence and Uniqueness of ASE-I Scheme Solution

Lemma 1. [17] Set $\rho >0$ , if matrix A is a non-negative real matrix, the matrix

${\left(I+\rho A\right)}^{-1}$ exists, and ${\Vert {\left(I+\rho A\right)}^{-1}\Vert }_2\le 1$ .

Lemma 2. The matrix $G_1$ and $G_2$ of ASE-I scheme are non-negative real matrix.

Proof. We only need to prove $G_1+G_1^{\text{T}}$ and $G_2+G_2^{\text{T}}$ are non-negative matrices. Obviously,

$G_{l^{*}}^{\left(i\right)}+{\left(G_{l^{*}}^{\left(i\right)}\right)}^{\text{T}}={\left[\begin{array}{ccccc}2& -2& & & \\ -2& 4& -2& & \\ & \ddots & \ddots & \ddots & \\ & & -2& 4& -2\\ & & & -2& 2\end{array}\right]}_{l^{*}\times l^{*}}{l}^{*}=l\text{or}\left(l+2\right)$

are diagonally dominant matrices and their diagonal elements are non-negative

real numbers, therefore $G_{l^{*}}^{\left(i\right)}+{\left(G_{l^{*}}^{\left(i\right)}\right)}^{\text{T}}$ are non-negative matrices. Similarly we can get ${\tilde{G}}_{l+1}^{\left(1\right)}+{\left({\tilde{G}}_{l+1}^{\left(1\right)}\right)}^{\text{T}}$ and ${\tilde{G}}_{l+1}^{\left(\frac{B+1}{2}\right)}+{\left({\tilde{G}}_{l+1}^{\left(\frac{B+1}{2}\right)}\right)}^{\text{T}}$ are also non-negative matrices.

According to Equation (15), ${\left(I+r{G}_1\right)}^{-1}$ and exists combing with Lemma 1 and Lemma 2. Then we can get

Theorem 1. The solution of the ASE-I scheme for solving time fractional sub-diffusion equation is uniquely solvable.

2.4. Stability of ASE-I Scheme

Lemma 3. Set, if A matrix is non-negative real matrix, we can get.

Proof. Because of

thus.

The growth matrix of ASE-I scheme for time fractional sub-diffusion equation is. Let

we can easily obtain by Lemma 3.

Suppose that is the solution of ASE-I scheme， is the approximate solution of the scheme，the error satisfies

(17)

in which. Let are the eigenvalues of respectively and the two matrices have the same eigenvalues. Hence

, , we can get

For,.

For, ,

,

Suppose that when, then we also have

Summing up we have. Hence, the following theorem is obtained.

Theorem 2. The ASE-I scheme for time fractional sub-diffusion equation is unconditionally stable.

2.5. Convergence of ASE-I Scheme

Because of

where,

hence we can get

At the same time the truncation error of improved Saul’yev asymmetric schemes used in spatial discretization is

, respectively.

When the two schemes are constructed alternatively on different time layers,

the truncation error is. In the same way, the truncation error of the explicit-implicit scheme (5-6) is

as well. Due to the first term of and can be cut off by the error term of temporal discretization, therefore the truncation error of ASE-I scheme is.

Define in which is exact solution of the Equation (1) and. Using to substitute into the Equation (15), we can get

in which

, C is a constant.

Lemma 4., , is a positive constant.

Proof . Lemma 4 can be proved using mathematical induction.

When,.

When,

Suppose that, then we also have

In summary, we have.

Due to, thus

where and T is a positive constant.

Theorem 3. The ASE-I scheme for time fractional sub-diffusion equation is unconditionally convergent and there is a positive A satisfies .

3. ASI-E Parallel Difference Method

Imitating the method constructed ASE-I scheme, we give the ASI-E scheme for solving the time fractional sub-diffusion equation. The difference between the ASE-I and ASI-E scheme is that the use of implicit segment and explicit segment is different.

On the odd level, we arrange the computation according to the rule of “the implicit segment-the explicit segment-the implicit segment”, when it turns to the even level, the rule changes into “the explicit segment-the implicit segment-the explicit segment”. Thus we get the ASI-E difference scheme

(17)

in which the definition of and are the same as above. Due to the implicit scheme on the first layer is unconditionally stable and convergent, we imitate the analytical and proved method of the ASE-I scheme (15) from the second time layer, and get the following theorem.

Theorem 4. The ASI-E scheme for time fractional sub-diffusion equation is unconditionally stable and convergent, meanwhile there is a positive A satisfies.

4. Numerical Examples

In this section, we present numerical examples to demonstrate that the ASE-I scheme is a computational effective numerical method for time fractional sub-diffusion equation compared with the implicit scheme as well as give the convergence rate of the ASE-I scheme. Numerical experiments will be done in MatlabR2015b, based on the Intel Core i5-2400 CPU@2.20GHz.

Considering the following time fractional sub-diffusion equation

The initial condition is

At, we compare the solution of ASE-I scheme with the exact solution and the numerical solution using the implicit scheme. For the exact solution, the series in Equation (3) is truncated after 20 terms. We take, when calculating numerical solutions, the computed results are listed in Table 1.

As these can be seen from Table 1, the numerical solutions of ASE-I and ASI-E scheme are better close to the exact solution compared with the implicit numerical solution, and the result of ASE-I scheme is obviously better, thus we mainly focus on ASE-I scheme. The surface of ASE-I numerical solution shown in Figure 3 describes the complete diffusion process which instructs the variations in concentration at different times and spaces. In particular, we consider the decay process curve of fractional diffusion model at the space point under the case of α taking different values. From Figure 4, it can be seen that the speed of diffusion is getting faster as α approaches to the number “1”, and the

Figure 3. Numerical solution surface of ASE-I scheme.

Table 1. Comparison of exact solution and numerical solutions.

Figure 4. Decay process curve of diffusion model.

diffusion velocity of solute becomes more and more slow with the reduction of the diffusion concentrations, this is consistent with measurements of some actual diffusion processes.

Next, to better validate the stability and compare the accuracy of the ASE-I, we will analyze the change cure of the sum of relative error with time steps (SRET) and the distribution of the difference total energy (DTE) at space grid points. Taking the exact solution as the control solution, we let the numerical solution of the scheme as the perturbation solutions. The definition of SRET and DTE are as follows:

From Figure 5, the SRET of ASE-I scheme is less than 5. The relative error is a little big in the first few steps, and decreases rapidly with the time step, thus we can know that the ASE-I scheme of the time fractional sub-diffusion equation is stable.

The values of ASE-I scheme’s DTE are between 0 and 0.0075 from Figure 6, and the calculation error of the diffusion concentration is getting reductive with the passage of space, this also can demonstrate that the ASE-I scheme of time fractional sub-diffusion equation is very close to the exact solution. The values of DTE appear to fluctuate near the grids 16, 32, 48 64, and its maximum values appear near the grids 16 and 32. Fortunately these grids are the “inter boundary point” of the ASE-I scheme, i.e. a couple of Saul’yev scheme alternatively applied in different temporal level. At the same time, the explicit-implicit scheme is applied in the “inter point” for the ASE-I scheme. So it is normal to conclusion that the values of DTE in “inter boundary point” are little bigger than that in “inter point”. Moreover we compare the DTE of the implicit scheme and ASE-I

Figure 5. Change curve of SRET.

Figure 6. Distribution of DTE at space points.

scheme. As shown in Figure 6, there is not much difference between the two schemes and the error of the ASE-I scheme is slightly smaller than the other one. Comprehensively considering the ASE-I scheme can be more effective to solve the time fractional sub-diffusion equation.

A test example will be performed to illustrate the convergence order of the ASE-I scheme. Denote [22]

Thus the numerical results are presented as follows.

Table 2 gives the computational errors with different temporal step sizes using the fractional order. We can see that the numerical accuracy in temporal direction is approximately order and compared with the implicit scheme, it has the higher accuracy. At the same time, we compute the numerical accuracy in spatial direction. Taking for Implicit and ASE-I scheme and let. From Table 3, we can see that the numerical accuracy in spatial direction is second-order for Implicit and ASE-I scheme, therefore the experimental results are basically consistent with the theoretical analysis.

At last we select and as the spatial grid number and the temporal grid number. In terms of computation time in Table 4, the computational efficiency (CPU time) of the ASE-I (ASI-E) scheme has big advantage compared with implicit scheme. With the increase of the grid number, the computation times of the implicit scheme rapidly grow up, and which the ASE-I (ASI-E) scheme’s has a lower growth rate by comparison. The computation time of the ASE-I (ASI-E) scheme can save nearly 75% compared with the implicit scheme and the Sp (Speedup) is approximately 4.38 and 4.16 respectively. Comprehensively considering the computing efficiency and the computing accuracy, the ASE-I scheme can be more effective to solve the time fractional sub-diffusion equation. When the long time course is calculated, the parallel computing advantages of ASE-I scheme will be more obvious.

Table 2. Numerical errors and convergence of ASE-I scheme in temporal direction (h = 0.02).

Table 3. Numerical errors and convergence of ASE-I scheme in spatial direction ().

Table 4. Comparison of the three difference schemes’ CPU time.

5. Conclusion

For the time fractional sub-diffusion equation, this paper constructs the ASE-I and ASI-E difference schemes with unconditional stability and convergence. Numerical experiments verify the theoretical analyses and show that the proposed scheme is of excellent computational accuracy and obvious parallel properties. The ASE-I (ASI-E) scheme given by this paper can be extended to solve other fractional diffusion models and the parallel computing advantages of the ASE-I (ASI-E) scheme will be more obvious for the long time course and the high dimensional fractional diffusion equation. But the application of the ASE-I (ASI-E) scheme in multidimensional fractional differential equations remains to be further studied.

Acknowledgements

The research is supported by the National Natural Science Foundation of China (Grant No. 11371135) and the Fundamental Research Funds of the Central Universities (Grant No. 2018MS168).

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Chen, W., Sun, H.G. and Li, X.C. (2010) Fractional Derivative Modeling of Mechanics and Engineering Problems. Science Press, Beijing. (In Chinese)
[2]	Uchaikin, V.V. (2013) Fractional Derivatives for Physicist and Engineers, Volume Ⅱ: Applications. Springer, New York.
[3]	Kai, D. (2010) The Analysis of Fractional Differential Equations. Journal of Mathematical Analysis and Applications, 265, 229-248.
[4]	Fulger, D.., Scalas, E. and Germano, G. (2008) Efficient Monte Carlo Simulation of Uncoupled Continuous-Time Random Walks and Approximate Solution of the Fractional Diffusion Equation. Physical Review E, 77, Article ID: 021122.
[5]	Guo, B.L., Pu, X.K. and Huang, F.H. (2015) Fractional Partial Differential Equations and their Numerical Solutions. Science Press, Beijing.
[6]	Sun, Z.Z. and Gao, G.H. (2015) Finite Difference Method for Fractional Differential Equations. Science Press, Beijing. (In Chinese)
[7]	Chen, W. and Sun, H.G. (2017) Fractional Differential Equations and Statistical Models for Anomalous Diffusion. Science Press, Beijing. (In Chinese)
[8]	Sweilam, N.H., Moharram, H., Moniem, N.K.A. and Ahmed, S. (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
[9]	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
[10]	Liu, Q.X., Gu, Y.T., Liu, F.W. and Zhuang, P.H. (2011) An Implicit RBF Meshless Approach for Time Fractional Diffusion Equations. Computational Mechanics, 48, 1-12. https://doi.org/10.1007/s00466-011-0573-x
[11]	Liu, F.W. and Zhuang, P.H. (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
[12]	Yuste, S.B. (2004) 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
[13]	Zhang, P. and Pu, H. (2017) The Error Analysis of Crank-Nicolson-Type Difference Scheme for Fractional Sub-Diffusion Equation with Spatially Variable Coefficient. Boundary Value Problems, 15, 1-19.
[14]	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
[15]	Chi, X.B., Wang, Y.W. and Wang, Y. (2015) Parallel Computing and Implementation Technology. Science Press, Beijing. (In Chinese)
[16]	Bjorstad, P. and Luskin, M. (2015) Parallel Solution of Partial Differential Equations. Science Press, Beijing.
[17]	Zhang, B., Yuan, G., Liu, X. and Chen, J. (1994) Parallel Finite Difference Methods for Partial Differential Equations. Science Press, Beijing. (In Chinese)
[18]	Yuan, G. and Hang, X. (2010) Conservative Parallel Schemes for Diffusion Equations. Chinese Journal of Computational Physics, 27, 475-491.
[19]	Liu, W. (2012) Actual Combat Matlab Parallel Programming. Beihang University Press, Beijing. (In Chinese)
[20]	Kai, D. (2011) An Efficient Parallel Algorithm for the Numerical Solution of Fractional Differential Equations. Fractional Calculus and Applied Analysis, 14, 475-490. https://doi.org/10.2478/s13540-011-0029-1
[21]	Gong, C., Bao, W. and Tang, G. (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
[22]	Zhang, Y., Sun, Z. and Wu, H. (2011) Error Estimates of Crank-Nicolson-Type Difference Schemes for the Sub-Diffusion Equation. Siam Journal on Numerical Analysis, 49, 2302-2322. https://doi.org/10.1137/100812707