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

Authors: Hengfei Ding; Changpin Li
Journal: Abstract and Applied Analysis
 (2013)

---

## 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.


---

## Body

## 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:(1)CD0,tαu(x,t)=Kα∂2u(x,t)∂x2-Cαu(x,t)+f(x,t),0≤t≤T,0<x<L.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:
(2)u(x,0)=0,∂u(x,0)∂t=0,0<x<L,(3)u(0,t)=φ1(t),0≤t≤T,(4)u(L,t)=φ2(t),0≤t≤T,
where f(x,t), φ1(x), and φ2(x) are sufficiently smooth functions and Kα>0 and Cα>0 are diffusion and reaction coefficients, respectively. CD0,tαu(x,t) 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 Section2, 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 operatorD0,t-α of order α>0 is defined in the following form:
(5)D0,t-αu(x,t)=1Γ(α)∫0t‍u(x,y)dy(t-y)1-α,
where Γ(·) is the Euler gamma function.Definition 2.
The left Caputo derivative operatorCD0,tα of order α>0 is defined in the following form:
(6)CD0,tαu(x,t)={1Γ(m-α)∫0t‍∂mu(x,y)∂ymdy(t-y)α-m+1,m-1<α<m∈ℤ+,∂mu(x,t)∂tm,α=m.Definition 3.
The left Riemann-Liouville derivative operatorRLD0,tα of order α>0 is defined by
(7)RLD0,tαu(x,t)={1Γ(m-α)(∂∂t)m∫0t‍u(x,y)dy(t-y)α-m+1,m-1<α<m∈ℤ+,∂m-1u(x,t)∂tm-1,α=m-1.Lemma 4.
Supposem-1<α<m∈ℤ+, for any t, if u(x,t)∈Cm-1[0,T] and ∂um-1(x,t)/∂tm-1 is integrable on [0,T] with respect to t, and then
(8)CD0,tαu(x,t)=RLD0,tαu(x,t)-∑j=0m-1‍∂ju(x,0)∂tjtj-αΓ(j+1-α).Lemma 5 (see [31]).
The eigenvalues of the tridiagonal matrixX with order n-1(9)X=(bacba⋱⋱⋱cbacb)
are given by
(10)ζι=b+2acacos(πιn),ι=1,2,…,n-1.Letxl=lh(l=0,1,…,M) and tk=kτ(k=0,1,…,N), where h=L/M and τ=T/N are the uniform spatial and temporal mesh sizes, respectively, and M, N are two positive integers.We consider a uniform meshΩ with nodal points xl on [0,L], for each segment [xl,xl+1](l=0,1,…,M-1) and point tk(k=0,1,…,N) we define the following hyperbolic-trigonometric spline function 𝒮l(x,tk) which has different form with the mixed spline function in [27]:
(11)𝒮l(x,tk)=alsinh(ϑ(x-xl))+blcosh(ϑ(x-xl))+clsin(ϑ(x-xl))+dl,x∈[xl,xl+1],l=0,1,…,M-1,k=0,1,…,N,whereal,bl,cl, and dl 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
(12)𝒮l(xl,tk)=u(xl,tk),𝒮l(xl+1,tk)=u(xl+1,tk),𝒮l′′(xl,tk)=Q(xl,tk),𝒮l′′(xl+1,tk)=Q(xl+1,tk).Through the straightforward calculations, we obtain the following expressions:(13)al=(1-2cosh(ϑh))Q(xl,tk)+Q(xl+1,tk)2ϑ2sinh(ϑh)+u(xl+1,tk)-u(xl,tk)2sinh(ϑh),bl=Q(xl,tk)ϑ2,cl=Q(xl,tk)-Q(xl+1,tk)2ϑ2sin(ϑh)+u(xl+1,tk)-u(xl,tk)2sin(ϑh),dl=u(xl,tk)-Q(xl,tk)ϑ2,in whichl=0,1,…,M-1,k=0,1,…,N.From (11), one has
(14)𝒮l-1(x,tk)=al-1sinh(ϑ(x-xl-1))+bl-1cosh(ϑ(x-xl-1))+cl-1sin(ϑ(x-xl-1))+dl-1,x∈[xl-1,xl],l=1,…,M,k=0,1,…,N.In view of the continuity of the first-order derivative at the common nodes(xl,u(xl,tk)), that is, 𝒮l-1′(xl,tk)=𝒮l′(xl,tk), we get the following relation:
(15)λ1Q(xl-1,tk)+λ2Q(xl,tk)+λ3Q(xl+1,tk)=1h2[μ1u(xi-1,tk)+μ2u(xl,tk)+μ3u(xl+1,tk)],l=1,…,M-1,k=0,1,…,N,
where
(16)λ1=1θ2(cosh(θ)(2cosh(θ)-1)sinh(θ)-cos(θ)sin(θ)-2sinh(θ)),λ2=1θ2(1-3cosh(θ)sinh(θ)+1+cos(θ)sin(θ)),λ3=1θ2(1sinh(θ)-1sin(θ)),μ1=-(cosh(θ)sinh(θ)+cos(θ)sin(θ)),μ2=1+cosh(θ)sinh(θ)+1+cos(θ)sin(θ),μ3=-(1sinh(θ)+1sin(θ)),andθ=ϑh.Obviously, from the expressions ofμ1,μ2, and μ3, we find that μ1+μ2+μ3=0.Whenϑ→0, then θ→0; the spline function given by (15) reduces to the ordinary cubic spline function:
(17)h2[Q(xl-1,tk)+4Q(xl,tk)+Q(xl+1,tk)]=6[u(xl-1,tk)-2u(xl,tk)+u(xl+1,tk)],l=1,…,M-1,k=0,1,…,N.If we use∂2u(xl-1,tk)/∂x2,∂2u(xl,tk)/∂x2; and ∂2u(xl+1,tk)/∂x2 to take the place of Q(xl-1,tk),Q(xl,tk); and Q(xl+1,tk) in (15), respectively, then one gets
(18)λ1∂2u(xl-1,tk)∂x2+λ2∂2u(xl,tk)∂x2+λ3∂2u(xl+1,tk)∂x2=1h2[μ1u(xl-1,tk)+μ2u(xl,tk)+μ3u(xl+1,tk)]+𝒪(hp),wherep>0.Using Taylor’s series for (18) at (xl,tk), we get the truncation error as follows:
(19)Eik={1h(μ1-μ3)∂∂x+12[2(λ1+λ2+λ3)-(μ1+μ3)]∂2∂x2+h6[6(λ3-λ1)-(μ3-μ1)]∂3∂x3+h224[12(λ1+λ3)-(μ1+μ3)]∂4∂x4+h3120[20(λ3-λ1)-(μ3-μ1)]∂5∂x5+h4720[30(λ1+λ3)-(μ1+μ3)]∂6∂x6}×u(xl,tk)+𝒪(hp).From the previous expression, we find that we can obtain different difference schemes by choosing different parametersλ1,λ2,λ3,μ1, and μ2. In this paper, we take λ1=λ3,λ2=10λ1, and μ1=μ3=12λ1, and then p=4. Under this assumption, (18) can be rewritten as
(20)∂2u(xl,tk)∂x2=1h212δx2u(xl,tk)12+δx2+𝒪(h4),
which is the same as in [32], where δx2 is the second-order central difference operator with respect to x.


## 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 Lemma4, (1) can be written as
(21)RLD0,tαu(x,t)=Kα∂2u(x,t)∂x2-Cαu(x,t)+f(x,t)+∑j=01‍∂ju(x,0)∂tjtj-αΓ(j+1-α).
Applying the initial value conditions (2) yields
(22)RLD0,tαu(x,t)=Kα∂2u(x,t)∂x2-Cαu(x,t)+f(x,t).The Riemann-Liouville derivative (with homogeneous initial value conditions) can be approximated by the following scheme:(23)RLD0,tαu(x,t)=1τα∑j=0[t/τ]‍ϖq,j(α)u(x,t-jτ)+𝒪(τq).There exist various choices of the generating functions which can lead to different approximation orderq.Letϖq(z,α) be the generating function with coefficients ϖq,j(α); that is,
(24)ϖq(z,α)=∑j=0∞‍ϖq,j(α)zj.If we take the generating function in the form(25)ϖ1(z,α)=(1-z)α,
that is,
(26)(1-z)α=∑j=0∞‍ϖ1,j(α)zj,
then we can obtain the following first-order scheme:
(27)RLD0,tαu(x,t)=1τα∑j=0[t/τ]‍ϖ1,j(α)u(x,t-jτ)+𝒪(τ),
in which
(28)ϖ1,0(α)=1,ϖ1,j(α)=(-1)j(αj)=(-1)jΓ(1+α)Γ(j+1)Γ(1+α-j),j≥1,
and  these coefficients have the following recursive relations:(29)ϖ1,0(α)=1,ϖ1,j(α)=(1-1+αj)ϖ1,j-1(α),j≥1.If we take the generating function as(30)ϖ2(z,α)=(32-2z+12z2)α,
that is,
(31)(32-2z+12z2)α=∑j=0∞‍ϖ2,j(α)zj,
which leads to order q=2 [30], then(32)DRL0,tαu(x,t)=1τα∑j=0[t/τ]ϖ2,j(α)u(x,t-jτ)+𝒪(τ2).Here we introduce two methods to determine the corresponding coefficientsϖ2,j(α) 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,j(α).By using(33)(32-2z+12z2)α=(32)α(1-z)α(1-13z)αand (26), one gets
(34)(32-2z+12z2)α=(32)α(1-z)α(1-13z)α=(32)α[∑j=0∞‍(-1)j(αj)zj][∑j=0∞‍(-13)j(αj)zj]=(32)α∑m=0∞‍∑n=0∞‍(-1)m(-13)n(αm)(αn)zm+n=(32)α∑j=0∞‍[∑m=0j‍(-1)m(-13)j-m(αm)(αj-m)]zj=∑j=0∞‍[(32)α∑m=0j‍(-1)j(13)m(αm)(αj-m)]zj=∑j=0∞‍[(32)α∑m=0j‍(13)mϖ1,m(α)ϖ1,j-m(α)]zj.Comparing the previous equation with (31), one obtains
(35)ϖ2,j(α)=(32)α∑m=0j‍(13)mϖ1,m(α)ϖ1,j-m(α),=(-1)j(32)α×∑m=0j‍(13)m×Γ2(1+α)m!(j-m)!Γ(1+α-m)Γ(1+α+m-j),which has the following recursive relation:(36)ϖ2,0(α)=(32)α,ϖ2,j+1(α)=(32)α∑m=0j+1‍(13)mϖ1,m(α)ϖ1,j+1-m(α)=(32)α∑m=0j‍(13)mϖ1,m(α)ϖ1,j+1-m(α)+(32)α(13)j+1ϖ1,j+1(α)=(32)α∑m=0j‍(13)m(1-1+αj+1-m)ϖ1,m(α)ϖ1,j-m(α)+(32)α(13)j+1ϖ1,j+1(α)=ϖ2,j(α)-(1+α)(32)α(13)j∑m=0j‍3mm+1ϖ1,m(α)ϖ1,j-m(α)+(32)α(13)j+1(1-1+αj+1)ϖ1,j(α),j≥0.
In summary, we have the following lemma.Lemma 6.
The coefficientsϖ2,j(α)(j=0,1,…)(1<α<2) satisfy
(37)(1)ϖ2,0(α)=(32)α>0,ϖ2,1(α)=-43α(32)α<0,ϖ2,2(α)=α(8α-5)9(32)α>0,ϖ2,3(α)=4α(1-α)(8α-7)81(32)α<0,ϖ2,4(α)=α(α-1)(64α2-176α+123)486(32)α>0,ϖ2,5(α)=2α(α-1)(2-α)(64α2-208α+183)3645(32)α>0,⋮(2)∑j=0∞‍ϖ2,j(α)=0,∑j=0∞‍(-1)jϖ2,j(α)=4α.Proof.
(
1
) In view of (36), we easily obtain the expressions of the coefficients ϖ2,j(α)(j=0,1,…).
(
2
) In (31), if we take z≡1, which leads to ∑j=0∞ϖ2,j(α)=0. Similarly, choosing z≡-1 and substitution it into (31) yield ∑j=0∞(-1)jϖ2,j(α)=4α. This completes the proof.Finally, letulk be the numerical approximation of u(xl,tk); substituting the expansions (20) and (32) in (22) and removing the higher order term, one obtains Numerical Algorithm I as follows:
(38)(12κ1-κ2-ϖ2,0(α))ul-1k-2(12κ1+5κ2+5ϖ2,0(α))ulk+(12κ1-κ2-ϖ2,0(α))ul+1k=ϖ2,1(α)(ul-1k-1+10ulk-1+ul+1k-1)+∑j=2k‍ϖ2,j(α)(ul-1k-j+10ulk-j+ul+1k-j)-τα(fl-1k+10flk+fl+1k),
for l=1,…,M-1,k=1,…,N. The initial and boundary value conditions can be discretized by
(39)ul0=0,l=0,1,…,M,1τδtul01+(1/2)δt=1τul1-ul01+(1/2)δt=0,l=0,1,…,M,u0k=φ1(kτ),uMk=φ2(kτ),k=0,1,…,N,whereflk=f(xl,tk),κ1=Kατα/h2, and κ2=Cατα.Next, we study the solvability and the local truncation error.For convenience, denote(40)U0=(0,…,0︸M-1)T,Uk=(u1k,u2k,…,uM-1k)T,k=1,2,…,N,Fk=(f1k,f2k,…,fM-1k)T,k=0,1,…,N,Λ1=(1011101⋱⋱⋱1101110),Λ2=(-211-21⋱⋱⋱1-211-2).Then we can give the following matrix form of the difference scheme (38):
(41)(12κ1Λ2-κ2Λ1-A^0)Uk=A^1Uk-1+∑j=2k‍A^jUk-j-ταΛ1Fk+Y^,
where
(42)A^j=ϖ2,j(α)Λ1,j=0,1,…,k,Y^=(Ψ^1,0,…,0︸M-3,Ψ^M-1)T,Ψ^1=-(12κ1-κ2-ϖ2,0(α))φ1(tk)+∑j=1k‍ϖ2,j(α)φ1(tk-j)-ταf(x0,tk),Ψ^M-1=-(12κ1-κ2-ϖ2,0(α))φ2(tk)+∑j=1k‍ϖ2,j(α)φ2(tk-j)-ταf(xM,tk).Theorem 7.
The difference equation (38) is uniquely solvable.Proof.
From Lemma5, we know that the eigenvalues of the tridiagonal matrix (12κ1Λ2-κ2Λ1-A^0) are
(43)ζ^l=-48κ1sin2(πl2M)-4(κ2+ϖ2,0(α))[3-sin2(πl2M)],l=1,2,…,M-1.
Due to κ1, κ2, and ϖ2,0(α)>0, hence,
(44)det(12κ1Λ2-κ2Λ1-A^0)=∏‍l=1M-1ζ^l≠0.
That is, the matrix (12κ1Λ2-κ2Λ1-A^0) 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 𝒪(τ2+h4).Proof.
We now define the local truncation error as(45)R^lk=1τα∑j=0k‍ϖ2,j(α)u(xl,tk-j)-(Kαh212δx2u(xl,tk)12+δx2)-Cαu(xl,tk)-f(xl,tk).
Applying (20), (22), and (32), we have
(46)R^lk=[1τα∑j=0k‍ϖ2,j(α)u(xl,tk-j)-RLD0,tαu(xl,tk)]+Kα[∂2u(xl,tk)∂x2-(1h212δx2u(xl,tk)12+δx2)]=𝒪(τ2+h4).
This finishes the proof.In the following, we discuss the stability of Numerical AlgorithmI by the fractional Fourier method [16].LetUlk be the approximate solution of (38) and define
(47)ɛlk=ulk-Ulk,l=1,2,…,M-1,k=0,1,…,N,ɛk=(ɛ1k,ɛ2k,…,ɛM-1k)T.So, we can easily obtain the following round-off error equation:(48)(12κ1-κ2-ϖ2,0(α))ɛl-1k-2(12κ1+5κ2+5ϖ2,0(α))ɛlk+(12κ1-κ2-ϖ2,0(α))ɛl+1k=ϖ2,1(α)(ɛl-1k-1+10ɛlk-1+ɛl-1k-1)+∑j=2k‍ϖ2,j(2-α)(ɛl-1k-j+10ɛlk-j+ɛl+1k-j).We suppose that the solution of (48) has the following form:
(49)ɛlk=ξkexp(iϕlh),
where ϕ is any of the real spatial wave numbers supported by the lattice. Substituting (49) into (48) gives
(50)[12κ1sin2((1/2)ϕh)3-sin2((1/2)ϕh)+κ2+ϖ2,0(α)]ξk=-∑j=1k‍ϖ2,j(α)ξk-j.
If we write
(51)ξk+1=ηξk
and assume that η≡η(ϕ) is independent of time, then we obtain the following expression of the amplification factor:
(52)12κ1sin2((1/2)ϕh)3-sin2((1/2)ϕh)+κ2=-∑j=0k‍ϖ2,j(α)η-j.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
(53)12κ1sin2((1/2)ϕh)3-sin2((1/2)ϕh)+κ2≥-∑j=0k‍(-1)jϖ2,j(α)=S^α,k.
Although Sα,k depends on k, we can estimate it by its limit value
(54)S^α=limk→∞S^α,k=-∑j=0∞‍(-1)jϖ2,j(α)=-4α.
In this case, the stability condition becomes
(55)12κ1sin2((1/2)ϕh)3-sin2((1/2)ϕh)+κ2≥S^α.Note thatS^α 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.


### 3.2.  Numerical Algorithm II

Differentiating (1) with order 2-α in the sense of Riemann-Liouville yields
(56)RLD0,t2-α(CD0,tαu(x,t))=RLD0,t2-α(Kα∂2u(x,t)∂x2-Cαu(x,t)+f(x,t)).
Due to
(57)CD0,tαu(x,t)=D0,t-(2-α)(∂2u(x,t)∂t2),
it immediately follows that
(58)RLD0,t2-α(CD0,tαu(x,t))=(RLD0,t2-αD0,t-(2-α))(∂2u(x,t)∂t2)=∂2u(x,t)∂t2.Combing (56) with (58) gives
(59)∂2u(x,t)∂t2=RLD0,t2-α(Kα∂2u(x,t)∂x2-Cαu(x,t))+g(x,t),whereg(x,t)=RLD0,t2-αf(x,t).Next, for the second-order time derivative∂2u(xl,tk)/∂t2, we approximate it by the usual three-point central difference formula
(60)∂2u(xl,tk)∂t2=u(xl,tk-1)-2u(xl,tk)+u(xl,tk+1)τ2+𝒪(τ2).Letulk be the numerical approximation of u(xl,tk). Substituting the expansions (20), (32), and (60) in (59) and removing the higher-order terms, one obtains Numerical Algorithm II as follows:
(61)ul-1k+1+10ulk+1+ul+1k+1=[2+ϖ2,0(2-α)(12κ1-κ2)]ul-1k+2[10-ϖ2,0(2-α)(12κ1+5κ2)]ulk+[2+ϖ2,0(2-α)(12κ1-κ2)]ul+1k-ul-1k-1-10ulk-1-ul+1k-1+12κ1∑j=1k‍ϖ2,j(2-α)[ul-1k-j-2ulk-j+ul+1k-j]-κ2∑j=1k‍ϖ2,j(2-α)[ul-1k-j+10ulk-j+ul+1k-j]+τ2[gl-1k+10glk+gl+1k],
where l=1,…,M-1,k=1,…,N-1. The initial and boundary value conditions can be discretized by
(62)ul0=0,l=0,1,…,M,1τδtul01+(1/2)δt=1τul1-ul01+(1/2)δt=0,l=0,1,…,M,u0k=φ1(kτ),uMk=φ2(kτ),k=0,1,…,N,whereflk=f(xl,tk),κ1=Kατα/h2, and κ2=Cατα.Similar to Numerical Algorithm I, we discuss the solvability and local truncation error of Numerical Algorithm II.Denote(63)Gk=(g1k,g2k,…,gM-1k)T,k=0,1,…,N.
Then the matrix form of Numerical Algorithm II (61) is written as
(64)Λ1Uk+1=(2Λ1+A~0)Uk-Λ1Uk-1+∑j=1k‍A~jUk-j+τ2Λ1Gk+Y~,k=1,…,N-1,
where
(65)A~j=ϖ2,j(2-α)(12κ1Λ2-κ2Λ1),j=0,1,…,k,Y~=(Ψ~1,0,…,0︸M-3,Ψ~M-1)T,Ψ~1=-φ1(tk+1)+2φ1(tk)-φ1(tk-1)+(12κ1-κ2)∑j=0k‍ϖ2,j(2-α)φ1(tk-j)+τ2g(x0,tk),Ψ~M-1=-φ2(tk+1)+2φ2(tk)-φ2(tk-1)+(12κ1-κ2)∑j=0k‍ϖ2,j(2-α)φ2(tk-j)+τ2g(xM,tk).Theorem 9.
The difference equation (61) is uniquely solvable.Proof.
Obviously, the eigenvalues of the tridiagonal matrixΛ1 are given by
(66)ζ~l=12-4sin2(πl2M),l=1,2,…,M-1.
Hence,(67)det(Λ1)=∏‍l=1M-1ζ~l>0.
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(68)1τ2-α∑j=0k‍ϖ2,j(2-α)=(kτ)α-2Γ(α-1)+𝒪(τ2),k=1,…,N.Proof.
Letu(x,t)≡1 and then from Definition 3 we have
(69)RLD0,t2-αu(x,t)|t=tk=1Γ(α-1)∂∂t∫0t‍u(x,s)(t-s)2-αds|t=tk=(kτ)α-2Γ(α-1).
It follows from (32) that
(70)RLD0,t2-αu(x,t)|t=tk=1τ2-α∑j=0k‍ϖ2,j(2-α)u(x,tk-j)+𝒪(τ2)=1τ2-α∑j=0k‍ϖ2,j(2-α)+𝒪(τ2).
Combining (69) and (77) leads to
(71)1τ2-α∑j=0k‍ϖ2,j(2-α)=(kτ)α-2Γ(α-1)+𝒪(τ2),k=1,…,N.
This ends the proof.Theorem 11.
The local truncation error of difference scheme (61) is 𝒪(τ2+h4).Proof.
We now define(72)R~lk=u(xl,tk-1)-2u(xl,tk)+u(xl,tk+1)τ2-Kατ2-α∑j=0k‍ϖ2,j(2-α)(1h212δx2u(xl,tk-j)12+δx2)-Cατ2-α∑j=0k‍ϖ2,j(2-α)u(xl,tk-j)-g(xl,tk).
By using (20), (32), (59), and (60), one gets
(73)R~lk=u(xl,tk-1)-2u(xl,tk)+u(xl,tk+1)τ2-Kατ2-α∑j=0k‍ϖ2,j(2-α)(1h212δx2u(xl,tk-l)12+δx2)+Cατ2-α∑j=0k‍ϖ2,j(2-α)u(xl,tk-j)-∂2u(xl,tk)∂t2+KαRLD0,t2-α(∂2u(xl,tk)∂x2)-CαRLD0,t2-αu(xl,tk)=[u(xl,tk-1)-2u(xl,tk)+u(xl,tk+1)τ2-∂2u(xl,tk)∂t2]+Kατ2-α∑j=0k‍ϖ2,j(2-α)(∂2u(xl,tk-l)∂x2-1h212δx2u(xl,tk-l)12+δx2)+Kα[∑j=0k‍ϖ2,j(2-α)(∂2u(xl,tk-l)∂x2)RLD0,t2-α(∂2u(xl,tk)∂x2)-1τ2-α∑j=0k‍ϖ2,j(2-α)(∂2u(xl,tk-l)∂x2)]+Cα[1τ2-α∑j=0k‍ϖ2,j(2-α)u(xl,tk-j)-RLD0,t2-αu(xl,tk)]=𝒪(τ2)+𝒪(τ2h4)+(kτ)α-2Γ(α-1)𝒪(h4).
Noticing that1≤k≤N, then
(74)(1T)2-α≤(kτ)α-2≤(NT)2-α,
so we can get
(75)R~lk=𝒪(τ2+h4).
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:(76)ɛl-1k+1+10ɛlk+1+ɛl+1k+1=[2+ϖ2,0(2-α)(12κ1-κ2)]ɛl-1k+2[10-ϖ2,0(2-α)(12κ1+5κ2)]ɛlk+[2+ϖ2,0(2-α)(12κ1-κ2)]ɛl+1k-ɛl-1k-1-10ɛlk-1-ɛl+1k-1+12κ1∑j=1k‍ϖ2,j(2-α)[ɛl-1k-j-2ɛlk-j+ɛl+1k-j]-κ2∑j=1k‍ϖ2,j(2-α)[ɛl-1k-j+10ɛlk-j+ɛl+1k-j],k=1,…,N-1.
Substituting (49) into (76) yields
(77)[3-sin2(12ϕh)]ξk+1={2[3-sin2(12ϕh)]-12κ1ϖ2,0(2-α)sin2(12ϕh)-κ2ϖ2,0(2-α)[3-sin2(12ϕh)]}ξk-[3-sin2(12ϕh)]ξk-1-κ2[3-sin2(12ϕh)]×∑j=1k‍ϖ2,j(2-α)ξk-j-12κ1sin2(12ϕh)×∑j=1k‍ϖ2,j(2-α)ξk-j.
From (51) and (77), we have
(78)[3-sin2(12ϕh)]η={2[3-sin2(12ϕh)]-12κ1ϖ2,0(2-α)sin2(12ϕh)-κ2ϖ2,0(2-α)[3-sin2(12ϕh)]}-[3-sin2(12ϕh)]η-1-κ2[3-sin2(12ϕh)]∑j=1k‍ϖ2,j(2-α)η-j-12κ1sin2(12ϕh)∑j=1k‍ϖ2,j(2-α)η-j.Considering the extreme valueη=-1, we can obtain the following stability condition from (78):
(79)sin2(12ϕh)≤12-3κ24+(12κ1-κ2)∑j=0k(-1)jϖ2,j(2-α)≡S~α,k.For large enoughk, we can estimate S~α,k by its limit value
(80)S~α=limk→∞S~α,k=12-3κ24+(12κ1-κ2)∑j=0∞(-1)jϖ2,j(2-α)=12-3κ24+(12κ1-κ2)42-α.It follows that(81)sin2(12ϕh)≤12-3κ24+(12κ1-κ2)42-α=S~α.Because the sine function is bounded by1, from (81) one finds that Numerical Algorithm II is to be stable if S~α≥1, that is,(82)12-3κ24+(12κ1-κ2)42-α≥1.Furthermore, stability condition (82) can be rewritten as
(83)3h2(4-Cατα)h2+41-ατα(12Kα-Cαh2)≥4,whereKα>0 and Cα>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 Theorems8 and 11, we know that Numerical Algorithm I (see (38)) and Numerical Algorithm II (see (61)) are all consistent with the local truncation errors R^lk=𝒪(τ2+h4) and R~lk=𝒪(τ2+h4) because R^lk,R~lk→0 as τ,h→0. Therefore, according to Lax’s Equivalence Theorem [33], Numerical Algorithm I unconditionally converges at the same order 𝒪(τ2+h4), but Numerical Algorithm II converges at the same order 𝒪(τ2+h4) under condition (83).


## 3.1. Numerical Algorithm I

Based on Lemma4, (1) can be written as
(21)RLD0,tαu(x,t)=Kα∂2u(x,t)∂x2-Cαu(x,t)+f(x,t)+∑j=01‍∂ju(x,0)∂tjtj-αΓ(j+1-α).
Applying the initial value conditions (2) yields
(22)RLD0,tαu(x,t)=Kα∂2u(x,t)∂x2-Cαu(x,t)+f(x,t).The Riemann-Liouville derivative (with homogeneous initial value conditions) can be approximated by the following scheme:(23)RLD0,tαu(x,t)=1τα∑j=0[t/τ]‍ϖq,j(α)u(x,t-jτ)+𝒪(τq).There exist various choices of the generating functions which can lead to different approximation orderq.Letϖq(z,α) be the generating function with coefficients ϖq,j(α); that is,
(24)ϖq(z,α)=∑j=0∞‍ϖq,j(α)zj.If we take the generating function in the form(25)ϖ1(z,α)=(1-z)α,
that is,
(26)(1-z)α=∑j=0∞‍ϖ1,j(α)zj,
then we can obtain the following first-order scheme:
(27)RLD0,tαu(x,t)=1τα∑j=0[t/τ]‍ϖ1,j(α)u(x,t-jτ)+𝒪(τ),
in which
(28)ϖ1,0(α)=1,ϖ1,j(α)=(-1)j(αj)=(-1)jΓ(1+α)Γ(j+1)Γ(1+α-j),j≥1,
and  these coefficients have the following recursive relations:(29)ϖ1,0(α)=1,ϖ1,j(α)=(1-1+αj)ϖ1,j-1(α),j≥1.If we take the generating function as(30)ϖ2(z,α)=(32-2z+12z2)α,
that is,
(31)(32-2z+12z2)α=∑j=0∞‍ϖ2,j(α)zj,
which leads to order q=2 [30], then(32)DRL0,tαu(x,t)=1τα∑j=0[t/τ]ϖ2,j(α)u(x,t-jτ)+𝒪(τ2).Here we introduce two methods to determine the corresponding coefficientsϖ2,j(α) 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,j(α).By using(33)(32-2z+12z2)α=(32)α(1-z)α(1-13z)αand (26), one gets
(34)(32-2z+12z2)α=(32)α(1-z)α(1-13z)α=(32)α[∑j=0∞‍(-1)j(αj)zj][∑j=0∞‍(-13)j(αj)zj]=(32)α∑m=0∞‍∑n=0∞‍(-1)m(-13)n(αm)(αn)zm+n=(32)α∑j=0∞‍[∑m=0j‍(-1)m(-13)j-m(αm)(αj-m)]zj=∑j=0∞‍[(32)α∑m=0j‍(-1)j(13)m(αm)(αj-m)]zj=∑j=0∞‍[(32)α∑m=0j‍(13)mϖ1,m(α)ϖ1,j-m(α)]zj.Comparing the previous equation with (31), one obtains
(35)ϖ2,j(α)=(32)α∑m=0j‍(13)mϖ1,m(α)ϖ1,j-m(α),=(-1)j(32)α×∑m=0j‍(13)m×Γ2(1+α)m!(j-m)!Γ(1+α-m)Γ(1+α+m-j),which has the following recursive relation:(36)ϖ2,0(α)=(32)α,ϖ2,j+1(α)=(32)α∑m=0j+1‍(13)mϖ1,m(α)ϖ1,j+1-m(α)=(32)α∑m=0j‍(13)mϖ1,m(α)ϖ1,j+1-m(α)+(32)α(13)j+1ϖ1,j+1(α)=(32)α∑m=0j‍(13)m(1-1+αj+1-m)ϖ1,m(α)ϖ1,j-m(α)+(32)α(13)j+1ϖ1,j+1(α)=ϖ2,j(α)-(1+α)(32)α(13)j∑m=0j‍3mm+1ϖ1,m(α)ϖ1,j-m(α)+(32)α(13)j+1(1-1+αj+1)ϖ1,j(α),j≥0.
In summary, we have the following lemma.Lemma 6.
The coefficientsϖ2,j(α)(j=0,1,…)(1<α<2) satisfy
(37)(1)ϖ2,0(α)=(32)α>0,ϖ2,1(α)=-43α(32)α<0,ϖ2,2(α)=α(8α-5)9(32)α>0,ϖ2,3(α)=4α(1-α)(8α-7)81(32)α<0,ϖ2,4(α)=α(α-1)(64α2-176α+123)486(32)α>0,ϖ2,5(α)=2α(α-1)(2-α)(64α2-208α+183)3645(32)α>0,⋮(2)∑j=0∞‍ϖ2,j(α)=0,∑j=0∞‍(-1)jϖ2,j(α)=4α.Proof.
(
1
) In view of (36), we easily obtain the expressions of the coefficients ϖ2,j(α)(j=0,1,…).
(
2
) In (31), if we take z≡1, which leads to ∑j=0∞ϖ2,j(α)=0. Similarly, choosing z≡-1 and substitution it into (31) yield ∑j=0∞(-1)jϖ2,j(α)=4α. This completes the proof.Finally, letulk be the numerical approximation of u(xl,tk); substituting the expansions (20) and (32) in (22) and removing the higher order term, one obtains Numerical Algorithm I as follows:
(38)(12κ1-κ2-ϖ2,0(α))ul-1k-2(12κ1+5κ2+5ϖ2,0(α))ulk+(12κ1-κ2-ϖ2,0(α))ul+1k=ϖ2,1(α)(ul-1k-1+10ulk-1+ul+1k-1)+∑j=2k‍ϖ2,j(α)(ul-1k-j+10ulk-j+ul+1k-j)-τα(fl-1k+10flk+fl+1k),
for l=1,…,M-1,k=1,…,N. The initial and boundary value conditions can be discretized by
(39)ul0=0,l=0,1,…,M,1τδtul01+(1/2)δt=1τul1-ul01+(1/2)δt=0,l=0,1,…,M,u0k=φ1(kτ),uMk=φ2(kτ),k=0,1,…,N,whereflk=f(xl,tk),κ1=Kατα/h2, and κ2=Cατα.Next, we study the solvability and the local truncation error.For convenience, denote(40)U0=(0,…,0︸M-1)T,Uk=(u1k,u2k,…,uM-1k)T,k=1,2,…,N,Fk=(f1k,f2k,…,fM-1k)T,k=0,1,…,N,Λ1=(1011101⋱⋱⋱1101110),Λ2=(-211-21⋱⋱⋱1-211-2).Then we can give the following matrix form of the difference scheme (38):
(41)(12κ1Λ2-κ2Λ1-A^0)Uk=A^1Uk-1+∑j=2k‍A^jUk-j-ταΛ1Fk+Y^,
where
(42)A^j=ϖ2,j(α)Λ1,j=0,1,…,k,Y^=(Ψ^1,0,…,0︸M-3,Ψ^M-1)T,Ψ^1=-(12κ1-κ2-ϖ2,0(α))φ1(tk)+∑j=1k‍ϖ2,j(α)φ1(tk-j)-ταf(x0,tk),Ψ^M-1=-(12κ1-κ2-ϖ2,0(α))φ2(tk)+∑j=1k‍ϖ2,j(α)φ2(tk-j)-ταf(xM,tk).Theorem 7.
The difference equation (38) is uniquely solvable.Proof.
From Lemma5, we know that the eigenvalues of the tridiagonal matrix (12κ1Λ2-κ2Λ1-A^0) are
(43)ζ^l=-48κ1sin2(πl2M)-4(κ2+ϖ2,0(α))[3-sin2(πl2M)],l=1,2,…,M-1.
Due to κ1, κ2, and ϖ2,0(α)>0, hence,
(44)det(12κ1Λ2-κ2Λ1-A^0)=∏‍l=1M-1ζ^l≠0.
That is, the matrix (12κ1Λ2-κ2Λ1-A^0) 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 𝒪(τ2+h4).Proof.
We now define the local truncation error as(45)R^lk=1τα∑j=0k‍ϖ2,j(α)u(xl,tk-j)-(Kαh212δx2u(xl,tk)12+δx2)-Cαu(xl,tk)-f(xl,tk).
Applying (20), (22), and (32), we have
(46)R^lk=[1τα∑j=0k‍ϖ2,j(α)u(xl,tk-j)-RLD0,tαu(xl,tk)]+Kα[∂2u(xl,tk)∂x2-(1h212δx2u(xl,tk)12+δx2)]=𝒪(τ2+h4).
This finishes the proof.In the following, we discuss the stability of Numerical AlgorithmI by the fractional Fourier method [16].LetUlk be the approximate solution of (38) and define
(47)ɛlk=ulk-Ulk,l=1,2,…,M-1,k=0,1,…,N,ɛk=(ɛ1k,ɛ2k,…,ɛM-1k)T.So, we can easily obtain the following round-off error equation:(48)(12κ1-κ2-ϖ2,0(α))ɛl-1k-2(12κ1+5κ2+5ϖ2,0(α))ɛlk+(12κ1-κ2-ϖ2,0(α))ɛl+1k=ϖ2,1(α)(ɛl-1k-1+10ɛlk-1+ɛl-1k-1)+∑j=2k‍ϖ2,j(2-α)(ɛl-1k-j+10ɛlk-j+ɛl+1k-j).We suppose that the solution of (48) has the following form:
(49)ɛlk=ξkexp(iϕlh),
where ϕ is any of the real spatial wave numbers supported by the lattice. Substituting (49) into (48) gives
(50)[12κ1sin2((1/2)ϕh)3-sin2((1/2)ϕh)+κ2+ϖ2,0(α)]ξk=-∑j=1k‍ϖ2,j(α)ξk-j.
If we write
(51)ξk+1=ηξk
and assume that η≡η(ϕ) is independent of time, then we obtain the following expression of the amplification factor:
(52)12κ1sin2((1/2)ϕh)3-sin2((1/2)ϕh)+κ2=-∑j=0k‍ϖ2,j(α)η-j.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
(53)12κ1sin2((1/2)ϕh)3-sin2((1/2)ϕh)+κ2≥-∑j=0k‍(-1)jϖ2,j(α)=S^α,k.
Although Sα,k depends on k, we can estimate it by its limit value
(54)S^α=limk→∞S^α,k=-∑j=0∞‍(-1)jϖ2,j(α)=-4α.
In this case, the stability condition becomes
(55)12κ1sin2((1/2)ϕh)3-sin2((1/2)ϕh)+κ2≥S^α.Note thatS^α 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.


## 3.2.  Numerical Algorithm II

Differentiating (1) with order 2-α in the sense of Riemann-Liouville yields
(56)RLD0,t2-α(CD0,tαu(x,t))=RLD0,t2-α(Kα∂2u(x,t)∂x2-Cαu(x,t)+f(x,t)).
Due to
(57)CD0,tαu(x,t)=D0,t-(2-α)(∂2u(x,t)∂t2),
it immediately follows that
(58)RLD0,t2-α(CD0,tαu(x,t))=(RLD0,t2-αD0,t-(2-α))(∂2u(x,t)∂t2)=∂2u(x,t)∂t2.Combing (56) with (58) gives
(59)∂2u(x,t)∂t2=RLD0,t2-α(Kα∂2u(x,t)∂x2-Cαu(x,t))+g(x,t),whereg(x,t)=RLD0,t2-αf(x,t).Next, for the second-order time derivative∂2u(xl,tk)/∂t2, we approximate it by the usual three-point central difference formula
(60)∂2u(xl,tk)∂t2=u(xl,tk-1)-2u(xl,tk)+u(xl,tk+1)τ2+𝒪(τ2).Letulk be the numerical approximation of u(xl,tk). Substituting the expansions (20), (32), and (60) in (59) and removing the higher-order terms, one obtains Numerical Algorithm II as follows:
(61)ul-1k+1+10ulk+1+ul+1k+1=[2+ϖ2,0(2-α)(12κ1-κ2)]ul-1k+2[10-ϖ2,0(2-α)(12κ1+5κ2)]ulk+[2+ϖ2,0(2-α)(12κ1-κ2)]ul+1k-ul-1k-1-10ulk-1-ul+1k-1+12κ1∑j=1k‍ϖ2,j(2-α)[ul-1k-j-2ulk-j+ul+1k-j]-κ2∑j=1k‍ϖ2,j(2-α)[ul-1k-j+10ulk-j+ul+1k-j]+τ2[gl-1k+10glk+gl+1k],
where l=1,…,M-1,k=1,…,N-1. The initial and boundary value conditions can be discretized by
(62)ul0=0,l=0,1,…,M,1τδtul01+(1/2)δt=1τul1-ul01+(1/2)δt=0,l=0,1,…,M,u0k=φ1(kτ),uMk=φ2(kτ),k=0,1,…,N,whereflk=f(xl,tk),κ1=Kατα/h2, and κ2=Cατα.Similar to Numerical Algorithm I, we discuss the solvability and local truncation error of Numerical Algorithm II.Denote(63)Gk=(g1k,g2k,…,gM-1k)T,k=0,1,…,N.
Then the matrix form of Numerical Algorithm II (61) is written as
(64)Λ1Uk+1=(2Λ1+A~0)Uk-Λ1Uk-1+∑j=1k‍A~jUk-j+τ2Λ1Gk+Y~,k=1,…,N-1,
where
(65)A~j=ϖ2,j(2-α)(12κ1Λ2-κ2Λ1),j=0,1,…,k,Y~=(Ψ~1,0,…,0︸M-3,Ψ~M-1)T,Ψ~1=-φ1(tk+1)+2φ1(tk)-φ1(tk-1)+(12κ1-κ2)∑j=0k‍ϖ2,j(2-α)φ1(tk-j)+τ2g(x0,tk),Ψ~M-1=-φ2(tk+1)+2φ2(tk)-φ2(tk-1)+(12κ1-κ2)∑j=0k‍ϖ2,j(2-α)φ2(tk-j)+τ2g(xM,tk).Theorem 9.
The difference equation (61) is uniquely solvable.Proof.
Obviously, the eigenvalues of the tridiagonal matrixΛ1 are given by
(66)ζ~l=12-4sin2(πl2M),l=1,2,…,M-1.
Hence,(67)det(Λ1)=∏‍l=1M-1ζ~l>0.
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(68)1τ2-α∑j=0k‍ϖ2,j(2-α)=(kτ)α-2Γ(α-1)+𝒪(τ2),k=1,…,N.Proof.
Letu(x,t)≡1 and then from Definition 3 we have
(69)RLD0,t2-αu(x,t)|t=tk=1Γ(α-1)∂∂t∫0t‍u(x,s)(t-s)2-αds|t=tk=(kτ)α-2Γ(α-1).
It follows from (32) that
(70)RLD0,t2-αu(x,t)|t=tk=1τ2-α∑j=0k‍ϖ2,j(2-α)u(x,tk-j)+𝒪(τ2)=1τ2-α∑j=0k‍ϖ2,j(2-α)+𝒪(τ2).
Combining (69) and (77) leads to
(71)1τ2-α∑j=0k‍ϖ2,j(2-α)=(kτ)α-2Γ(α-1)+𝒪(τ2),k=1,…,N.
This ends the proof.Theorem 11.
The local truncation error of difference scheme (61) is 𝒪(τ2+h4).Proof.
We now define(72)R~lk=u(xl,tk-1)-2u(xl,tk)+u(xl,tk+1)τ2-Kατ2-α∑j=0k‍ϖ2,j(2-α)(1h212δx2u(xl,tk-j)12+δx2)-Cατ2-α∑j=0k‍ϖ2,j(2-α)u(xl,tk-j)-g(xl,tk).
By using (20), (32), (59), and (60), one gets
(73)R~lk=u(xl,tk-1)-2u(xl,tk)+u(xl,tk+1)τ2-Kατ2-α∑j=0k‍ϖ2,j(2-α)(1h212δx2u(xl,tk-l)12+δx2)+Cατ2-α∑j=0k‍ϖ2,j(2-α)u(xl,tk-j)-∂2u(xl,tk)∂t2+KαRLD0,t2-α(∂2u(xl,tk)∂x2)-CαRLD0,t2-αu(xl,tk)=[u(xl,tk-1)-2u(xl,tk)+u(xl,tk+1)τ2-∂2u(xl,tk)∂t2]+Kατ2-α∑j=0k‍ϖ2,j(2-α)(∂2u(xl,tk-l)∂x2-1h212δx2u(xl,tk-l)12+δx2)+Kα[∑j=0k‍ϖ2,j(2-α)(∂2u(xl,tk-l)∂x2)RLD0,t2-α(∂2u(xl,tk)∂x2)-1τ2-α∑j=0k‍ϖ2,j(2-α)(∂2u(xl,tk-l)∂x2)]+Cα[1τ2-α∑j=0k‍ϖ2,j(2-α)u(xl,tk-j)-RLD0,t2-αu(xl,tk)]=𝒪(τ2)+𝒪(τ2h4)+(kτ)α-2Γ(α-1)𝒪(h4).
Noticing that1≤k≤N, then
(74)(1T)2-α≤(kτ)α-2≤(NT)2-α,
so we can get
(75)R~lk=𝒪(τ2+h4).
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:(76)ɛl-1k+1+10ɛlk+1+ɛl+1k+1=[2+ϖ2,0(2-α)(12κ1-κ2)]ɛl-1k+2[10-ϖ2,0(2-α)(12κ1+5κ2)]ɛlk+[2+ϖ2,0(2-α)(12κ1-κ2)]ɛl+1k-ɛl-1k-1-10ɛlk-1-ɛl+1k-1+12κ1∑j=1k‍ϖ2,j(2-α)[ɛl-1k-j-2ɛlk-j+ɛl+1k-j]-κ2∑j=1k‍ϖ2,j(2-α)[ɛl-1k-j+10ɛlk-j+ɛl+1k-j],k=1,…,N-1.
Substituting (49) into (76) yields
(77)[3-sin2(12ϕh)]ξk+1={2[3-sin2(12ϕh)]-12κ1ϖ2,0(2-α)sin2(12ϕh)-κ2ϖ2,0(2-α)[3-sin2(12ϕh)]}ξk-[3-sin2(12ϕh)]ξk-1-κ2[3-sin2(12ϕh)]×∑j=1k‍ϖ2,j(2-α)ξk-j-12κ1sin2(12ϕh)×∑j=1k‍ϖ2,j(2-α)ξk-j.
From (51) and (77), we have
(78)[3-sin2(12ϕh)]η={2[3-sin2(12ϕh)]-12κ1ϖ2,0(2-α)sin2(12ϕh)-κ2ϖ2,0(2-α)[3-sin2(12ϕh)]}-[3-sin2(12ϕh)]η-1-κ2[3-sin2(12ϕh)]∑j=1k‍ϖ2,j(2-α)η-j-12κ1sin2(12ϕh)∑j=1k‍ϖ2,j(2-α)η-j.Considering the extreme valueη=-1, we can obtain the following stability condition from (78):
(79)sin2(12ϕh)≤12-3κ24+(12κ1-κ2)∑j=0k(-1)jϖ2,j(2-α)≡S~α,k.For large enoughk, we can estimate S~α,k by its limit value
(80)S~α=limk→∞S~α,k=12-3κ24+(12κ1-κ2)∑j=0∞(-1)jϖ2,j(2-α)=12-3κ24+(12κ1-κ2)42-α.It follows that(81)sin2(12ϕh)≤12-3κ24+(12κ1-κ2)42-α=S~α.Because the sine function is bounded by1, from (81) one finds that Numerical Algorithm II is to be stable if S~α≥1, that is,(82)12-3κ24+(12κ1-κ2)42-α≥1.Furthermore, stability condition (82) can be rewritten as
(83)3h2(4-Cατα)h2+41-ατα(12Kα-Cαh2)≥4,whereKα>0 and Cα>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 Theorems8 and 11, we know that Numerical Algorithm I (see (38)) and Numerical Algorithm II (see (61)) are all consistent with the local truncation errors R^lk=𝒪(τ2+h4) and R~lk=𝒪(τ2+h4) because R^lk,R~lk→0 as τ,h→0. Therefore, according to Lax’s Equivalence Theorem [33], Numerical Algorithm I unconditionally converges at the same order 𝒪(τ2+h4), but Numerical Algorithm II converges at the same order 𝒪(τ2+h4) 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:(84)CD0,tαu(x,t)=∂2u(x,t)∂x2-u(x,t)+f(x,t),0≤t≤1,0<x<2,
where f(x,t)=α3[(6t3-α/Γ(4-α))+t3(1+(π2/4))]sin((π/2)x).
The analytical solution of this equation isu(x,t)=(αt)3sin((π/2)x). The initial conditions and the boundary conditions are satisfied with the exact solution u(x,t) 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 h in Figures 1 and 2, respectively, and we observe that the numerical solution is in line with the analytical solution and the velocity u 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)).
α
The  maximum error
Temporal  convergence order
Spatial  convergence order
1.1
h
=
1
4
,
τ
=
1
2
8.057755
e
-
002
—
—
h
=
1
8
,
τ
=
1
8
9.573435
e
-
003
1.5366
3.0733
h
=
1
16
,
τ
=
1
32
6.103333
e
-
004
1.9857
3.9714
h
=
1
32
,
τ
=
1
128
3.827561
e
-
005
1.9976
3.9951
h
=
1
64
,
τ
=
1
512
2.394875
e
-
006
1.9992
3.9984
1.3
h
=
1
4
,
τ
=
1
2
8.079213
e
-
002
—
—
h
=
1
8
,
τ
=
1
8
1.170294
e
-
002
1.3937
2.7873
h
=
1
16
,
τ
=
1
32
7.818506
e
-
004
1.9519
3.9038
h
=
1
32
,
τ
=
1
128
5.015017
e
-
005
1.9813
3.9626
h
=
1
64
,
τ
=
1
512
3.166977
e
-
006
1.9925
3.9851
1.5
h
=
1
4
,
τ
=
1
2
7.117932
e
-
002
—
—
h
=
1
8
,
τ
=
1
8
1.300290
e
-
002
1.2263
2.4526
h
=
1
16
,
τ
=
1
32
9.519179
e
-
004
1.8859
3.7719
h
=
1
32
,
τ
=
1
128
6.385188
e
-
005
1.9490
3.8980
h
=
1
64
,
τ
=
1
512
4.135300
e
-
006
1.9743
3.9487
1.7
h
=
1
4
,
τ
=
1
2
5.069482
e
-
002
—
—
h
=
1
8
,
τ
=
1
8
1.157045
e
-
002
1.0657
2.1314
h
=
1
16
,
τ
=
1
32
9.766015
e
-
004
1.7833
3.5665
h
=
1
32
,
τ
=
1
128
7.054045
e
-
005
1.8956
3.7912
h
=
1
64
,
τ
=
1
512
4.819461
e
-
006
1.9358
3.8715
1.9
h
=
1
4
,
τ
=
1
2
1.938949
e
-
002
—
—
h
=
1
8
,
τ
=
1
8
4.643245
e
-
003
1.0310
2.0621
h
=
1
16
,
τ
=
1
32
5.122387
e
-
004
1.5901
3.1802
h
=
1
32
,
τ
=
1
128
4.251159
e
-
005
1.7954
3.5909
h
=
1
64
,
τ
=
1
512
3.231450
e
-
006
1.8588
3.7176Figure 1
Dependence of the numerical solution for different values ofα by Numerical Algorithm I (see (38)) with τ=0.002,h=0.0025.Figure 2
Dependence of the numerical solution for different values ofα by Numerical Algorithm I (see (38)) with τ=0.0001,h=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,h=0.01.Next, we use Numerical Algorithm II (61). In this situation, we take (60), the equivalence equation of (1), where g(x,t)=6tα3[1+(tα/Γ(2+α))(1+(π2/4))]sin((π/2)x).Figures4 and 5 display that the velocity u increases with the increase of order α for different τ and h. 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)).
α
The maximum error
The temporal convergence order
1.1
h
=
1
5
,
τ
=
1
50
5.030598
e
-
004
—
h
=
1
5
,
τ
=
1
100
1.398353
e
-
004
1.8470
h=15,τ=1200
4.521910
e
-
005
1.6287
h
=
1
5
,
τ
=
1
400
2.132246
e
-
005
1.0846
1.2
h
=
1
6
,
τ
=
1
50
4.445704
e
-
004
—
h
=
1
6
,
τ
=
1
100
1.173968
e
-
004
1.9210
h
=
1
6
,
τ
=
1
200
3.399818
e
-
005
1.7879
h
=
1
6
,
τ
=
1
400
1.295545
e
-
005
1.3919
1.3
h
=
1
5
,
τ
=
1
50
3.816458
e
-
004
—
h
=
1
7
,
τ
=
1
100
9.909362
e
-
005
1.9454
h
=
1
7
,
τ
=
1
200
2.713495
e
-
005
1.8686
h
=
1
7
,
τ
=
1
400
8.986675
e
-
006
1.5943
1.4
h
=
1
5
,
τ
=
1
50
3.193648
e
-
004
—
h
=
1
8
,
τ
=
1
100
8.224534
e
-
005
1.9572
h
=
1
8
,
τ
=
1
200
2.186799
e
-
005
1.9111
h
=
1
8
,
τ
=
1
400
6.641262
e
-
006
1.7193
1.5
h
=
1
10
,
τ
=
1
50
2.574223
e
-
004
—
h
=
1
10
,
τ
=
1
100
6.580906
e
-
005
1.9678
h
=
1
10
,
τ
=
1
200
1.700153
e
-
005
1.9526
h
=
1
10
,
τ
=
1
400
4.689779
e
-
006
1.8581
1.6
h
=
1
15
,
τ
=
1
50
1.972538
e
-
004
—
h
=
1
15
,
τ
=
1
100
5.020982
e
-
005
1.9740
h
=
1
15
,
τ
=
1
200
1.272752
e
-
005
1.9800
h
=
1
15
,
τ
=
1
400
3.268637
e
-
006
1.9612
1.7
h
=
1
20
,
τ
=
1
50
1.402871
e
-
004
—
h
=
1
20
,
τ
=
1
100
3.570742
e
-
005
1.9741
h
=
1
20
,
τ
=
1
200
9.023412
e
-
006
1.9845
h
=
1
20
,
τ
=
1
400
2.285957
e
-
006
1.9809
1.8
h
=
1
20
,
τ
=
1
50
8.770513
e
-
005
—
h
=
1
20
,
τ
=
1
100
2.235014
e
-
005
1.9724
h
=
1
20
,
τ
=
1
200
5.655983
e
-
006
1.9824
h
=
1
20
,
τ
=
1
400
1.438349
e
-
006
1.9754
1.9
h
=
1
20
,
τ
=
1
50
4.067501
e
-
005
—
h
=
1
20
,
τ
=
1
100
1.038309
e
-
005
1.9699
h
=
1
20
,
τ
=
1
200
2.636242
e
-
006
1.9777
h
=
1
20
,
τ
=
1
400
6.778414
e
-
007
1.9595Figure 4
Dependence of the numerical solution for different values ofα by Numerical Algorithm II (see (61)) with τ=0.002,h=0.0025.Figure 5
Dependence of the numerical solution for different values ofα by Numerical Algorithm II (see (61)) with τ=0.0001,h=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,h=0.002.Finally, we checked the stability condition given in (83) in the following way. Firstly, we choose h=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 h=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 h=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 h=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 𝒪(τ2+h4), which are the best results till now. Finally, the numerical results given in this paper show the effectiveness of the derived numerical algorithms.



---

*Source: 493406-2013-07-21.xml*