Numerical Algorithms for the Fractional Diffusion-Wave Equation with Reaction Term

and Applied Analysis 3 where a l , b l , c l , and d l are constants and θ is an arbitrary parameter which can improve the accuracy of the numerical method. To obtain the expressions of the coefficients in (11), we denote S l (x l , t k ) = u (x l , t k ) , S l (x l+1 , t k ) = u (x l+1 , t k ) , S 󸀠󸀠 l (xl, tk) = Q (xl, tk) , S 󸀠󸀠 l (xl+1, tk) = Q (xl+1, tk) . (12) Through the straightforward calculations, we obtain the following expressions: a l = (1 − 2 cosh (θh)) Q (xl, tk) + Q (xl+1, tk) 2θ sinh (θh) + u (x l+1 , t k ) − u (x l , t k ) 2 sinh (θh) , b l = Q (xl, tk) θ , c l = Q (x l , t k ) − Q (x l+1 , t k ) 2θ sin (θh) + u (x l+1 , t k ) − u (x l , t k ) 2 sin (θh) , d l = u (xl, tk) − Q (x l , t k ) θ , (13) in which l = 0, 1, . . . ,M − 1, k = 0, 1, . . . , N. From (11), one has


Introduction
The realm of fractional differential equations has drawn increasing attention and interest due to their important applications in biology, physics, chemistry, biochemistry, medicine, and engineering [1][2][3][4][5][6][7][8].Generally speaking, the analytical solutions of most fractional differential equations are not easy, even impossible, to obtain, so seeking numerical solutions of these equations becomes more and more important.Till now, there have been some studies in this respect.For example, Cui [9] proposed a compact finite difference method for the fractional diffusion equation.Chen et al. [10] analyzed an implicit difference scheme for the subdiffusion equation and proved the unconditional stability and convergence.Li et al. [11] constructed some novel methods for the fractional calculus and fractional ordinary differential equation.Liu et al. [12] considered the numerical method and analytical technique for the modified anomalous subdiffusion equation with a nonlinear source term.Sousa [13] gave three finite difference schemes for a fractional advection diffusion problem.In [14], she furthermore presented an improved difference scheme to enhance the spatial accuracy.H. Wang and K. Wang [15] developed a fast finite difference method for the fractional diffusion equations.Yuste and Acedo [16] derived an explicit finite difference scheme and gave a new von Neumann-type stability analysis for the fractional diffusion equations.Roop and his coworkers pioneered in considering the fractional differential equation by using the finite element method in [17][18][19], where the temporal derivative is the classical derivative and the spatial derivative is the fractional derivative.Li et al. [20] used the finite difference/finite element mixed method to numerically solve the nonlinear subdiffusion and superdiffusion partial differential equation, where both temporal and spatial derivatives are fractional derivatives.Zhang et al. [21] obtained a numerical method for the symmetric space fractional partial differential equations by using Galerkin finite element method and a backward difference technique.In [22], Zheng et al. established a fully discrete approximation for the nonlinear spatial fractional Fokker-Planck equation, where the discontinuous Galerkin finite element approach was utilized in time domain and the Galerkin finite element approach was utilized in spatial domain.Yang et al. used the matrix transform method for the Riesz space fractional diffusion and advection-dispersion equations and the time and space fractional Fokker-Planck equations [23,24], respectively.Recently, the spectral method was used to approximate the fractional calculus and the fractional differential equations [25,26].Very recently, Ding and Li [27] proposed two kinds of finite difference schemes for the reaction-subdiffusion equation.

Preliminary Knowledge
In the present section, we introduce some important definitions and lemmas which are used later on.The following definitions and lemmas can be found in [30].
Definition 1.The left Riemann-Liouville integral operator  − 0, of order  > 0 is defined in the following form: where Γ(⋅) is the Euler gamma function.
In view of the continuity of the first-order derivative at the common nodes (  , (  ,   )), that is, S  −1 (  ,   ) = S   (  ,   ), we get the following relation: where and  = ℎ.

Numerical Algorithms
In the present section, we introduce two schemes for the fractional diffusion-wave equation ( 1) together with the homogeneous initial value conditions (2) and boundary values conditions ( 3) and (4).

Numerical Algorithm I. Based on Lemma 4, (1) can be written as
Applying the initial value conditions (2) yields The Riemann-Liouville derivative (with homogeneous initial value conditions) can be approximated by the following scheme: There exist various choices of the generating functions which can lead to different approximation order .
Let   (, ) be the generating function with coefficients  ()  , ; that is, If we take the generating function in the form that is, then we can obtain the following first-order scheme: in which and these coefficients have the following recursive relations: If we take the generating function as that is, which leads to order  = 2 [30], then Here we introduce two methods to determine the corresponding coefficients  ()  2, of the formula (32).One of the most general methods is to use the fast Fourier transform to calculate them [30].
In the following, we give another simple but interesting method to obtain the coefficients  () 2, .By using Abstract and Applied Analysis 5 and ( 26), one gets ) ) Comparing the previous equation with (31), one obtains ) ) which has the following recursive relation: ) ) ) ) ) ) In summary, we have the following lemma.
Finally, let    be the numerical approximation of (  ,   ); substituting the expansions (20) and (32) in (22) and removing the higher order term, one obtains Numerical Algorithm I as follows: for  = 1, . . .,  − 1,  = 1, . . ., .The initial and boundary value conditions can be discretized by where ) , Then we can give the following matrix form of the difference scheme (38): where Theorem 7. The difference equation (38) is uniquely solvable.
Theorem 8.The local truncation error of the difference scheme Proof.We now define the local truncation error as Applying ( 20), (22), and (32), we have This finishes the proof.
In the following, we discuss the stability of Numerical Algorithm I by the fractional Fourier method [16].
Let    be the approximate solution of (38) and define    =    −    ,  = 1, 2, . . .,  − 1,  = 0, 1, . . ., , So, we can easily obtain the following round-off error equation: We suppose that the solution of (48) has the following form: where  is any of the real spatial wave numbers supported by the lattice.Substituting (49) into (48) gives If we write and assume that  ≡ () is independent of time, then we obtain the following expression of the amplification factor: When || > 1 for some , then the temporal factor of the solution grows to infinity according to (51); the numerical method is unstable.Considering the extreme value  = −1, we find that the numerical method is stable when Although  , depends on , we can estimate it by its limit value In this case, the stability condition becomes Note that Ŝ is always negative, and then (55) holds for all  1 > 0 and  2 > 0. That is to say, Numerical Algorithm I (see (38)) is unconditionally stable.
Then the matrix form of Numerical Algorithm II (61) is written as where Therefore, the solution of (61) exists and is uniquely solvable.
In order to get the local truncation error analysis, we need the following lemma.Proof.We now define Abstract and Applied Analysis 9 By using ( 20), ( 32), (59), and (60), one gets Noticing that 1 ≤  ≤ , then so we can get This completes the proof.
Next, we study the stability of Numerical Algorithm II (see (61)).
As before, we can easily obtain the following round-off error equation: Substituting ( 49) into (76) yields From ( 51) and (77), we have Considering the extreme value  = −1, we can obtain the following stability condition from (78):  For large enough , we can estimate S, by its limit value Because the sine function is bounded by 1, from (81) one finds that Numerical Algorithm II is to be stable if S ≥ 1, that is, Furthermore, stability condition (82) can be rewritten as where   > 0 and   > 0 are diffusion and reaction coefficients, respectively.In addition, the numerical check of the validity of the stability condition (83) will be discussed in Section 4.
The analytical solution of this equation is (, ) = () 3 sin((/2)).The initial conditions and the boundary conditions are satisfied with the exact solution (, ) given previously.
On one hand, we use Numerical Algorithm I (see (38)) to solve the previous equation.At first, we give the comparisons of the analytical and numerical solutions for different order, , , and ℎ in Figures 1 and 2, respectively, and we observe that the numerical solution is in line with the analytical solution and the velocity  increases with the increase of order .
Meanwhile, the errors are shown in Figure 3 with  = 1.5.
Finally, the temporal and spatial convergence orders are listed in Table 1.
Figures 4 and 5 display that the velocity  increases with the increase of order  for different  and ℎ. Figure 6 presents the error between the analytical solution and numerical solution of (60).Table 2 lists the temporal convergence orders.
Finally, we checked the stability condition given in (83) in the following way.Firstly, we choose ℎ = 0.04,  = 0.005, and  = 1.6, which satisfy stability condition (83), the comparison result of the numerical solution and analytical solution is shown in Figure 7. Secondly, we choose ℎ = 0.01,  = 0.01, and  = 1.6, and then these parameters do not satisfy stability condition (83); the comparison result of the numerical solution and analytical solution is shown in Figure 8. From Figures 7 and 8, we declare that stability condition (83) is valid.

Conclusion
In this paper, we construct two difference schemes for solving the fractional diffusion-wave equation with reaction term.It is proved that Numerical Algorithm I (38) is unconditionally stable and Numerical Algorithm II (61) is conditionally stable by using the fractional Fourier method.The local truncation errors of two difference schemes are both O( 2 +ℎ 4 ), which are the best results till now.Finally, the numerical results given in this paper show the effectiveness of the derived numerical algorithms.

Figure 3 :
Figure 3: The error of the analytical solution and the numerical solution for  = 1.5 by Numerical Algorithm I (see (38)) with  = 0.0025, ℎ = 0.01.

Table 1 :
The maximum error and temporal and spatial convergence orders by Numerical Algorithm I (see (38)).

Table 2 :
The maximum error and temporal convergence order by Numerical Scheme II (see (61)).