# An Efficient Parallel Approximate Algorithm for Solving Time Fractional Reaction-Diffusion Equations

Authors: Xiaozhong Yang; Lifei Wu
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

In this paper, we construct pure alternative segment explicit-implicit (PASE-I) and implicit-explicit (PASI-E) difference algorithms for time fractional reaction-diffusion equations (FRDEs). They are a kind of difference schemes with intrinsic parallelism and based on classical explicit scheme and classical implicit scheme combined with alternating segment technology. The existence and uniqueness analysis of solutions of the parallel difference schemes are given. Both the theoretical proof and the numerical experiment show that PASE-I and PASI-E schemes are unconditionally stable and convergent with second-order spatial accuracy and2−α order time accuracy. Compared with implicit scheme and E-I (I-E) scheme, the computational efficiency of PASE-I and PASI-E schemes is greatly improved. PASE-I and PASI-E schemes have obvious parallel computing properties, which shows that the difference schemes with intrinsic parallelism in this paper are feasible to solve the time FRDEs.


---

## Body

## 1. Introduction

Fractional differential equations can be used to describe some physical phenomena more accurately than the classical integer order differential equations. The time fractional reaction-diffusion equations (FRDEs) play an important role in dynamical systems of physics [1, 2], bioinformatics [3, 4], image processing [5], and other research areas [6, 7]. With the further application of FRDEs, the numerical solutions of FRDEs have become an urgent research work. The analytical solutions of FRDEs are difficult to give explicitly, even the analytical solution of linear fractional differential equation contains special functions, such as Mittag-Leffler function. The series corresponding to these functions converge very slowly; therefore, it is quite difficult to calculate these special functions in practical application. It makes the development of efficient numerical algorithms for FRDEs particularly important [8–10].In recent years, many scholars have studied the numerical algorithm of the fractional diffusion model [11–13]. The finite difference method is still dominant among the existing algorithms now. Tadjeran et al. [14] examined a second-order accurate numerical method in time and in space to solve a class of initial-boundary value fractional diffusion equations with variable coefficients. Gao and Sun [15] derived a compact finite difference scheme for the subdiffusion equation, which is fourth-order accuracy compact approximate for the second-order space derivative. For time fractional subdiffusion equation with Dirichelt boundary value conditions, Luo et al. [16] established a novel collocation method via taking quadratic spline polynomialsas basic functions with Oτ4+h4 order accuracy, when the solution υx,t has four-order continual derivative with respects to x and t. Moreover, the accuracy of the proposed method does not depend on the order α. Luo et al. [17] proposed a novel numerical approximate method for the Caputo fractional derivative by employing the piecewise linear and quadratic Lagrange interpolation functions, which convergence order is Oτ3−α with α∈0,1. Yaseen et al. [18] proposed a finite difference scheme with a cubic trigonometric B-spline function for time fractional diffusion-wave equation with reaction term. The scheme is unconditional stabile and convergent. For a class of one-dimensional and two-dimensional time FRDE with variable coefficients and time drift term, Zhao and Gu [19] presented an implicit finite difference scheme based on the L2−1σ formula. The unconditional stability and convergence of this scheme are proved rigorously by the discrete energy method, and the optimal convergence order in the L2−norm is Oτ2+h2. For semilinear parabolic singularly perturbed systems of reaction-diffusion type, Clavero and Jorge [20] constructed a numerical method, which combined the central finite difference schemes to discretize in space and a linearized fractional implicit Euler method together with a splitting by components technique to integrate in time. The computational cost of the algorithm is considerably less than that of classical schemes. For three kinds of typical Caputo-type PDEs (including Caputo-type reaction-diffusion equation, reaction-diffusion-wave equation, and cable equation), Li and Wang [21, 22] proposed the local discontinuous Galerkin finite element methods. They studied the existence, uniqueness, and regularity of the solutions to these equations.However, the existing implicit serial difference schemes are of relatively high computational complexity and relatively low computation efficiency. The reason is the fractional differential operators have nonlocal property [4, 11]. So, its numerical methods have full coefficient matrices, which require the storage of ON2 and computational cost of ON3. Even with high-performance computers, it is difficult to simulate long time or large computational domain. How to construct an effectively numerical algorithm with high precision and low computational complexity becomes an urgent problem.In order to overcome the difficulty of numerical algorithms of fractional order differential equations which required large amount of computational cost and storage, there are two main types of solutions in the existing literature. One is the fast algorithm which needs to combine the characteristics of the algorithm for reducing computational complexity of the algorithms. The another is the parallel algorithm which makes full use of multicore and cluster devices. Parallel algorithms have been successfully applied in integer order differential equations [23, 24].The research achievements of fast algorithms of fractional differential equations are as follows. Wang et al. [25] proposed a fast algorithm for solving the spatial fractional diffusion equation based on the special structure of difference scheme. Lu et al. [26] proposed a fast numerical difference algorithm for time fractional subdiffusion equation. For time space fractional convection-diffusion equations, Gu et al. [27] presented a new unconditionally stable implicit difference method, which converges with the second-order accuracy in both time and space variables. Then, they solved the Toeplitz-like systems of linear equations by Krylov subspace solvers with suitable circulant preconditioners. Each time level of these methods reduces the memory requirement of the proposed implicit difference scheme from ON2 to ON and the computational complexity from ON3 to ONlogN in each iterative step. Jiang et al. [28] have developed a fast algorithm for the evaluation of the Caputo fractional derivative 0<α<1. The algorithm relies on an efficient sum-of-exponentials approximation for the convolution kernel t−1−K with the absolute error over the interval △t,T. The algorithm has nearly optimal complexity in both CPU time and storage. Zhao et al. [29] discussed a block lower triangular Toeplitz system arising from the time space fractional diffusion equation and applied the preconditioned biconjugate gradient stabilized method and the flexible general minimal residual method to obtain the solutions of the linear system. They develpoed a block bidiagonal Toeplitz preconditioner for the block lower triangular Toeplitz system, whose storage is of ON with N being the spatial grid number.With the rapid development of multicore and cluster technology, parallel algorithms have become one of the mainstream technologies to improve computing efficiency. For integer order parabolic equation, Evans and Abdullab [30] put forward the idea of group explicit (GE), and designed the alternating group explicit (AGE) scheme, which not only guarantees the stability of numerical calculation, but also has good parallel property. Inspired by the AGE method, Zhang et al. [31] proposed the idea of constructing alternating segment implicit scheme by using Saul’yev asymmetric scheme, then established a variety of alternating segment explicit-implicit (ASE-I) and alternating segment pure implicit parallel difference formulas by appropriate use of alternating technique. The methods have been well applied to numerical solving integer partial differential equations (PDEs) [32, 33]. The existing numerical methods cannot be applied to numerically solve fractional order PDEs directly and even produce completely different numerical analysis. How to extend the existing parallel difference method of integer order PDEs to fractional order PDEs is a great challenge to the computational mathematics (physics) fields.The research achievements of parallel algorithms for fractional differential equations are as follows. Diethelm [34] implements parallel computation for Adams-Bashforth-Moulton method of the second order of the fractional derivative and discusses the accuracy of the parallel algorithm. Gong et al. [35, 36] parallelizes the explicit difference scheme of the time fractional diffusion equation and implicit scheme of the Riesz type space FRDE. The core of parallelization is the product of matrix and vector, the addition vector and vector. Sweilam et al. [37] constructed a class of parallel Crank-Nicolson difference schemes for time fractional parabolic equations. The preconditioned conjugate gradient method was applied to solve discrete algebraic equations in parallel. Wang et al. [38] studies the parallel algorithm of the implicit difference scheme for Caputo-type FRDE. Compared with the classical schemes, the parallel algorithm has some improvements in computational efficiency. Under the principle of minimizing communication, reasonably allocating computing tasks and minimizing the change of the original serial difference format, the serial algorithm is parallelized. Biala and Khaliq [39] developed a time stepping scheme for solving nonlinear time space fractional PDEs. The schemes are implemented in parallel using the distributed (MPI), shared memory systems (OpenMP), and a combination of both (Hybrid). They discussed the advantages of the parallel algorithms over the sequential version. Fu and Wang [40] developed a parallel finite difference scheme for space-time fractional PDEs. By exploring the structure of the stiffness matrices, they developed a matrix-free preconditioned fast Krylov subspace iterative method with significantly reduced computational work and memory requirement. Most of the parallel algorithms of fractional PDEs are from the perspective of numerical algebra.For a long time, most of the existing parallel schemes are conditionally stable or unconditionally stable but with only first order accuracy in space. In order to obtain more accurate and unconditionally stable parallel difference schemes, we do not study parallel algorithms from the point of numerical algebra, but based on the parallelization of traditional difference schemes. When faced with long time history problem, large computational domain problem or high-dimensional problem, the amount of computation and storage of serial algorithms will increase exponentially [31]. At this point, considering the calculation accuracy and efficiency, the parallel algorithm with second-order spatial accuracy will be the good choice. In this paper, for timeFRDEs, the pure alternative segment explicit-implicit (PASE-I) and implicit-explicit (PASI-E) parallel difference schemes are obtained by processing the explicit and implicit difference schemes with alternating segment technique. Comparing with the existing difference schemes, this kind of schemes can greatly improve the efficiency of numerical simulation of fractional order models.The rest of the paper is organized as follows. In Section2, we structure the PASE-I parallel difference scheme of time FRDE. In Section 3, we prove the uniqueness, stability, and convergence of PASE-I scheme and PASI-E difference scheme. In Section 4, we present some numerical experiments to demonstrate the accuracy and efficiency of the PASE-I and PASI-E methods. Finally, some concluding remarks are presented in Section 5.


## 2. Parallel Intrinsic Difference Schemes for Time Fractional Reaction-Diffusion Equation

### 2.1. Time Fractional Reaction-Diffusion Equation

In this section, we consider the time FRDE [1, 8]:(1)∂αux,t∂tα=∂2ux,t∂x2−pux,t+fx,t,0<x<L,0<t≤T,ux,0=u0x,0≤x≤L,u0,t=uL,t=0,0≤t≤T,where u0x, fx,t are given functions, p is a nonnegative constant, 0<α<1, ∂αux,t/∂tα is a fractional derivative of Captuo type:(2)∂αux,t∂tα=1Γ1−α∫0t∂ux,τ∂τdτt−τα.Whenfx,t=0, by using the finite sinusoidal transformation and Laplace transformation, the analytic solution of equation (1) can be obtained [1, 8]:(3)ux,t=2L∑n=1∞Eα−p+a2n2tαsinanx∫0Lu0rsinanxdr,where Eαz=∑k=0∞zα/Γαk+1, a=π/L is a special function of type Mittag-Leffler, and(4)E1−z=exp−z,E1/2z=∑k=0∞zkΓk/2+1=∑k=0∞zkΓk!=expz2erfc−z,where erfcz is error function, erfcz=2/π∫z∞exp−t2dt.


### 2.2. Construction of PASE-I Parallel Difference Scheme

We will divide the computational domain into mesh by taking the space steph and the time step k. The mesh points are xi,tn, here h=L/M, xi=ih,i=0,1,2,…,M; k=T/N, tn=nk,n=0,1,2,…,N; M,N are positive integers. Let the value of Vik be approximated to the value of uxi,tk.The time fractional derivative term can be approximated by the following format:(5)∂αuxi,tn+1∂tα=1Γ1−α∑j=0n∫jkj+1k∂Vx,ξ∂ξdξtn+1−ξα≈1Γ1−α∑j=0nVxi,tj+1−Vxi,tjk∫jkj+1kdξtn+1−ξα=k−αΓ2−α∑j=0nljVxi,tn+1−j−Vxi,tn−j.LetD0Ctαuxi,tn+1=k−α/Γ2−α∑j=0nljVxi,tn+1−j−Vxi,tn−j, lj=j+11−α−j1−α.In order to construct PASE-I scheme of equation (1), the classical explicit and implicit schemes of equation (1) are given. Let c=kαΓ2−α,b=pc, a=c/h2.(1)
The classical explicit scheme(6)D0Ctαuxi,tn+1=1h2Vi−1n−2Vin+Vi+1n−pVin+fin.Simplify,(7)∑j=0nljVin+1−j−Vin−j=aVi−1n−2Vin+Vi+1n−bVin+cfin.(2)
The classical implicit scheme(8)D0Ctαuxi,tn+1=1h2Vi−1n+1−2Vin+1+Vi+1n+1−pVin+1+fin+1.Simplify,(9)∑j=0nljVin+1−j−Vin−j=aVi−1n+1−2Vin+1+Vi+1n+1−bVin+1+cfin+1.Whenn=0,(10)Vi1=aVi−10+1−2a−bVi0+aVi+10+cfi0,−aVi−11+1+2a+bVi1−aVi+11=Vi0−cfi1.Whenn>0,(11)Vin+1=aVi−1n+1−2a−l1−bVin+aVi+1n+∑j=1n−1lj−lj+1Vin−j+lnVi0+cfin,−aVi−1n+1+1+2a+bVin+1−aVi+1n+1=1−l1Vin+∑j=1n−1lj−lj+1Vin−j+lnVi0+cfin+1.The classical explicit scheme (7) has ideal parallelism and is very suitable for parallel computing. However, it is conditionally stable, especially in multidimensional problems its time step is restricted severely. The classical implicit scheme (9) is unconditionally stable, and usually discrete into tridiagonal systems. The tridiagonal systems are solved inON complexity using the Thomas algorithm. It is inconvenient to obtain the numerical results directly and quickly. We will apply the classical explicit scheme (7) and implicit scheme (9) combining with alternating segmentation technology to construct PASE-I scheme.We construct PASE-I scheme as follows. LetM−1=Bl (B≥3,l≥3 are integer numbers and B is an odd number). On the same time level, the space grid points are divided into B segments, which are recorded in the order of S1,S2,…,SB. On the odd time level, each segment is calculated according to the rule of “classical explicit-classical implicit-classical explicit” from left to right; on the even time level, each segment is calculated according to the order of “classical implicit-classical explicit-classical implicit” from left to right. For each segment, the explicit and implicit schemes are applied alternately in time. A piecewise sketch of grid points is shown in Figure 1. In Figure 1, ◯ denotes the classical explicit format at the grid point, and □ denotes classical implicit schemes, respectively. The internal boundary of each implicit segment depended on the value of the first or last point of the adjacent explicit segment.Figure 1
Point diagram of PASE-I scheme (B = 5).Therefore, the PASE-I scheme of equation (1) is constructed as follows:(12)I+aG1¯V1=I−aG2¯V0+g0+l0V0+cf0,I+aG2¯Vn+1=c1I−aG1¯Vn+gn+∑j=1n−1cj+1Vn−j+lnV0+cfn,n=1,3,…,I+aG1¯Vn+2=c1I−aG2¯Vn+1+gn+1+∑j=1ncj+1Vn+1−j+ln+1V0+cfn+1,where cj=lj−1−ljj=1,2,…,N,Vk=V1k,V2k,…,VM−1kTk=0,1,…,N,fn=f1n,f2n,…,fM−1nT, gn=aV0n,0,…,aVMnTn=0,1,…,N,Ql is l×l order zero matrix.(13)G1¯=QlWlQlWl⋱QlWlQlM−1×M−1,G2¯=WlQlWlQl⋱WlQlWlM−1×M−1,Wl=2+ba−1−12+ba−1⋱⋱⋱−12+ba−1−12+bal×l.


## 2.1. Time Fractional Reaction-Diffusion Equation

In this section, we consider the time FRDE [1, 8]:(1)∂αux,t∂tα=∂2ux,t∂x2−pux,t+fx,t,0<x<L,0<t≤T,ux,0=u0x,0≤x≤L,u0,t=uL,t=0,0≤t≤T,where u0x, fx,t are given functions, p is a nonnegative constant, 0<α<1, ∂αux,t/∂tα is a fractional derivative of Captuo type:(2)∂αux,t∂tα=1Γ1−α∫0t∂ux,τ∂τdτt−τα.Whenfx,t=0, by using the finite sinusoidal transformation and Laplace transformation, the analytic solution of equation (1) can be obtained [1, 8]:(3)ux,t=2L∑n=1∞Eα−p+a2n2tαsinanx∫0Lu0rsinanxdr,where Eαz=∑k=0∞zα/Γαk+1, a=π/L is a special function of type Mittag-Leffler, and(4)E1−z=exp−z,E1/2z=∑k=0∞zkΓk/2+1=∑k=0∞zkΓk!=expz2erfc−z,where erfcz is error function, erfcz=2/π∫z∞exp−t2dt.


## 2.2. Construction of PASE-I Parallel Difference Scheme

We will divide the computational domain into mesh by taking the space steph and the time step k. The mesh points are xi,tn, here h=L/M, xi=ih,i=0,1,2,…,M; k=T/N, tn=nk,n=0,1,2,…,N; M,N are positive integers. Let the value of Vik be approximated to the value of uxi,tk.The time fractional derivative term can be approximated by the following format:(5)∂αuxi,tn+1∂tα=1Γ1−α∑j=0n∫jkj+1k∂Vx,ξ∂ξdξtn+1−ξα≈1Γ1−α∑j=0nVxi,tj+1−Vxi,tjk∫jkj+1kdξtn+1−ξα=k−αΓ2−α∑j=0nljVxi,tn+1−j−Vxi,tn−j.LetD0Ctαuxi,tn+1=k−α/Γ2−α∑j=0nljVxi,tn+1−j−Vxi,tn−j, lj=j+11−α−j1−α.In order to construct PASE-I scheme of equation (1), the classical explicit and implicit schemes of equation (1) are given. Let c=kαΓ2−α,b=pc, a=c/h2.(1)
The classical explicit scheme(6)D0Ctαuxi,tn+1=1h2Vi−1n−2Vin+Vi+1n−pVin+fin.Simplify,(7)∑j=0nljVin+1−j−Vin−j=aVi−1n−2Vin+Vi+1n−bVin+cfin.(2)
The classical implicit scheme(8)D0Ctαuxi,tn+1=1h2Vi−1n+1−2Vin+1+Vi+1n+1−pVin+1+fin+1.Simplify,(9)∑j=0nljVin+1−j−Vin−j=aVi−1n+1−2Vin+1+Vi+1n+1−bVin+1+cfin+1.Whenn=0,(10)Vi1=aVi−10+1−2a−bVi0+aVi+10+cfi0,−aVi−11+1+2a+bVi1−aVi+11=Vi0−cfi1.Whenn>0,(11)Vin+1=aVi−1n+1−2a−l1−bVin+aVi+1n+∑j=1n−1lj−lj+1Vin−j+lnVi0+cfin,−aVi−1n+1+1+2a+bVin+1−aVi+1n+1=1−l1Vin+∑j=1n−1lj−lj+1Vin−j+lnVi0+cfin+1.The classical explicit scheme (7) has ideal parallelism and is very suitable for parallel computing. However, it is conditionally stable, especially in multidimensional problems its time step is restricted severely. The classical implicit scheme (9) is unconditionally stable, and usually discrete into tridiagonal systems. The tridiagonal systems are solved inON complexity using the Thomas algorithm. It is inconvenient to obtain the numerical results directly and quickly. We will apply the classical explicit scheme (7) and implicit scheme (9) combining with alternating segmentation technology to construct PASE-I scheme.We construct PASE-I scheme as follows. LetM−1=Bl (B≥3,l≥3 are integer numbers and B is an odd number). On the same time level, the space grid points are divided into B segments, which are recorded in the order of S1,S2,…,SB. On the odd time level, each segment is calculated according to the rule of “classical explicit-classical implicit-classical explicit” from left to right; on the even time level, each segment is calculated according to the order of “classical implicit-classical explicit-classical implicit” from left to right. For each segment, the explicit and implicit schemes are applied alternately in time. A piecewise sketch of grid points is shown in Figure 1. In Figure 1, ◯ denotes the classical explicit format at the grid point, and □ denotes classical implicit schemes, respectively. The internal boundary of each implicit segment depended on the value of the first or last point of the adjacent explicit segment.Figure 1
Point diagram of PASE-I scheme (B = 5).Therefore, the PASE-I scheme of equation (1) is constructed as follows:(12)I+aG1¯V1=I−aG2¯V0+g0+l0V0+cf0,I+aG2¯Vn+1=c1I−aG1¯Vn+gn+∑j=1n−1cj+1Vn−j+lnV0+cfn,n=1,3,…,I+aG1¯Vn+2=c1I−aG2¯Vn+1+gn+1+∑j=1ncj+1Vn+1−j+ln+1V0+cfn+1,where cj=lj−1−ljj=1,2,…,N,Vk=V1k,V2k,…,VM−1kTk=0,1,…,N,fn=f1n,f2n,…,fM−1nT, gn=aV0n,0,…,aVMnTn=0,1,…,N,Ql is l×l order zero matrix.(13)G1¯=QlWlQlWl⋱QlWlQlM−1×M−1,G2¯=WlQlWlQl⋱WlQlWlM−1×M−1,Wl=2+ba−1−12+ba−1⋱⋱⋱−12+ba−1−12+bal×l.


## 3. Numerical Analysis of PASE-I Difference Scheme

### 3.1. Uniqueness of PASE-I Difference Scheme’s Solution

Theorem 1.
The pure alternative segment explicit-implicit (PASE-I) difference scheme (12) for the time fractional reaction-diffusion equation is uniquely solvable.Proof.
From the definition ofG1¯ and G2¯, we can know that G1¯ and G2¯ are strictly diagonally dominant matrices, and the main diagonal elements of the two matrices are positive real numbers. I+aG1¯ and I+aG2¯ are nonsingular matrices. From that I+aG2¯ is a nonsingular matrix, the second equation of PASE-I scheme (12) is uniquely solvable. From that I+aG1¯ is a nonsingular matrix, the first and third equations of PASE-I scheme (12) are uniquely solvable. Therefore, Theorem 1 can be obtained.


### 3.2. Stability of PASE-I Difference Scheme

Using the property of the functiongx=x1−αx≥1, the following conclusions are obtained [4, 9]:(14)1=l0>l1>l2>⋯>ln>0,ln⟶0n⟶∞,∑j=1ncj=1−ln,∑j=1∞cj=1.Lemma 1.
If the matrixC is a nonnegative real matrix, then there is an estimation formula σI−ρCI+ρC−12≤1 for any parameter ρ≥0,0≤σ≤1.Proof.
(15)σI−ρCI+ρC−122=maxφ∈Rn,φ≠0σI−ρCI+ρC−1φ,σI−ρCI+ρC−1φφ,φ.
Make a transformationψ=I+ρC−1φ, then(16)σI−ρCI+ρC−122=maxψ∈Rn,ψ≠0σI−ρCψ,σI−ρCψI+ρCψ,I+ρCψ=maxψ∈Rn,ψ≠0σ2ψ,ψ−2σρCψ,ψ+ρ2Cψ,Cψψ,ψ+2ρCψ,ψ+ρ2Cψ,Cψ≤1.
Therefore,σI−ρCI+ρC−12≤1. That is the proof of Lemma 1.
For stability, assume thatVin is the numerical solution of the PASE-I difference scheme and V¯in is the numerical solution of the PASE-I difference scheme from another initial value. Error εin=Vin−V¯in,ε0n=εMn=0,εk=ε1n,ε2n,…,εM−1n, The error εin is substituted into scheme (12) to obtain the following formula:(17)I+aG1¯ε1=I−aG2¯ε0+l0ε0,I+aG2¯εn+1=c1I−aG1¯εn+∑j=1n−1cj+1εn−j+lnε0,n=1,3,…,I+aG1¯εn+2=c1I−aG2¯εn+1+∑j=1ncj+1εn+1−j+ln+1ε0.
Letλ1,λ2λ1,λ2≥0 is any eigenvalue of the matrix aG1¯, and meet c1−λ1 takes the maximum value and 1/1+λ2 takes the minimum value. λ3,λ4λ3,λ4≥0 is any eigenvalue of the matrix aG2¯, and meet c1−λ3 takes the maximum value and 1/1+λ4 takes the minimum value. Because the minimum eigenvalues of the matrices aG1¯ and aG2¯ are both 0, which means that λ2=λ4=0. Therefore, 1/1+λ2=1/1+λ4=1. From Lemma 1, let Tj=c1I−aGj¯I+aGj¯−1,j=1,2, we can know(18)T12=c1I−aG1¯I+aG1¯−12≤c1I−aG1¯2I+aG1¯−12=c1−λ11+λ2≤1.
Whenλ1<c1, we have c1−λ1+l1/1+λ2=c1−λ1+l1/1+λ2=1−λ1/1+λ2<1. When λ1≥c1, let λ1=r1c1, here r1≥1. From c1−λ1/1+λ2≤1, we can get r1≤2 and c1≤λ1≤2c1. c1−λ1+l1/1+λ2=λ1−c1+l1/1+λ2≤c1+l1/1+λ2≤1.
Whenm=1,(19)ε1=I+aG1¯−1I−aG2¯ε0,ε12=I+aG1¯−1I−aG2¯ε02≤1−λ31+λ2ε02≤ε02.
Whenm=2,(20)ε2=I+aG2¯−1c1I−aG1¯ε1+l1ε0,ε22=I+aG2¯−1c1I−aG1¯ε1+l1ε02≤I+aG2¯−1c1I−aG1¯ε12+I+aG2¯−1l1ε02≤c1−λ1+l11+λ4ε02=c1−λ1+l11+λ2ε02≤ε02.
Suppose whenm≤2n, the equation εm2≤ε02 is established. When m=2n+1, we have(21)ε2n+1=I+aG1¯−1c1I−aG2¯ε2n+∑j=12n−1cj+1ε2n+1−j+l2nε0,ε2n+12=I+aG1¯−1c1I−aG2¯ε2n+∑j=12n−1cj+1ε2n+1−j+l2nε02≤I+aG1¯−12c1I−aG2¯ε2n+∑j=12n−1cj+1ε2n+1−j+l2nε02≤I+aG1¯−12c1I−aG2¯2ε02+∑j=12n−1cj+1ε02+l2nε02=c1−λ3+∑j=12n−1cj+1+l2n1+λ2ε02=c1−λ3+l11+λ4ε02.
Whenλ3<c1, we have c1−λ3+l1/1+λ4=c1−λ3+l1/1+λ4=1−λ3/1+λ4<1. When λ3≥c1, let λ3=r1c1, here r1≥1; from c1−λ3/1+λ4≤1, we can get r1≤2, c1≤λ3≤2c1; and c1−λ3+l1/1+λ4=λ3−c1+l1/1+λ4≤c1+l1/1+λ4≤1. Therefore, ε2n+12≤ε02.
Whenm=2n+2,(22)ε2n+2=I+aG2¯−1c1I−aG1¯ε2n+2+∑j=12ncj+1ε2n+1−j+l2n+1ε0,ε2n+22=I+aG2¯−1c1I−aG1¯ε2n+2+∑j=12ncj+1ε2n+1−j+l2n+1ε02≤I+aG2¯−12c1I−aG1¯ε2n+2+∑j=12ncj+1ε2n+1−j+l2n+1ε02≤I+aG2¯−12c1I−aG1¯2ε02+∑j=12ncj+1ε02+l2n+1ε02=c1−λ1+∑j=12ncj+1+l2n+11+λ4ε02=c1−λ1+l11+λ2ε02≤ε02.
In the same way, we can getε2n+22≤ε02.
To sum up,εn2≤ε02,n=1,2,3,…, and the following theorem can be obtained.Theorem 2.
The pure alternative segment explicit-implicit (PASE-I) difference scheme (12) for time fractional reaction-diffusion equation is unconditionally stable.


### 3.3. Convergence of PASE-I Difference Scheme

First, the accuracy of PASE-I scheme is analyzed. The PASE-I scheme is alternatively applied explicit and implicit segments, not only in the direction of space, but also in the direction of time. Examination of grid pointsxi,t2n+1 and xi,t2n+2. The classical explicit scheme is(23)D0Ctαuxi,t2n+1=1h2Vi−12n−2Vi2n+Vi+12n−pVi2n+fi2n,(24)Vi2n+1=aVi−12n+1−2a−l1−bVi2n+aVi+12n+∑j=12n−1lj−lj+1Vi2n−j+l2nVi0+cfi2n.The classical implicit scheme is(25)D0Ctαuxi,t2n+2=1h2Vi−12n+2−2Vi2n+2+Vi+12n+2−pVi2n+2+fi2n+2,(26)−aVi−12n+2+1+2a+bVi2n+2−aVi+12n+2=1−l1Vi2n+1+∑j=12nlj−lj+1Vi2n+1−j+l2n+1Vi0+cfi+12n+2.We expand (24) and (26) at two points xi,t2n+1 and xi,t2n+2 for Taylor, respectively. The truncation error is recorded as T1k,h,≤2n+11−αCvk1−α/Γ2−αT2k,h. It is known that the precision of D0Ctαuxi,t2n+1 is 2−α order [4, 9]. Let CV=max0≤t≤t2n+1∂2Vx,t/∂t2, we have(27)∂α+1Vxi,t2n+1∂tα+1=1Γ1−α∫0t2n+1∂2Vxi,τ∂τ2dτt2n+1−τα=1Γ1−α∑j=02n∫jkj+1kdτt2n+1−τα≤1Γ1−αmax0≤t≤t2n+1∂2Vx,t∂t2∑j=02n∫jkj+1kdτt2n+1−τα=CVΓ1−α∑j=02n∫jkj+1kdτt2n+1−τα=CVk1−α1−αΓ1−α∑j=02n2n+1−j1−α−2n−j1−α≤2n+11−αCVk1−αΓ2−α,T1k,h=∂αVxi,t2n+1∂tα−1h2Vi−12n−2Vi2n+Vi+12n+pVi2n=∂αVx,t∂tα−Vxx+pV+kVxxt−k22Vxxtt−h212Vxxxx−pkVt+pk22Vtt+Ok2−α+h2,T2k,h=∂αVxi,t2n+2∂tα−1h2Vi−12n+2−2Vi2n+2+Vi+12n+2+pVi2n+2=∂αVxi,t2n+1∂tα+k∂α+1Vxi,t2n+1∂tα+1−Vxx+pV−kVxxt−k22Vxxtt−h212Vxxxx+pkVt+pk22Vtt+Ok2−α+h2.The two termskVxxt,pkVt of T1k,h have the same form as two terms kVxxt,pkVt of T2k,h, but the opposite of the symbol. When the explicit format and implicit format are applied alternately, the errors of the two terms will be offset. And k∂α+1Vxi,t2n+1/∂tα+1≤Ck2−α,C>0. When applying explicit and implicit formats alternately, the accuracy is 2 order in space direction and 2−α in time direction. Similarly, the same results can be obtained when implicit and explicit formats are applied alternately. Therefore, the truncation of PASE-I scheme is Rin=Ok2−α+h2.Whenm=1, the first layer is calculated applying the explicit format.(28)Vi1=aVi−10+1−2a−bVi0+aVi+10+cfi0.Whenm=2, the second layer is calculated applying the implicit format.(29)Vi1=aVi0+1−2a−bVi0+aVi+10+cfi0,(30)−aVi−12+1+2a+bVi2−aVi+12=c1Vi1+l1Vi0+cfi2.Substitute (29) into (30) to eliminate Vi1:(31)−aVi−12+1+2a+bVi2−aVi+12=c1aVi−10+c11−2a−b+l1Vi0+c1aVi+10+cc1fi0+cfi2.Whenm>2, we substitute the explicit format (24) into the implicit format (26) to eliminate Vi2n+1:(32)−aVi−12n+2+1+2a+bVi2n+2−aVi+12n+2=c1aVi−12n+c1c1−2a−bVi2n+c1aVi+12n+c2Vi2n+∑j=12n−1c1cj+1+cj+2Vi2n−j+c1l2n+l2n+1Vi0+cc1fi2n+cfi2n+2.For convenience, assume thatV˜in is the exact solution of time FRDE, and Vin is the numerical solution of the difference scheme, Ein=Vin−V˜in. For n=0,1,2,…,N, (28), (31) and (32) we have(33)Ei1=aEi−10+1−2a−bEi0+aEi+10+kαRi1,−aEi−12+1+2a+bEi2−aEi+12=c1aEi−10+c11−2a−b+l1Ei0+c1aEi+10+kαRi2,−aEi+12n+2+1+2a+bEi2n+2−aEi−12n+2=c1aEi−12n+c1c1−2a−bEi2n+c1aEi+12n+c2Ei2n+∑j=12n−1c1cj+1+cj+2Ei2n−j+c1l2n+l2n+1Ei0+Ri2n+2,where Rin=Ok2−α+h2,E0n=EMn=Ei0=0,n=0,1,…,N,i=0,1,…,M.Define two grid functions as follows:(34)Enx=Eln,x∈xl−1/2,xl+1/2,l=1,2,…,M−1,0,x∈0,h2∪L−h2,L,Rnx=Rln,x∈xl−1/2,xl+1/2,l=1,2,…,M−1,0,x∈0,h2∪L−h2,L,where Rn=R1n,R2n,…,RM−1nT,En=E1n,E2n,…,EM−1nT,n=0,1,2,…,N. We do Fourier expansion on Enx and Rnx:(35)Enx=∑j=−∞+∞Ejnei2jπx/L,i=−1,n=0,1,…,N,j=0,1,…,M,Rnx=∑j=−∞+∞Rjnei2jπx/L,i=−1,n=0,1,…,N,j=0,1,…,M,where Ejn=1/L∫0LEnei2jπx/Ldx,Rjn=1/L∫0LRnei2jπx/Ldx,j=±0,±1,…Define normEnx=∑j=1M−1hEjn2, Rnx=∑j=1M−1hRjn2.We have the following formula:(36)Enx2=∑j=1M−1hEjn2=L∑j=−∞+∞Ejn2,n=0,1,…,N,(37)Rnx2=∑j=1M−1hRjn2=L∑j=−∞+∞Rjn2,n=0,1,…,N.SetEln=eneiλlh,Rln=rneiλlhλ=2jπ/L, and substitute (33), we have(38)e1=kαr1,−ae−iλh+1+2a+b−aeiλhe2=kαr2,−ae−iλh+1+2a+b−aeiλhe2n+2=c1ae−iλh+c1c1−2a−b+c1aeiλh+c2e2n+∑j=12n−1c1cj+1+cj+2e2n−j+kαr2n+2.Becausesin2λh/2=−1/4e−iλh−2+eiλh, (38) can be converted to(39)e1=kαr1,1+b+4asin2λh2e2=kαr2,1+b+4asin2λh2e2n+2=c12+c2−c1b−4c1asin2λh2e2n+∑j=12n−1c1cj+1+cj+2e2n−j+kαr2n+2.Lemma 2.
Letek be the solution of (39), ∃C1>0 makes e2n+2≤C1l2n+1−1kαr1,n=0,1,…Proof.
FromRjn=Ok2−α+h2,j=0,1,…,M,n=0,1,…,N, we know ∃C2>0, makes(40)Rjn≤C2k2−α+h2.
According toRnx2=∑j=1M−1hRjn2=L∑j=−∞+∞Rjn2,n=0,1,…,N, we have(41)Rn≤C2Lk2−α+h2.
FromRln=rneiλlhλ=2jπ/L and (37), we can know(42)rn≡Rjn≤C3Rj1≡C3r1,n=2,3,…,N,where C3 is a positive number. From (39), we have(43)e1=kαr1=l0−1kαr1,e2=kαr21+b+4asin2λh/2≤C3r1=C3l1−1r1.
Supposeek≤C1lk−1−1kαr1k≤2n,C1=max1,C3, we have(44)1+b+4asin2λh2e2n+2=c12+c2−c1b−4c1asin2λh2e2n+∑j=12n−1c1cj+1+cj+2e2n−j+kαr2n+2≤c1b+4c1asin2λh2e2n+∑j=12nc1cj+cj+1e2n+1−j+C3kαr1≤c1b+4c1asin2λh2e2n+∑j=12nc1cj+cj+1C1l2n−j−1kαr1+C3kαr1≤c1b+4c1asin2λh2C1l2n−1−1kαr1+∑j=12nc1cj+cj+1C1l2n−j−1kαr1+C3kαr1≤c1b+4c1asin2λh2C1l2n+1−1kαr1+∑j=12nc1cj+cj+1C1l2n+1−1kαr1+C3kαr1≤c1b+4c1asin2λh2+c1+∑j=12ncj+1+l2n+1−c1l2nC1l2n+1−1r1≤1+c1b+4c1asin2λh2C1l2n+1−1kαr1≤1+b+4asin2λh2C1l2n+1−1kαr1.
Therefore, we havee2n+2≤C1l2n+1−1kαr1. Lemma 2 proof is complete.
Becausenαkα≤Tαn=1,2,…,N,limn⟶∞ln−1/nα=limn⟶∞n−α/n+11−α−n1−α=limn⟶∞n−1/1+1/n1−α−1=1/1−α, we have(45)e2n+2≤C1l2n+1−1kαr1≤C1l2n+1−12n+1−α2n+1αkαr1≤C11−α2n+1−αkαr1≤C11−αTαr1.
From (36), (37) and (41) we can get(46)En≤C11−αTαR1≤C11−αTαC2Lk2−α+h2=Ck2−α+h2.
To sum up, the following theorem can be obtained.Theorem 3.
The pure alternative segmented explicit-implicit (PASE-I) scheme (12) of the time fractional reaction-diffusion equation is unconditionally convergent, and satisfying(47)Vin−V˜in≤Ck2−α+h2,i=1,2,…,M,n=1,2,…,N.here C is a positive number.


### 3.4. PASI-E Scheme of Time Fractional Reaction-Diffusion Equation

Similar to the construction of PASE-I scheme, we give the PASI-E scheme of FRDE.When solving the odd time layer, the calculation is performed according to the rules of “classical implicit-classical explicit-classical implicit”. On the even time layer, the rule of “classical explicit-classical implicit-classical explicit” is applied. Then the PASI-E scheme for solving the time FRDE is obtained as follows:(48)I+aG2¯V1=I−aG1¯V0+g1+l0V0+cf1,I+aG1¯Vn+1=c1I−aG2¯Vn+gn+1+∑j=1n−1cj+1Vn−j+lnV0+cfn+1,n=1,3,…,I+aG2¯Vn+2=c1I−aG1¯Vn+1+gn+1+∑j=1ncj+1Vn+1−j+ln+1V0+cfn+1.The implicit scheme of the first time layer is unconditionally stable and convergent. Starting from the second time layer, with reference to the proof analysis of the PASE-I scheme, there are the following theorems.Theorem 4.
The pure alternative segmented implicit-explicit (PASI-E) scheme (48) of the time fractional reaction-diffusion equation is uniquely solvable, unconditionally stable and convergent, satisfying(49)Vin−V˜in≤Ck2−α+h2,i=1,2,…,M,n=1,2,…,N.here C is a positive number.


## 3.1. Uniqueness of PASE-I Difference Scheme’s Solution

Theorem 1.
The pure alternative segment explicit-implicit (PASE-I) difference scheme (12) for the time fractional reaction-diffusion equation is uniquely solvable.Proof.
From the definition ofG1¯ and G2¯, we can know that G1¯ and G2¯ are strictly diagonally dominant matrices, and the main diagonal elements of the two matrices are positive real numbers. I+aG1¯ and I+aG2¯ are nonsingular matrices. From that I+aG2¯ is a nonsingular matrix, the second equation of PASE-I scheme (12) is uniquely solvable. From that I+aG1¯ is a nonsingular matrix, the first and third equations of PASE-I scheme (12) are uniquely solvable. Therefore, Theorem 1 can be obtained.


## 3.2. Stability of PASE-I Difference Scheme

Using the property of the functiongx=x1−αx≥1, the following conclusions are obtained [4, 9]:(14)1=l0>l1>l2>⋯>ln>0,ln⟶0n⟶∞,∑j=1ncj=1−ln,∑j=1∞cj=1.Lemma 1.
If the matrixC is a nonnegative real matrix, then there is an estimation formula σI−ρCI+ρC−12≤1 for any parameter ρ≥0,0≤σ≤1.Proof.
(15)σI−ρCI+ρC−122=maxφ∈Rn,φ≠0σI−ρCI+ρC−1φ,σI−ρCI+ρC−1φφ,φ.
Make a transformationψ=I+ρC−1φ, then(16)σI−ρCI+ρC−122=maxψ∈Rn,ψ≠0σI−ρCψ,σI−ρCψI+ρCψ,I+ρCψ=maxψ∈Rn,ψ≠0σ2ψ,ψ−2σρCψ,ψ+ρ2Cψ,Cψψ,ψ+2ρCψ,ψ+ρ2Cψ,Cψ≤1.
Therefore,σI−ρCI+ρC−12≤1. That is the proof of Lemma 1.
For stability, assume thatVin is the numerical solution of the PASE-I difference scheme and V¯in is the numerical solution of the PASE-I difference scheme from another initial value. Error εin=Vin−V¯in,ε0n=εMn=0,εk=ε1n,ε2n,…,εM−1n, The error εin is substituted into scheme (12) to obtain the following formula:(17)I+aG1¯ε1=I−aG2¯ε0+l0ε0,I+aG2¯εn+1=c1I−aG1¯εn+∑j=1n−1cj+1εn−j+lnε0,n=1,3,…,I+aG1¯εn+2=c1I−aG2¯εn+1+∑j=1ncj+1εn+1−j+ln+1ε0.
Letλ1,λ2λ1,λ2≥0 is any eigenvalue of the matrix aG1¯, and meet c1−λ1 takes the maximum value and 1/1+λ2 takes the minimum value. λ3,λ4λ3,λ4≥0 is any eigenvalue of the matrix aG2¯, and meet c1−λ3 takes the maximum value and 1/1+λ4 takes the minimum value. Because the minimum eigenvalues of the matrices aG1¯ and aG2¯ are both 0, which means that λ2=λ4=0. Therefore, 1/1+λ2=1/1+λ4=1. From Lemma 1, let Tj=c1I−aGj¯I+aGj¯−1,j=1,2, we can know(18)T12=c1I−aG1¯I+aG1¯−12≤c1I−aG1¯2I+aG1¯−12=c1−λ11+λ2≤1.
Whenλ1<c1, we have c1−λ1+l1/1+λ2=c1−λ1+l1/1+λ2=1−λ1/1+λ2<1. When λ1≥c1, let λ1=r1c1, here r1≥1. From c1−λ1/1+λ2≤1, we can get r1≤2 and c1≤λ1≤2c1. c1−λ1+l1/1+λ2=λ1−c1+l1/1+λ2≤c1+l1/1+λ2≤1.
Whenm=1,(19)ε1=I+aG1¯−1I−aG2¯ε0,ε12=I+aG1¯−1I−aG2¯ε02≤1−λ31+λ2ε02≤ε02.
Whenm=2,(20)ε2=I+aG2¯−1c1I−aG1¯ε1+l1ε0,ε22=I+aG2¯−1c1I−aG1¯ε1+l1ε02≤I+aG2¯−1c1I−aG1¯ε12+I+aG2¯−1l1ε02≤c1−λ1+l11+λ4ε02=c1−λ1+l11+λ2ε02≤ε02.
Suppose whenm≤2n, the equation εm2≤ε02 is established. When m=2n+1, we have(21)ε2n+1=I+aG1¯−1c1I−aG2¯ε2n+∑j=12n−1cj+1ε2n+1−j+l2nε0,ε2n+12=I+aG1¯−1c1I−aG2¯ε2n+∑j=12n−1cj+1ε2n+1−j+l2nε02≤I+aG1¯−12c1I−aG2¯ε2n+∑j=12n−1cj+1ε2n+1−j+l2nε02≤I+aG1¯−12c1I−aG2¯2ε02+∑j=12n−1cj+1ε02+l2nε02=c1−λ3+∑j=12n−1cj+1+l2n1+λ2ε02=c1−λ3+l11+λ4ε02.
Whenλ3<c1, we have c1−λ3+l1/1+λ4=c1−λ3+l1/1+λ4=1−λ3/1+λ4<1. When λ3≥c1, let λ3=r1c1, here r1≥1; from c1−λ3/1+λ4≤1, we can get r1≤2, c1≤λ3≤2c1; and c1−λ3+l1/1+λ4=λ3−c1+l1/1+λ4≤c1+l1/1+λ4≤1. Therefore, ε2n+12≤ε02.
Whenm=2n+2,(22)ε2n+2=I+aG2¯−1c1I−aG1¯ε2n+2+∑j=12ncj+1ε2n+1−j+l2n+1ε0,ε2n+22=I+aG2¯−1c1I−aG1¯ε2n+2+∑j=12ncj+1ε2n+1−j+l2n+1ε02≤I+aG2¯−12c1I−aG1¯ε2n+2+∑j=12ncj+1ε2n+1−j+l2n+1ε02≤I+aG2¯−12c1I−aG1¯2ε02+∑j=12ncj+1ε02+l2n+1ε02=c1−λ1+∑j=12ncj+1+l2n+11+λ4ε02=c1−λ1+l11+λ2ε02≤ε02.
In the same way, we can getε2n+22≤ε02.
To sum up,εn2≤ε02,n=1,2,3,…, and the following theorem can be obtained.Theorem 2.
The pure alternative segment explicit-implicit (PASE-I) difference scheme (12) for time fractional reaction-diffusion equation is unconditionally stable.


## 3.3. Convergence of PASE-I Difference Scheme

First, the accuracy of PASE-I scheme is analyzed. The PASE-I scheme is alternatively applied explicit and implicit segments, not only in the direction of space, but also in the direction of time. Examination of grid pointsxi,t2n+1 and xi,t2n+2. The classical explicit scheme is(23)D0Ctαuxi,t2n+1=1h2Vi−12n−2Vi2n+Vi+12n−pVi2n+fi2n,(24)Vi2n+1=aVi−12n+1−2a−l1−bVi2n+aVi+12n+∑j=12n−1lj−lj+1Vi2n−j+l2nVi0+cfi2n.The classical implicit scheme is(25)D0Ctαuxi,t2n+2=1h2Vi−12n+2−2Vi2n+2+Vi+12n+2−pVi2n+2+fi2n+2,(26)−aVi−12n+2+1+2a+bVi2n+2−aVi+12n+2=1−l1Vi2n+1+∑j=12nlj−lj+1Vi2n+1−j+l2n+1Vi0+cfi+12n+2.We expand (24) and (26) at two points xi,t2n+1 and xi,t2n+2 for Taylor, respectively. The truncation error is recorded as T1k,h,≤2n+11−αCvk1−α/Γ2−αT2k,h. It is known that the precision of D0Ctαuxi,t2n+1 is 2−α order [4, 9]. Let CV=max0≤t≤t2n+1∂2Vx,t/∂t2, we have(27)∂α+1Vxi,t2n+1∂tα+1=1Γ1−α∫0t2n+1∂2Vxi,τ∂τ2dτt2n+1−τα=1Γ1−α∑j=02n∫jkj+1kdτt2n+1−τα≤1Γ1−αmax0≤t≤t2n+1∂2Vx,t∂t2∑j=02n∫jkj+1kdτt2n+1−τα=CVΓ1−α∑j=02n∫jkj+1kdτt2n+1−τα=CVk1−α1−αΓ1−α∑j=02n2n+1−j1−α−2n−j1−α≤2n+11−αCVk1−αΓ2−α,T1k,h=∂αVxi,t2n+1∂tα−1h2Vi−12n−2Vi2n+Vi+12n+pVi2n=∂αVx,t∂tα−Vxx+pV+kVxxt−k22Vxxtt−h212Vxxxx−pkVt+pk22Vtt+Ok2−α+h2,T2k,h=∂αVxi,t2n+2∂tα−1h2Vi−12n+2−2Vi2n+2+Vi+12n+2+pVi2n+2=∂αVxi,t2n+1∂tα+k∂α+1Vxi,t2n+1∂tα+1−Vxx+pV−kVxxt−k22Vxxtt−h212Vxxxx+pkVt+pk22Vtt+Ok2−α+h2.The two termskVxxt,pkVt of T1k,h have the same form as two terms kVxxt,pkVt of T2k,h, but the opposite of the symbol. When the explicit format and implicit format are applied alternately, the errors of the two terms will be offset. And k∂α+1Vxi,t2n+1/∂tα+1≤Ck2−α,C>0. When applying explicit and implicit formats alternately, the accuracy is 2 order in space direction and 2−α in time direction. Similarly, the same results can be obtained when implicit and explicit formats are applied alternately. Therefore, the truncation of PASE-I scheme is Rin=Ok2−α+h2.Whenm=1, the first layer is calculated applying the explicit format.(28)Vi1=aVi−10+1−2a−bVi0+aVi+10+cfi0.Whenm=2, the second layer is calculated applying the implicit format.(29)Vi1=aVi0+1−2a−bVi0+aVi+10+cfi0,(30)−aVi−12+1+2a+bVi2−aVi+12=c1Vi1+l1Vi0+cfi2.Substitute (29) into (30) to eliminate Vi1:(31)−aVi−12+1+2a+bVi2−aVi+12=c1aVi−10+c11−2a−b+l1Vi0+c1aVi+10+cc1fi0+cfi2.Whenm>2, we substitute the explicit format (24) into the implicit format (26) to eliminate Vi2n+1:(32)−aVi−12n+2+1+2a+bVi2n+2−aVi+12n+2=c1aVi−12n+c1c1−2a−bVi2n+c1aVi+12n+c2Vi2n+∑j=12n−1c1cj+1+cj+2Vi2n−j+c1l2n+l2n+1Vi0+cc1fi2n+cfi2n+2.For convenience, assume thatV˜in is the exact solution of time FRDE, and Vin is the numerical solution of the difference scheme, Ein=Vin−V˜in. For n=0,1,2,…,N, (28), (31) and (32) we have(33)Ei1=aEi−10+1−2a−bEi0+aEi+10+kαRi1,−aEi−12+1+2a+bEi2−aEi+12=c1aEi−10+c11−2a−b+l1Ei0+c1aEi+10+kαRi2,−aEi+12n+2+1+2a+bEi2n+2−aEi−12n+2=c1aEi−12n+c1c1−2a−bEi2n+c1aEi+12n+c2Ei2n+∑j=12n−1c1cj+1+cj+2Ei2n−j+c1l2n+l2n+1Ei0+Ri2n+2,where Rin=Ok2−α+h2,E0n=EMn=Ei0=0,n=0,1,…,N,i=0,1,…,M.Define two grid functions as follows:(34)Enx=Eln,x∈xl−1/2,xl+1/2,l=1,2,…,M−1,0,x∈0,h2∪L−h2,L,Rnx=Rln,x∈xl−1/2,xl+1/2,l=1,2,…,M−1,0,x∈0,h2∪L−h2,L,where Rn=R1n,R2n,…,RM−1nT,En=E1n,E2n,…,EM−1nT,n=0,1,2,…,N. We do Fourier expansion on Enx and Rnx:(35)Enx=∑j=−∞+∞Ejnei2jπx/L,i=−1,n=0,1,…,N,j=0,1,…,M,Rnx=∑j=−∞+∞Rjnei2jπx/L,i=−1,n=0,1,…,N,j=0,1,…,M,where Ejn=1/L∫0LEnei2jπx/Ldx,Rjn=1/L∫0LRnei2jπx/Ldx,j=±0,±1,…Define normEnx=∑j=1M−1hEjn2, Rnx=∑j=1M−1hRjn2.We have the following formula:(36)Enx2=∑j=1M−1hEjn2=L∑j=−∞+∞Ejn2,n=0,1,…,N,(37)Rnx2=∑j=1M−1hRjn2=L∑j=−∞+∞Rjn2,n=0,1,…,N.SetEln=eneiλlh,Rln=rneiλlhλ=2jπ/L, and substitute (33), we have(38)e1=kαr1,−ae−iλh+1+2a+b−aeiλhe2=kαr2,−ae−iλh+1+2a+b−aeiλhe2n+2=c1ae−iλh+c1c1−2a−b+c1aeiλh+c2e2n+∑j=12n−1c1cj+1+cj+2e2n−j+kαr2n+2.Becausesin2λh/2=−1/4e−iλh−2+eiλh, (38) can be converted to(39)e1=kαr1,1+b+4asin2λh2e2=kαr2,1+b+4asin2λh2e2n+2=c12+c2−c1b−4c1asin2λh2e2n+∑j=12n−1c1cj+1+cj+2e2n−j+kαr2n+2.Lemma 2.
Letek be the solution of (39), ∃C1>0 makes e2n+2≤C1l2n+1−1kαr1,n=0,1,…Proof.
FromRjn=Ok2−α+h2,j=0,1,…,M,n=0,1,…,N, we know ∃C2>0, makes(40)Rjn≤C2k2−α+h2.
According toRnx2=∑j=1M−1hRjn2=L∑j=−∞+∞Rjn2,n=0,1,…,N, we have(41)Rn≤C2Lk2−α+h2.
FromRln=rneiλlhλ=2jπ/L and (37), we can know(42)rn≡Rjn≤C3Rj1≡C3r1,n=2,3,…,N,where C3 is a positive number. From (39), we have(43)e1=kαr1=l0−1kαr1,e2=kαr21+b+4asin2λh/2≤C3r1=C3l1−1r1.
Supposeek≤C1lk−1−1kαr1k≤2n,C1=max1,C3, we have(44)1+b+4asin2λh2e2n+2=c12+c2−c1b−4c1asin2λh2e2n+∑j=12n−1c1cj+1+cj+2e2n−j+kαr2n+2≤c1b+4c1asin2λh2e2n+∑j=12nc1cj+cj+1e2n+1−j+C3kαr1≤c1b+4c1asin2λh2e2n+∑j=12nc1cj+cj+1C1l2n−j−1kαr1+C3kαr1≤c1b+4c1asin2λh2C1l2n−1−1kαr1+∑j=12nc1cj+cj+1C1l2n−j−1kαr1+C3kαr1≤c1b+4c1asin2λh2C1l2n+1−1kαr1+∑j=12nc1cj+cj+1C1l2n+1−1kαr1+C3kαr1≤c1b+4c1asin2λh2+c1+∑j=12ncj+1+l2n+1−c1l2nC1l2n+1−1r1≤1+c1b+4c1asin2λh2C1l2n+1−1kαr1≤1+b+4asin2λh2C1l2n+1−1kαr1.
Therefore, we havee2n+2≤C1l2n+1−1kαr1. Lemma 2 proof is complete.
Becausenαkα≤Tαn=1,2,…,N,limn⟶∞ln−1/nα=limn⟶∞n−α/n+11−α−n1−α=limn⟶∞n−1/1+1/n1−α−1=1/1−α, we have(45)e2n+2≤C1l2n+1−1kαr1≤C1l2n+1−12n+1−α2n+1αkαr1≤C11−α2n+1−αkαr1≤C11−αTαr1.
From (36), (37) and (41) we can get(46)En≤C11−αTαR1≤C11−αTαC2Lk2−α+h2=Ck2−α+h2.
To sum up, the following theorem can be obtained.Theorem 3.
The pure alternative segmented explicit-implicit (PASE-I) scheme (12) of the time fractional reaction-diffusion equation is unconditionally convergent, and satisfying(47)Vin−V˜in≤Ck2−α+h2,i=1,2,…,M,n=1,2,…,N.here C is a positive number.


## 3.4. PASI-E Scheme of Time Fractional Reaction-Diffusion Equation

Similar to the construction of PASE-I scheme, we give the PASI-E scheme of FRDE.When solving the odd time layer, the calculation is performed according to the rules of “classical implicit-classical explicit-classical implicit”. On the even time layer, the rule of “classical explicit-classical implicit-classical explicit” is applied. Then the PASI-E scheme for solving the time FRDE is obtained as follows:(48)I+aG2¯V1=I−aG1¯V0+g1+l0V0+cf1,I+aG1¯Vn+1=c1I−aG2¯Vn+gn+1+∑j=1n−1cj+1Vn−j+lnV0+cfn+1,n=1,3,…,I+aG2¯Vn+2=c1I−aG1¯Vn+1+gn+1+∑j=1ncj+1Vn+1−j+ln+1V0+cfn+1.The implicit scheme of the first time layer is unconditionally stable and convergent. Starting from the second time layer, with reference to the proof analysis of the PASE-I scheme, there are the following theorems.Theorem 4.
The pure alternative segmented implicit-explicit (PASI-E) scheme (48) of the time fractional reaction-diffusion equation is uniquely solvable, unconditionally stable and convergent, satisfying(49)Vin−V˜in≤Ck2−α+h2,i=1,2,…,M,n=1,2,…,N.here C is a positive number.


## 4. Numerical Experiments

The numerical examples are based on Intel Core i5-2400 CPU@3.10GHz, which is a quad-core CPU. By using the Parallel Computing Toolbox of Matlab R2014a to run the PASE-I and PASI-E parallel algorithms and complete numerical experiments.Example 1.
In equation (1) take p=1, L=2, T=1, ux,0=0, and fx,t=2/Γ3−αx2−xt2−α+t2x2−x+2t2. Then, the exact solution of the equation is ux,t=t2x2−x [35, 41].
TakeM=81,N=1000,α=0.7,B=5,l=M−1/B=16. Figures 2–4 give the surfaces of exact solution, the numerical solutions’ surface of PASE-I and PASI-E schemes, respectively. It can be seen that they are of the same shape and the solution surfaces PASE-I and PASI-E formats are smooth. Therefore, the PASE-I and PASI-E schemes are high-precision differential methods for solving time FRDE.
TakeM=81,N=1000,α=0.7. At t=0.4, we compare the numerical solutions of the PASE-I and PASI-E schemes with exact solution, the numerical solutions of implicit scheme, E-I and I-E schemes [41]. The calculation results are shown in Table 1. It can be seen that the numerical solutions of PASE-I and PASI-E schemes are very close to the exact solution.
In order to verify the stability of PASE-I and PASI-E methods, the variation of relative error over time is investigated. The Sum of Relative Error for every Time level (SRET) is defined as(50)SRETj=∑i=1MV˜ij−VijV˜ij.
As shown in Figure5, for M=81,N=1000,α=0.7, it can be seen that the SRET of PASE-I and PASI-E solutions is relatively large when the number of time layers is small. But with the increase of time layers, SRET decreases rapidly and is bounded. Therefore, Figure 5 shows that the PASE-I and PASI-E methods for time FRDE are computationally stable.
At the same time, Figure5 shows that the SRET of the numerical solution of PASE-I (PASI-E) format is large near the initial moment, which is similar to other formats. Here we will encrypt the grid near the initial moment to improve the accuracy of PASE-I (PASI-E) format. Let t=1, M=81,N=1000. The grids of the interval 0,0.2 are encrypted by a factor of 2 to 5, respectively. From Figure 6, we can know that with the encryption of the initial time grid, the SRET of the initial moment decreases obviously. When the local encryption grid is 3 times, the SRET of initial time decreases about from 50 to 20. The local encryption technology can improve the numerical simulation error.
As can be seen from Table2, the computation time (CPU time) of PASE-I or PASI-E scheme increases linearly with the multiples of local grid refinement. We can choose the appropriate grid encryption multiples according to the accuracy and time requirements of the numerical simulation process. We can also see that the computation time of the two parallel difference schemes is very close.
In order to better verify the accuracy of PASE-I and PASI-E methods, the distribution of errors on spatial grid points is investigated below. The Difference Total Energy (DTE) is defined as(51)DTEi=12∑j=1NV˜ij−Vij2.
From Figure7, we can see that when M=81,N=1000,α=0.7, the DTE of PASE-I and PASI-E schemes are between 0 and 1.4e-04 on the whole spatial grid points. In other words, the numerical solutions of PASE-I and PASI-E schemes are very close to the exact solution. Figure 7 shows that that the PASE-I and PASI-E methods have good accuracy for solving time FRDE.
Because the numerical solutions of the PASE-I scheme and PASI-E scheme are very close, we only consider PASE-I scheme. The numerical experiments of the order of time convergence and the order of space convergence are given below for the implicit scheme, E-I scheme, and PASE-I scheme. Define the normL2 of error, the order of space convergence Order1, and the order of time convergence Order2 as follows [26, 42]:
(52)E2h,k=max0≤n≤NV˜in−Vin2,Order1=log2E22h,kE2h,k,Order2=log2E2h,2kE2h,k.
First, we verify the order of spatial direction convergence order of the PASE-I scheme. SetM=21,41,81,161,321,N≈M2/4k≈h2, α was taken as 0.3, 0.5 and 0.7 respectively. Table 3 indicates that the order of spatial direction convergence of the PASE-I scheme is approximately second order. Secondly, we verify and calculate the time direction convergence order of PASE-I scheme. Set N=400,800,1600,3200,6400 and space step h=2/81M=81, α was taken as 0.3, 0.5, and 0.7 respectively. We can see that the time direction convergence order of PASE-I scheme is approximately 2−α order from Table 4. It is shown that the PASE-I scheme in this paper is equivalent to the existing implicit and E-I scheme in spatial and temporal accuracy, and the results of numerical experiments are consistent with theoretical analysis.Figure 2
The surface of exact solution.Figure 3
The solution surface of PASE-I scheme.Figure 4
The solution surface of PASI-E scheme.Table 1
Comparison of numerical solutions and exact solution.α=0.7,t=0.4.
x0.40.81.21.6Exact solution0.1024000.1536000.1536000.102400Implicit scheme0.1023960.1537500.1537500.102396E-I scheme [41]0.1025260.1533430.1533430.102526I-E scheme [41]0.1025220.1535800.1535800.102522PASE-I scheme0.1023390.1532880.1532880.102339PASI-E scheme0.1024740.1532620.1536130.102356Figure 5
SRET of PASE-I and PASI-E schemes’ solutions.Figure 6
SRET of PASE-I and PASI-E schemes when local grid refinement. (a) Local grid refinement 2 times. (b) Local grid refinement 3 times. (c) Local grid refinement 4 times. (d) Local grid refinement 5 times.
(a)(b)(c)(d)Table 2
The CPU time of PASE-I and PASI-E schemes after grid local refinement (unite: second).
Grid local refinement1 time2 times3 times4 times5 timesCPU time of PASE-I2.542.633.083.784.51CPU time of PASI-E2.402.643.103.774.55Figure 7
The distribution of DTE of two schemes in spatial grid points.Table 3
The order of spatial convergence of implicit, E-I, and PASE-I schemesk≈h2.
αMImplicit schemeE-I scheme [41]PASE-I schemeE2h,kOrder1E2h,kOrder1E2h,kOrder10.3217.748165e − 022.563728e − 021.645750e − 02411.982738e − 021.9663606.165582e − 032.0559343.591584e − 032.196053814.988066e − 031.9909421.496886e − 032.0422728.037972e − 042.15971716111.249155e − 031.9975283.666129e − 042.0296351.829138e − 042.1356673213.124381e − 041.9993109.038107e − 052.0201654.203903e − 052.1213630.5217.686890e − 022.349958e − 021.540929e − 01411.974082e − 021.9612185.775474e − 032.0246233.959968e − 021.960240814.976324e − 031.9880301.429413e − 032.0145179.648636e − 032.0370921611.247601e − 031.9959243.554208e − 042.0078222.376859e − 032.0212693213.122356e − 041.9984498.860591e − 052.0040545.895655e − 042.0113320.7217.521184e − 022.255208e − 021.194496e − 01411.942855e − 021.9527815.639231e − 031.9996902.855516e − 022.064579814.920092e − 031.9814211.411534e − 031.9982356.730069e − 032.0850581611.237740e − 031.9909773.534096e − 041.9978511.613612e − 032.0603283213.105362e − 041.9948758.846820e − 052.0059513.912510e − 042.044127Table 4
The order of time convergence of implicit, E-I, and PASE-I schemesh=1/40.
αNImplicit schemee − I scheme [41]PASe − I schemeE2h,kOrder2E2h,kOrder2E2h,kOrder20.34002.550088e − 06—9.449930e − 04—3.078795e − 03—8007.988661e − 071.6745222.857286e − 041.7256589.740270e − 041.66033216002.493529e − 071.6797658.536685e − 051.7428973.027275e − 041.68594232007.761226e − 081.6838332.525755e − 051.7569629.422965e − 051.68376764002.410372e − 081.6870287.407991e − 061.7695602.968082e − 051.6666500.54001.445964e − 05—4.636766e − 04—8.550850e − 03—8005.130327e − 061.4949091.664891e − 041.4776913.122569e − 031.45333416001.818307e − 061.4964555.762628e − 051.5306291.151259e − 031.43952132006.439744e − 071.4975201.925467e − 051.5815194.302524e − 041.41995764002.279541e − 071.4982606.216571e − 061.6310171.636565e − 041.3945120.74006.390559e − 05—6.535001e − 04—1.748703e − 02—8002.600370e − 051.2972252.704299e − 041.2729337.230957e − 031.27402816001.057243e − 051.2984101.114194e − 041.2792542.969345e − 031.28404232004.296485e − 061.2990784.533637e − 051.2972601.216933e − 031.28689564001.745569e − 061.2994591.655564e − 051.4533464.988235e − 041.286649Example 2.
In equation (1) take p=0, L=2π, T=1, ux,0=0 and fx,t=tα+Γ1+α+μ/Γ1+μtμsinx. The exact solution is ux,t=tα+μsinx [35, 43].
In the Example2, we consider the time FRDE under the condition of weak regularity. We take M=41,α=0.8,μ=0.1and0.5 for Example 2. For Table 5, when α+μ>2−α, the orders of time convergence of PASE-I and PASE-I schemes are approximately 2−α order. When α+μ<2−α, the orders of time convergence of PASE-I and PASE-I schemes are generally lower than 2−α order. From Table 5, we can see that the PASE-I and PASI-E schemes have satisfactory robustness for solving time FRDEs.
Next, we consider the performance of PASE-I and PASI-E schemes in the case ofα=0.8,μ=0.1. We take M=41,N=100. Figure 8 gives the absolute error of implicit scheme, E-I scheme, and PASE-I and PASI-E schemes, respectively. The absolute error is less than 8×10−3. Therefore, the PASE-I and PASI-E schemes can numerically solve the time FRDEs effectively.
In order to better compare the computational efficiency of PASE-I scheme and implicit and E-I schemes, the speedupSp is defined. We define the speedup as Sp=T1/Tp (T1 is the CPU time of implicit or E-I scheme, Tp is the CPU time of PASE-I parallel scheme) [44]. The speedup of PASE-I scheme and implicit and E-I schemes are Sp1,Sp2, respectively. Set B=5 and the time layer as N=500, and select the number of space grids M=101,301,501,701,901,1101,1301. The calculation results are shown in Table 6.
From Table6, we can see that the CPU time of three kinds of difference schemes for time FRDE increases with the increase of space grid number. The speedup Sp1,Sp2 are greater than 1 and gradually increase. The CPU time of E-I scheme is about 60% of the implicit scheme. When number of grid points M≤100, the PASE-I scheme takes about the same CPU time as E-I scheme. The reason is the influence of data communication among different workers. Compared with the implicit scheme or E-I scheme, the PASE-I scheme has more and more advantages in computing time with the increase in the number of spatial grids. We can see the PASE-I has obvious advantages in computation time. So the PASE-I parallel scheme of time FRDE has higher computational efficiency than the serial implicit scheme, E-I scheme, especially for numerical simulation of long time history or large computational domain.Table 5
Order of time convergence of PASE-I and PASE-I schemesα=0.8.
μNPASE-I schemePASE-I schemeE2h,kOrder2E2h,kOrder20.11005.932488e − 3—5.924955e − 3—2003.200570e − 30.8903083.194804e − 30.8910764001.749330e − 30.8715261.745987e − 30.8716848009.950064e − 40.8140259.932580e − 40.81380216006.077651e − 40.7111916.068948e − 40.7107210.51002.229026e − 3—2.244422e − 3—2009.960646e − 41.1621029.993384e − 31.1672994004.357703e − 41.1926714.364595e − 41.1951258001.929037e − 41.1756861.930400e − 41.17694716008.913430e − 51.1138288.915908e − 51.114446Figure 8
Absolute error of the four kinds of difference schemes. (a) Absolute error of implicit scheme. (b) Absolute error of E−I scheme. (c) Absolute error of PASE−I scheme. (d) Absolute error of PASI−E scheme.
(a)(b)(c)(d)Table 6
Comparison of CPU time and speedup among implicit, E-I, and PASE-I schemes (unit: second).
Number of grid points10130150170190111011301Implicit scheme0.812.224.889.4018.4029.9139.90E-I scheme [41]0.691.553.736.7712.1820.4727.99PASE-I scheme0.590.740.911.281.983.434.44Sp11.373.005.367.349.298.728.98Sp21.162.094.095.286.155.966.30


## 5. Conclusion

The PASE-I and PASI-E parallel difference methods for solving time FRDE are unconditionally stable and convergent. Both schemes are second order in space precision and2−α order in time precision. At the same time, numerical experiments show that the PASE-I and PASI-E schemes have better computational accuracy and computational efficiency than implicit and E-I difference schemes. In particular, when the space points are large enough, PASE-I and PASI-E methods have obvious localization features in computing and communication, which are suitable for large-scale parallel computing systems with distributed storage. Compared with implicit scheme and E-I (I-E) scheme, the computational efficiency of PASE-I (PASI-E) scheme is greatly improved. PASE-I and PASI-E schemes are efficient algorithms for solving time FRDEs.



---

*Source: 4524387-2020-08-26.xml*