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

Abstract

Two numerical algorithms are derived to compute the fractional diffusion-wave equation with a reaction term. Firstly, using the relations between Caputo and Riemann-Liouville derivatives, we get two equivalent forms of the original equation, where we approximate Riemann-Liouville derivative by a second-order difference scheme. Secondly, for second-order derivative in space dimension, we construct a fourth-order difference scheme in terms of the hyperbolic-trigonometric spline function. The stability analysis of the derived numerical methods is given by means of the fractional Fourier method. Finally, an illustrative example is presented to show that the numerical results are in good agreement with the theoretical analysis.

1. 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–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–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.

As far as we know, there are some numerical methods for the subdiffusion equations [9, 10, 12, 16]; however, there seem to be very limited ones for the fractional diffusion-wave (seldom called “superdiffusion”) equation. The possible reason is that there may be no better method to analyze the stability and convergence. And in the studies available, the accuracy of the temporal direction for the subdiffusion and fractional diffusion-wave equations is often less than second order. This motivates us to construct some higher order methods.

In this paper, we study the fractional diffusion-wave equation with a reaction term in the following form:_𝐶𝐷^𝛼_0,𝑡𝑢(𝑥,𝑡)=𝐾_𝛼𝜕²𝑢(𝑥,𝑡)𝜕𝑥²−𝐶_𝛼𝑢(𝑥,𝑡)+𝑓(𝑥,𝑡),0≤𝑡≤𝑇,0<𝑥<𝐿. (1)

Roughly speaking, the previous fractional differential equation is obtained from the classical diffusion-wave equation by replacing the second-order time derivative by a fractional derivative of order 𝛼(1<𝛼<2). It has been proposed to deal with disordered media to comb structures, dielectrics and semiconductors, and viscoelastic problems, for example, in the description of the propagation of stress waves in viscoelastic solids; see [28, 29] and references therein. In particular, in order to carry out the numerical comparisons, we study the following initial boundary values: 𝑢()()𝑥,0=0,𝜕𝑢𝑥,0𝜕𝑡=0,0<𝑥<𝐿,𝑢(0,𝑡)=𝜑₁(𝑡),0≤𝑡≤𝑇,𝑢(𝐿,𝑡)=𝜑₂(𝑡),0≤𝑡≤𝑇, (2)(3)(4) where 𝑓(𝑥,𝑡), 𝜑₁(𝑥), and 𝜑₂(𝑥) are sufficiently smooth functions and 𝐾_𝛼>0 and 𝐶_𝛼>0 are diffusion and reaction coefficients, respectively. _𝐶𝐷^𝛼_0,𝑡𝑢(𝑥,𝑡) is the Caputo derivative of order 𝛼(1<𝛼<2), which will be introduced in the following.

The remainder of this paper is organized as follows. In Section 2, we present preliminary knowledge, which is necessary for our study. In Section 3, we develop two numerical methods. Meanwhile, we investigate the solvability and stability by the matrix method and fractional Fourier method, respectively. In Section 4, a numerical experiment is carried out which supports the theoretical analysis. Finally, the paper concludes with a summary in the last section.

2. 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: 𝐷^−𝛼_0,𝑡𝑢()=1𝑥,𝑡Γ(𝛼)^𝑡₀𝑢𝑥,𝑦𝑑𝑦𝑡−𝑦^1−𝛼, (5) where Γ(⋅) is the Euler gamma function.

Definition 2. The left Caputo derivative operator _𝐶𝐷^𝛼_0,𝑡 of order 𝛼>0 is defined in the following form: _𝐶𝐷^𝛼_0,𝑡𝑢()=1𝑥,𝑡Γ()∫𝑚−𝛼^𝑡0𝜕^𝑚𝑢𝑥,𝑦𝜕𝑦^𝑚𝑑𝑦𝑡−𝑦^{𝛼−𝑚+1},𝑚−1<𝛼<𝑚∈ℤ⁺,𝜕^𝑚𝑢(𝑥,𝑡)𝜕𝑡^𝑚,𝛼=𝑚. (6)

Definition 3. The left Riemann-Liouville derivative operator _RL𝐷^𝛼_0,𝑡 of order 𝛼>0 is defined by _RL𝐷^𝛼_0,𝑡𝑢()=1𝑥,𝑡Γ()𝜕𝑚−𝛼𝜕𝑡^𝑚∫^𝑡₀𝑢𝑥,𝑦𝑑𝑦𝑡−𝑦^{𝛼−𝑚+1},𝑚−1<𝛼<𝑚∈ℤ⁺,𝜕^𝑚−1𝑢(𝑥,𝑡)𝜕𝑡^𝑚−1,𝛼=𝑚−1. (7)

Lemma 4. Suppose 𝑚−1<𝛼<𝑚∈ℤ⁺, for any 𝑡, if 𝑢(𝑥,𝑡)∈𝐶^𝑚−1[0,𝑇] and 𝜕𝑢^𝑚−1(𝑥,𝑡)/𝜕𝑡^𝑚−1 is integrable on [0,𝑇] with respect to 𝑡, and then _𝐶𝐷^𝛼_0,𝑡𝑢(𝑥,𝑡)=_𝑅𝐿𝐷^𝛼_0,𝑡𝑢()−𝑥,𝑡_{𝑚−1𝑗=0}𝜕^𝑗𝑢(𝑥,0)𝜕𝑡^𝑗𝑡^𝑗−𝛼Γ.𝑗+1−𝛼 (8)

Lemma 5 (see [31]). The eigenvalues of the tridiagonal matrix 𝑋 with order 𝑛−1𝑋=𝑏𝑎𝑐𝑏𝑎⋱⋱⋱𝑐𝑏𝑎𝑐𝑏 (9) are given by 𝜁_𝜄=𝑏+2𝑎𝑐𝑎cos𝜋𝜄𝑛,𝜄=1,2,…,𝑛−1. (10)

Let 𝑥_𝑙=𝑙ℎ(𝑙=0,1,…,𝑀) and 𝑡_𝑘=𝑘𝜏(𝑘=0,1,…,𝑁), where ℎ=𝐿/𝑀 and 𝜏=𝑇/𝑁 are the uniform spatial and temporal mesh sizes, respectively, and 𝑀, 𝑁 are two positive integers.

We consider a uniform mesh Ω with nodal points 𝑥_𝑙 on [0,𝐿], for each segment [𝑥_𝑙,𝑥_𝑙+1](𝑙=0,1,…,𝑀−1) and point 𝑡_𝑘(𝑘=0,1,…,𝑁) we define the following hyperbolic-trigonometric spline function 𝒮_𝑙(𝑥,𝑡_𝑘) which has different form with the mixed spline function in [27]: 𝒮_𝑙𝑥,𝑡_𝑘=𝑎_𝑙sinh𝜗𝑥−𝑥_𝑙+𝑏_𝑙cosh𝜗𝑥−𝑥_𝑙+𝑐_𝑙sin𝜗𝑥−𝑥_𝑙+𝑑_𝑙,𝑥∈𝑥_𝑙,𝑥_𝑙+1,𝑙=0,1,…,𝑀−1,𝑘=0,1,…,𝑁, (11)

where 𝑎_𝑙,𝑏_𝑙,𝑐_𝑙, and 𝑑_𝑙 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 𝒮_𝑙𝑥_𝑙,𝑡_𝑘=𝑢𝑥_𝑙,𝑡_𝑘,𝒮_𝑙𝑥_𝑙+1,𝑡_𝑘=𝑢𝑥_𝑙+1,𝑡_𝑘,𝒮_𝑙^′′𝑥_𝑙,𝑡_𝑘=Q𝑥_𝑙,𝑡_𝑘,𝒮_𝑙^′′𝑥_𝑙+1,𝑡_𝑘=𝑄𝑥_𝑙+1,𝑡_𝑘. (12)

Through the straightforward calculations, we obtain the following expressions: 𝑎_𝑙=(1−2cosh(𝜗ℎ))𝑄𝑥_𝑙,𝑡_𝑘+𝑄𝑥_𝑙+1,𝑡_𝑘2𝜗²()+𝑢𝑥sinh𝜗ℎ_𝑙+1,𝑡_𝑘−𝑢𝑥_𝑙,𝑡_𝑘(),𝑏2sinh𝜗ℎ_𝑙=𝑄𝑥_𝑙,𝑡_𝑘𝜗²,𝑐_𝑙=𝑄𝑥_𝑙,𝑡_𝑘−𝑄𝑥_𝑙+1,𝑡_𝑘2𝜗²()+𝑢𝑥sin𝜗ℎ_𝑙+1,𝑡_𝑘−𝑢𝑥_𝑙,𝑡_𝑘(),𝑑2sin𝜗ℎ_𝑙=𝑢𝑥_𝑙,𝑡_𝑘−𝑄𝑥_𝑙,𝑡_𝑘𝜗², (13)

in which 𝑙=0,1,…,𝑀−1,𝑘=0,1,…,𝑁.

From (11), one has 𝒮_𝑙−1𝑥,𝑡_𝑘=𝑎_𝑙−1sinh𝜗𝑥−𝑥_𝑙−1+𝑏_𝑙−1cosh𝜗𝑥−𝑥_𝑙−1+𝑐_𝑙−1sin𝜗𝑥−𝑥_𝑙−1+𝑑_𝑙−1,𝑥∈𝑥_𝑙−1,𝑥_𝑙,𝑙=1,…,𝑀,𝑘=0,1,…,𝑁. (14)

In view of the continuity of the first-order derivative at the common nodes (𝑥_𝑙,𝑢(𝑥_𝑙,𝑡_𝑘)), that is, 𝒮^′_𝑙−1(𝑥_𝑙,𝑡_𝑘)=𝒮^′_𝑙(𝑥_𝑙,𝑡_𝑘), we get the following relation: 𝜆₁𝑄𝑥_𝑙−1,𝑡_𝑘+𝜆₂𝑄𝑥_𝑙,𝑡_𝑘+𝜆₃𝑄𝑥_𝑙+1,𝑡_𝑘=1ℎ²𝜇₁𝑢𝑥_𝑖−1,𝑡_𝑘+𝜇₂𝑢𝑥_𝑙,𝑡_𝑘+𝜇₃𝑢𝑥_𝑙+1,𝑡_𝑘,𝑙=1,…,𝑀−1,𝑘=0,1,…,𝑁, (15) where 𝜆₁=1𝜃²cosh(𝜃)(2cosh(𝜃)−1)(𝜃)−(𝜃)sinhcos(𝜃)(𝜃),𝜆sin−2sinh₂=1𝜃²1−3cosh(𝜃)(𝜃)+(𝜃)sinh1+cos(𝜃),𝜆sin₃=1𝜃²1(𝜃)−1sinh(𝜃),𝜇sin₁(𝜃)=−cosh(𝜃)+(𝜃)sinhcos(𝜃),𝜇sin₂=1+cosh(𝜃)(𝜃)+(𝜃)sinh1+cos(𝜃),𝜇sin₃1=−(𝜃)+1sinh(𝜃),sin (16)

and 𝜃=𝜗ℎ.

Obviously, from the expressions of 𝜇₁,𝜇₂, and 𝜇₃, we find that 𝜇₁+𝜇₂+𝜇₃=0.

When 𝜗→0, then 𝜃→0; the spline function given by (15) reduces to the ordinary cubic spline function: ℎ²𝑄𝑥_𝑙−1,𝑡_𝑘+4𝑄𝑥_𝑙,𝑡_𝑘+𝑄𝑥_𝑙+1,𝑡_𝑘=6𝑢𝑥_𝑙−1,𝑡_𝑘−2𝑢𝑥_𝑙,𝑡_𝑘+𝑢𝑥_𝑙+1,𝑡_𝑘,𝑙=1,…,𝑀−1,𝑘=0,1,…,𝑁. (17)

If we use 𝜕²𝑢(𝑥_𝑙−1,𝑡_𝑘)/𝜕𝑥²,𝜕²𝑢(𝑥_𝑙,𝑡_𝑘)/𝜕𝑥²; and 𝜕²𝑢(𝑥_𝑙+1,𝑡_𝑘)/𝜕𝑥² to take the place of 𝑄(𝑥_𝑙−1,𝑡_𝑘),𝑄(𝑥_𝑙,𝑡_𝑘); and 𝑄(𝑥_𝑙+1,𝑡_𝑘) in (15), respectively, then one gets 𝜆₁𝜕²𝑢𝑥_𝑙−1,𝑡_𝑘𝜕𝑥²+𝜆₂𝜕²𝑢𝑥_𝑙,𝑡_𝑘𝜕𝑥²+𝜆₃𝜕²𝑢𝑥_𝑙+1,𝑡_𝑘𝜕𝑥²=1ℎ²𝜇₁𝑢𝑥_𝑙−1,𝑡_𝑘+𝜇₂𝑢𝑥_𝑙,𝑡_𝑘+𝜇₃𝑢𝑥_𝑙+1,𝑡_𝑘+𝒪ℎ^𝑝, (18)

where 𝑝>0.

Using Taylor’s series for (18) at (𝑥_𝑙,𝑡_𝑘), we get the truncation error as follows: 𝐸^𝑘_𝑖=1ℎ𝜇₁−𝜇₃𝜕+1𝜕𝑥22𝜆₁+𝜆₂+𝜆₃−𝜇₁+𝜇₃𝜕²𝜕𝑥²+ℎ66𝜆₃−𝜆₁−𝜇₃−𝜇₁𝜕³𝜕𝑥³+ℎ²𝜆2412₁+𝜆₃−𝜇₁+𝜇₃𝜕⁴𝜕𝑥⁴+ℎ³𝜆12020₃−𝜆₁−𝜇₃−𝜇₁𝜕⁵𝜕𝑥⁵+ℎ⁴𝜆72030₁+𝜆₃−𝜇₁+𝜇₃𝜕⁶𝜕𝑥⁶×𝑢𝑥_𝑙,𝑡_𝑘+𝒪ℎ^𝑝. (19)

From the previous expression, we find that we can obtain different difference schemes by choosing different parameters 𝜆₁,𝜆₂,𝜆₃,𝜇₁, and 𝜇₂. In this paper, we take 𝜆₁=𝜆₃,𝜆₂=10𝜆₁, and 𝜇₁=𝜇₃=12𝜆₁, and then 𝑝=4. Under this assumption, (18) can be rewritten as 𝜕²𝑢𝑥_𝑙,𝑡_𝑘𝜕𝑥²=1ℎ²12𝛿²_𝑥𝑢𝑥_𝑙,𝑡_𝑘12+𝛿²_𝑥+𝒪ℎ⁴, (20) which is the same as in [32], where 𝛿²_𝑥 is the second-order central difference operator with respect to 𝑥.

3. 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).

3.1. Numerical Algorithm I

Based on Lemma 4, (1) can be written as _RL𝐷^𝛼_0,𝑡𝑢(𝑥,𝑡)=𝐾_𝛼𝜕²𝑢(𝑥,𝑡)𝜕𝑥²−𝐶_𝛼𝑢()()+𝑥,𝑡+𝑓𝑥,𝑡_1𝑗=0𝜕^𝑗𝑢(𝑥,0)𝜕𝑡^𝑗𝑡^𝑗−𝛼Γ.𝑗+1−𝛼 (21) Applying the initial value conditions (2) yields _RL𝐷^𝛼_0,𝑡𝑢(𝑥,𝑡)=𝐾_𝛼𝜕²𝑢(𝑥,𝑡)𝜕𝑥²−𝐶_𝛼𝑢(𝑥,𝑡)+𝑓(𝑥,𝑡). (22)

The Riemann-Liouville derivative (with homogeneous initial value conditions) can be approximated by the following scheme: _RL𝐷^𝛼_0,𝑡𝑢()=1𝑥,𝑡𝜏^{𝛼[𝑡/𝜏]}_𝑗=0𝜛^(𝛼)_𝑞,𝑗𝑢𝑥,𝑡−𝑗𝜏+𝒪𝜏^𝑞. (23)

There exist various choices of the generating functions which can lead to different approximation order 𝑞.

Let 𝜛_𝑞(𝑧,𝛼) be the generating function with coefficients 𝜛^(𝛼)_𝑞,𝑗; that is, 𝜛_𝑞(𝑧,𝛼)=^∞_𝑗=0𝜛^(𝛼)_𝑞,𝑗𝑧^𝑗. (24)

If we take the generating function in the form 𝜛₁(𝑧,𝛼)=(1−𝑧)^𝛼, (25) that is, (1−𝑧)^𝛼=^∞_𝑗=0𝜛^(𝛼)_1,𝑗𝑧^𝑗, (26) then we can obtain the following first-order scheme: _RL𝐷^𝛼_0,𝑡𝑢()=1𝑥,𝑡𝜏^{𝛼[𝑡/𝜏]}_𝑗=0𝜛^(𝛼)_1,𝑗𝑢𝑥,𝑡−𝑗𝜏+𝒪(𝜏), (27) in which 𝜛^(𝛼)_1,0𝜛=1,_(𝛼)1,𝑗=(−1)^𝑗𝛼𝑗=(−1)^𝑗Γ(1+𝛼)Γ𝑗+1Γ1+𝛼−𝑗,𝑗≥1, (28) and these coefficients have the following recursive relations:𝜛^(𝛼)_1,0=1,𝜛^(𝛼)_1,𝑗=1−1+𝛼𝑗𝜛^(𝛼)_1,𝑗−1,𝑗≥1. (29)

If we take the generating function as 𝜛₂()=3𝑧,𝛼21−2𝑧+2𝑧^2𝛼, (30) that is, 321−2𝑧+2𝑧^2𝛼=^∞_𝑗=0𝜛^(𝛼)_2,𝑗𝑧^𝑗, (31) which leads to order 𝑞=2 [30], then_RL𝐷^𝛼_0,𝑡𝑢()=1𝑥,𝑡𝜏^{𝛼[𝑡/𝜏]}_𝑗=0𝜛^(𝛼)_2,𝑗𝑢𝑥,𝑡−𝑗𝜏+𝒪𝜏². (32)

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 321−2𝑧+2𝑧^2𝛼=32^𝛼(1−𝑧)^𝛼11−3𝑧^𝛼 (33)

and (26), one gets 321−2𝑧+2𝑧^2𝛼=32^𝛼(1−𝑧)^𝛼11−3𝑧^𝛼=32^𝛼∞_𝑗=0(−1)^𝑗𝛼𝑗𝑧^𝑗∞_𝑗=0−13^𝑗𝛼𝑗𝑧^𝑗=32^𝛼∞∞_𝑚=0𝑛=0(−1)^𝑚−13^𝑛𝛼𝑚𝛼𝑛𝑧^𝑚+𝑛=32^𝛼∞_𝑗=0^𝑗_𝑚=0(−1)^𝑚−13^𝑗−𝑚𝛼𝑚𝛼𝑧𝑗−𝑚^𝑗=^∞_𝑗=032^𝛼𝑗_𝑚=0(−1)^𝑗13^𝑚𝛼𝑚𝛼𝑧𝑗−𝑚^𝑗=^∞_𝑗=032^𝛼𝑗_𝑚=013^𝑚𝜛^(𝛼)_1,𝑚𝜛^(𝛼)_{1,𝑗−𝑚}𝑧^𝑗. (34)

Comparing the previous equation with (31), one obtains 𝜛^(𝛼)_2,𝑗=32^𝛼𝑗_𝑚=013^𝑚𝜛^(𝛼)_1,𝑚𝜛^(𝛼)_{1,𝑗−𝑚},=(−1)^𝑗32^𝛼×^𝑗_𝑚=013^𝑚×Γ²(1+𝛼)()Γ,𝑚!𝑗−𝑚!Γ1+𝛼−𝑚1+𝛼+𝑚−𝑗 (35)

which has the following recursive relation: 𝜛^(𝛼)_2,0=32^𝛼,𝜛^(𝛼)_2,𝑗+1=32^𝛼𝑗+1_𝑚=013^𝑚𝜛^(𝛼)_1,𝑚𝜛^(𝛼)_{1,𝑗+1−𝑚}=32^𝛼𝑗_𝑚=013^𝑚𝜛^(𝛼)_1,𝑚𝜛^(𝛼)_{1,𝑗+1−𝑚}+32^𝛼13^𝑗+1𝜛^(𝛼)_1,𝑗+1=32^𝛼𝑗_𝑚=013^𝑚1−1+𝛼𝜛𝑗+1−𝑚^(𝛼)1,𝑚𝜛^(𝛼)_{1,𝑗−𝑚}+32^𝛼13^𝑗+1𝜛^(𝛼)_1,𝑗+1=𝜛^(𝛼)_2,𝑗−()31+𝛼2^𝛼13^𝑗𝑗_𝑚=03^𝑚𝜛𝑚+1^(𝛼)1,𝑚𝜛^(𝛼)_{1,𝑗−𝑚}+32^𝛼13^𝑗+11−1+𝛼𝜛𝑗+1^(𝛼)1,𝑗,𝑗≥0. (36) In summary, we have the following lemma.

Lemma 6. The coefficients 𝜛^(𝛼)_2,𝑗(𝑗=0,1,…)(1<𝛼<2) satisfy (1)𝜛^(𝛼)_2,0=32^𝛼𝜛>0,_(𝛼)2,14=−3𝛼32^𝛼𝜛<0,_(𝛼)2,2=𝛼(8𝛼−5)932^𝛼𝜛>0,_(𝛼)2,3=4𝛼(1−𝛼)(8𝛼−7)3812^𝛼𝜛<0,_(𝛼)2,4=𝛼(𝛼−1)64𝛼²−176𝛼+12334862^𝛼𝜛>0,_(𝛼)2,5=2𝛼(𝛼−1)(2−𝛼)64𝛼²−208𝛼+183336452^𝛼⋮(2)>0,_∞𝑗=0𝜛^(𝛼)_2,𝑗=0,^∞_𝑗=0(−1)^𝑗𝜛^(𝛼)_2,𝑗=4^𝛼. (37)

Proof. (1) In view of (36), we easily obtain the expressions of the coefficients 𝜛^(𝛼)_2,𝑗(𝑗=0,1,…).
(2) In (31), if we take 𝑧≡1, which leads to ∑^∞_𝑗=0𝜛^(𝛼)_2,𝑗=0. Similarly, choosing 𝑧≡−1 and substitution it into (31) yield ∑^∞_𝑗=0(−1)^𝑗𝜛^(𝛼)_2,𝑗=4^𝛼. This completes the proof.

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: 12𝜅₁−𝜅₂−𝜛^(𝛼)_2,0𝑢^𝑘_𝑙−1−212𝜅₁+5𝜅₂+5𝜛^(𝛼)_2,0𝑢^𝑘_𝑙+12𝜅₁−𝜅₂−𝜛^(𝛼)_2,0𝑢^𝑘_𝑙+1=𝜛^(𝛼)_2,1𝑢^𝑘−1_𝑙−1+10𝑢_𝑙^𝑘−1+𝑢^𝑘−1_𝑙+1+^𝑘_𝑗=2𝜛^(𝛼)_2,𝑗𝑢^𝑘−𝑗_𝑙−1+10𝑢_𝑙^𝑘−𝑗+𝑢^𝑘−𝑗_𝑙+1−𝜏^𝛼𝑓^𝑘_𝑙−1+10𝑓^𝑘_𝑙+𝑓^𝑘_𝑙+1, (38) for 𝑙=1,…,𝑀−1,𝑘=1,…,𝑁. The initial and boundary value conditions can be discretized by 𝑢⁰_𝑙1=0,𝑙=0,1,…,𝑀,𝜏𝛿_𝑡𝑢⁰_𝑙1+(1/2)𝛿_𝑡=1𝜏𝑢¹_𝑙−𝑢⁰_𝑙1+(1/2)𝛿_𝑡𝑢=0,𝑙=0,1,…,𝑀,_𝑘0=𝜑₁(𝑘𝜏),𝑢^𝑘_𝑀=𝜑₂(𝑘𝜏),𝑘=0,1,…,𝑁, (39)

where 𝑓^𝑘_𝑙=𝑓(𝑥_𝑙,𝑡_𝑘),𝜅₁=𝐾_𝛼𝜏^𝛼/ℎ², and 𝜅₂=𝐶_𝛼𝜏^𝛼.

Next, we study the solvability and the local truncation error.

For convenience, denote 𝐔⁰=0,…,0_𝑀−1^𝑇,𝐔^𝑘=𝑢^𝑘₁,𝑢^𝑘₂,…,𝑢^𝑘_𝑀−1^𝑇𝐅,𝑘=1,2,…,𝑁,_𝑘=𝑓^𝑘₁,𝑓^𝑘₂,…,𝑓^𝑘_𝑀−1^𝑇Λ,𝑘=0,1,…,𝑁,₁=,Λ1011101⋱⋱⋱1101110₂=.−211−21⋱⋱⋱1−211−2 (40)

Then we can give the following matrix form of the difference scheme (38): 12𝜅₁Λ₂−𝜅₂Λ₁−𝐴₀𝐔^𝑘=𝐴₁𝐔^𝑘−1+^𝑘_𝑗=2𝐴_𝑗𝐔^𝑘−𝑗−𝜏^𝛼Λ₁𝐅^𝑘+𝐘, (41) where 𝐴_𝑗=𝜛^(𝛼)_2,𝑗Λ₁Ψ,𝑗=0,1,…,𝑘,𝐘=₁,0,…,0_𝑀−3,Ψ_𝑀−1^𝑇,Ψ₁=−12𝜅₁−𝜅₂−𝜛^(𝛼)_2,0𝜑₁𝑡_𝑘+^𝑘_𝑗=1𝜛^(𝛼)_2,𝑗𝜑₁𝑡_𝑘−𝑗−𝜏^𝛼𝑓𝑥₀,𝑡_𝑘,Ψ_𝑀−1=−12𝜅₁−𝜅₂−𝜛^(𝛼)_2,0𝜑₂𝑡_𝑘+^𝑘_𝑗=1𝜛^(𝛼)_2,𝑗𝜑₂𝑡_𝑘−𝑗−𝜏^𝛼𝑓𝑥_𝑀,𝑡_𝑘. (42)

Theorem 7. The difference equation (38) is uniquely solvable.

Proof. From Lemma 5, we know that the eigenvalues of the tridiagonal matrix (12𝜅₁Λ₂−𝜅₂Λ₁−𝐴₀) are 𝜁̂_𝑙=−48𝜅₁sin²𝜋𝑙𝜅2𝑀−4₂+𝜛^(𝛼)_2,03−sin²𝜋𝑙,2𝑀𝑙=1,2,…,𝑀−1. (43) Due to 𝜅₁, 𝜅₂, and 𝜛^(𝛼)_2,0>0, hence, det12𝜅₁Λ₂−𝜅₂Λ₁−𝐴₀=^𝑀−1_𝑙=1𝜁̂_𝑙≠0. (44) That is, the matrix (12𝜅₁Λ₂−𝜅₂Λ₁−𝐴₀) is nonsingular. Thus (41) is uniquely solvable, and so is (38). The proof is finished.

Theorem 8. The local truncation error of the difference scheme (38) is 𝒪(𝜏²+ℎ⁴).

Proof. We now define the local truncation error as 𝑅^𝑘_𝑙=1𝜏^𝛼𝑘_𝑗=0𝜛^(𝛼)_2,𝑗𝑢𝑥_𝑙,𝑡_𝑘−𝑗−𝐾_𝛼ℎ²12𝛿²_𝑥𝑢𝑥_𝑙,𝑡_𝑘12+𝛿²_𝑥−𝐶_𝛼𝑢𝑥_𝑙,𝑡_𝑘−𝑓𝑥_𝑙,𝑡_𝑘. (45) Applying (20), (22), and (32), we have 𝑅^𝑘_𝑙=1𝜏^𝛼𝑘_𝑗=0𝜛^(𝛼)_2,𝑗𝑢𝑥_𝑙,𝑡_𝑘−𝑗−_RL𝐷^𝛼_0,𝑡𝑢𝑥_𝑙,𝑡_𝑘+𝐾_𝛼𝜕²𝑢𝑥_𝑙,𝑡_𝑘𝜕𝑥²−1ℎ²12𝛿²_𝑥𝑢𝑥_𝑙,𝑡_𝑘12+𝛿²_𝑥=𝒪𝜏²+ℎ⁴. (46)
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,…,𝑁,_𝑘=𝜀^𝑘₁,𝜀^𝑘₂,…,𝜀^𝑘_𝑀−1^𝑇. (47)

So, we can easily obtain the following round-off error equation: 12𝜅₁−𝜅₂−𝜛^(𝛼)_2,0𝜀^𝑘_𝑙−1−212𝜅₁+5𝜅₂+5𝜛^(𝛼)_2,0𝜀^𝑘_𝑙+12𝜅₁−𝜅₂−𝜛^(𝛼)_2,0𝜀^𝑘_𝑙+1=𝜛^(𝛼)_2,1𝜀^𝑘−1_𝑙−1+10𝜀_𝑙^𝑘−1+𝜀^𝑘−1_𝑙−1+^𝑘_𝑗=2𝜛^(2−𝛼)_2,𝑗𝜀^𝑘−𝑗_𝑙−1+10𝜀_𝑙^𝑘−𝑗+𝜀^𝑘−𝑗_𝑙+1. (48)

We suppose that the solution of (48) has the following form: 𝜀^𝑘_𝑙=𝜉_𝑘exp𝑖𝜙𝑙ℎ, (49) where 𝜙 is any of the real spatial wave numbers supported by the lattice. Substituting (49) into (48) gives 12𝜅₁sin²(1/2)𝜙ℎ3−sin²(1/2)𝜙ℎ+𝜅₂+𝜛^(𝛼)_2,0𝜉_𝑘=−^𝑘_𝑗=1𝜛^(𝛼)_2,𝑗𝜉_𝑘−𝑗. (50) If we write 𝜉_𝑘+1=𝜂𝜉_𝑘 (51) and assume that 𝜂≡𝜂(𝜙) is independent of time, then we obtain the following expression of the amplification factor: 12𝜅₁sin²(1/2)𝜙ℎ3−sin²(1/2)𝜙ℎ+𝜅₂=−^𝑘_𝑗=0𝜛^(𝛼)_2,𝑗𝜂^−𝑗. (52)

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 12𝜅₁sin²(1/2)𝜙ℎ3−sin²(1/2)𝜙ℎ+𝜅₂≥−^𝑘_𝑗=0(−1)^𝑗𝜛^(𝛼)_2,𝑗=𝑆_𝛼,𝑘. (53) Although 𝑆_𝛼,𝑘 depends on 𝑘, we can estimate it by its limit value 𝑆_𝛼=lim_𝑘→∞𝑆_𝛼,𝑘=−^∞_𝑗=0(−1)^𝑗𝜛^(𝛼)_2,𝑗=−4^𝛼. (54) In this case, the stability condition becomes 12𝜅₁sin²(1/2)𝜙ℎ3−sin²(1/2)𝜙ℎ+𝜅₂≥𝑆_𝛼. (55)

Note that 𝑆_𝛼 is always negative, and then (55) holds for all 𝜅₁>0 and 𝜅₂>0. That is to say, Numerical Algorithm I (see (38)) is unconditionally stable.

3.2. Numerical Algorithm II

Differentiating (1) with order 2−𝛼 in the sense of Riemann-Liouville yields _RL𝐷^2−𝛼_0,𝑡𝐶𝐷^𝛼_0,𝑡𝑢()=𝑥,𝑡_RL𝐷^2−𝛼_0,𝑡𝐾_𝛼𝜕²𝑢(𝑥,𝑡)𝜕𝑥²−𝐶_𝛼𝑢(𝑥,𝑡)+𝑓(𝑥,𝑡). (56) Due to _𝐶𝐷^𝛼_0,𝑡𝑢(𝑥,𝑡)=𝐷^{−(2−𝛼)}_0,𝑡𝜕²𝑢(𝑥,𝑡)𝜕𝑡², (57) it immediately follows that _RL𝐷^2−𝛼_0,𝑡𝐶𝐷^𝛼_0,𝑡𝑢(𝑥,𝑡)=_RL𝐷^2−𝛼_0,𝑡𝐷^{−(2−𝛼)}_0,𝑡𝜕²𝑢(𝑥,𝑡)𝜕𝑡²=𝜕²𝑢(𝑥,𝑡)𝜕𝑡². (58)

Combing (56) with (58) gives 𝜕²𝑢(𝑥,𝑡)𝜕𝑡²=_RL𝐷^2−𝛼_0,𝑡𝐾_𝛼𝜕²𝑢(𝑥,𝑡)𝜕𝑥²−𝐶_𝛼𝑢(𝑥,𝑡)+𝑔(𝑥,𝑡), (59)

where 𝑔(𝑥,𝑡)=_RL𝐷^2−𝛼_0,𝑡𝑓(𝑥,𝑡).

Next, for the second-order time derivative 𝜕²𝑢(𝑥_𝑙,𝑡_𝑘)/𝜕𝑡², we approximate it by the usual three-point central difference formula 𝜕²𝑢𝑥_𝑙,𝑡_𝑘𝜕𝑡²=𝑢𝑥_𝑙,𝑡_𝑘−1−2𝑢𝑥_𝑙,𝑡_𝑘+𝑢𝑥_𝑙,𝑡_𝑘+1𝜏²+𝒪𝜏². (60)

Let 𝑢^𝑘_𝑙 be the numerical approximation of 𝑢(𝑥_𝑙,𝑡_𝑘). Substituting the expansions (20), (32), and (60) in (59) and removing the higher-order terms, one obtains Numerical Algorithm II as follows: 𝑢^𝑘+1_𝑙−1+10𝑢_𝑙^𝑘+1+𝑢^𝑘+1_𝑙+1=2+𝜛^(2−𝛼)_2,012𝜅₁−𝜅₂𝑢^𝑘_𝑙−1+210−𝜛^(2−𝛼)_2,012𝜅₁+5𝜅₂𝑢^𝑘_𝑙+2+𝜛^(2−𝛼)_2,012𝜅₁−𝜅₂𝑢^𝑘_𝑙+1−𝑢^𝑘−1_𝑙−1−10𝑢_𝑙^𝑘−1−𝑢^𝑘−1_𝑙+1+12𝜅₁^𝑘_𝑗=1𝜛^(2−𝛼)_2,𝑗𝑢^𝑘−𝑗_𝑙−1−2𝑢_𝑙^𝑘−𝑗+𝑢^𝑘−𝑗_𝑙+1−𝜅₂^𝑘_𝑗=1𝜛^(2−𝛼)_2,𝑗𝑢^𝑘−𝑗_𝑙−1+10𝑢_𝑙^𝑘−𝑗+𝑢^𝑘−𝑗_𝑙+1+𝜏²𝑔^𝑘_𝑙−1+10𝑔^𝑘_𝑙+𝑔^𝑘_𝑙+1, (61) where 𝑙=1,…,𝑀−1,𝑘=1,…,𝑁−1. The initial and boundary value conditions can be discretized by 𝑢⁰_𝑙1=0,𝑙=0,1,…,𝑀,𝜏𝛿_𝑡𝑢⁰_𝑙1+(1/2)𝛿_𝑡=1𝜏𝑢¹_𝑙−𝑢⁰_𝑙1+(1/2)𝛿_𝑡𝑢=0,𝑙=0,1,…,𝑀,_𝑘0=𝜑₁(𝑘𝜏),𝑢^𝑘_𝑀=𝜑₂(𝑘𝜏),𝑘=0,1,…,𝑁, (62)

where 𝑓^𝑘_𝑙=𝑓(𝑥_𝑙,𝑡_𝑘),𝜅₁=𝐾_𝛼𝜏^𝛼/ℎ², and 𝜅₂=𝐶_𝛼𝜏^𝛼.

Similar to Numerical Algorithm I, we discuss the solvability and local truncation error of Numerical Algorithm II.

Denote 𝐆^𝑘=𝑔^𝑘₁,𝑔^𝑘₂,…,𝑔^𝑘_𝑀−1^𝑇,𝑘=0,1,…,𝑁. (63) Then the matrix form of Numerical Algorithm II (61) is written as Λ₁𝐔^𝑘+1=2Λ₁+𝐴₀𝐔^𝑘−Λ₁𝐔^𝑘−1+^𝑘_𝑗=1𝐴_𝑗𝐔^𝑘−𝑗+𝜏²Λ₁𝐆^𝑘+𝐘,𝑘=1,…,𝑁−1, (64) where 𝐴_𝑗=𝜛^(2−𝛼)_2,𝑗12𝜅₁Λ₂−𝜅₂Λ₁Ψ,𝑗=0,1,…,𝑘,𝐘=₁,0,…,0_𝑀−3,Ψ_𝑀−1^𝑇,Ψ₁=−𝜑₁𝑡_𝑘+1+2𝜑₁𝑡_𝑘−𝜑₁𝑡_𝑘−1+12𝜅₁−𝜅₂^𝑘_𝑗=0𝜛^(2−𝛼)_2,𝑗𝜑₁𝑡_𝑘−𝑗+𝜏²𝑔𝑥₀,𝑡_𝑘,Ψ_𝑀−1=−𝜑₂𝑡_𝑘+1+2𝜑₂𝑡_𝑘−𝜑₂𝑡_𝑘−1+12𝜅₁−𝜅₂^𝑘_𝑗=0𝜛^(2−𝛼)_2,𝑗𝜑₂𝑡_𝑘−𝑗+𝜏²𝑔𝑥_𝑀,𝑡_𝑘. (65)

Theorem 9. The difference equation (61) is uniquely solvable.

Proof. Obviously, the eigenvalues of the tridiagonal matrix Λ₁ are given by 𝜁̃_𝑙=12−4sin²𝜋𝑙2𝑀,𝑙=1,2,…,𝑀−1. (66)
Hence, detΛ₁=^𝑀−1_𝑙=1𝜁̃_𝑙>0. (67)
Therefore, the solution of (61) exists and is uniquely solvable.

In order to get the local truncation error analysis, we need the following lemma.

Lemma 10. Consider 1𝜏^𝑘2−𝛼_𝑗=0𝜛^(2−𝛼)_2,𝑗=(𝑘𝜏)^𝛼−2Γ()𝜏𝛼−1+𝒪²,𝑘=1,…,𝑁. (68)

Proof. Let 𝑢(𝑥,𝑡)≡1 and then from Definition 3 we have _RL𝐷^2−𝛼_0,𝑡𝑢(𝑥,𝑡)_{𝑡=𝑡𝑘}=1Γ()𝜕𝛼−1𝜕𝑡^𝑡₀𝑢(𝑥,𝑠)(𝑡−𝑠)^2−𝛼𝑑𝑠_{𝑡=𝑡𝑘}=(𝑘𝜏)^𝛼−2Γ().𝛼−1 (69) It follows from (32) that _RL𝐷^2−𝛼_0,𝑡𝑢(𝑥,𝑡)_{𝑡=𝑡𝑘}=1𝜏^𝑘2−𝛼_𝑗=0𝜛^(2−𝛼)_2,𝑗𝑢𝑥,𝑡_𝑘−𝑗+𝒪𝜏²=1𝜏^𝑘2−𝛼_𝑗=0𝜛^(2−𝛼)_2,𝑗+𝒪𝜏². (70) Combining (69) and (77) leads to 1𝜏^𝑘2−𝛼_𝑗=0𝜛^(2−𝛼)_2,𝑗=(𝑘𝜏)^𝛼−2Γ()𝜏𝛼−1+𝒪²,𝑘=1,…,𝑁. (71) This ends the proof.

Theorem 11. The local truncation error of difference scheme (61) is 𝒪(𝜏²+ℎ⁴).

Proof. We now define 𝑅^𝑘_𝑙=𝑢𝑥_𝑙,𝑡_𝑘−1−2𝑢𝑥_𝑙,𝑡_𝑘+𝑢𝑥_𝑙,𝑡_𝑘+1𝜏²−𝐾_𝛼𝜏^𝑘2−𝛼_𝑗=0𝜛^(2−𝛼)_2,𝑗1ℎ²12𝛿²_𝑥𝑢𝑥_𝑙,𝑡_𝑘−𝑗12+𝛿²_𝑥−𝐶_𝛼𝜏^𝑘2−𝛼_𝑗=0𝜛^(2−𝛼)_2,𝑗𝑢𝑥_𝑙,𝑡_𝑘−𝑗−𝑔𝑥_𝑙,𝑡_𝑘. (72)
By using (20), (32), (59), and (60), one gets 𝑅^𝑘_𝑙=𝑢𝑥_𝑙,𝑡_𝑘−1−2𝑢𝑥_𝑙,𝑡_𝑘+𝑢𝑥_𝑙,𝑡_𝑘+1𝜏²−𝐾_𝛼𝜏^𝑘2−𝛼_𝑗=0𝜛^(2−𝛼)_2,𝑗1ℎ²12𝛿²_𝑥𝑢𝑥_𝑙,𝑡_𝑘−𝑙12+𝛿²_𝑥+𝐶_𝛼𝜏^𝑘2−𝛼_𝑗=0𝜛^(2−𝛼)_2,𝑗𝑢𝑥_𝑙,𝑡_𝑘−𝑗−𝜕²𝑢𝑥_𝑙,𝑡_𝑘𝜕𝑡²+𝐾_𝛼RL𝐷^2−𝛼_0,𝑡𝜕²𝑢𝑥_𝑙,𝑡_𝑘𝜕𝑥²−𝐶_𝛼RL𝐷^2−𝛼_0,𝑡𝑢𝑥_𝑙,𝑡_𝑘=𝑢𝑥_𝑙,𝑡_𝑘−1−2𝑢𝑥_𝑙,𝑡_𝑘+𝑢𝑥_𝑙,𝑡_𝑘+1𝜏²−𝜕²𝑢𝑥_𝑙,𝑡_𝑘𝜕𝑡²+𝐾_𝛼𝜏^𝑘2−𝛼_𝑗=0𝜛^(2−𝛼)_2,𝑗𝜕²𝑢𝑥_𝑙,𝑡_𝑘−𝑙𝜕𝑥²−1ℎ²12𝛿²_𝑥𝑢𝑥_𝑙,𝑡_𝑘−𝑙12+𝛿²_𝑥+𝐾_𝛼RL𝐷^2−𝛼_0,𝑡𝜕²𝑢𝑥_𝑙,𝑡_𝑘𝜕𝑥²−1𝜏^𝑘2−𝛼_𝑗=0𝜛^(2−𝛼)_2,𝑗𝜕²𝑢𝑥_𝑙,𝑡_𝑘−𝑙𝜕𝑥²+𝐶_𝛼1𝜏^𝑘2−𝛼_𝑗=0𝜛^(2−𝛼)_2,𝑗𝑢𝑥_𝑙,𝑡_𝑘−𝑗−_RL𝐷^2−𝛼_0,𝑡𝑢𝑥_𝑙,𝑡_𝑘=𝒪𝜏²+𝒪𝜏²ℎ⁴+(𝑘𝜏)^𝛼−2Γ()𝒪ℎ𝛼−1⁴. (73)
Noticing that 1≤𝑘≤𝑁, then 1𝑇^2−𝛼≤(𝑘𝜏)^𝛼−2≤𝑁𝑇^2−𝛼, (74) so we can get 𝑅^𝑘_𝑙=𝒪𝜏²+ℎ⁴. (75) 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: 𝜀^𝑘+1_𝑙−1+10𝜀_𝑙^𝑘+1+𝜀^𝑘+1_𝑙+1=2+𝜛^(2−𝛼)_2,012𝜅₁−𝜅₂𝜀^𝑘_𝑙−1+210−𝜛^(2−𝛼)_2,012𝜅₁+5𝜅₂𝜀^𝑘_𝑙+2+𝜛^(2−𝛼)_2,012𝜅₁−𝜅₂𝜀^𝑘_𝑙+1−𝜀^𝑘−1_𝑙−1−10𝜀_𝑙^𝑘−1−𝜀^𝑘−1_𝑙+1+12𝜅₁^𝑘_𝑗=1𝜛^(2−𝛼)_2,𝑗𝜀^𝑘−𝑗_𝑙−1−2𝜀_𝑙^𝑘−𝑗+𝜀^𝑘−𝑗_𝑙+1−𝜅₂^𝑘_𝑗=1𝜛^(2−𝛼)_2,𝑗𝜀^𝑘−𝑗_𝑙−1+10𝜀_𝑙^𝑘−𝑗+𝜀^𝑘−𝑗_𝑙+1,𝑘=1,…,𝑁−1. (76) Substituting (49) into (76) yields 3−sin²12𝜙ℎ𝜉_𝑘+1=23−sin²12𝜙ℎ−12𝜅₁𝜛^(2−𝛼)_2,0sin²12𝜙ℎ−𝜅₂𝜛^(2−𝛼)_2,03−sin²12𝜙ℎ𝜉_𝑘−3−sin²12𝜙ℎ𝜉_𝑘−1−𝜅₂3−sin²12×𝜙ℎ_𝑘𝑗=1𝜛^(2−𝛼)_2,𝑗𝜉_𝑘−𝑗−12𝜅₁sin²12×𝜙ℎ_𝑘𝑗=1𝜛^(2−𝛼)_2,𝑗𝜉_𝑘−𝑗. (77) From (51) and (77), we have 3−sin²12𝜙ℎ𝜂=23−sin²12𝜙ℎ−12𝜅₁𝜛^(2−𝛼)_2,0sin²12𝜙ℎ−𝜅₂𝜛^(2−𝛼)_2,03−sin²12−𝜙ℎ3−sin²12𝜙ℎ𝜂⁻¹−𝜅₂3−sin²12𝜙ℎ^𝑘_𝑗=1𝜛^(2−𝛼)_2,𝑗𝜂^−𝑗−12𝜅₁sin²12𝜙ℎ^𝑘_𝑗=1𝜛^(2−𝛼)_2,𝑗𝜂^−𝑗. (78)

Considering the extreme value 𝜂=−1, we can obtain the following stability condition from (78): sin²12𝜙ℎ≤12−3𝜅₂4+12𝜅₁−𝜅₂∑^𝑘_𝑗=0(−1)^𝑗𝜛^(2−𝛼)_2,𝑗≡𝑆_𝛼,𝑘. (79)

For large enough 𝑘, we can estimate 𝑆_𝛼,𝑘 by its limit value 𝑆_𝛼=lim_𝑘→∞𝑆_𝛼,𝑘=12−3𝜅₂4+12𝜅₁−𝜅₂∑^∞_𝑗=0(−1)^𝑗𝜛^(2−𝛼)_2,𝑗=12−3𝜅₂4+12𝜅₁−𝜅₂4^2−𝛼. (80)

It follows that sin²12𝜙ℎ≤12−3𝜅₂4+12𝜅₁−𝜅₂4^2−𝛼=𝑆_𝛼. (81)

Because the sine function is bounded by 1, from (81) one finds that Numerical Algorithm II is to be stable if 𝑆_𝛼≥1, that is,12−3𝜅₂4+12𝜅₁−𝜅₂4^2−𝛼≥1. (82)

Furthermore, stability condition (82) can be rewritten as 3ℎ²4−𝐶_𝛼𝜏^𝛼ℎ²+4^1−𝛼𝜏^𝛼12𝐾_𝛼−𝐶_𝛼ℎ²≥4, (83)

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.

Moreover, from Theorems 8 and 11, we know that Numerical Algorithm I (see (38)) and Numerical Algorithm II (see (61)) are all consistent with the local truncation errors 𝑅^𝑘_𝑙=𝒪(𝜏²+ℎ⁴) and 𝑅^𝑘_𝑙=𝒪(𝜏²+ℎ⁴) because 𝑅^𝑘_𝑙,𝑅^𝑘_𝑙→0 as 𝜏,ℎ→0. Therefore, according to Lax’s Equivalence Theorem [33], Numerical Algorithm I unconditionally converges at the same order 𝒪(𝜏²+ℎ⁴), but Numerical Algorithm II converges at the same order 𝒪(𝜏²+ℎ⁴) under condition (83).

4. Numerical Example

In this section, we give a numerical example to demonstrate the efficiency of the derived numerical algorithms.

Example 1. Consider the following equation: _𝐶𝐷^𝛼_0,𝑡𝑢()=𝜕𝑥,𝑡²𝑢(𝑥,𝑡)𝜕𝑥²−𝑢(𝑥,𝑡)+𝑓(𝑥,𝑡),0≤𝑡≤1,0<𝑥<2, (84) where 𝑓(𝑥,𝑡)=𝛼³[(6𝑡^3−𝛼/Γ(4−𝛼))+𝑡³(1+(𝜋²/4))]sin((𝜋/2)𝑥).
The analytical solution of this equation is 𝑢(𝑥,𝑡)=(𝛼𝑡)³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.

Table 1 

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

Figure 1 

Dependence of the numerical solution for different values of 𝛼 by Numerical Algorithm I (see (38)) with 𝜏=0.002,ℎ=0.0025.

Figure 2 

Dependence of the numerical solution for different values of 𝛼 by Numerical Algorithm I (see (38)) with 𝜏=0.0001,ℎ=0.05.

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.

Next, we use Numerical Algorithm II (61). In this situation, we take (60), the equivalence equation of (1), where 𝑔(𝑥,𝑡)=6𝑡𝛼³[1+(𝑡^𝛼/Γ(2+𝛼))(1+(𝜋²/4))]sin((𝜋/2)𝑥).

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.

Table 2 

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

Figure 4 

Dependence of the numerical solution for different values of 𝛼 by Numerical Algorithm II (see (61)) with 𝜏=0.002,ℎ=0.0025.

Figure 5 

Dependence of the numerical solution for different values of 𝛼 by Numerical Algorithm II (see (61)) with 𝜏=0.0001,ℎ=0.05.

Figure 6 

The error of the analytical solution and the numerical solution for 𝛼=1.5 by Numerical Algorithm II (see (61)) with 𝜏=0.04,ℎ=0.002.

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.

Figure 7 

Comparison of the analytical solution and numerical solution for 𝛼=1.6 by Numerical Algorithm II (see (61)) with ℎ=0.04,𝜏=0.005, which satisfy the stability (83).

Figure 8 

Comparison of the analytical solution and numerical solution for 𝛼=1.6 by Numerical Algorithm II (see (61)) with ℎ=0.01,𝜏=0.01, which do not satisfy the stability (83).

5. 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 𝒪(𝜏²+ℎ⁴), which are the best results till now. Finally, the numerical results given in this paper show the effectiveness of the derived numerical algorithms.

Acknowledgments

The work was partially supported by the Key Program of Shanghai Municipal Education Commission under Grant no. 12ZZ084 and the grant of “The First-Class Discipline of Universities in Shanghai.”