# Fully Discrete Local Discontinuous Galerkin Approximation for Time-Space Fractional Subdiffusion/Superdiffusion Equations

Authors: Meilan Qiu; Liquan Mei; Dewang Li
Journal: Advances in Mathematical Physics
 (2017)

---

## Abstract

We focus on developing the finite difference (i.e., backward Euler difference or second-order central difference)/local discontinuous Galerkin finite element mixed method to construct and analyze a kind of efficient, accurate, flexible, numerical schemes for approximately solving time-space fractional subdiffusion/superdiffusion equations. Discretizing the time Caputo fractional derivative by using the backward Euler difference for the derivative parameter (0<α<1) or second-order central difference method for (1<α<2), combined with local discontinuous Galerkin method to approximate the spatial derivative which is defined by a fractional Laplacian operator, two high-accuracy fully discrete local discontinuous Galerkin (LDG) schemes of the time-space fractional subdiffusion/superdiffusion equations are proposed, respectively. Through the mathematical induction method, we show the concrete analysis for the stability and the convergence under the L2 norm of the LDG schemes. Several numerical experiments are presented to validate the proposed model and demonstrate the convergence rate of numerical schemes. The numerical experiment results show that the fully discrete local discontinuous Galerkin (LDG) methods are efficient and powerful for solving fractional partial differential equations.


---

## Body

## 1. Introduction

During the last few decades, fractional calculus emerges as a natural description for a broad range of nonclassical phenomena in the applied sciences and engineering. For example, anomalous diffusion is one of the most important concepts in modern physics [1–3] and is presented in extremely diverse engineering fields such as vibration and acoustic dissipation in soft matter [4], polymer dynamics [5], turbulence flow [6], glassy and porous media [7], charges transport in amorphous semiconductor [8], chemistry and biochemistry [9], contaminant transport in groundwater flow [10], and chaotic dynamics of classical conservative systems [11].To date, many solution techniques for fractional differential equations have exploited the properties of the Fourier and Laplace transforms of the operators to confirm a classical solution. Typically, Gorenflo et al. in [12] used the method of Laplace transform to obtain the scale-invariant solution of the time fractional diffusion-wave equation according to the Wright function. Starting from the Fourier-Laplace representation, Mainardi et al. in [13] derived the explicit expression of the Green function (i.e., the fundamental solution) for the Cauchy problem with respect to its scaling and similarity properties. Since fractional order differential operators possess their own properties, (i) they are nonlocal operators and (ii) the adjoint operator itself of fractional differential operator is a nonnegative operator. Therefore, most fractional differential equations do not have exact solutions. Thus, some potent and accurate numerical methods for these equations are developed. For example, a finite difference method was applied to construct the numerical approximation for the time fractional diffusion equation and introduce an implicit difference scheme for the space-time fractional diffusion equation was studied in [14–17]. Ervin and Roop in [18] considered the stationary fractional advection dispersion equation by using the standard finite element method. Zhao et al. adopted the similar method to discuss the numerical solutions of the multiterm time-space Riesz fractional advection-diffusion equations in their paper [19] and Deng in [20] presented the finite element approximation for the space-time fractional Fokker-Planck equation. Lin and Xu in [21] and Li and Xu in [22] proposed a spectral method for time or space fractional diffusion equations based on a weak formulation and the detailed error analysis was carried out. A new matrix approach to solve fractional partial differential equations numerically has been developed by Podlubny in [23]. Li et al. in [24] used the fractional difference method in time and the finite element method in space, for the time-space fractional order subdiffusion and superdiffusion equations in which they presented the error estimates for the full discrete schemes. And for the subdiffusion linear problem, they obtain the approximation of order Δt2-α+Δt-αh2 in theoretical; for the superdiffusion linear problem, they got the approximation of order Δt2+Δt1-αh2. However, their numerical experiments cannot be consistent with the corresponding theoretical results very well, particularly for the superdiffusion problems.Moreover, high order accuracy numerical methods for solving fractional differential equations are still very limited. Discontinuous Galerkin method is famous for extreme flexibility and high-accuracy properties [25]. By applying this method, differential equations can be solved locally and accurately while we only need to assemble and solve a large linear system. Deng and Hesthaven in [26] studied the convergence rate analysis of local discontinuous Galerkin method (LDG) for solving spatial fractional diffusion equation, which was the first attempt to solve fractional derivative equations by using the LDG method; at the same time, the authors in this paper speculated that the LDG method has significant potential in solving fractional calculus problems. After that lots of researchers use the LDG method to analyze the time discretization of fractional differential equations [27, 28]. Recently, Xu and Hesthaven [29] obtained a consistent and high-accuracy numerical scheme for fractional convection-diffusion equations with a space fractional Laplacian operator of order α(1<α<2) through exploiting LDG method. They introduced a new technique by rewriting the fractional Laplacian operator as a composite of first-order derivatives and a fractional integral; then they converted the fractional convection-diffusion equation into a system of low order equations. Optimal order of convergence O(hk+1) for the space fractional diffusion problem and suboptimal order of convergence O(hk+1/2) for the general spatial fractional convection-diffusion problem are established, respectively.In this paper, motivated by the research work of [24], we intend to construct and analyze the numerical schemes for solving the time-space fractional subdiffusion equation as follows: (1)∂αux,t∂tα=ϵ--△β/2ux,t,x∈Ω,t∈0,T,ux,0=u0x,x∈Ω,and the time-space fractional superdiffusion equation as follows:(2)∂αux,t∂tα=ϵ--△β/2ux,t,x∈Ω,t∈0,T,ux,0=ϕx,x∈Ω,utx,0=ψx,x∈Ω,where Ω=(a,b), T>0, 1<β<2 denotes the order number of space fractional derivatives, 0<α<1 describes the fractional subdiffusion equation (1), 1<α<2 describes the fractional superdiffusion equation (2), ϵ is a real positive parameter, and u0,ϕ(x),ψ(x)∈H01(Ω) are the given smooth functions.For the sake of simplicity, we just consider the periodic boundary conditions in this paper. Notice that this assumption is not essential; our method can be directly extended to the problems with aperiodic boundary condition.The time fractional derivative in (1) uses the Caputo fractional partial derivative of order α, defined as [30, 31] (3)∂αux,t∂tα=1Γ1-α∫0t∂ux,s∂sdst-sα,0<α<1,∂ux,t∂t,α=1,and the time fractional Caputo derivative in (2) defined as [24](4)∂αux,t∂tα=1Γ2-α∫0tt-s1-α∂2ux,s∂s2ds,1<α<2,∂2ux,t∂t2,α=2,where Γ(·) is the Gamma function.The space fractional terms in (1) and (2) are defined through the fractional Laplacian operator (-(-△)β/2)u(x,t), 1<β<2, which is also known as the Riesz fractional derivative ∂β/∂|x|βu(x,t) under the assumption that ux,t,u′x,t,…,u(n-1)(x,t) vanishes at the end point of an infinite domain (i.e.,x=∓∞) [32].The spatial fractional differential operator is also an important tool to describe the anomalous diffusion phenomenon; when1<β<2, it represents a Lévy β-stable flight [33]; when β→2, it models a Brownian diffusion process.The rest of this paper is organized as follows. In the next section, we review the appropriate functional spaces and some important definitions and properties of fractional derivatives. By using the technique which is proposed in [29], that is, through converting the space fractional operator of order β into a composite of first-order derivative and a fractional integral of order 2-β, the fully discrete LDG schemes for the time-space fractional subdiffusion equation and the time-space fractional superdiffusion equation are constructed, respectively, in Section 3. Sections 4 and 5 provide the detailed stability and L2 error analysis for the two fully discrete schemes, respectively. Several examples are presented to test our theoretical results in Section 6. Section 7 contains the concluding remarks.


## 2. Definitions and Auxiliary Results

In this section, we present a few definitions and recall some properties of fractional calculus [34, 35].Definition 1.
For functionu(x,t) given in the domain Ω=(a,b) (a,b can be finite or infinite), the expressions(5)​aRIxαux,t=1Γα∫axx-ξα-1uξ,tdξ,x>a,​xRIbαux,t=1Γα∫xbξ-xα-1uξ,tdξ,x<bare called left- and right-side Riemann-Liouville fractional integrals of order α (0≤n-1<α<n, n being an integer), respectively. Γ(·) is a Gamma function.Definition 2.
The left and right of Riemann-Liouville fractional derivatives of orderα(n-1<α<n) in the domain Ω=(a,b) (a,b can be finite or infinite) can be derived by Riemann-Liouville fractional integrals as follows:(6)​aRDxαux,t=1Γn-αdndxn∫axx-ξn-α-1uξ,tdξ,x>a,(7)​xRDbαux,t=-1nΓn-αdndxn∫xbξ-xn-α-1uξ,tdξ,x<b,where Γ(·) is a Gamma function. When α is an integer, Dα=dα=dxα.Definition 3 (see [34]).
The fractional Riesz derivatives of orderα(n-1<α≤n) in the domain Ω=(a,b) (a,b can be finite or infinite) are defined as follows:(8)∂α∂xαux,t=-Cα​aDxα+​xDbαux,t,where Cα=1/2cos⁡πα/2,α≠1, aDxα, xDbα denote the left and the right Riemann-Liouville fractional derivatives of order α, respectively.Definition 4.
For a functionu(x,t) defined on the infinite domain (-∞<x<∞),t∈(0,T], the following equality holds [34]: (9)--△α/2ux,t=-12cos⁡πα/2​-∞Dxα+​xD∞αux,t=∂α∂xαux,t.Remark 5.
The definition of the fractional Laplacian operator uses the Fourier transform on an infinite domain, with a natural extension to include finite domains when the functionu(x,t) is subject to homogeneous Dirichlet boundary conditions. In this paper, the developments and analysis are based on this definition.Lemma 6.
Suppose that∂j/∂xu(x,t)=0 for x→∓∞, ∀0≤j≤n, (n-1<α<n), and then the relationship between Riemann-Liouville fractional integrals and Riemann-Liouville fractional derivatives satisfies [30](10)​aRDxαux,t=Dn​aIxn-αux,t=​aIxn-αDnux,t,​xRDbαux,t=-Dn​xIbn-αux,t=​xIbn-α-Dnux,t.In order to carry out the analysis of the fully discrete LDG scheme, we need the following properties of fractional integrals and derivatives and some lemmas which were collected from [29].Lemma 7 (linearity).
One has(11)​aIxαλfx+μgx=λ​aIxαfx+μ​aIxαgx,​aDxαλfx+μgx=λ​aDxαfx+μ​aDxαgx.Lemma 8 (semigroup property).
(12)
​
a
I
x
α
+
β
f
x
=
​
a
I
x
α
​
a
I
x
β
f
x
=
​
a
I
x
β
​
a
I
x
α
f
x
.By equalities (10), for any continuous function u(x,t) with limx→±∞uj(x,t)=0,(0≤j≤2) and 1<β<2, we can obtain the following result [29]: (13)--△β/2ux,t=-12cos⁡πβ/2∂2∂x2​-∞Ix2-βux,t+​xI∞2-βux,t=-12cos⁡πβ/2∂∂x​-∞Ix2-β∂ux,t∂x+​xI∞2-β∂ux,t∂x=∂∂x-​-∞Ix2-β+​xI∞2-β2cos⁡πβ/2∂ux,t∂x.In the above equation, if β<0, the fractional Laplacian operator becomes the fractional integral operator. Similar to [29], since 1<β<2, then -1<β-2<0, and we define (14)△β-2/2=-​-∞Ix2-β+​xI∞2-β2cos⁡πβ/2.Lemma 9.
When1<β<2, we can rewrite the space fractional Laplacian operator in a composite form of first-order derivative and a fractional integral of order 2-β as follows [29]:(15)--△β/2ux,t=∂∂x△β-2/2∂ux,t∂x.Definition 10.
We define the left- and right-fractional space norms as [18](16)uJLαR=uJLαR2+uL2R21/2,uJRαR=uJRαR2+uL2R21/2,where |u|JαR=DαuL2(R) is the seminorm of fractional space and JLα(R), JRα(R) refer to the closure of C0∞(R) with respect to ·JLα(R) and ·JRα(R), respectively.Lemma 11 (adjoint property).
One has(17)​-∞Ixαu,uR=u,​xI∞αuR.Lemma 12 (see [18]).
One has(18)​-∞Ixαu,​xI∞αuR=cos⁡απuJL-αR2=cos⁡απuJR-αR2.Lemma 13 (see [19]).
Let0<θ<2,θ≠1. Then, for any w,v∈H0θ/2(0,X), then (19)​0RDxθw,vL20,X=​0RDxθ/2w,​tRDXθ/2vL20,X.From the above definitions and Lemma8, together with the properties of Riemann-Liouville fractional integrals, the following lemmas hold.Lemma 14.
For any0<s<1, the fractional integral satisfies the following property: (20)△-su,uR=uJL-αR2=uJR-αR2.Proof.
Since0<s<1, by Lemmas 9, 11, 12, and 13, we obtain (21)△-su,uR=​-∞Ix2su+​xI∞2su2cos⁡sπ,uR=12cos⁡sπ​-∞Ix2su,uR+​xI∞2su,uR=12cos⁡sπ​-∞Ixsu,​xI∞suR+​xI∞su,​-∞IxsuR=12cos⁡sπ2cos⁡sπuJL-sR2=12cos⁡sπ2cos⁡sπuJR-sR2.Theorem 15 (see [26]).
If-α2<-α1<0, then JL,0-α1(Ω) (or JR,0-α1(Ω) or JS,0-α1(Ω)) is embedded into JL,0-α2(Ω) (or JR,0-α2(Ω) or JS,0-α2(Ω)), and L2(Ω) is embedded into both of them.Lemma 16 (fractional Poincaré-Friedrichs inequality [18]).
One has(22)uL2Ω≤CuJL,0μΩ,u∈JL,0μΩ,μ∈R.uL2Ω≤CuJR,0μΩ,u∈JR,0μΩ,μ∈R.Lemma 17.
Suppose the fractional integral is defined on[0,b] and let g(y)=f(b-y). Then (23)​xIbαfx=Iyα0gy.Proof.
One has(24)​xIbαfx=∫xbt-xα-1ftdt=t=b-ξ-∫b-x0b-ξ-xα-1fb-ξdξ=∫0b-xb-x-ξα-1fb-ξdξ=y=b-x∫0yy-ξα-1fb-ξdξ=fb-ξ=gξ∫0yy-ξα-1gξdξ=Iyα0gy.Lemma 18 (see [29]).
Supposeu(x) is a smooth function defined on Ω⊂R. Ωh is discretization of the domain with interval width h, and uh(x) is an approximation of u in Phk. For all j, uh(x)∈Ij is a polynomial of degree up to order k, and (u,v)Ij=(uh,v)Ij∀v∈Pk. k is the degree of the polynomial. Then for -1<α≤0, we obtain (25)△β/2ux-△β/2uhxL2Ω≤Chk+1,where C is a constant independent of h.Assume the following mesh to cover the computational domainΩ=[a,b], consisting of cells Ij=[xj-1/2,xj+1/2] for j=1,…,N, where a=x1/2<x3/2<⋯<xN+1/2=b, denoting the cell lengths △xj=xj+1/2-xj-1/2, 1≤j≤N, h=max1≤j≤N△xj. Define Vhk as the space of polynomials of the degree up to k in each cell Ij; that is,(26)Vhk=v:v∈PkIj,x∈Ij,j=1,2,…,N.To prepare for the analysis of the error estimates, we need two projections in one dimension[a,b], denoted by Q, that is, for each j,(27)∫IjQωx-ωxυxdx=0,∀υ∈PkIj,and special projection P±, that is, for each j,(28)∫IjP+ωx-ωxυxdx=0,∀υ∈Pk-1Ij,P+ωxj-1/2+=ωxj-1/2,∫IjP-ωx-ωxυxdx=0,∀υ∈Pk-1Ij,P-ωxj+1/2-=ωxj+1/2.The two projections satisfy the following approximation inequality [36]:(29)ωe+hωe∞+h1/2ωeτh≤Chk+1,where ωe=Pω(x)-ω(x), or ωe=P±ω(x)-ω(x). The positive constant C, solely depending on ω, is independent of h. τh denotes the set of boundary points of all elements Ij.In the present paper, we useC to denote a positive constant which may have a different value in each occurrence. The usual notation of norms in Sobolev spaces will be used. Let the scalar inner product on L2(D) be denoted by (·,·)D and the associated norm by ·D. If D=Ω, we drop D.


## 3. Numerical Schemes

In this section, we present the fully discrete local discontinuous Galerkin (LDG) schemes for solving the time-space fractional subdiffusion/superdiffusion problems, respectively.Letτ=T/M be the time mesh size, M is a positive integer, let tn=nτ,n=0,1,…,M be mesh point, and let tn-1/2=(tn+tn-1)/2, n=0,1,…,M be mid-mesh points. Here and below an unmarked norm ∥·∥ refers to the usual L2 norm.


### 3.1. Subdiffusion Case

In the rest of this section, we will follow the time discrete method of [24] to show the time discretization for the time fractional derivative ∂αu(x,t)/∂tα of subdiffusion equation.For the case of subdiffusion(0<α<1) equation (1), we adopt the backward Euler method to numerically approximate the time fractional derivative ∂αu(x,t)/∂tα at tn as follows [21, 27]:(30)∂αux,tn∂tα=1Γ1-α∑i=0n-1∫titi+1∂ux,s∂sdstn-sα=1Γ1-α∑i=0n-1ux,ti+1-ux,tiτ∫iτi+1τdstn-sα+rnx=τ1-αΓ2-α∑i=0n-1ux,tn-i-ux,tn-i-1τi+11-α-i1-α+rnx=τ1-αΓ2-α∑i=0n-1biux,tn-i-ux,tn-i-1τ+rnx.Here bi=(i+1)1-α-i1-α satisfy the following properties: (31)1=b0>b1>b2>⋯>bn>0,bn⟶0n⟶∞,∑i=1nbi-1-bi+bn=1.rn(x) is the truncation error and satisfies the following inequality [21]:(32)rnx≤CuΓ2-ατ2-α,where Cu is a constant only depending on u. We rewrite (1) as a first-order system(33)∂αux,t∂tα-ϵ∂∂xq=0,q-△β-2/2p=0,p-ϵ∂∂xux,t=0.And with (30), we give the weak form of system (33) at tn as follows:(34)∫Ωux,tnvdx+α0ϵ∫Ωqx,tnvxdx-∑j=1Nqx,tnvj+1/2--qx,tnvj-1/2+=∑i=1n-1bi-1-bi∫Ωux,tn-ivdx+bn-1∫Ωux,t0vdx-α0∫Ωrnxvdx,∫Ωqx,tnwdx-∫Ω△β-2/2px,tnwdx=0,∫Ωpx,tnzdx+ϵ∫Ωux,tnzxdx-∑j=1Nux,tnzj+1/2--ux,tnzj-1/2+=0,∫Ωux,t0vdx-∫Ωu0xvdx=0,where α0=ταΓ(2-α).Letuhn,phn,qhn∈Vhk be the approximation of u(·,tn),p(·,tn),q(·,tn), respectively. We define a fully discrete local discontinuous Galerkin (LDG) scheme as follows: Find uhn,phn,qhn∈Vhk, such that for all v,w,z∈Vhk(35)∫Ωuhnvdx+α0ϵ∫Ωqhnvxdx-∑j=1Nqhn^vj+1/2--qhn^vj-1/2+=∑i=1n-1bi-1-bi∫Ωuhn-ivdx+bn-1∫Ωuh0vdx,∫Ωqhnwdx-∫Ω△β-2/2phnwdx=0,∫Ωphnzdx+ϵ∫Ωuhnzxdx-∑j=1Nuhn^zj+1/2--uhn^zj-1/2+=0,∫Ωuh0vdx-∫Ωu0xvdx=0,where α0=ταΓ(2-α).The “hat” terms in (35) in the cell boundary terms from integration by parts are the so-called “numerical fluxes” [37], which are single valued functions defined on the edges and should be designed based on different guiding principles for different PDEs to ensure the stability. Here we take the simple choices such that(36)uhn^=uhn-,qhn^=qhn+.We remark that the choice for the fluxes (36) is not unique. In fact the crucial part is taking uhn and qhn from opposite sides. We know the truncation error is rn(x) from (30).


### 3.2. Superdiffusion Case

Similar to the subdiffusion case, we adopt the center difference scheme to estimate the time fractional derivative∂αu(x,t)/∂tα of superdiffusion case (2) at tn as follows [24]:(37)∂αux,tn∂tα=1Γ2-α∫0tn∂2ux,s∂s2dstn-sα-1=1Γ2-α∑i=0n-1∫titi+1∂2ux,s∂s2dstn-sα-1=1Γ2-α∑i=0n-1ux,ti+1-2ux,ti+ux,ti-1τ2∫iτi+1τtn-s1-αds+Rnx=τ2-αΓ3-α∑i=0n-1ux,tn-i-2ux,tn-i-1+ux,tn-i-2τ2i+12-α-i2-α+Rnx=τ2-αΓ3-α∑i=0n-1νiux,tn-i-ux,tn-i-1τ+Rnx,where νi=(i+1)2-α-i2-α,i=0,1,…,n-1, satisfy the following properties: (38)1=ν0>ν1>ν2>⋯>νn>0,νn⟶0n⟶∞,∑i=1nνi-1-νi+νn=1.Rn(x) is the truncation error; that is, (39)Rnx=1Γ2-α∑i=0n-1∫titi+1tn-s1-α∂2ux,s∂s2-∂2u∂t2x,ti-τ212∂4u∂t4x,ti+Oτ4ds=1Γ2-α∑i=0n-1∫titi+1tn-s1-αs-ti∂3u∂t3x,ti+Oτ2ds,which satisfies the following estimate:(40)Rnx≤CτΓ2-αmaxt0≤t≤tn∂3u∂t3x,ti∑i=0n-1∫titi+1tn-s1-αds+Oτ2≤CuT3-αΓ2-ατ+Oτ2,where Cu is the constant which only depends on u.Analogous to the subdiffusion equation, we also can rewrite (2) as the first-order system (33). From the time discretization equation (37), we obtain the weak formulation of (2) at tn:(41)∫Ωux,tnvdx+α1ϵ∫Ωqx,tnvxdx-∑j=1Nqx,tnvj+1/2--qx,tnvj-1/2+=∑i=1n-2-νi-1+2νi-νi+1∫Ωux,tn-i-1vdx+2ν0-ν1∫Ωux,tn-1vdx+2νn-1-νn-2∫Ωux,t0vdx-νn-1∫Ωux,t-1vdx-α1∫ΩRnxvdx,∫Ωqx,tnwdx-∫Ω△β-2/2px,tnwdx=0,∫Ωpx,tnzdx+ϵ∫Ωux,tnzxdx-∑j=1Nux,tnzj+1/2--ux,tnzj-1/2+=0,∫Ωux,t0vdx-∫Ωϕxvdx=0,∫Ωutx,t0vdx-∫Ωψxvdx=0,where α1=ταΓ(3-α).Letuhn,phn,qhn∈Vhk be the approximation of u(·,tn),p(·,tn),q(·,tn), respectively. We define a fully discrete local discontinuous Galerkin (LDG) scheme as follows: find uhn,phn,qhn∈Vhk, such that for all v,w,z∈Vhk(42)∫Ωuhnvdx+α1ϵ∫Ωqhnvxdx-∑j=1Nqhn^vj+1/2--qhn^vj-1/2+=∑i=1n-2-νi-1+2νi-νi+1∫Ωuhn-i-1vdx+2ν0-ν1∫Ωuhn-1vdx+2νn-1-νn-2∫Ωuh0vdx-νn-1∫Ωuh-1vdx,∫Ωqhnwdx-∫Ω△β-2/2phnwdx=0,∫Ωphnzdx+ϵ∫Ωuhnzxdx-∑j=1Nuhn^zj+1/2--uhn^zj-1/2+=0,∫Ωuh0vdx-∫Ωϕxvdx=0,∫Ωuh0tvdx-∫Ωψxvdx=0.Remarks 1.
Originally, we assume thatu has compact support and restrict the problem to the bounded domain Ω. As a consequence, we impose homogeneous Dirichlet boundary conditions for u∉Ω.
First, we define the fluxes at the boundary asu^N+(1/2)=u(b,t), q^N+(1/2)=qN+(1/2)-+μ/h[uN+(1/2)] for the right boundary and u^1/2=u(a,t), q^1/2=q1/2-+μ/h[u1/2] for the left boundary [29], where μ is a positive constant and [·] is the jump term operator of function u.
Next, we define the boundary term operatorL as follows:(43)Lv,w,z=ϵua,tz1/2+-ϵμhub,tvN+1/2--ϵub,tzN+1/2-=0.


## 3.1. Subdiffusion Case

In the rest of this section, we will follow the time discrete method of [24] to show the time discretization for the time fractional derivative ∂αu(x,t)/∂tα of subdiffusion equation.For the case of subdiffusion(0<α<1) equation (1), we adopt the backward Euler method to numerically approximate the time fractional derivative ∂αu(x,t)/∂tα at tn as follows [21, 27]:(30)∂αux,tn∂tα=1Γ1-α∑i=0n-1∫titi+1∂ux,s∂sdstn-sα=1Γ1-α∑i=0n-1ux,ti+1-ux,tiτ∫iτi+1τdstn-sα+rnx=τ1-αΓ2-α∑i=0n-1ux,tn-i-ux,tn-i-1τi+11-α-i1-α+rnx=τ1-αΓ2-α∑i=0n-1biux,tn-i-ux,tn-i-1τ+rnx.Here bi=(i+1)1-α-i1-α satisfy the following properties: (31)1=b0>b1>b2>⋯>bn>0,bn⟶0n⟶∞,∑i=1nbi-1-bi+bn=1.rn(x) is the truncation error and satisfies the following inequality [21]:(32)rnx≤CuΓ2-ατ2-α,where Cu is a constant only depending on u. We rewrite (1) as a first-order system(33)∂αux,t∂tα-ϵ∂∂xq=0,q-△β-2/2p=0,p-ϵ∂∂xux,t=0.And with (30), we give the weak form of system (33) at tn as follows:(34)∫Ωux,tnvdx+α0ϵ∫Ωqx,tnvxdx-∑j=1Nqx,tnvj+1/2--qx,tnvj-1/2+=∑i=1n-1bi-1-bi∫Ωux,tn-ivdx+bn-1∫Ωux,t0vdx-α0∫Ωrnxvdx,∫Ωqx,tnwdx-∫Ω△β-2/2px,tnwdx=0,∫Ωpx,tnzdx+ϵ∫Ωux,tnzxdx-∑j=1Nux,tnzj+1/2--ux,tnzj-1/2+=0,∫Ωux,t0vdx-∫Ωu0xvdx=0,where α0=ταΓ(2-α).Letuhn,phn,qhn∈Vhk be the approximation of u(·,tn),p(·,tn),q(·,tn), respectively. We define a fully discrete local discontinuous Galerkin (LDG) scheme as follows: Find uhn,phn,qhn∈Vhk, such that for all v,w,z∈Vhk(35)∫Ωuhnvdx+α0ϵ∫Ωqhnvxdx-∑j=1Nqhn^vj+1/2--qhn^vj-1/2+=∑i=1n-1bi-1-bi∫Ωuhn-ivdx+bn-1∫Ωuh0vdx,∫Ωqhnwdx-∫Ω△β-2/2phnwdx=0,∫Ωphnzdx+ϵ∫Ωuhnzxdx-∑j=1Nuhn^zj+1/2--uhn^zj-1/2+=0,∫Ωuh0vdx-∫Ωu0xvdx=0,where α0=ταΓ(2-α).The “hat” terms in (35) in the cell boundary terms from integration by parts are the so-called “numerical fluxes” [37], which are single valued functions defined on the edges and should be designed based on different guiding principles for different PDEs to ensure the stability. Here we take the simple choices such that(36)uhn^=uhn-,qhn^=qhn+.We remark that the choice for the fluxes (36) is not unique. In fact the crucial part is taking uhn and qhn from opposite sides. We know the truncation error is rn(x) from (30).


## 3.2. Superdiffusion Case

Similar to the subdiffusion case, we adopt the center difference scheme to estimate the time fractional derivative∂αu(x,t)/∂tα of superdiffusion case (2) at tn as follows [24]:(37)∂αux,tn∂tα=1Γ2-α∫0tn∂2ux,s∂s2dstn-sα-1=1Γ2-α∑i=0n-1∫titi+1∂2ux,s∂s2dstn-sα-1=1Γ2-α∑i=0n-1ux,ti+1-2ux,ti+ux,ti-1τ2∫iτi+1τtn-s1-αds+Rnx=τ2-αΓ3-α∑i=0n-1ux,tn-i-2ux,tn-i-1+ux,tn-i-2τ2i+12-α-i2-α+Rnx=τ2-αΓ3-α∑i=0n-1νiux,tn-i-ux,tn-i-1τ+Rnx,where νi=(i+1)2-α-i2-α,i=0,1,…,n-1, satisfy the following properties: (38)1=ν0>ν1>ν2>⋯>νn>0,νn⟶0n⟶∞,∑i=1nνi-1-νi+νn=1.Rn(x) is the truncation error; that is, (39)Rnx=1Γ2-α∑i=0n-1∫titi+1tn-s1-α∂2ux,s∂s2-∂2u∂t2x,ti-τ212∂4u∂t4x,ti+Oτ4ds=1Γ2-α∑i=0n-1∫titi+1tn-s1-αs-ti∂3u∂t3x,ti+Oτ2ds,which satisfies the following estimate:(40)Rnx≤CτΓ2-αmaxt0≤t≤tn∂3u∂t3x,ti∑i=0n-1∫titi+1tn-s1-αds+Oτ2≤CuT3-αΓ2-ατ+Oτ2,where Cu is the constant which only depends on u.Analogous to the subdiffusion equation, we also can rewrite (2) as the first-order system (33). From the time discretization equation (37), we obtain the weak formulation of (2) at tn:(41)∫Ωux,tnvdx+α1ϵ∫Ωqx,tnvxdx-∑j=1Nqx,tnvj+1/2--qx,tnvj-1/2+=∑i=1n-2-νi-1+2νi-νi+1∫Ωux,tn-i-1vdx+2ν0-ν1∫Ωux,tn-1vdx+2νn-1-νn-2∫Ωux,t0vdx-νn-1∫Ωux,t-1vdx-α1∫ΩRnxvdx,∫Ωqx,tnwdx-∫Ω△β-2/2px,tnwdx=0,∫Ωpx,tnzdx+ϵ∫Ωux,tnzxdx-∑j=1Nux,tnzj+1/2--ux,tnzj-1/2+=0,∫Ωux,t0vdx-∫Ωϕxvdx=0,∫Ωutx,t0vdx-∫Ωψxvdx=0,where α1=ταΓ(3-α).Letuhn,phn,qhn∈Vhk be the approximation of u(·,tn),p(·,tn),q(·,tn), respectively. We define a fully discrete local discontinuous Galerkin (LDG) scheme as follows: find uhn,phn,qhn∈Vhk, such that for all v,w,z∈Vhk(42)∫Ωuhnvdx+α1ϵ∫Ωqhnvxdx-∑j=1Nqhn^vj+1/2--qhn^vj-1/2+=∑i=1n-2-νi-1+2νi-νi+1∫Ωuhn-i-1vdx+2ν0-ν1∫Ωuhn-1vdx+2νn-1-νn-2∫Ωuh0vdx-νn-1∫Ωuh-1vdx,∫Ωqhnwdx-∫Ω△β-2/2phnwdx=0,∫Ωphnzdx+ϵ∫Ωuhnzxdx-∑j=1Nuhn^zj+1/2--uhn^zj-1/2+=0,∫Ωuh0vdx-∫Ωϕxvdx=0,∫Ωuh0tvdx-∫Ωψxvdx=0.Remarks 1.
Originally, we assume thatu has compact support and restrict the problem to the bounded domain Ω. As a consequence, we impose homogeneous Dirichlet boundary conditions for u∉Ω.
First, we define the fluxes at the boundary asu^N+(1/2)=u(b,t), q^N+(1/2)=qN+(1/2)-+μ/h[uN+(1/2)] for the right boundary and u^1/2=u(a,t), q^1/2=q1/2-+μ/h[u1/2] for the left boundary [29], where μ is a positive constant and [·] is the jump term operator of function u.
Next, we define the boundary term operatorL as follows:(43)Lv,w,z=ϵua,tz1/2+-ϵμhub,tvN+1/2--ϵub,tzN+1/2-=0.


## 4. Stability

In this subsection, we present the stability analysis for the fully discrete LDG schemes (35) and (42) correspondingly.Theorem 19.
For periodic or compactly supported boundary conditions, the fully discrete LDG scheme (35) is unconditionally stable, and the numerical solution uhn satisfies(44)uhn≤uh0.Using the fluxes (36) and the fluxes at the boundaries combining equality (43), after a simple rearrangement of terms, the fully LDG scheme (35) becomes(45)∫Ωuhnvdx+α0ϵ∫Ωqhnvxdx+ϵ∫Ωuhnzxdx-∫Ω△β-2/2phnwdx+∫Ωqhnwdx+∫Ωphnzdx+α0ϵ∑j=1N-1qhnj+1/2+vj+1/2+ϵ∑j=1N-1uhnj+1/2-zj+1/2+α0ϵqhn1/2+v1/2+-qhnN+1/2-vN+1/2-+ϵμhuhnN+1/2-2=∑i=1n-1bi-1-bi∫Ωuhn-ivdx+bn-1∫Ωuh0vdx.Now, we use the mathematical induction to proof Theorem19.Proof.
First, considern=1; the LDG scheme (35) is(46)∫Ωuh1vdx+α0ϵ∫Ωqh1vxdx+ϵ∫Ωuh1zxdx-∫Ω△β-2/2ph1wdx+∫Ωqh1wdx+∫Ωph1zdx+α0ϵ∑j=1N-1qh1j+1/2+vj+1/2+ϵ∑j=1N-1uh1j+1/2-zj+1/2+α0ϵqh11/2+v1/2+-α0ϵqh1N+1/2-vN+1/2-+ϵμhuh1N+1/2-2=∫Ωuh0vdx.Taking the test functions v=uh1,w=-α0ph1,z=α0qh1, we get(47)∫Ωuh1uh1dx+α0ϵ∫Ωqh1uh1xdx+α0ϵ∫Ωuh1qh1xdx+α0∫Ω△β-2/2ph1ph1dx-α0∫Ωqh1ph1dx+α0∫Ωph1qh1dx+α0ϵ∑j=1N-1qh1j+1/2+uh1j+1/2+α0ϵ∑j=1N-1uh1j+1/2-qh1j+1/2+α0ϵqh11/2+uh11/2+-α0ϵqh1N+1/2-uh1N+1/2-+α0ϵμhuh1N+1/2-2=∫Ωuh0uh1dx.Considering the integration by parts formulation: (qh1,(uh1)x)Ij+(uh1,(qh1)x)Ij=uh1qh1xj-(1/2)+xj+(1/2)-, we obtain the interface condition, (48)α0ϵ∫Ωqh1uh1xdx+α0ϵ∫Ωuh1qh1xdx+α0ϵ∑j=1N-1qh1j+1/2+uh1j+1/2+α0ϵ∑j=1N-1uh1j+1/2-qh1j+1/2=α0ϵ∑j=1Nuh1qh1xj-1/2+xj+1/2-+α0ϵ∑j=1N-1qh1j+1/2+uh1j+1/2+α0ϵ∑j=1N-1uh1j+1/2-qh1j+1/2=-α0ϵqh11/2+uh11/2++α0ϵqh1N+1/2-uh1N+1/2-.Thus, we obtain(49)uh12+α0△β-2/2ph1,ph1+α0ϵμhuh1N+1/2-2=∫Ωuh0uh1dx.From the fact that (50)∫Ωuh0uh1dx≤12uh12+12uh02,we have (51)uh12+α0△β-2/2ph1,ph1+α0ϵμhuh1N+1/2-2≤uh02.Combining the boundary condition and Lemma 14, we get(52)uh1≤uh0.
Now we assume that the following inequality holds:(53)uhm≤uh0,m=1,2,…,P.We need to prove uhP+1≤uh0. Let n=P+1 and take the test functions v=uhP+1,w=-α0phP+1,z=α0qhP+1 in scheme (35); we obtain (54)uhP+12+α0△β-2/2phP+1,phP+1+α0ϵμhuhP+1N+1/2-2=∑i=1Pbi-1-bi∫ΩuhP+1-iuhP+1dx+bP∫Ωuh0uhP+1dx≤∑i=1Pbi-1-biuhP+1-iuhP+1+bPuh0uhP+1≤∑i=1Pbi-1-biuh0uhP+1+bPuh0uhP+1=∑i=1Pbi-1-bi+bPuh0uhP+1≤12uh02+12uhP+12.Calling Lemma 14 again and considering the boundary condition, we have(55)uhP+1≤uh0.Theorem 20.
For a sufficiently small step size0<τ<T, the fully discrete LDG scheme (42) is unconditional stable and the numerical solution uhn satisfies(56)uhn≤Cϕx+τψx,n=0,1,…,N.Proof.
Similar to the subdiffusion case, using the inequalityu-1≤C(ϕ(x)+τψ(x)) and adopting the same methods and techniques as Theorem 19, we can obtain the conclusions of Theorem 20. For brevity, we ignore the details in proof process here.


## 5. Error Estimates

In this section, we discuss the error estimates for the fully discrete LDG schemes (35), (42) respectively.Theorem 21.
Letu(x,tn) be the exact solution of the problem (1), which is sufficiently smooth with bounded derivatives, and uhn be the numerical solution of the fully discrete LDG scheme (35); then the following error estimates hold:
When0<α<1,(57)ux,tn-uhn≤CuTα1-ατ-αhk+1+τ2-α+τ-α/2hk+1+hk+1.When α→1,(58)ux,tn-uhn≤CuTτ-1hk+1+τ+τ-1/2hk+1+hk+1,where Cu is a constant only depending on u.Proof.
Denote(59)eun=ux,tn-uhn=P-eun-P-ux,tn-ux,tn,epn=px,tn-phn=Qepn-Qpx,tn-px,tn,eqn=qx,tn-qhn=P+eqn-P+qx,tn-qx,tn.Subtracting equation (35) from equation (34), with the fluxes on boundaries, we obtain the error equation:(60)∫Ωeunvdx+α0ϵ∫Ωeqnvxdx-∑j=1Neqn+vj+1/2--eqn+vj-1/2+-∑i=1n-1bi-1-bi∫Ωeun-ivdx-bn-1∫Ωeu0vdx+α0∫Ωrnxvdx+∫Ωeqnwdx-∫Ω△β-2/2epnwdx+∫Ωepnzdx+ϵ∫Ωeunzxdx-∑j=1Neun-zj+1/2--eun-zj-1/2+=0.
Using equations (59) and (43), after a simple rearrangement of terms, the error equation (60) can be rewritten as(61)∫ΩP-eunvdx+α0ϵ∫ΩP+eqnvxdx+ϵ∫ΩP-eunzxdx-∫Ω△β-2/2Qepnwdx+∫ΩP+eqnwdx+∫ΩQepnzdx+α0ϵ∑j=1N-1P+eqnj+1/2+vj+1/2+ϵ∑j=1N-1P-eunj+1/2-zj+1/2+α0ϵP+eqn1/2+v1/2+-α0ϵP+eqnN+1/2-vN+1/2-+ϵμhP-eunN+1/2-vN+1/2-=∑i=1n-1bi-1-bi∫ΩP-eun-ivdx+bn-1∫ΩP-eu0vdx-α0∫Ωrnxvdx+∫ΩP-ux,tn-ux,tnvdx+α0ϵ∫ΩP+qx,tn-qx,tnvxdx+ϵ∫ΩP-ux,tn-ux,tnzxdx-∫Ω△β-2/2Qpx,tn-px,tnwdx+∫ΩP+qx,tn-qx,tnwdx+∫ΩQpx,tn-px,tnzdx-∑i=1n-1bi-1-bi∫ΩP-ux,tn-i-ux,tn-ivdx-bn-1∫ΩP-ux,t0-ux,t0vdx+α0ϵ∑j=1N-1P+qx,tn-qx,tnj+1/2+vj+1/2+ϵ∑j=1N-1P-ux,tn-ux,tnj+1/2-zj+1/2+α0ϵP+qx,tn-qx,tn1/2+v1/2+-α0ϵP+qx,tn-qx,tnN+1/2-vN+1/2-+ϵμhP-ux,tn-ux,tnN+1/2-vN+1/2-.Taking the test functions v=P-eun, w=-α0Qepn, z=α0P+eqn in (61), using the properties (27)-(28), then the following equality holds:(62)P-eun2+α0△β-2/2Qepn,Qepn+α0ϵμhP-eunN+1/22=∑i=1n-1bi-1-bi∫ΩP-eun-iP-eundx+bn-1∫ΩP-eu0P-eundx-α0∫ΩrnxP-eundx+∫ΩP-ux,tn-ux,tnP-eundx+α0ϵ∫ΩP+qx,tn-qx,tnP-eunxdx+α0ϵ∫ΩP-ux,tn-ux,tnP+eqnxdx+α0∫Ω△β-2/2Qpx,tn-px,tnQepndx-α0∫ΩP+qx,tn-qx,tnQepndx+α0∫ΩQpx,tn-px,tnP+eqndx-∑i=1n-1bi-1-bi∫ΩP-ux,tn-i-ux,tn-iP-eundx-bn-1∫ΩP-ux,t0-ux,t0P-eundx+α0ϵ∑j=1N-1P+qx,tn-qx,tnj+1/2+P-eunj+1/2+α0ϵ∑j=1N-1P-ux,tn-ux,tnj+1/2-P+eqnj+1/2+α0ϵP+qx,tn-qx,tn1/2+P-eun1/2+-α0ϵP+qx,tn-qx,tnN+1/2-P-eunN+1/2-+ϵμhP-ux,tn-ux,tnN+1/2-P-eun-N+1/2.That is,(63)P-eun2+α0△β-2/2Qepn,Qepn+α0ϵμhP-eunN+1/22=∑i=1n-1bi-1-bi∫ΩP-eun-iP-eundx+bn-1∫ΩP-eu0P-eundx-α0∫ΩrnxP-eundx+∫ΩP-ux,tn-ux,tnP-eundx+α0△β-2/2Qpx,tn-px,tn-P+qx,tn-qx,tn,Qepn-α0ϵP+qx,tn-qx,tn-N+1/2P-eunN+1/2--∑i=1n-1bi-1-bi∫ΩP-ux,tn-i-ux,tn-iP-eundx-bn-1∫ΩP-ux,t0-ux,t0P-eundx.
We prove Theorem21 by mathematical induction.
First, we consider the case whenn=1, (63) becomes(64)P-eu12+α0△β-2/2Qep1,Qep1+α0ϵμhP-eu1N+1/22=∫ΩP-eu0P-eu1dx-α0∫Ωr1xP-eu1dx+∫ΩP-ux,t1-ux,t1P-eu1dx+α0△β-2/2Qpx,t1-px,t1-P+qx,t1-qx,t1,Qep1-α0ϵP+qx,t1-qx,t1-N+1/2P-eu1N+1/2--∫ΩP-ux,t0-ux,t0P-eu1dx.Notice the fact that P-eu0≤Chk+1,ab≤ϵa2+1/4ϵb2.
Together with property (29) and Lemma 18, we can get (65)△β-2/2Qpx,t1-px,t1-P+qx,t1-qx,t1=△β-2/2Qpx,t1-px,t1-△β-2/2P+px,t1-px,t1≤Chk+1.Combining this with Definition 4 and (59), using Young’s inequality we obtain(66)P-eu12+α0△β-2/2Qep1,Qep1+α0ϵμhP-eu1N+1/22≤P-eu0+α0r1x+P-ux,t1-ux,t1+P-ux,t0-ux,t0P-eu1+Ch2k+2τα+α0ϵCQep12+α0ϵμhP-eu1N+1/22.Using Young’s inequality again in (66), we get(67)P-eu12+α0△β-2/2Qep1,Qep1≤Chk+1+τ22+ϵP-eu12+Ch2k+2τα+α0ϵCQep12.If we choose ϵ≤min⁡1,1/C and recall the fractional Poincaré-Friedrichs Lemma 16, we have(68)P-eu1≤b0-1Chk+1+τ2+τα/2hk+1,where C is a constant depending on u,T,α.
Next, we suppose that the following inequality holds:(69)P-eum≤bm-1-1Chk+1+τ2+τα/2hk+1,m=1,2,…,l.When n=l+1, from (63) and Young’s inequality, we can obtain (70)P-eul+12+α0△β-2/2Qepl+1,Qepl+1+α0ϵμhP-eul+1N+1/22=∑i=1lbi-1-bi∫ΩP-eul+1-iP-eul+1dx+bl∫ΩP-eu0P-eul+1dx-α0∫Ωrl+1xP-eul+1dx+∫ΩP-ux,tl+1-ux,tl+1P-eul+1dx-α0ϵP+qx,tl+1-qx,tl+1-N+1/2P-eul+1N+1/2-+α0△β-2/2Qpx,tl+1-px,tl+1-P+qx,tl+1-qx,tl+1,Qepl+1-∑i=1lbi-1-bi∫ΩP-ux,tl+1-i-ux,tl+1-iP-eul+1dx-bl∫ΩP-ux,t0-ux,t0P-eul+1dx.≤∑i=1lbi-1-biP-eul+1-i+blP-eu0+α0rl+1x+P-ux,tl+1-ux,tl+1+∑i=1lbi-1-biP-ux,tl+1-i-ux,tl+1-i+blP-ux,t0-ux,t0×P-eul+1+Ch2k+2τα+α0ϵCQepl+12+α0ϵμhP-eul+1N+1/22≤∑i=1lbi-1-bibl-i-1Chk+1+τ2+τα/2hk+1+blP-eu0+α0rl+1x+P-ux,tl+1-ux,tl+1+∑i=1lbi-1-bi+blChk+1×P-eul+1+Ch2k+2τα+α0ϵCQepl+12+α0ϵμhP-eul+1N+1/22≤∑i=1lbi-1-bibl-i-1Chk+1+τ2+τα/2hk+1+Chk+1+τ2P-eul+1+Ch2k+2τα+α0ϵCQepl+12+α0ϵμhP-eul+1N+1/22.Since(71)bi-1-1<bi-1,τα/2hk+1>0,then(72)P-eul+12+α0△β-2/2Qepl+1,Qepl+1+α0ϵμhP-eul+1N+1/22≤∑i=1lbi-1-bibl-1Chk+1+τ2+τα/2hk+1+blbl-1Chk+1+τ2+τα/2hk+1×P-eul+1+Ch2k+2τα+α0ϵCQepl+12+α0ϵμhP-eul+1N+1/22≤∑i=1lbi-1-bi+blbl-1Chk+1+τ2+τα/2hk+1P-eul+1+Ch2k+2τα+α0ϵCQepl+12+α0ϵμhP-eul+1N+1/22.Using Definition (7) into above inequalities and combining Young’s inequalities, we have(73)P-eul+12+α0△β-2/2Qepl+1,Qepl+1≤bl-1Chk+1+τ2+τα/2hk+12+ϵP-eul+12+Ch2k+2τα+α0ϵCQepl+12.Analogously, if we choose ϵ small enough and recalling the fractional Poincaré-Friedrichs Lemma 16 again, we obtain(74)P-eul+1≤bl-1Chk+1+τ2+τα/2hk+1,where C is a constant related to u,T,α.
According to the principle of mathematical induction, we get(75)P-eun≤bn-1-1Chk+1+τ2+τα/2hk+1.
By [21, 27], we know that n-αbn-1-1 is monotonly increasing and tends to 1/1-α,n→∞; thus, (76)P-eun≤bn-1-1Chk+1+τ2+τα/2hk+1≤nαn-αbn-1-1Chk+1+τ2+τα/2hk+1≤Tτα11-αChk+1+τ2+τα/2hk+1≤CuTα1-ατ-αhk+1+τ2-α+τ-α/2hk+1,where Cu is a constant only depending on u.
Whenα→1, 1/1-α→∞, it is meaningless for (76), so as for the case of α→1, we should give out the estimate again.
We prove the following estimate by using the mathematical induction:(77)P-eun≤nChk+1+τ2+τα/2hk+1.Obviously, when n=1, (77) is workable clearly.
Whenn=l+1, by (63), we obtain (78)P-eul+12+α0△β-2/2Qepl+1,Qepl+1+α0ϵμhP-eul+1N+1/22≤∑i=1lbi-1-biP-eul+1-i+blP-eu0+α0rl+1x+P-ux,tl+1-ux,tl+1+∑i=1lbi-1-biP-ux,tl+1-i-ux,tl+1-i+blP-ux,t0-ux,t0×P-eul+1+Ch2k+2τα+α0ϵCQepl+12+α0ϵμhP-eul+1N+1/22≤∑i=1lbi-1-bil+1-iChk+1+τ2+τα/2hk+1+Chk+1+τ2P-eul+1+Ch2k+2τα+α0ϵCQepl+12+α0ϵμhP-eul+1N+1/22≤∑i=1lbi-1-bil+1-i+1l+1l+1Chk+1+τ2+τα/2hk+1P-eul+1+Ch2k+2τα+α0ϵCQepl+12+α0ϵμhP-eul+1N+1/22.Using Young’s inequality into the last inequality above, and choosing ϵ≤min⁡1,1/C, then combining the property of bi and Lemma 16 (Poincaré-Friedrichs inequality), we obtain (79)P-eul+12≤∑i=1lbi-1-bil+1-i+1l+1l+1Chk+1+τ2+τα/2hk+12+Ch2k+2τα.P-eul+1≤∑i=1lbi-1-bil+1-i+1l+1l+1Chk+1+τ2+τα/2hk+1≤∑i=1lbi-1-bil+1l+1l+1Chk+1+τ2+τα/2hk+1≤l+1Chk+1+τ2+τα/2hk+1.By mathematical induction, we get inequality (77).
Sincenτ≤T, n=0,1,…,M, that is, n≤T/τ, hence, when α→1, inequality (77) becomes(80)P-eun≤nChk+1+τ2+τα/2hk+1≤TτChk+1+τ2+τα/2hk+1≤TCuτ-1hk+1+τ+τ-1/2hk+1,where Cu is a constant only depending on u.
Combined with inequalities (76) and (80), and according to the triangle inequality and the standard approximation property (29), coupling with the error associated with the projection error, we can prove the results of Theorem 21.Theorem 22.
Letu(x,tn) be the exact solution of problem (2), which is sufficiently smooth with bounded derivatives, and uhn be the numerical solution of the fully discrete LDG scheme (42); then the following error estimates hold.
When1<α<2,(81)ux,tn-uhn≤CuTα-12-ατ1-αhk+1+τ2+τ1-α/2hk+1+hk+1;when α→2,(82)ux,tn-uhn≤CuTτ-1hk+1+τ2+hk+1,where Cu is a constant only depending on u.Proof.
Subtracting equation (42) from equation (41), and combining the fluxes on boundaries, we obtain the error equation as follows:(83)∫Ωeunvdx+α1ϵ∫Ωeqnvxdx-∑j=1Neqn+vj+1/2--eqn+vj-1/2+-∑i=1n-2-νi-1+2νi-νi+1∫Ωeun-i-1vdx-2ν0-ν1∫Ωeun-1vdx-2νn-1-νn-2∫Ωeu0vdx+νn-1∫Ωeu-1vdx+α1∫ΩRnxvdx+∫Ωeqnwdx-∫Ω△β-2/2epnwdx+∫Ωepnzdx+ϵ∫Ωeunzxdx-∑j=1Neun-zj+1/2--eun-zj-1/2+=0.
Similar to the error estimate of subdiffusion equation, first we can prove the following error inequality:(84)eu-1≤Chk+1+τ2+τhk+1,where C is a constant related to u,T,α.
In fact,(85)eu-1=ux,t-1-uh-1=ux,t-1-u-1+u-1-uh-1.ux,t-1-u-1=ux,t-1-ux,t0+τutx,t0=τ22!uttx,t0+Oτ3,u-1-uh-1=ux,t0-τutx,t0-ux,t0-τutx,t0h=ux,t0-uh0-τutx,t0-uth0≤Chk+1+τhk+1.Substituting the final results of the last two equalities into the first equality, we get the error inequality (84).
Next, we adopt the same methods and techniques as Theorem21 to obtain (86)P-eun≤νn-1-1Chk+1+τ1+α+τα/2hk+1.
Now we consider the function,(87)fx=x1-αx2-α-x-12-α,x≥1.Differentiating the function f(x) with respect to x, we have(88)f′x=1-αx-αx2-α-x-12-α-2-αx1-αx1-α-x-11-αx2-α-x-12-α2.When x>1, the denominator of equality (88) is strictly greater than zero, so we only need to consider the condition of the molecular, because (89)1-αx-αx2-α-x-12-α-2-αx1-αx1-α-x-11-α=1-αx2-2α-1-αx-αx-12-α-2-αx2-2α+2-αx1-αx-11-α=-x2-2α+x-αx-12-α-2-αx-αx-12-α+2-αx1-αx-11-α=-x2-2α+x-αx-12-α+2-αx-αx-11-αx-x-1=-x2-2α+x-αx-12-α+2-αx-αx-11-α=x-αx-12-α-x2-α+2-αx-αx-11-α=x-αx-12-α-x2-α+x-11-α+1-αx-αx-11-α=x1-αx-11-α-x1-α+1-αx-αx-11-α=xxαx-1x-1α-xxα+1-α1xαx-1x-1α=1x2αx-1αxx-1xα-x2x-1α+1-αx-1xα.When 1<α<2 and x>1, always x2α(x-1)α>0, so we only need to estimate the positive or negative property of the molecular of the last equality above. (90)xx-1xα-x2x-1α+1-αx-1xα=x-1xαx+1-α-x2x-1α=x-1xαx+1-α-x2-αx-1α-1.For x>1, (x-1)xα>0 always holds.
For any givenx, we define the function as follows:(91)gα=x+1-α-x2-αx-1α-1.Differentiating the function g(α) with respect to α, we obtain(92)g′α=-1-x2-αx-1α-1lnx-lnx-1.Let g′(α)=0; there is only one zero root of the derived function g′(α) on interval (1,2): (93)α=logx-1/x1-x/x2lnx/x-1.Since g(1)=g(2)=0, for any given x, the function g(α) with respect to α can only increase first and then decrease or decrease first and then increase on interval (1,2).
Now if we can prove that the value of any point on interval(1,2) is positive, then the function g(α) must be increasing first and then decreasing. For example, we choose the parameter α=3/2 and estimate the positive or negative property of the following equality: g(3/2)=(x-1/2)-x1/2(x-1)1/2. Since x-1/2/x1/2(x-1)1/2=[(x-1/2)2/x(x-1)]1/2>1, we get g(3/2)>0; hence, the function g(α) with respect to α can only increase first and then decrease on interval (1,2). That is, for any given x, when 1<α<2, we have g(α)>0. When x=1, (x-1)2xα-x2(x-1)α+(1-α)(x-1)xα=0.
From the above discussion we obtainf′(x)≥0,x∈[1,+∞); that is, the function f(x) is a monotone increasing function on interval [1,+∞).
Considering the function limitlimx→∞f(x), by Hospital’s rule, we have (94)limx→∞⁡fx=limx→∞⁡x1-αx2-α-x-12-α=limx→∞⁡1-αx-α2-αx1-α-x-11-α=limx→∞⁡1-α2-αx-11-x/x-1α-1=limx→∞⁡12-αx-2x-1-2x/x-1α-2=limx→∞⁡12-αxx-1-2xx-12-α=limx→∞⁡12-αxx-1-α=limx→∞⁡12-α1-1xα=12-α.Hence, (95)fx≤12-α,x≥1.Therefore, the sequence n1-ανn-1-1 monotone increasing tends to 1/2-α,n→∞. We obtain(96)P-eun≤νn-1-1Chk+1+τ1+α+τα/2hk+1≤nα-1n1-ανn-1-1Chk+1+τ1+α+τα/2hk+1≤Tτα-112-αChk+1+τ1+α+τα/2hk+1≤CuTα-12-ατ1-αhk+1+τ2+τ1-α/2hk+1,where Cu is a constant only depending on u.
Whenα→2, 1/2-α→∞, estimate (96) is meaningless, so for the case of α→2, we need to estimate again.
Analogous to the subdiffusion case, we can prove the following estimate still using the mathematical induction:(97)P-eun≤nChk+1+τ1+α+τα/2hk+1.Because nτ≤T, n=0,1,…,M, that is, n≤T/τ, as a result, when α→2, the inequality (97) becomes(98)P-eun≤nChk+1+τ1+α+τα/2hk+1≤TτChk+1+τ1+α+τα/2hk+1≤TCuτ-1hk+1+τ2+hk+1,where Cu is a constant only depending on u.
Combining inequality (96) with inequality (98), according to the triangle inequality and standard approximation Theorem (29), together with the property of projection operator, we obtain the results of Theorem 22.


## 6. Numerical Examples

In this section, we present some numerical examples to demonstrate the effectiveness of the finite difference/LDG approximation for solving the time-space fractional subdiffusion/superdiffusion equations. We check the convergence behavior of numerical solutions with respect to the time step sizeτ and the space step size h. By choosing the appropriate time step size, we avoid the contamination of the temporal error compared with the spatial error.


### 6.1. Implementation

As examples, we consider the time-space fractional subdiffusion/superdiffusion equations as follows:(99)∂αux,t∂tα-ϵ--△β/2ux,t=fx,t,inx∈Ω,t∈0,T,ux,0=u0x,inx∈Ω,where Ω=[a,b], T>0, 0<α<1 for the subdiffusion case (1), or 1<α<2 for the superdiffusion case (2), and 1<β<2, ϵ>0.For brevity, we only derive the linear system obtained from the fully discrete LDG scheme (35). For scheme (42), we can obtain the corresponding linear system similarly. We consider the case of discontinuous piecewise polynomials basis function sequences {ϕi(x)}i=1k in space Vkh. By the definition of the space Vkh, we know that, for all v∈Vkh, v(x)=∑i=1N‍viϕi(x) holds, where vi=v(xi).By the LDG scheme (35), we obtain the fully discrete scheme of (99) as follows:(100)∫Ωuhnvdx+α0ϵ∫Ωqhnvxdx-∑j=1Nqhn^vj+1/2--qhn^vj-1/2+=∑i=1n-1bi-1-bi∫Ωuhn-ivdx+bn-1∫Ωuh0vdx+∫Ωfx,tvdx,∫Ωqhnwdx-∫Ω△β-2/2phnwdx=0,∫Ωphnzdx+ϵ∫Ωuhnzxdx-∑j=1Nuhn^zj+1/2--uhn^zj-1/2+=0,∫Ωuh0vdx-∫Ωu0xvdx=0,where α0=ταΓ(2-α).In terms of the basis{ϕi(x)}i=1N, we denote approximate solution functions by uhn=∑i=1N‍uj(t)ϕj(x), phn=∑i=1N‍pj(t)ϕj(x), and qhn=∑i=1N‍qj(t)ϕj(x) and let test functions v=∑i=1N‍viϕi(x),vi=v(xi), w=∑i=1N‍wiϕi(x),wi=w(xi), and z=∑i=1N‍ziϕi(x),zi=z(xi). Putting these expressions into the LDG scheme (100) and taking the flux in (36), we get the linear system of the scheme in jth element at tn as follows:(101)M0α0ϵA-H11+H120RMϵA-H13+H14M0UjnPjnQjn=Fnv+∑i=1n-1bi-1-biUn-iv+bn-1U0v,where M={(ϕmj(x),ϕkj(x))}j=1,…,N is the mass matrix, m,k=0,1,… denote the polynomials order, and (102)A=ϕmjx,ϕkj′xj=1,…,N,H11=ϕmjxj-1/2,ϕkjxj+1/2j=1,…,N,H12=ϕmjxj-1/2,ϕkjxj-1/2j=1,…,N,H13=ϕmjxj+1/2,ϕkjxj+1/2j=1,…,N,H14=ϕmj-1xj-1/2,ϕkjxj-1/2j=1,…,N.R is the quadrature of the fractional integral operator △(β-2)/2phn and test function w in jth element. From [29], we know that (103)△β-2/2phn=-​aIx2-βphn+xIb2-βphn2cos⁡βπ/2,where(104)​aIx2-βphnx,t=1Γ2-β∫axx-ξ1-βphnξ,tdξ,​xIb2-βphnx,t=1Γ2-β∫xbξ-x1-βphnξ,tdξ,and Fn=(f1,f2,…,fn)Tis a vector valued function, fj(tn)=f(xj,tn). Then, by solving the linear system (101), we can obtain the solution Un.


### 6.2. Numerical Results

Example 1.
Consider the fractional time-space subdiffusion equation (99) in domain [0,1]×[0,0.5] which satisfies the following initial boundary conditions: (105)ux,0=x21-x2,x∈0,1.In order to obtain the exact solution u(x,t)=(2t-1)2x2(1-x)2, we choose the source term as (106)fx,t=x21-x28t2-αΓ3-α-4t1-αΓ2-α+t-αΓ1-α+2t-122cos⁡πβ+1/224x3-βΓ4-β-12x2-βΓ3-β+2x1-βΓ2-β+241-x3-βΓ4-β-121-x2-βΓ3-β+21-x1-βΓ2-β.Choosing the fluxes in (36) and adopting the linear piecewise polynomials, (99) is solved numerically. We investigate the numerical solutions with respect to different fractional order parameters α and β; we select the appropriate time step size τ such that the time discretization errors are negligible as compared with the spatial errors. Then choose space step size sequence as h=1/2i,(i=2,…,6); Tables 1 and 2 list out the L2 errors and the corresponding convergence order, respectively, fixing the value of the one parameter β (or α), with the values change of the other one parameter α (or β). As can be seen from the data in the two tables, our schemes can achieve the optimal convergence O(hk+1) in spatial direction; this matches the theoretical results of Theorem 21.Table 1
Example1. The L2 errors and convergence rate in space direction. β=1.6, T=0.5.
h
α=0.1
Rate
α=0.5
rate
α=0.9
Rate
1/4
3.3592×10-2
—
2.7826
×
1
0
-
2
—
1.4821×10-2
—
1/8
9.1511
×
1
0
-
3
1.8761
6.9632
×
1
0
-
3
1.9986
3.7341×10-3
1.9888
1/16
2.3499
×
1
0
-
3
1.9613
1.7411
×
1
0
-
3
1.9997
9.3353×10-4
2.0000
1/32
5.9437×10-4
1.9832
4.3333
×
1
0
-
4
2.0065
2.3241×10-4
2.0060
1/64
1.4875
×
1
0
-
4
1.9984
1.0748
×
1
0
-
4
2.0113
5.7422×10-5
2.0170Table 2
Example1. The L2 errors and convergence rate in space direction. α=0.6, T=0.5.
h
β=1.1
Rate
β=1.5
rate
β=1.9
Rate
1/4
1.4938
×
1
0
-
2
—
5.5316
×
1
0
-
2
—
4.6063
×
1
0
-
2
—
1/8
4.0147
×
1
0
-
3
1.8956
1.3845
×
1
0
-
2
1.9983
1.1679
×
1
0
-
2
1.9796
1/16
1.0400
×
1
0
-
3
1.9487
3.4615
×
1
0
-
3
1.9999
2.9465
×
1
0
-
3
1.9869
1/32
2.6052
×
1
0
-
4
1.9971
8.6539
×
1
0
-
4
2.0000
7.5009
×
1
0
-
4
1.9739
1/64
6.5141
×
1
0
-
5
1.9998
2.1589
×
1
0
-
4
2.0030
1.9156
×
1
0
-
4
1.9692The numerical example results listed in Tables3 and 4 show that when we choose the appropriate space step size h and the time step size sequence τ=1/2i,(i=2,…,6), fixing β (or α), and changing α (or β), the temporal accuracy of the time-space fractional subdiffusion equation under the L2 norm can converge to O(τ2-α) for 0<α<1, and close to 1 when α→1, the results agree with the theoretical computation results in Theorem 21.Table 3
Example1. The L2 errors and convergence rate in time direction. β=1.6, T=0.5.
τ
α=0.2
Rate
α=0.6
Rate
α
=
0.999
Rate
1/4
4.8064
×
1
0
-
3
—
5.5078
×
1
0
-
3
—
5.4127
×
1
0
-
3
—
1/8
1.4486
×
1
0
-
3
1.7302
2.1182
×
1
0
-
3
1.3786
2.8977
×
1
0
-
3
0.9014
1/16
4.2962
×
1
0
-
4
1.7536
8.0864
×
1
0
-
4
1.3893
1.5246
×
1
0
-
3
0.9265
1/32
1.2320
×
1
0
-
4
1.8020
3.0576
×
1
0
-
4
1.4031
7.6910
×
1
0
-
4
0.9872
1/64
3.5291
×
1
0
-
5
1.8037
1.1501
×
1
0
-
4
1.4106
3.8372
×
1
0
-
4
1.0031Table 4
Example1. The L2 errors and convergence rate in time direction. α=0.5, T=0.5.
τ
β=1.05
Rate
β=1.5
Rate
β=1.99
Rate
1/4
6.7348
×
1
0
-
3
—
6.8961
×
1
0
-
3
—
5.9468
×
1
0
-
3
—
1/8
2.5767
×
1
0
-
3
1.3861
2.6180
×
1
0
-
3
1.3973
2.2038
×
1
0
-
3
1.4321
1/16
9.6141
×
1
0
-
4
1.4223
9.8328
×
1
0
-
4
1.4128
7.9725
×
1
0
-
4
1.4669
1/32
3.4618
×
1
0
-
4
1.4736
3.6092
×
1
0
-
4
1.4459
2.8224
×
1
0
-
4
1.4981
1/64
1.2263
×
1
0
-
4
1.4972
1.2211
×
1
0
-
4
1.4943
9.9691
×
1
0
-
5
1.5014Example 2.
Consider the time-space fractional subdiffusion(0<α<1) equation (99); for t∈[0,1],x∈[0,1], we choose the exact solution as u(x,t)=sin⁡(2πt)x2(1-x)2 and ϵ=0.05.
Similar to Example1, we give out the L2 errors and convergence rate with fixing α (or β) and changing β or α. We also select the appropriate time step size τ, which guarantee that the time discretization errors are negligible as compared with the spatial errors. Choosing the space step size sequence as h=1/2i,(i=2,…,6), we obtain the optimal convergence rate in space direction (as shown in Tables 5 and 6, this corresponds to the theoretical results as illustrated in Theorem 21). Tables 7 and 8 listed out the temporal accuracy; here we choose the appropriate space step size h and the time step size sequence as τ=1/2i,(i=2,…,6). The data in Table 7 display that, for 0<α<1, the convergence order in time direction can be up to O(τ2-α) and close to 1 order for α→1, which are consistent with Theorem 21.Table 5
Example2. The L2 errors and convergence rate in space direction. β=1.6, T=1.
h
α=0.2
Rate
α=0.6
Rate
α=0.99
Rate
1/4
3.1592
×
1
0
-
2
—
1.9328
×
1
0
-
2
—
2.1532
×
1
0
-
2
—
1/8
9.5228
×
1
0
-
3
1.7301
4.8649
×
1
0
-
3
1.9902
5.4460
×
1
0
-
3
1.9832
1/16
2.4728
×
1
0
-
3
1.9452
1.2174
×
1
0
-
3
1.9985
1.3693
×
1
0
-
3
1.9917
1/32
1.2521
×
1
0
-
3
1.9818
3.0317
×
1
0
-
4
2.0057
3.3469
×
1
0
-
4
2.0326
1/64
3.1377
×
1
0
-
4
1.9966
7.4936
×
1
0
-
5
2.0164
8.3644
×
1
0
-
5
2.0005Table 6
Example2. The L2 errors and convergence rate in space direction. α=0.6, T=1.
h
β
=
1.05
Rate
β
=
1.5
Rate
β
=
1.98
Rate
1/4
3.2624
×
1
0
-
2
—
1.6352
×
1
0
-
2
—
2.3573
×
1
0
-
2
—
1/8
8.8793
×
1
0
-
3
1.8782
4.1930
×
1
0
-
3
1.9961
5.9182
×
1
0
-
3
1.9939
1/16
2.3194
×
1
0
-
3
1.9367
1.0570
×
1
0
-
3
1.9879
1.4854
×
1
0
-
3
1.9943
1/32
5.8251
×
1
0
-
4
1.9934
2.6453
×
1
0
-
4
1.9986
3.7137
×
1
0
-
4
1.9999
1/64
1.4571
×
1
0
-
4
1.9991
6.5808
×
1
0
-
5
2.0071
9.2703
×
1
0
-
5
2.0022Table 7
Example2. The L2 errors and convergence rate in time direction. β=1.6, T=1.
τ
α=0.2
Rate
α=0.6
Rate
α=0.999
Rate
1/4
6.3463
×
1
0
-
3
—
5.9462
×
1
0
-
3
—
5.8016
×
1
0
-
3
—
1/8
1.9707
×
1
0
-
3
1.6872
2.4095
×
1
0
-
3
1.3032
3.3481
×
1
0
-
3
0.7931
1/16
6.0938
×
1
0
-
4
1.6933
9.3755
×
1
0
-
4
1.3618
1.8392
×
1
0
-
3
0.8642
1/32
1.7484
×
1
0
-
4
1.8013
3.5496
×
1
0
-
4
1.4012
9.2200
×
1
0
-
4
0.9963
1/64
5.0206
×
1
0
-
5
1.8001
1.3445
×
1
0
-
4
1.4006
4.4695
×
1
0
-
4
0.9877Table 8
Example2. The L2 errors and convergence rate in time direction. α=0.5, T=1.
τ
β
=
1.05
Rate
β
=
1.5
Rate
β
=
1.99
Rate
1/4
6.7819
×
1
0
-
3
—
6.2487
×
1
0
-
3
—
5.1918
×
1
0
-
3
—
1/8
2.6331
×
1
0
-
3
1.3649
2.3006
×
1
0
-
3
1.4415
1.9810
×
1
0
-
3
1.3900
1/16
9.9694
×
1
0
-
4
1.4012
8.2688
×
1
0
-
4
1.4763
7.2403
×
1
0
-
4
1.4521
1/32
3.5724
×
1
0
-
4
1.4806
2.9434
×
1
0
-
4
1.4902
2.5607
×
1
0
-
4
1.4995
1/64
1.2685
×
1
0
-
4
1.4937
1.0393
×
1
0
-
4
1.5018
9.1235
×
1
0
-
5
1.4889Example 3.
In this example, we consider the case for1<α<2; that is, (99) is a time-space fractional superdiffusion model. We assume that the equation satisfies the following initial boundary conditions in the domain [0,1]×[0,0.5]:(107)ux,0=x21-x2,x∈0,1,utx,0=-4x21-x2,x∈0,1.Here we choose ϵ=1, and the source term is (108)fx,t=x21-x28t2-αΓ3-α+t-αΓ1-α+2t-122cos⁡πβ+1/224x3-βΓ4-β-12x2-βΓ3-β+2x1-βΓ2-β+241-x3-βΓ4-β-121-x2-βΓ3-β+21-x1-βΓ2-β.Then the time-space fractional superdiffusion obtains the same exact solution as Example 1: u(x,t)=(2t-1)2x2(1-x)2.Analogous to the case of fractional subdiffusion equation, we consider the errors under theL2 norm and convergence orders of the LDG numerical solutions for the original equation with different values of the fractional derivative parameters α and β. First, we choose the appropriate time step size τ such that the time discretization errors are negligible as compared with the spatial errors. Then choose the space step size sequence as h=1/2i,(i=2,…,6), Tables 9 and 10 list out the errors and the convergence rate when we fix β (or α), and change α (or β). As can be seen from the table, we obtain the optimal spatial convergence rate O(hk+1), which agrees with the theoretical results proved in Theorem 22. The numerical results in Tables 11 and 12 show that when h=C′τ (C′ is a positive constant), we choose the time step size sequence τ=1/2i,(i=2,⋯,6), fix β (or α), and then change α (or β); the temporal accuracy of the time-space fractional superdiffusion can reach O(τ3-α) for 1<α<2 and be close to 1 when α→2. The numerical results are consistent with the theoretical results proved in Theorem 22.Table 9
Example3. The L2 errors and convergence rate in space direction. β=1.6, T=0.5.
h
α
=
1.02
Rate
α
=
1.5
Rate
α
=
1.99
Rate
1/4
2.9816
×
1
0
-
2
—
3.1419
×
1
0
-
2
—
3.8907
×
1
0
-
2
—
1/8
8.3473
×
1
0
-
3
1.8367
9.0427
×
1
0
-
3
1.7968
1.1709
×
1
0
-
2
1.7324
1/16
2.1928
×
1
0
-
3
1.9285
2.4998
×
1
0
-
3
1.8549
3.3840
×
1
0
-
3
1.7908
1/32
5.9729
×
1
0
-
4
1.8763
6.5214
×
1
0
-
4
1.9386
9.2694
×
1
0
-
4
1.8682
1/64
1.5433
×
1
0
-
4
1.9524
1.6745
×
1
0
-
4
1.9614
2.5002
×
1
0
-
4
1.8904Table 10
Example3. The L2 errors and convergence rate in space direction. α=1.5, T=0.5.
h
β
=
1.1
Rate
β
=
1.5
Rate
β
=
1.9
Rate
1/4
3.9209
×
1
0
-
2
—
4.2305
×
1
0
-
2
—
4.6712
×
1
0
-
2
—
1/8
1.0526
×
1
0
-
2
1.8972
1.2116
×
1
0
-
2
1.8035
1.2515
×
1
0
-
2
1.9011
1/16
2.8406
×
1
0
-
3
1.8897
3.2656
×
1
0
-
3
1.8915
3.2087
×
1
0
-
3
1.9636
1/32
7.5049
×
1
0
-
4
1.9203
8.2774
×
1
0
-
4
1.9801
8.0486
×
1
0
-
4
1.9952
1/64
1.9265
×
1
0
-
4
1.9618
2.0789
×
1
0
-
4
1.9933
2.0139
×
1
0
-
4
1.9987Table 11
Example3. The L2 errors and convergence rate in time direction. h=C′τ, C′=0.5, and β=1.6.
τ
α
=
1.05
Rate
α
=
1.6
Rate
α
=
1.999
Rate
1/4
5.6293
×
1
0
-
2
—
4.9876
×
1
0
-
2
—
4.3574
×
1
0
-
2
—
1/8
2.0868
×
1
0
-
2
1.4316
2.2149
×
1
0
-
2
1.2309
3.3742
×
1
0
-
2
0.3689
1/16
7.2314
×
1
0
-
3
1.5290
8.7045
×
1
0
-
3
1.2876
2.3331
×
1
0
-
2
0.5323
1/32
2.3828
×
1
0
-
3
1.6016
3.3831
×
1
0
-
3
1.3634
1.4645
×
1
0
-
2
0.6718
1/64
7.1811
×
1
0
-
4
1.7304
1.2871
×
1
0
-
3
1.3942
8.8796
×
1
0
-
3
0.7219Table 12
Example3. The L2 errors and convergence rate in time direction. h=C′τ, C′=0.5, and α=1.5.
τ
β
=
1.1
Rate
β
=
1.5
Rate
β
=
1.99
Rate
1/4
6.0219
×
1
0
-
2
—
4.8892
×
1
0
-
2
—
3.6905
×
1
0
-
2
—
1/8
2.3996
×
1
0
-
2
1.3274
1.8938
×
1
0
-
2
1.3691
1.3836
×
1
0
-
2
1.4153
1/16
9.7672
×
1
0
-
3
1.2968
6.9743
×
1
0
-
3
1.4412
4.8223
×
1
0
-
3
1.5207
1/32
3.8297
×
1
0
-
3
1.3507
2.4718
×
1
0
-
3
1.4649
1.7062
×
1
0
-
3
1.4989
1/64
1.4505
×
1
0
-
3
1.4006
8.4813
×
1
0
-
4
1.5432
5.8790
×
1
0
-
4
1.5372


## 6.1. Implementation

As examples, we consider the time-space fractional subdiffusion/superdiffusion equations as follows:(99)∂αux,t∂tα-ϵ--△β/2ux,t=fx,t,inx∈Ω,t∈0,T,ux,0=u0x,inx∈Ω,where Ω=[a,b], T>0, 0<α<1 for the subdiffusion case (1), or 1<α<2 for the superdiffusion case (2), and 1<β<2, ϵ>0.For brevity, we only derive the linear system obtained from the fully discrete LDG scheme (35). For scheme (42), we can obtain the corresponding linear system similarly. We consider the case of discontinuous piecewise polynomials basis function sequences {ϕi(x)}i=1k in space Vkh. By the definition of the space Vkh, we know that, for all v∈Vkh, v(x)=∑i=1N‍viϕi(x) holds, where vi=v(xi).By the LDG scheme (35), we obtain the fully discrete scheme of (99) as follows:(100)∫Ωuhnvdx+α0ϵ∫Ωqhnvxdx-∑j=1Nqhn^vj+1/2--qhn^vj-1/2+=∑i=1n-1bi-1-bi∫Ωuhn-ivdx+bn-1∫Ωuh0vdx+∫Ωfx,tvdx,∫Ωqhnwdx-∫Ω△β-2/2phnwdx=0,∫Ωphnzdx+ϵ∫Ωuhnzxdx-∑j=1Nuhn^zj+1/2--uhn^zj-1/2+=0,∫Ωuh0vdx-∫Ωu0xvdx=0,where α0=ταΓ(2-α).In terms of the basis{ϕi(x)}i=1N, we denote approximate solution functions by uhn=∑i=1N‍uj(t)ϕj(x), phn=∑i=1N‍pj(t)ϕj(x), and qhn=∑i=1N‍qj(t)ϕj(x) and let test functions v=∑i=1N‍viϕi(x),vi=v(xi), w=∑i=1N‍wiϕi(x),wi=w(xi), and z=∑i=1N‍ziϕi(x),zi=z(xi). Putting these expressions into the LDG scheme (100) and taking the flux in (36), we get the linear system of the scheme in jth element at tn as follows:(101)M0α0ϵA-H11+H120RMϵA-H13+H14M0UjnPjnQjn=Fnv+∑i=1n-1bi-1-biUn-iv+bn-1U0v,where M={(ϕmj(x),ϕkj(x))}j=1,…,N is the mass matrix, m,k=0,1,… denote the polynomials order, and (102)A=ϕmjx,ϕkj′xj=1,…,N,H11=ϕmjxj-1/2,ϕkjxj+1/2j=1,…,N,H12=ϕmjxj-1/2,ϕkjxj-1/2j=1,…,N,H13=ϕmjxj+1/2,ϕkjxj+1/2j=1,…,N,H14=ϕmj-1xj-1/2,ϕkjxj-1/2j=1,…,N.R is the quadrature of the fractional integral operator △(β-2)/2phn and test function w in jth element. From [29], we know that (103)△β-2/2phn=-​aIx2-βphn+xIb2-βphn2cos⁡βπ/2,where(104)​aIx2-βphnx,t=1Γ2-β∫axx-ξ1-βphnξ,tdξ,​xIb2-βphnx,t=1Γ2-β∫xbξ-x1-βphnξ,tdξ,and Fn=(f1,f2,…,fn)Tis a vector valued function, fj(tn)=f(xj,tn). Then, by solving the linear system (101), we can obtain the solution Un.


## 6.2. Numerical Results

Example 1.
Consider the fractional time-space subdiffusion equation (99) in domain [0,1]×[0,0.5] which satisfies the following initial boundary conditions: (105)ux,0=x21-x2,x∈0,1.In order to obtain the exact solution u(x,t)=(2t-1)2x2(1-x)2, we choose the source term as (106)fx,t=x21-x28t2-αΓ3-α-4t1-αΓ2-α+t-αΓ1-α+2t-122cos⁡πβ+1/224x3-βΓ4-β-12x2-βΓ3-β+2x1-βΓ2-β+241-x3-βΓ4-β-121-x2-βΓ3-β+21-x1-βΓ2-β.Choosing the fluxes in (36) and adopting the linear piecewise polynomials, (99) is solved numerically. We investigate the numerical solutions with respect to different fractional order parameters α and β; we select the appropriate time step size τ such that the time discretization errors are negligible as compared with the spatial errors. Then choose space step size sequence as h=1/2i,(i=2,…,6); Tables 1 and 2 list out the L2 errors and the corresponding convergence order, respectively, fixing the value of the one parameter β (or α), with the values change of the other one parameter α (or β). As can be seen from the data in the two tables, our schemes can achieve the optimal convergence O(hk+1) in spatial direction; this matches the theoretical results of Theorem 21.Table 1
Example1. The L2 errors and convergence rate in space direction. β=1.6, T=0.5.
h
α=0.1
Rate
α=0.5
rate
α=0.9
Rate
1/4
3.3592×10-2
—
2.7826
×
1
0
-
2
—
1.4821×10-2
—
1/8
9.1511
×
1
0
-
3
1.8761
6.9632
×
1
0
-
3
1.9986
3.7341×10-3
1.9888
1/16
2.3499
×
1
0
-
3
1.9613
1.7411
×
1
0
-
3
1.9997
9.3353×10-4
2.0000
1/32
5.9437×10-4
1.9832
4.3333
×
1
0
-
4
2.0065
2.3241×10-4
2.0060
1/64
1.4875
×
1
0
-
4
1.9984
1.0748
×
1
0
-
4
2.0113
5.7422×10-5
2.0170Table 2
Example1. The L2 errors and convergence rate in space direction. α=0.6, T=0.5.
h
β=1.1
Rate
β=1.5
rate
β=1.9
Rate
1/4
1.4938
×
1
0
-
2
—
5.5316
×
1
0
-
2
—
4.6063
×
1
0
-
2
—
1/8
4.0147
×
1
0
-
3
1.8956
1.3845
×
1
0
-
2
1.9983
1.1679
×
1
0
-
2
1.9796
1/16
1.0400
×
1
0
-
3
1.9487
3.4615
×
1
0
-
3
1.9999
2.9465
×
1
0
-
3
1.9869
1/32
2.6052
×
1
0
-
4
1.9971
8.6539
×
1
0
-
4
2.0000
7.5009
×
1
0
-
4
1.9739
1/64
6.5141
×
1
0
-
5
1.9998
2.1589
×
1
0
-
4
2.0030
1.9156
×
1
0
-
4
1.9692The numerical example results listed in Tables3 and 4 show that when we choose the appropriate space step size h and the time step size sequence τ=1/2i,(i=2,…,6), fixing β (or α), and changing α (or β), the temporal accuracy of the time-space fractional subdiffusion equation under the L2 norm can converge to O(τ2-α) for 0<α<1, and close to 1 when α→1, the results agree with the theoretical computation results in Theorem 21.Table 3
Example1. The L2 errors and convergence rate in time direction. β=1.6, T=0.5.
τ
α=0.2
Rate
α=0.6
Rate
α
=
0.999
Rate
1/4
4.8064
×
1
0
-
3
—
5.5078
×
1
0
-
3
—
5.4127
×
1
0
-
3
—
1/8
1.4486
×
1
0
-
3
1.7302
2.1182
×
1
0
-
3
1.3786
2.8977
×
1
0
-
3
0.9014
1/16
4.2962
×
1
0
-
4
1.7536
8.0864
×
1
0
-
4
1.3893
1.5246
×
1
0
-
3
0.9265
1/32
1.2320
×
1
0
-
4
1.8020
3.0576
×
1
0
-
4
1.4031
7.6910
×
1
0
-
4
0.9872
1/64
3.5291
×
1
0
-
5
1.8037
1.1501
×
1
0
-
4
1.4106
3.8372
×
1
0
-
4
1.0031Table 4
Example1. The L2 errors and convergence rate in time direction. α=0.5, T=0.5.
τ
β=1.05
Rate
β=1.5
Rate
β=1.99
Rate
1/4
6.7348
×
1
0
-
3
—
6.8961
×
1
0
-
3
—
5.9468
×
1
0
-
3
—
1/8
2.5767
×
1
0
-
3
1.3861
2.6180
×
1
0
-
3
1.3973
2.2038
×
1
0
-
3
1.4321
1/16
9.6141
×
1
0
-
4
1.4223
9.8328
×
1
0
-
4
1.4128
7.9725
×
1
0
-
4
1.4669
1/32
3.4618
×
1
0
-
4
1.4736
3.6092
×
1
0
-
4
1.4459
2.8224
×
1
0
-
4
1.4981
1/64
1.2263
×
1
0
-
4
1.4972
1.2211
×
1
0
-
4
1.4943
9.9691
×
1
0
-
5
1.5014Example 2.
Consider the time-space fractional subdiffusion(0<α<1) equation (99); for t∈[0,1],x∈[0,1], we choose the exact solution as u(x,t)=sin⁡(2πt)x2(1-x)2 and ϵ=0.05.
Similar to Example1, we give out the L2 errors and convergence rate with fixing α (or β) and changing β or α. We also select the appropriate time step size τ, which guarantee that the time discretization errors are negligible as compared with the spatial errors. Choosing the space step size sequence as h=1/2i,(i=2,…,6), we obtain the optimal convergence rate in space direction (as shown in Tables 5 and 6, this corresponds to the theoretical results as illustrated in Theorem 21). Tables 7 and 8 listed out the temporal accuracy; here we choose the appropriate space step size h and the time step size sequence as τ=1/2i,(i=2,…,6). The data in Table 7 display that, for 0<α<1, the convergence order in time direction can be up to O(τ2-α) and close to 1 order for α→1, which are consistent with Theorem 21.Table 5
Example2. The L2 errors and convergence rate in space direction. β=1.6, T=1.
h
α=0.2
Rate
α=0.6
Rate
α=0.99
Rate
1/4
3.1592
×
1
0
-
2
—
1.9328
×
1
0
-
2
—
2.1532
×
1
0
-
2
—
1/8
9.5228
×
1
0
-
3
1.7301
4.8649
×
1
0
-
3
1.9902
5.4460
×
1
0
-
3
1.9832
1/16
2.4728
×
1
0
-
3
1.9452
1.2174
×
1
0
-
3
1.9985
1.3693
×
1
0
-
3
1.9917
1/32
1.2521
×
1
0
-
3
1.9818
3.0317
×
1
0
-
4
2.0057
3.3469
×
1
0
-
4
2.0326
1/64
3.1377
×
1
0
-
4
1.9966
7.4936
×
1
0
-
5
2.0164
8.3644
×
1
0
-
5
2.0005Table 6
Example2. The L2 errors and convergence rate in space direction. α=0.6, T=1.
h
β
=
1.05
Rate
β
=
1.5
Rate
β
=
1.98
Rate
1/4
3.2624
×
1
0
-
2
—
1.6352
×
1
0
-
2
—
2.3573
×
1
0
-
2
—
1/8
8.8793
×
1
0
-
3
1.8782
4.1930
×
1
0
-
3
1.9961
5.9182
×
1
0
-
3
1.9939
1/16
2.3194
×
1
0
-
3
1.9367
1.0570
×
1
0
-
3
1.9879
1.4854
×
1
0
-
3
1.9943
1/32
5.8251
×
1
0
-
4
1.9934
2.6453
×
1
0
-
4
1.9986
3.7137
×
1
0
-
4
1.9999
1/64
1.4571
×
1
0
-
4
1.9991
6.5808
×
1
0
-
5
2.0071
9.2703
×
1
0
-
5
2.0022Table 7
Example2. The L2 errors and convergence rate in time direction. β=1.6, T=1.
τ
α=0.2
Rate
α=0.6
Rate
α=0.999
Rate
1/4
6.3463
×
1
0
-
3
—
5.9462
×
1
0
-
3
—
5.8016
×
1
0
-
3
—
1/8
1.9707
×
1
0
-
3
1.6872
2.4095
×
1
0
-
3
1.3032
3.3481
×
1
0
-
3
0.7931
1/16
6.0938
×
1
0
-
4
1.6933
9.3755
×
1
0
-
4
1.3618
1.8392
×
1
0
-
3
0.8642
1/32
1.7484
×
1
0
-
4
1.8013
3.5496
×
1
0
-
4
1.4012
9.2200
×
1
0
-
4
0.9963
1/64
5.0206
×
1
0
-
5
1.8001
1.3445
×
1
0
-
4
1.4006
4.4695
×
1
0
-
4
0.9877Table 8
Example2. The L2 errors and convergence rate in time direction. α=0.5, T=1.
τ
β
=
1.05
Rate
β
=
1.5
Rate
β
=
1.99
Rate
1/4
6.7819
×
1
0
-
3
—
6.2487
×
1
0
-
3
—
5.1918
×
1
0
-
3
—
1/8
2.6331
×
1
0
-
3
1.3649
2.3006
×
1
0
-
3
1.4415
1.9810
×
1
0
-
3
1.3900
1/16
9.9694
×
1
0
-
4
1.4012
8.2688
×
1
0
-
4
1.4763
7.2403
×
1
0
-
4
1.4521
1/32
3.5724
×
1
0
-
4
1.4806
2.9434
×
1
0
-
4
1.4902
2.5607
×
1
0
-
4
1.4995
1/64
1.2685
×
1
0
-
4
1.4937
1.0393
×
1
0
-
4
1.5018
9.1235
×
1
0
-
5
1.4889Example 3.
In this example, we consider the case for1<α<2; that is, (99) is a time-space fractional superdiffusion model. We assume that the equation satisfies the following initial boundary conditions in the domain [0,1]×[0,0.5]:(107)ux,0=x21-x2,x∈0,1,utx,0=-4x21-x2,x∈0,1.Here we choose ϵ=1, and the source term is (108)fx,t=x21-x28t2-αΓ3-α+t-αΓ1-α+2t-122cos⁡πβ+1/224x3-βΓ4-β-12x2-βΓ3-β+2x1-βΓ2-β+241-x3-βΓ4-β-121-x2-βΓ3-β+21-x1-βΓ2-β.Then the time-space fractional superdiffusion obtains the same exact solution as Example 1: u(x,t)=(2t-1)2x2(1-x)2.Analogous to the case of fractional subdiffusion equation, we consider the errors under theL2 norm and convergence orders of the LDG numerical solutions for the original equation with different values of the fractional derivative parameters α and β. First, we choose the appropriate time step size τ such that the time discretization errors are negligible as compared with the spatial errors. Then choose the space step size sequence as h=1/2i,(i=2,…,6), Tables 9 and 10 list out the errors and the convergence rate when we fix β (or α), and change α (or β). As can be seen from the table, we obtain the optimal spatial convergence rate O(hk+1), which agrees with the theoretical results proved in Theorem 22. The numerical results in Tables 11 and 12 show that when h=C′τ (C′ is a positive constant), we choose the time step size sequence τ=1/2i,(i=2,⋯,6), fix β (or α), and then change α (or β); the temporal accuracy of the time-space fractional superdiffusion can reach O(τ3-α) for 1<α<2 and be close to 1 when α→2. The numerical results are consistent with the theoretical results proved in Theorem 22.Table 9
Example3. The L2 errors and convergence rate in space direction. β=1.6, T=0.5.
h
α
=
1.02
Rate
α
=
1.5
Rate
α
=
1.99
Rate
1/4
2.9816
×
1
0
-
2
—
3.1419
×
1
0
-
2
—
3.8907
×
1
0
-
2
—
1/8
8.3473
×
1
0
-
3
1.8367
9.0427
×
1
0
-
3
1.7968
1.1709
×
1
0
-
2
1.7324
1/16
2.1928
×
1
0
-
3
1.9285
2.4998
×
1
0
-
3
1.8549
3.3840
×
1
0
-
3
1.7908
1/32
5.9729
×
1
0
-
4
1.8763
6.5214
×
1
0
-
4
1.9386
9.2694
×
1
0
-
4
1.8682
1/64
1.5433
×
1
0
-
4
1.9524
1.6745
×
1
0
-
4
1.9614
2.5002
×
1
0
-
4
1.8904Table 10
Example3. The L2 errors and convergence rate in space direction. α=1.5, T=0.5.
h
β
=
1.1
Rate
β
=
1.5
Rate
β
=
1.9
Rate
1/4
3.9209
×
1
0
-
2
—
4.2305
×
1
0
-
2
—
4.6712
×
1
0
-
2
—
1/8
1.0526
×
1
0
-
2
1.8972
1.2116
×
1
0
-
2
1.8035
1.2515
×
1
0
-
2
1.9011
1/16
2.8406
×
1
0
-
3
1.8897
3.2656
×
1
0
-
3
1.8915
3.2087
×
1
0
-
3
1.9636
1/32
7.5049
×
1
0
-
4
1.9203
8.2774
×
1
0
-
4
1.9801
8.0486
×
1
0
-
4
1.9952
1/64
1.9265
×
1
0
-
4
1.9618
2.0789
×
1
0
-
4
1.9933
2.0139
×
1
0
-
4
1.9987Table 11
Example3. The L2 errors and convergence rate in time direction. h=C′τ, C′=0.5, and β=1.6.
τ
α
=
1.05
Rate
α
=
1.6
Rate
α
=
1.999
Rate
1/4
5.6293
×
1
0
-
2
—
4.9876
×
1
0
-
2
—
4.3574
×
1
0
-
2
—
1/8
2.0868
×
1
0
-
2
1.4316
2.2149
×
1
0
-
2
1.2309
3.3742
×
1
0
-
2
0.3689
1/16
7.2314
×
1
0
-
3
1.5290
8.7045
×
1
0
-
3
1.2876
2.3331
×
1
0
-
2
0.5323
1/32
2.3828
×
1
0
-
3
1.6016
3.3831
×
1
0
-
3
1.3634
1.4645
×
1
0
-
2
0.6718
1/64
7.1811
×
1
0
-
4
1.7304
1.2871
×
1
0
-
3
1.3942
8.8796
×
1
0
-
3
0.7219Table 12
Example3. The L2 errors and convergence rate in time direction. h=C′τ, C′=0.5, and α=1.5.
τ
β
=
1.1
Rate
β
=
1.5
Rate
β
=
1.99
Rate
1/4
6.0219
×
1
0
-
2
—
4.8892
×
1
0
-
2
—
3.6905
×
1
0
-
2
—
1/8
2.3996
×
1
0
-
2
1.3274
1.8938
×
1
0
-
2
1.3691
1.3836
×
1
0
-
2
1.4153
1/16
9.7672
×
1
0
-
3
1.2968
6.9743
×
1
0
-
3
1.4412
4.8223
×
1
0
-
3
1.5207
1/32
3.8297
×
1
0
-
3
1.3507
2.4718
×
1
0
-
3
1.4649
1.7062
×
1
0
-
3
1.4989
1/64
1.4505
×
1
0
-
3
1.4006
8.4813
×
1
0
-
4
1.5432
5.8790
×
1
0
-
4
1.5372


## 7. Conclusions

In this paper, we proposed the fully discrete local discontinuous Galerkin scheme for solving the time-space fractional subdiffusion/superdiffusion equations, respectively. Through carefully choosing the numerical flux on boundary terms and making a detailed theoretical analysis, we proved that our numerical schemes are unconditionally stable with convergence under theL2 norm. From the numerical experiment results, we can observe that a convergence rate O(τ2-α+hk+1) is achieved for 0<α<1, and for 1<α<2, when the space step size h and the time step size τ satisfy the relationship h=C′τ (C′ is a constant), the convergence rate of our scheme can approach O(τ3-α+hk+1). The numerical results show that the fully discrete local discontinuous Galerkin (LDG) approximation is effective and powerful for solving time-space fractional subdiffusion/superdiffusion equations.



---

*Source: 4961797-2017-03-16.xml*