# Local Discontinuous Galerkin Method for Nonlinear Time-Space Fractional Subdiffusion/Superdiffusion Equations

Authors: Meilan Qiu; Dewang Li; Yanyun Wu
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

Fractional partial differential equations with time-space fractional derivatives describe some important physical phenomena. For example, the subdiffusion equation (time order0<α<1) is more suitable to describe the phenomena of charge carrier transport in amorphous semiconductors, nuclear magnetic resonance (NMR) diffusometry in percolative, Rouse, or reptation dynamics in polymeric systems, the diffusion of a scalar tracer in an array of convection rolls, or the dynamics of a bead in a polymeric network, and so on. However, the superdiffusion case (1<α<2) is more accurate to depict the special domains of rotating flows, collective slip diffusion on solid surfaces, layered velocity fields, Richardson turbulent diffusion, bulk-surface exchange controlled dynamics in porous glasses, the transport in micelle systems and heterogeneous rocks, quantum optics, single molecule spectroscopy, the transport in turbulent plasma, bacterial motion, and even for the flight of an albatross (for more physical applications of fractional sub-super diffusion equations, one can see Metzler and Klafter in 2000). In this work, we establish two fully discrete numerical schemes for solving a class of nonlinear time-space fractional subdiffusion/superdiffusion equations by using backward Euler difference 1<α<2 or second-order central difference 1<α<2/local discontinuous Galerkin finite element mixed method. By introducing the mathematical induction method, we show the concrete analysis for the stability and the convergence rate under the L2 norm of the two LDG schemes. In the end, we adopt several numerical experiments to validate the proposed model and demonstrate the features of the two numerical schemes, such as the optimal convergence rate in space direction is close to Ohk+1. The convergence rate in time direction can arrive at Oτ2−α when the fractional derivative is 0<α<1. If the fractional derivative parameter is 1<α<2 and we choose the relationship as h=C′τ (h denotes the space step size, C′ is a constant, and τ is the time step size), then the time convergence rate can reach to Oτ3−α. The experiment results illustrate that the proposed method is effective in solving nonlinear time-space fractional subdiffusion/superdiffusion equations.


---

## Body

## 1. Introduction

Over the past several decades, fractional differential equations attract more and more scholar’s attention due to their wide applications in science and engineering [1–4]. Because of the nonlocality and hereditariness, as well as other properties of fractional derivatives, researchers find that fractional partial differential equations (FPDEs) can be more accurately and effectively used to depict the complex systems with memory, hereditary, and long-range interactions than classical partial differential equations. Among them, the space-time fractional diffusion equation we mean is an evolution equation; they imply for the flux fractional Fick’s law which accounts for spatial and temporal nonlocality. The fundamental solution (for the Cauchy problem) of this type of fractional diffusion equation can be used to interpret a probability density evolving in time of a peculiar self-similar stochastic process, and we also view it as a generalized diffusion process [5]. Especially, fractional subdiffusion/superdiffusion equations (FSEs) play important roles in describing a special type of anomalous diffusion process, and thus become very popular for many real applications [6, 7].It is of great importance to seek efficient methods to solve fractional subdiffusion/superdiffusion equations (FSEs). However, generally speaking, it is difficult to obtain the analytic solutions for most of the fractional differential equations. Moreover, most nonlinear fractional subdiffusion/superdiffusion equations (FSEs) are not solvable by analytic methods. This motivates researchers to resort to efficient and reliable numerical methods for solving these fractional subdiffusion/superdiffusion equations and then understand the behaviors of them.In recent years, a large number of researchers focus on various numerical methods for solving FSEs. For example, Zeng et al. [8] used the fast finite difference method for solving the high-dimensional time-fractional subdiffusion equation. Gao and Sun [9] applied the L1 approximation for the time-fractional derivative and developed a compact finite difference scheme for the fractional subdiffusion equation by using the energy method; they prove the solvability and stability, as well as L∞ convergence. Zhang et al. presented a Crank–Nicolson-type difference scheme for solving the subdiffusion equation [10]. Wang and Vong [11] established a compact finite difference scheme for solving the modified anomalous fractional subdiffusion equation and obtained a second-order time accuracy convergence rate. Jin et al. [12] obtained a fully discrete scheme for the subdiffusion equation by developing essential initial corrections at the starting two steps for the Crank–Nicolson scheme in time and the Galerkin finite element method in space. Shen et al. in [13] applied the implicit and explicit difference methods to solve the space-time fractional convection-subdiffusion equations. Saadatmandi et al. [14] constructed a sinc-Legendre collocation method for solving the fractional convection-subdiffusion equation. Liu’s group [15] analyzed a radial basis function approximation scheme for solving the fractional mobile/immobile convection-subdiffusion equation. Zhang et al. [10] proposed the two-order Crank–Nicolson difference method for solving the subdiffusion equation, and they proved the discrete H1 norm convergence rate as well as the maximum norm error estimate. In [16], Jiang and Ma adopted the high-order finite element method to solve the time-fractional subdiffusion equation. Zeng et al. [17] proposed a finite difference/element mixed approach to solve the time-fractional subdiffusion problem. Hu et al. in [18], Hu et al. established the center box approach to solve the time-fractional subdiffusion system. In [19], Bhrawy et al. gave out the Legendre spectral collocation scheme for solving nonlinear fractional subdiffusion and reaction subdiffusion equations. Du et al. [20] developed a fully discrete local discontinuous Galerkin (LDG) method for solving nonlinear time-fractional fourth-order partial differential equations. Aboelenen in [21] constructed a local discontinuous Galerkin (LDG) method for solving the distributed-order time and Riesz space fractional convection-subdiffusion- and Schrödinger-type equations. Liu et al. [22] established a high-order local discontinuous Galerkin (LDG) approach combined with weighted and shifted Grünwald difference (WSGD) approximation for solving a Caputo-type time-fractional subdiffusion equation.However, to the best of our knowledge so far, the LDG method have not been considered in solving the nonlinear time-space fractional subdiffusion and/or superdiffusion equations. Therefore, our aim in this paper is to establish two fully discrete LDG schemes to solve the nonlinear time-space fractional subdiffusion0<α<1 equation:(1)∂αux,t∂tα+∂∂xfu=ε−−Δβ/2ux,t,x∈Ω=a,b,t∈0,T,ux,t=0,x∈∂Ω,t∈0,T,ux,0=u0x,x∈Ω=a,b,and the nonlinear time-space fractional superdiffusion 1<α<2 equation:(2)∂αux,t∂tα+∂∂xfu=ε−−Δβ/2ux,t,x∈Ω=a,b,t∈0,T,ux,t=0,x∈∂Ω,t∈0,T,ux,0=ϕx,x∈Ω=a,b,utx,0=ψx,x∈Ω=a,b,where Ω=a,b, T>0, 1<β<2, denotes the order number of space fractional derivatives, ϵ is a real positive parameter, u0,ϕx,ψx∈H01Ω are the given smooth functions, and the nonlinear term fu satisfies the Lipschitz continuous condition on the domain Ω.For simplicity, we only consider the periodic boundary conditions in this article. Notice that this assumption is not essential, and our method can be directly extended to the problems with aperiodic boundary conditions.The time fractional derivative in nonlinear time-space fractional subdiffusion equation (1) uses the following Caputo fractional derivative [2, 3]:(3)∂αux,t∂tα=1Γ1−α∫0t∂ux,s∂sdst−sα,0<α<1,∂ux,t∂t,α=1,and the time-fractional Caputo derivative in nonlinear time-space fractional superdiffusion equation (2) is defined as follows [23]:(4)∂αux,t∂tα=1Γ2−α∫0tt−s1−α∂2ux,s∂s2ds,1<α<2,∂2ux,t∂t2,α=2.The spatial fractional derivative in both equations (1) and (2) is an important tool to describe the anomalous diffusion phenomenon [24]; when 1<β<2, it represents a Lévy β-stable flight [25], and when β⟶2, it depicts a Brownian diffusion process. In this paper, we still use the definition of the fractional Laplacian operator −−Δβ/2ux,t, 1<β<2, which is also known as the Riesz fractional derivative ∂β/∂xβux,t, and we assume that the unknown function and all its derivatives ux,t,u′x,t,…,un−1x,t vanish at the end points on an infinite domain i.e.,x=∓∞ [26].The discontinuous Galerkin (DG) method is a classical finite element method using discontinuous, piecewise polynomials as the solution and the test spaces in the spatial direction. There are various DG methods suggested in the literature to solve diffusion problems; its generalization calls the local discontinuous Galerkin (LDG) methods which are introduced in [27] by Cockburn and Shu and further studied in [28, 29]. Recently, the LDG method is applied to solving the fractional partial differential equations [30–32]. In [33, 34], Xu et al. developed the LDG method for solving the time fractional diffusion equation and fractional convection-diffusion equations, respectively.The organization of this article is as follows. In the second part, we review the appropriate functional spaces and some important definitions and properties of fractional derivatives. Based on the technique introduced in [34], we rewrite the fractional Laplacian operator (with the order β) as a composite of an integer derivative with one order and a fractional integral with an order 2−β, and we, respectively, get the fully discrete LDG schemes for the nonlinear time-space fractional subdiffusion equation and the nonlinear time-space fractional superdiffusion equation in the third part. In the fourth and fifth part, we provide the concrete stability and L2 error analysis for the two fully discrete schemes, respectively. In the sixth part, we present several examples to verify the theoretical results and elaborate the concluding remarks in the end.


## 2. Preliminary Results

In order to prove the convergence rate of two fully discrete LDG numerical schemes better, we briefly list some relevant properties of fractional integrals and derivatives as well as several basic lemmas that are needed as follows; some of them have already been proved in [24]. Therefore, we ignore the detailed proof processes here.Lemma 1.
(linearity).(5)Iaxαλfx+μgx=λaIxαfx+μaIxαgx.Daxαλfx+μgx=λaDxαfx+μaDxαgx.Lemma 2.
(semigroup property).(6)Iaxα+βfx=IaxαIaxβfx=IaxβIaxαfx.These lemmas combine with definitions 2.1–2.3 in [24]; we can convert the fractional Laplacian derivative into a composite form of an integer derivative with one order and a fractional integral with an order 2−β1<β<2 [34]:(7)−−Δβ/2ux,t=∂∂xΔβ−2/2∂ux,t∂x.Lemma 3.
(adjoint property).(8)Iaxαu,uΩ=u,IxbαuΩ.Lemma 4.
(see [35]).(9)I−∞xαu,R=cosαπuJL−αR2=cosαπuJR−αR2.Lemma 5.
(see [36]). Let 0<θ<2,θ≠1. Then, for any w,v∈H0θ/20,X, then(10)D0Rxθw,vL20,X=D0Rxθ/2w,DtRXθ/2vL20,X.By composite form (7) of the fractional Laplacian operator and Lemmas 3–5, together with the properties of Riemann–Liouville fractional integrals, the following lemmas hold.Lemma 6.
For all real numbers,0<s<1, the fractional integral satisfies the property(11)Δ−su,uΩ=uJL−αΩ2=uJR−αΩ2.Lemma 7.
(fractional Poincaré–Friedrichs inequality, see [35]).(12)uL2Ω≤CuJL,0μΩ,u∈JL,0μΩ,μ∈R,uL2Ω≤CuJR,0μΩ,u∈JR,0μΩ,μ∈R.Lemma 8.
(see [34]). Suppose ux is a smooth function defined on Ω⊂R. Ωh is a discretization of the domain with interval width h, uhx is an approximation solution of u in polynomial function space Phk. For all j, uhx∈Ij is a polynomial of degree up to order k, and u,vIj=uh,vIj,∀v∈Pk, k is degree of polynomials. For all −1<α≤0, we have(13)Δβ/2ux−Δβ/2uhxL2Ω≤Chk+1,where C is a constant independent of h.To prepare the analysis of the error estimates, we need the following two projection operators in one dimensiona,b; denote them by Q, i.e., for each j,(14)∫IjQωx−ωxυxdx=0,∀υ∈PkIj,and special projection P−+, i.e., for any j,(15)∫IjP+ωx−ωxυxdx=0,∀υ∈Pk−1Ij,P+ωxj−1/2+=ωxj−1/2,(16)∫IjP−ωx−ωxυxdx=0,∀υ∈Pk−1Ij,P−ωxj+1/2−=ωxj+1/2.The above two projection operators satisfy the following approximation inequality [37]:(17)ωe+hωe∞+h1/2ωeτh≤Chk+1,where ωe=Pωx−ωx or ωe=P−+ωx−ωx, where the constant C is just depending on ω but independent on h and τh denotes the set of boundary points of all elements Ij.In this article, we useC to denote the positive constant which different values from each case and apply the usual notation of norms in Sobolev spaces. Let ⋅,⋅D denote the scalar inner product on L2D, and denote the corresponding norm by ⋅D. If D=Ω, we drop D.


## 3. Numerical Schemes

In this part, we present two fully discrete local discontinuous Galerkin (LDG) schemes for solving the nonlinear time-space fractional subdiffusion/superdiffusion problems, respectively.Letτ=T/M be the time mesh size; M is a positive integer, tn=nτ,n=0,1,…,M are mesh points, and tn−1/2=tn+tn−1/2, n=0,1,…,M are midmesh points. The unmarked norm ⋅ refers to the usual L2 norm in the whole text.


### 3.1. Subdiffusion Case

We use the backward Euler scheme to discrete the time fractional derivative∂αux,t/∂tα of nonlinear subdiffusion 0<α<1 case (1) at tn which is as follows [38, 39]:(18)∂αux,tn∂tα=τ1−αΓ2−α∑i=0n−1biux,tn−i−ux,tn−i−1τ+rnx.Here,bi=i+11−α−i1−α satisfies properties as follows:(19)1=b0>b1>b2>⋯>bn>0,bn⟶0n⟶∞,∑i=1nbi−1−bi+bn=1,where rnx is the truncation error and satisfies the inequality rnx≤Cτ2−α [38], constant C depends on u,T,α.Now, we rewrite equation (1) as a first-order system:(20)∂αux,t∂tα+∂∂xfu−ε∂∂xq=0,(21)q−Δβ−2/2p=0,(22)p−ε∂∂xux,t=0.Combining with equation (18), we obtain the weak formulation of system (21) at tn:(23)∫Ωux,tnvdx−α0∫Ωfux,tnvxdx+α0ε∫Ωqx,tnvxdx+α0∑j=1Nfux,tnvj+1/2−−fux,tnvj−1/2+−α0ε∑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 numerical solutions of u⋅,tn,p⋅,tn,q⋅,tn, respectively. We define the fully discrete local discontinuous Galerkin (LDG) scheme as follows: find uhn,phn,qhn∈Vhk such that, for all v,w,z∈Vhk, we have(24)∫Ωuhnvdx−α0∫Ωfuhnvxdx+α0ε∫Ωqhnvxdx+α0∑j=1Nf^uhnvj+1/2−−f^uhnvj−1/2+−α0ε∑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,where α0=ταΓ2−α.The “hat” terms in equation (18) in the cell boundary terms from integration by parts are the so-called “numerical fluxes” [40], 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 of the LDG numerical scheme. Here, for the sake of concise, we take the same alternating numerical fluxes on the boundary as in [24], i.e.,(25)uhn^=uhn−,qhn^=qhn+.We remark that the choice for the fluxes (25) is not unique; in fact, the crucial part is taking uhn and qhn from opposite sides. For the nonlinear term f^, it is applicable for taking any monotonic flux. According to equation (18), we obtain a truncation error of the above LDG scheme which is rnx.


### 3.2. Superdiffusion Case

In this part, we exploit a center difference scheme [24] to estimate the second-order time-fractional Caputo derivative ∂αux,t/∂tα in the nonlinear superdiffusion case 1<α<2; the discretization form at tn is as follows:(26)∂αux,tn∂tα=τ2−αΓ3−α∑i=0n−1νiux,tn−i−ux,tn−i−1τ+Rnx.Here,νi=i+12−α−i2−α,i=0,1,…,n−1, and truncation error satisfies the inequality Rnx≤CuT3−α/Γ2−ατ+τ2 [23], Cu is a constant which only depends on u.Analogous to the nonlinear subdiffusion equation, we also convert equation (2) into first-order system (21); combined with time discretization form equation (26), we obtain the weak formulation for nonlinear time-space fractional equation (2) at tn:(27)∫Ωux,tnvdx−α1∫Ωfux,tnvxdx+α1ε∫Ωqx,tnvxdx+α1∑j=1Nfux,tnvj+1/2−−fux,tnvj−1/2+−α1ε∑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 ∀v,w,z∈Vhk:(28)∫Ωuhnvdx−α1∫Ωfuhnvxdx+α1ε∫Ωqhnvxdx+α1∑j=1Nf^uhnvj+1/2−−f^uhnvj−1/2+−α1ε∑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−∫Ωu0xvdx=0,∫Ωuh0tvdx−∫Ωψxvdx=0,where α1=ταΓ3−α.The “hat” terms in equation (28) in the cell boundary terms from integration by parts are the so-called “numerical fluxes” [40], which are single-valued functions defined on the edges. By equation (27), we know that the truncation error of the fully discrete scheme (28) is Rnx.Remarks 1.
Originally, we assume thatu has compact support and restricts the problem to the bounded domain Ω. Thus, we impose homogeneous Dirichlet boundary conditions for u∉Ω.


## 3.1. Subdiffusion Case

We use the backward Euler scheme to discrete the time fractional derivative∂αux,t/∂tα of nonlinear subdiffusion 0<α<1 case (1) at tn which is as follows [38, 39]:(18)∂αux,tn∂tα=τ1−αΓ2−α∑i=0n−1biux,tn−i−ux,tn−i−1τ+rnx.Here,bi=i+11−α−i1−α satisfies properties as follows:(19)1=b0>b1>b2>⋯>bn>0,bn⟶0n⟶∞,∑i=1nbi−1−bi+bn=1,where rnx is the truncation error and satisfies the inequality rnx≤Cτ2−α [38], constant C depends on u,T,α.Now, we rewrite equation (1) as a first-order system:(20)∂αux,t∂tα+∂∂xfu−ε∂∂xq=0,(21)q−Δβ−2/2p=0,(22)p−ε∂∂xux,t=0.Combining with equation (18), we obtain the weak formulation of system (21) at tn:(23)∫Ωux,tnvdx−α0∫Ωfux,tnvxdx+α0ε∫Ωqx,tnvxdx+α0∑j=1Nfux,tnvj+1/2−−fux,tnvj−1/2+−α0ε∑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 numerical solutions of u⋅,tn,p⋅,tn,q⋅,tn, respectively. We define the fully discrete local discontinuous Galerkin (LDG) scheme as follows: find uhn,phn,qhn∈Vhk such that, for all v,w,z∈Vhk, we have(24)∫Ωuhnvdx−α0∫Ωfuhnvxdx+α0ε∫Ωqhnvxdx+α0∑j=1Nf^uhnvj+1/2−−f^uhnvj−1/2+−α0ε∑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,where α0=ταΓ2−α.The “hat” terms in equation (18) in the cell boundary terms from integration by parts are the so-called “numerical fluxes” [40], 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 of the LDG numerical scheme. Here, for the sake of concise, we take the same alternating numerical fluxes on the boundary as in [24], i.e.,(25)uhn^=uhn−,qhn^=qhn+.We remark that the choice for the fluxes (25) is not unique; in fact, the crucial part is taking uhn and qhn from opposite sides. For the nonlinear term f^, it is applicable for taking any monotonic flux. According to equation (18), we obtain a truncation error of the above LDG scheme which is rnx.


## 3.2. Superdiffusion Case

In this part, we exploit a center difference scheme [24] to estimate the second-order time-fractional Caputo derivative ∂αux,t/∂tα in the nonlinear superdiffusion case 1<α<2; the discretization form at tn is as follows:(26)∂αux,tn∂tα=τ2−αΓ3−α∑i=0n−1νiux,tn−i−ux,tn−i−1τ+Rnx.Here,νi=i+12−α−i2−α,i=0,1,…,n−1, and truncation error satisfies the inequality Rnx≤CuT3−α/Γ2−ατ+τ2 [23], Cu is a constant which only depends on u.Analogous to the nonlinear subdiffusion equation, we also convert equation (2) into first-order system (21); combined with time discretization form equation (26), we obtain the weak formulation for nonlinear time-space fractional equation (2) at tn:(27)∫Ωux,tnvdx−α1∫Ωfux,tnvxdx+α1ε∫Ωqx,tnvxdx+α1∑j=1Nfux,tnvj+1/2−−fux,tnvj−1/2+−α1ε∑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 ∀v,w,z∈Vhk:(28)∫Ωuhnvdx−α1∫Ωfuhnvxdx+α1ε∫Ωqhnvxdx+α1∑j=1Nf^uhnvj+1/2−−f^uhnvj−1/2+−α1ε∑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−∫Ωu0xvdx=0,∫Ωuh0tvdx−∫Ωψxvdx=0,where α1=ταΓ3−α.The “hat” terms in equation (28) in the cell boundary terms from integration by parts are the so-called “numerical fluxes” [40], which are single-valued functions defined on the edges. By equation (27), we know that the truncation error of the fully discrete scheme (28) is Rnx.Remarks 1.
Originally, we assume thatu has compact support and restricts the problem to the bounded domain Ω. Thus, we impose homogeneous Dirichlet boundary conditions for u∉Ω.


## 4. Stability

In this part, we prove the stability for the above two fully discrete LDG schemes (18) and (28).Theorem 1.
For periodic or compactly supported boundary conditions, the fully discrete LDG scheme (18) is unconditionally stable, and the numerical solution uhn satisfies(29)uhn≤uh0,n=1,2,…,M.Proof.
Firstly, in order to better analyze the fully discrete LDG numerical format, we introduce the following boundary term operatorℒ:(30)ℒv,w,z=εua,tz1/2+−εμhub,tvN+1/2−−εub,tzN+1/2−=0.
Using the fluxes (25) and the fluxes at the boundaries combining equality (30), after a simple rearrangement of terms, the fully LDG scheme (24) becomes(31)∫Ωuhnvdx−α0∫Ωfuhnvxdx+α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∑j=1N−1f^uhnj+1/2vj+1/2−α0f^uhn1/2v1/2+−f^uhnN+1/2vN+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 prove Theorem1 by using the mathematical induction. First, consider n=1; then, the LDG scheme (31) is(32)∫Ωuh1vdx−α0∫Ωfuh1vxdx+α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∑j=1N−1f^uh1j+1/2vj+1/2−α0f^uh11/2v1/2+−f^uh1N+1/2vN+1/2−+α0εqh11/2+v1/2+−qh1N+1/2−vN+1/2−+εμhuh1N+1/2−2=∫Ωuh0vdx.
DefineΦu=∫ufudu, and take the test functions as v=uh1,w=−α0ph1,z=α0qh1; plug them into equation (32), and we get(33)uh12+α0Δβ−2/2ph1,ph1+α0εμhuh1N+1/2−2+α0Φuh11/2−Φuh1N+1/2−f^uh11/2uh11/2++f^uh1N+1/2uh1N+1/2−+α0∑j=1N−1Φuh1j+1/2−f^uh1j+1/2=∫Ωuh0uh1dx,where we use the following results in the above equation:(34)∫Ωfuuxdx=∫Ωfudu=∫ΩΦ′udx=∑j=1NΦxxj−1/2+xj+1/2−=∑j=1NΦxj+1/2−−Φxj−1/2+=−∑j=1N−1Φuj+1/2−Φu1/2+ΦuN+1/2,α0ε∫Ωqh1uh1xdx+α0ε∫Ωuh1qh1xdx+α0ε∑j=1N−1qh1j+1/2+uh1j+1/2+α0ε∑j=1N−1uh1j+1/2−qh1j+1/2=α0ε∑j=1Nquxj−1/2+xj+1/2−+α0ε∑j=1N−1qj+1/2+uj+1/2+α0ε∑j=1N−1uj+1/2−qj+1/2=−α0εq1/2+u1/2++α0εqN+1/2−uN+1/2−−α0ε∑j=1N−1qj+1/2+uj+1/2−α0ε∑j=1N−1uj+1/2−qj+1/2+α0ε∑j=1N−1qj+1/2+uj+1/2+α0ε∑j=1N−1uj+1/2−qj+1/2=−α0εq1/2+u1/2+−qN+1/2−uN+1/2−.
From the properties of monotonic fluxes, we know that the functionf^u−,u+ is nondecreasing with its first argument and nonincreasing with its second argument. Hence, we have Φuj+1/2−fh^uhj+1/2>0, 1≤j≤N−1. Recall the orthogonality property of Galerkin and the following Cauchy–Schwartz inequality:(35)∫Ωuh0uh1dx≤12uh12+12uh02.
We can rewrite equation (33) as(36)uh12+α0Δβ−2/2ph1,ph1+α0ϵμhuh1N+1/2−2+α0Φuh11/2−Φuh1N+1/2−f^uh112uh112++f^uh1N+1/2uh1N+1/2−≤12uh12+12uh02.
By boundary conditions (periodic) and recalling Lemma6, we get(37)uh1≤uh0.
Next, we suppose that the following inequality holds:(38)uhm≤uh0,m=1,2,…,p.
Now, we need to proveuhp+1≤uh0. Let n=p+1, and take the test functions v=uhp+1,w=−α0php+1,andz=α0qhp+1 in equation (31):(39)uhp+12+α0Δβ−2/2php+1,php+1+α0εμhuhp+1N+1/2−2+α0Φuhp+11/2−Φuhp+1N+1/2−f^uhp+11/2uhp+11/2++f^uhp+1N+1/2uhp+1N+1/2−=∑i=1pbi−1−bi∫Ωuhp+1−iuhp+1dx+bp∫Ωuh0uhp+1dx.
Analogously, consider the properties of the monotone fluxf^ and boundary conditions again, and we get(40)uhp+12+α0Δβ−2/2php+1,php+1≤∑i=1pbi−1−biuhp+1−iuhp+1+bpuh0uhp+1≤∑i=1pbi−1−bi+bpuh0uhp+1≤uh0uhp+1.
Combining with Lemma6, we know that(41)α0Δβ−2/2php+1,php+1=α0php+1JLβ−2/2Ω2≥0.
Then, the last inequality in (40) gives(42)uhp+1≤uh0.Theorem 2.
For a sufficiently small step size,0<τ<T, the fully discrete LDG scheme (28) is unconditional stable, and the numerical solution uhn satisfies(43)uhn≤Cϕx+τψx,n=0,1,…,M.Proof.
Similar to the nonlinear subdiffusion case, by using the inequalityu−1≤Cϕx+τψx and adopting the same methods and techniques as in Theorem 1, we can obtain the conclusions of Theorem 2. For brevity, we ignore the details of proof process here.


## 5. Error Estimation

Theorem 3.
Letux,tn be the exact solution of nonlinear subdiffusion problem (1) which is sufficiently smooth with bounded derivatives, and assume that f∈C3; let uhn be numerical solutions of the LDG scheme (24), and we adopt eu=ux,tn−uhn to denote the corresponding numerical error; then, there hold the following error estimates:
When0<α<1,(44)ux,tn−uhn≤CuTα1−ατ−αhk+1+τ2−α+τ−α/2hk+1/2+hk+1,and when α⟶1,(45)ux,tn−uhn≤CuTτ−1hk+1+τ+τ−1/2hk+1/2+hk+1,where Cu is a constant which only depends on u.Proof.
Denote(46)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.
By subtracting equation (24) from equation (23) and combining with the fluxes on boundaries, we obtain the error equation:(47)∫Ωeunvdx+α0ε∫Ωeqnvxdx−∑j=1Neqn+vj+1/2−−eqn+vj−1/2+−α0∫Ωfux,tn−fuhnvxdx+α0∑j=1Nfux,tn−f^uhnvj+1/2−−α0∑j=1Nfux,tn−f^uhnvj−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,for all v,w,z∈Vhk.
Define(48)ℱf:u,uh:v=α0∫Ωfux,tn−fuhnvxdx−α0∑j=1Nfux,tn−f^uhnvj+1/2−+α0∑j=1Nfux,tn−f^uhnvj−1/2+.
Together with equations (46) and (30), after a simple rearrangement of terms, we can rewrite error equation (47) as(49)∫Ω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.
Next, we estimate the last term on the right-hand side of energy equation (49). We redefine it as in the following form [41]:(50)ℱf:u,uh:v=α0∫Ωfux,tn−fuhnvxdx+α0∑j=1Nfux,tn−fu¯hnvj+1/2+α0∑j=1Nfu¯hn−f^uhnvj+1/2,where u¯hn is an average value which is defined by u¯hn=1/2uhn++uhn−.
In order to obtain the estimation of the nonlinear term in equation (5), we need the following results:Lemma 9.
(see [42]). For each piecewise smooth function ω∈L2Ω, we define the following function on every cell boundary point:(51)κf^;ω≡κf^;ω−,ω+≜ω−1fω¯−f^ω,ifω≠0,12f′ω¯,ifω=0.Here,f^ω≡f^ω−,ω+ is a monotone numerical flux, which is consistent with the given flux function f. And κf^;ω is nonnegative and bounded for all ω−,ω+∈R. Moreover, we obtain(52)12f′ω¯≤κf^;ω+C∗ω,−18f″ω¯ω≤κf^;ω+C∗ω2,where the positive constant C∗ depends solely on the maximum of f″ and/or f‴.In order to estimate the nonlinear termfux,t, we need to assume that, for sufficiently small space step size h and for all k≥1, we prove that [41](53)eu=ux,t−uh≤h.Lemma 10.
For the operatorℱf:u,uh:v defined in equation (50), we have the following estimate:(54)ℱf:u,uh:v≤−14af^;uhv2τα+C+C∗v∞+h−1eu∞2vL2Ω2τα+C+C∗h−1eu∞2h2k+1τα.The proof of Lemma10 is similar to Xu and Shu’s work in [41]; here, we ignore the detail process of the proof.Take the test functions asv=P−eun,w=−α0Qepn,z=α0P+eqn, and plug them into equation (49); then, by using the properties (14)–(16), we get the following equality:(55)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,tn×Qepndx−α0∫ΩP+qx,tn−qx,tnQepndx+α0∫ΩQpx,tn−px,tn×P+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+ℱf:u,uh:v,i.e.,(56)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+ℱf:u,uh:v.Similarly, we analyze the error estimate by applying the mathematical induction method. First, we consider the case whenn=1, equation (56) becomes(57)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+ℱf:u,uh:v.Notice that(58)P−eu0≤Chk+1,ab≤εa2+14εb2.Combining with the standard approximation theorem (17), we obtain(59)Δβ−2/2Qpx,t1−px,t1−P+qx,t1−qx,t1=Δβ−2/2Qpx,t1−px,t1−△β−2/2P+px,t1−px,t1≤Chk+1.Combining Lemmas9 and 10 in [24] with Definition 2.4, equation (46), and Young’s inequality, we obtain(60)P−eu12+α0Δβ−2/2Qep1,Qep1+α0εμhP−eu1N+1/22≤P−eu0+α0r1x+P−ux,t1−ux,t1+P−ux,t0−ux,t0×P−eu1+Ch2k+2τα+α0εCQep12+α0εμhP−eu1N+1/22−14af^;uhv2τα+C+C∗v∞+h−1eu∞2vL2Ω2τα+C+C∗h−1eu∞2h2k+1τα.Thus,(61)P−eu12+α0Δβ−2/2Qep1,Qep1≤Chk+1+τ22+εP−eu12+Ch2k+2τα+α0εCQep12−14af^;uhv2τα+C+C∗v∞+h−1eu∞2vL2Ω2τα+C+C∗h−1eu∞2h2k+1τα.If we chooseε≤min1,1/C and recall the fractional Poincaré–Friedrichs, Lemmas 7 and 9 and inequality (53), the standard approximation theorem (17), and the positive property of af^;uh, we obtain(62)P−eu1≤b0−1Chk+1+τ2+τα/2hk+1/2,where C is a positive constant which depends on u,T,α.Next, we assume that the following inequality holds:(63)P−eum≤bm−1−1Chk+1+τ2+τα/2hk+1/2,m=1,2,…,l.Whenn=l+1, from equation (56) and Young’s inequality, we obtain(64)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Δβ−2/2Qpx,tl+1−px,tl+1−P+qx,tl+1−qx,tl+1,Qepl+1−α0εP+qx,tl+1−qx,tl+1−N+1/2P−eul+1N+1/2−−∑i=1lbi−1−bi∫ΩP−ux,tl+1−i−ux,tl+1−iP−eul+1dx−bl∫ΩP−ux,t0−ux,t0P−eul+1dx+ℱf:u,uh:v.≤∑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−14af^;uhv2τα+C+C∗v∞+h−1eu∞2vL2Ω2τα+C+C∗h−1eu∞2h2k+1τα.Because(65)bi−1−1<bi−1,τα/2hk+1>0,so(66)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−14af^;uhv2τα+C+C∗v∞+h−1eu∞2vL2Ω2τα+C+C∗h−1eu∞2h2k+1τα≤∑i=1lbi−1−bi+blbl−1Chk+1+τ2+τα/2hk+1×P−eul+1+Ch2k+2τα+α0εCQepl+12+α0εμhP−eul+1N+1/22−14af^;uhv2τα+C+C∗v∞+h−1eu∞2vL2Ω2τα+C+C∗h−1eu∞2h2k+1τα.Analogously, apply Young’s inequality to equation(66) and choose sufficiently small ϵ; by using the fractional Poincaré–Friedrichs, Lemmas 7 and 9 and inequality (53), the standard approximation theorem (17), and the positive property of af^;uh again, we obtain(67)P−eul+1≤bl−1Chk+1+τ2+τα/2hk+1/2,where C is a positive constant which depends on u,T,α.According to the mathematical induction method, we get(68)P−eun≤bn−1−1Chk+1+τ2+τα/2hk+1/2.From [39], we know that the term n−αbn−1−1 monotonic increase tends to 1/1−αn⟶∞; then,(69)P−eun≤bn−1−1Chk+1+τ2+τα/2hk+1/2≤nαn−αbn−1−1Chk+1+τ2+τα/2hk+1/2≤Tτα11−αChk+1+τ2+τα/2hk+1/2≤CuTα1−ατ−αhk+1+τ2−α+τ−α/2hk+1/2,where Cu is a constant which only depends on u.Whenα⟶1, 1/1−α⟶∞, estimation formula (69) is meaningless, so we should re-estimate the case for α⟶1.Again, we can prove the following estimator by using mathematical induction:(70)P−eun≤nChk+1+τ2+τα/2hk+1/2.The proof techniques are almost the same as the subdiffusion case in [24], and in order to save space, here, we ignore the concrete proof; thus, when α⟶1, we obtain(71)P−eun≤nChk+1+τ2+τα/2hk+1/2≤TCuτ−1hk+1+τ+τ−1/2hk+1/2,where Cu is a constant which only depends on u.According to inequalities (69) and (71), the triangle inequality, and the standard approximation theorem (17), coupling with the property of the projection operator, we can prove that the results of Theorem 3 hold.Theorem 4.
Letux,tn be the exact solution of nonlinear superdiffusion problem (2), which is sufficiently smooth with bounded derivatives, and assume that f∈C3, let uhn denote the numerical solution of the LDG scheme (28), and use eu=ux,tn−uhn to denote the corresponding numerical error; then, the following error estimates hold:
When1<α<2,(72)ux,tn−uhn≤CuTα−12−ατ1−αhk+1+τ3−α+τ1−α/2hk+1/2+hk+1,and when α⟶2,(73)ux,tn−uhn≤CuTτ−1hk+1+τ+hk+1/2+hk+1,where Cu is a constant which only depends on u.Proof.
Subtracting equation (28) from equation (27) and combining with the fluxes on boundaries, we obtain the following error equation:(74)∫Ωeunvdx+α1ε∫Ωeqnvxdx−∑j=1Neqn+vj+1/2−−eqn+vj−1/2+−α1∫Ωfux,tn−fuhnvxdx+α1∑j=1Nfux,tn−f^uhnvj+1/2−−α1∑j=1Nfux,tn−f^uhnvj−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 the nonlinear subdiffusion case, together with the following error estimator(75)eu−1≤Chk+1+τ2+τhk+1,where C is a constant related to u,T,α [24], we can prove that the following error inequality holds by using mathematical induction:(76)P−eun≤νn−1−1Chk+1+τ2+τα/2hk+1/2.
According to Qiu et al. [24], we know that n1−ανn−1−1 monotone increasing tends to 1/2−αn⟶∞, so(77)P−eun≤νn−1−1Chk+1+τ2+τα/2hk+1/2≤nα−1n1−ανn−1−1Chk+1+τ2+τα/2hk+1/2≤Tτα−112−αChk+1+τ2+τα/2hk+1/2≤CuTα−12−ατ1−αhk+1+τ3−α+τ1−α/2hk+1/2,where Cu is a constant which only depends on u.
Whenα⟶2 and 1/2−α⟶∞, the estimator (77) is insignificant; hence, we should re-estimate the case for α⟶2.
Similar to the nonlinear subdiffusion case, we still use the mathematical induction to prove the following estimator:(78)P−eun≤nChk+1+τ2+τα/2hk+1/2.
The proof techniques are almost the same as the superdiffusion case in [24]; for simplicity, here, we ignore the proof process; thus, when α⟶2, we obtain(79)P−eun≤nChk+1+τ2+τα/2hk+1/2≤TCuτ−1hk+1+τ+hk+1/2,where Cu is a constant which only depends on u.
According to inequalities (77) and (79), the triangle inequality, and the standard approximation theorem (17), coupling with the property of the projection operator, we can prove that the results of Theorem 4 hold.


## 6. Implementation Procedure

In this part, we give out several examples to demonstrate the effectiveness of the finite difference/LDG approximation for solving nonlinear time-space fractional subdiffusion/superdiffusion equations. At the same time, we propose the convergence behavior of numerical solutions with respect to the time step sizeτ and/or the space step size h, i.e., through carefully choosing the appropriate time step size, and let the Lax–Friedrichs numerical flux be a numerical flux function of the nonlinear terms f^; we observe that a convergence rate can arrive at Oτ2−α+hk+1 when 0<α<1 and approach to Oτ3−α+hk+1 when the derivative parameter is 1<α<2.


### 6.1. Numerical Examples

As examples, we consider the following nonlinear time-space fractional subdiffusion/superdiffusion equations:(80)∂αux,t∂tα+ux,t∂ux,t∂x−axε−−Δβ/2ux,t=hx,t,x∈Ω,t∈0,T,ux,0=u0x,x∈Ω,ux,t=0,x∈∂Ω,t∈0,T,where Ω=a,b, fu=1/2u2, and T>0; when 0<α<1, it is the nonlinear subdiffusion case, and when the derivative parameter is 1<α<2, it is the nonlinear superdiffusion case, and 1<β<2, ϵ>0.For brevity, we only derive the linear system obtained from the fully discrete LDG scheme (24). Similarly, we can get the corresponding linear system of scheme (28).We consider the case of discontinuous piecewise polynomial basis function sequencesϕixi=1k in the space Vkh. By the definition of the space Vkh, we know that, for all v∈Vkh, there holds vx=∑i=1Nviϕix, where vi=vxi.From the LDG scheme (24), we obtain the fully discrete scheme of equation (80) as follows:(81)∫Ωuhnvdx−α0∫Ωfuhnvxdx+α0∑j=1Nf^uhnvj+1/2−−f^uhnvj−1/2++α0εax∫Ωqhnvxdx−∑j=1Nqhn^vj+1/2−−qhn^vj−1/2+=∑i=1n−1bi−1−bi∫Ωuhn−ivdx+bn−1∫Ωuh0vdx+∫Ωhx,tvdx,∫Ωqhnwdx−∫ΩΔβ−2/2phnwdx=0,∫Ωphnzdx+εax∫Ωuhnzxdx−∑j=1Nuhn^zj+1/2−−uhn^zj−1/2+=0,∫Ωuh0vdx−∫Ωu0xvdx=0,where α0=ταΓ2−α.In terms of the basisϕixi=1N, we denote approximate solutions by uhn=∑i=1Nujtϕjx, phn=∑i=1Npjtϕjx, and qhn=∑i=1Nqjtϕjx, and take test functions as v=∑i=1Nviϕix,vi=vxi, w=∑i=1Nwiϕix,wi=wxi, and z=∑i=1Nziϕix,zi=zxi. Plug these expressions in the LDG scheme (81) and take flux equation (25); we get the linear system of the scheme in the jth element at tn as follows:(82)M+α0B0α0εA−H11+H120RMεA−H13+H14M0UjnPjnQjn=Fnv+∑i=1n−1bi−1−biUn−iv+bn−1U0v,where B represents the combined part of nonlinear terms; here, we select the Lax–Friedrichs flux as a numerical flux function of nonlinear term f^uhn, M=ϕmjx,ϕkjxj=1,…,N is the mass matrix, m,k=0,1,…, denote the polynomials order, andA=ϕmjx,ϕkj′xj=1,…,N,(83)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 a quadrature for fractional integral operator △β−2/2phn and test function w in the jth element. Recalling the works in [34], we know that(84)Δβ−2/2phn=−Iax2−βphn+Ixb2−βphn2cosβπ/2,where(85)Iax2−βphnx,t=1Γ2−β∫axx−ξ1−βphnξ,tdξ,Ixb2−βphnx,t=1Γ2−β∫xbξ−x1−βphnξ,tdξ,and Fn=h1,h2,…,hnT is a vector-valued function, hjtn=hxj,tn. Then, by solving linear system (82), we can solve the solution Un. This study uses MATLAB programs to execute the mathematical calculations involved in the proposed numerical method.


### 6.2. Numerical Results

Example 1.
Consider the nonlinear fractional time-space subdiffusion0<α<1 equation (80), and satisfy initial boundary conditions:(86)ux,0=x21−x,x∈0,1,u0,t=u1,t=0,t>0.
In order to validate the optimal convergence of the proposed method, we selectax=Γ2.8x/2,ε=1 and the source\term as(87)hx,t=2x21−xt1.2Γ2.2+0.2x1.81+t2.
Thus, equation (80) has an exact solution ux,t=x21−x1+t2.
By choosing the fluxes in equation (25) for the linear term and Lax–Friedrichs flux for the nonlinear term f^ and adopting the linear piecewise polynomials, then we numerically solve equation (80). We investigate the L2 errors and the corresponding convergence rate of the solutions with respect to different fractional-order parameters α and β. First, we select 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 1 and 2 list out the L2 errors and the corresponding convergence order, respectively, which is fixed the value of one parameter β (or α), with the values change of the other parameter α (or β). As can be seen from the data in the two tables, our schemes can achieve the optimal convergence Ohk+1 in spatial direction; this matches the theoretical results of Theorem 3.
The numerical data 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, fix the values of one derivative parameter β (or α), and then change the values of the other derivative constant α (or β), when 0<α<1, the accuracy in time direction of the nonlinear time-space fractional subdiffusion equation can achieve to Oτ2−α under the L2 norm and can reach to 1 when α⟶1; this experiment result agrees with the theoretical computation results in Theorem 3.
Although an order of convergenceOhk+1 is predicted by Theorem 3 in the spatial direction for any monotone flux of the nonlinear term, however, optimal order of superconvergence is observed when we choose a Lax–Friedrichs flux for the nonlinear term f^ and an alternating direction flux for the linear term. We do not understand the reason for this improved rate at present, but it may be associated with the “global dependency” of the fractional operator which exhibits optimal convergence.Table 1
The spatial convergence rate can achieve toOhk+1; here, β=1.6 and T=1.
hα=0.1Rateα=0.6Rateα=0.99Rate1/47.1268 × 10−2—6.3749 × 10−2—6.7888 × 10−2—1/81.9575 × 10−21.86421.6376 × 10−21.96081.8660 × 10−21.86321/165.2918 × 10−31.88724.1156 × 10−31.99244.9554 × 10−31.91291/321.3757 × 10−31.94361.0279 × 10−32.00131.2498 × 10−31.98731/643.4718 × 10−41.98642.5439 × 10−42.01473.1373 × 10−41.9941Table 2
The spatial convergence rate can achieve toOhk+1; here, α=0.3 and T=1.
hβ=1.01Rateβ=1.48Rateβ=1.97Rate1/45.6314 × 10−2—6.3326 × 10−2—4.7829 × 10−2—1/81.4588 × 10−21.94871.6475 × 10−21.94251.3025 × 10−21.87661/163.7356 × 10−31.96154.2249 × 10−31.96333.2786 × 10−31.99011/329.4929 × 10−41.98031.0590 × 10−31.99628.2262 × 10−41.99481/642.3761 × 10−41.99822.6413 × 10−42.00342.0523 × 10−42.0003Table 3
The temporal convergence rate can arrive atOτ2−α; β=1.1 and T=1.
τα=0.03Rateα=0.7Rateα=0.99Rate1/46.0622 × 10−3—6.3740 × 10−3—4.0562 × 10−3—1/81.7016 × 10−31.83292.7488 × 10−31.21342.2067 × 10−30.87821/164.5762 × 10−41.89471.1518 × 10−31.25491.1845 × 10−30.89761/321.1936 × 10−41.93884.6859 × 10−41.29756.1069 × 10−40.95581/643.0484 × 10−51.96921.9025 × 10−41.30043.0668 × 10−40.9973Table 4
The temporal convergence rate can arrive atOτ2−α; α=0.9 and T=1.
τβ=1.1Rateβ=1.6Rateβ=1.98Rate1/45.3762 × 10−2—4.1429 × 10−2—3.2679 × 10−2—1/82.6873 × 10−21.00042.0677 × 10−21.00261.6281 × 10−21.00511/161.3400 × 10−21.00291.0239 × 10−21.01398.0914 × 10−31.00881/326.3922 × 10−31.06794.8191 × 10−31.08733.9726 × 10−31.02631/642.9818 × 10−31.10012.2427 × 10−31.10351.8549 × 10−31.0987Example 2.
Consider the following nonlinear time-space fractional subdiffusion0<α<1 equation with Dirichlet boundary conditions:(88)∂αux,t∂tα+ux,t∂ux,t∂x+cxux,tε−−△β2ux,t,inx∈Ω,t∈0,T,ux,0=0,inx∈Ω,u−1,t=u1,t=0,inx∈∂Ω,t∈0,T.
Here,Ω=−1,1,T=1,cx=x2+1, ϵ=0.5; the equation has an exact solution ux,t=sintCx+12,2π. Here,(89)Cx+nγ,ω≜x+nγ∑j=0∞−1jωx+n2jΓγ+2j+1.
In order to complete the scheme, we choose an alternating direction flux for the linear term and Lax–Friedrichs flux for the nonlinear term. Tables5 and 6 reveal the L2 norm and the convergence rate, respectively, in which one of the parameters β (or α) is fixed and the other parameter α (or β) is changed, with an appropriate time step size τ. Analogously, if we choose an appropriate space step size h and fix one parameter β (or α) and then change the other constant α (or β), the results in Tables 7 and 8 show that when 0<α<1, the convergence rate in time direction of the nonlinear time-space fractional subdiffusion equation can converge to Oτ2−α under the L2 norm, and the rate can be close to 1 when α⟶1; this result agrees with the theoretical analysis conclusion in Theorem 3.Table 5
The spatial convergence rate can achieve toOhk+1; β=1.6 and T=1.
hα=0.1Rateα=0.5Rateα=0.9Rate1/41.7936 × 10−2—3.4408 × 10−2—5.3221 × 10−2—1/84.7607 × 10−31.91368.6408 × 10−31.99351.4560 × 10−21.86991/161.2281 × 10−31.95472.1630 × 10−31.99813.7084 × 10−31.97321/323.1192 × 10−41.97725.4001 × 10−42.00209.2814 × 10−41.99841/647.8041 × 10−51.99891.2619 × 10−42.09732.2497 × 10−42.0446Table 6
The spatial convergence rate can achieve toOhk+1; α=0.5 and T=1.
hβ=1.02Rateβ=1.4Rateβ=1.9Rate1/44.5016 × 10−2—5.1678 × 10−3—8.4526 × 10−3—1/81.1740 × 10−21.93891.3028 × 10−31.98792.1608 × 10−31.96781/163.0014 × 10−31.96783.2606 × 10−41.99845.4081 × 10−41.99841/327.2299 × 10−42.05367.7774 × 10−52.06781.3516 × 10−42.00041/641.6926 × 10−42.09471.8152 × 10−52.09913.9256 × 10−52.0361Table 7
The temporal convergence rate can arrive atOτ2−α; β=1.6 and T=1.
τα=0.1Rateα=0.7Rateα=0.97Rate1/46.0398 × 10−3—5.5331 × 10−3—8.9365 × 10−3—1/81.6957 × 10−31.83262.3184 × 10−31.25494.8416 × 10−30.88421/164.6458 × 10−41.86799.5726 × 10−41.27622.5280 × 10−30.93751/321.2436 × 10−41.90143.8866 × 10−41.30041.2411 × 10−31.02631/643.3314 × 10−51.90031.5764 × 10−41.30186.1737 × 10−41.0075Table 8
The temporal convergence rate can arrive atOτ2−α; α=0.1 and T=1.
τβ=1.2Rateβ=1.6Rateβ=1.99Rate1/43.2985 × 10−3—7.0233 × 10−3—6.4021 × 10−3—1/88.9844 × 10−41.87631.8923 × 10−31.89201.7360 × 10−31.88271/162.4326 × 10−41.88495.1141 × 10−41.88764.6633 × 10−41.89641/326.5050×10−51.90291.3624 × 10−41.90831.2490 × 10−41.90051/641.7311 × 10−51.90983.5915 × 10−51.92353.3368 × 10−51.9043Example 3.
In this example, we consider the nonlinear time-space fractional superdiffusion case1<α<2 with the nonhomogeneous boundary condition:(90)∂αux,t∂tα+ux,t∂ux,t∂x−axε−−Δβ/2ux,t=hx,t,x∈Ω,t∈0,T,ux,0=x3,utx,0=−4x3,x∈Ω,u0,t=0,u1,t=2t−12,x∈∂Ω,where Ω=0,1,T=1,1<α<2,1<β<2, ax=1,ϵ=1; we choose the exact solution of the equation as ux,t=2t−12x3, and hx,t is from equation (90). Similar to the nonlinear subdiffusion model, we also choose alternating direction flux for the linear term and Lax–Friedrichs flux for the nonlinear term and select the appropriate time step size τ such that the time discretization errors are negligible as compared to the spatial errors. The L2 error and convergence rate in spatial direction are indicated in Tables 9 and 10, in which one of the parameters β (or α) is fixed and the other parameter α (or β) is changed. As can be seen from the data in tables, both schemes have optimal convergence rate Ohk+1 in the spatial direction; this conclusion is consistent with the theoretical results predicted in Theorem 4. When h=C′τ (C′), we choose the time step size sequence as τ=1/2i,i=2,⋯,6 and fix β (or α), and then change α (or β); the results in Tables 11 and 12 show that when 1<α<2, the convergence rate in time direction for the time-space fractional superdiffusion can reach to Oτ3−α and tends to 1 when α⟶2. These results are consistent with the theoretical consequence proved in Theorem 4.Table 9
The spatial convergence rate can achieve toOhk+1; β=1.6 and T=1.
hα=1.05Rateα=1.6Rateα=1.999Rate1/43.8952 × 10−2—7.9493 × 10−2—5.0248 × 10−2—1/89.9577 × 10−31.96782.1840 × 10−21.86381.3661 × 10−21.87891/162.5203 × 10−31.98225.7575 × 10−31.92353.5435 × 10−31.94691/326.2991 × 10−42.00041.4499 × 10−31.98958.9173 × 10−41.99051/641.5363 × 10−42.03563.5356 × 10−42.03592.2245 × 10−42.0031Table 10
The spatial convergence rate can achieve toOhk+1; α=1.3 and T=1.
hβ=1.05Rateβ=1.5Rateβ=1.9Rate1/46.0285 × 10−2—4.9039 × 10−2—6.7648 × 10−2—1/81.6179 × 10−21.89761.3190 × 10−21.89441.7680 × 10−21.93591/164.2071 × 10−31.94333.4816 × 10−31.92174.5017 × 10−31.97361/321.0652 × 10−31.98168.8475 × 10−41.97641.1318 × 10−31.99181/642.6718 × 10−41.99532.2166 × 10−41.99692.6435 × 10−42.0981Table 11
The temporal convergence rate can reach toOτ3−α; β=1.5, C′=0.5,andT=1.
τα=1.1Rateα=1.3Rateα=1.99Rate1/43.7566 × 10−2—3.2989 × 10−2—3.0408 × 10−3—1/81.0324 × 10−21.86341.0920 × 10−21.59491.6243 × 10−30.90461/162.6767 × 10−31.94753.5105 × 10−31.63738.3049 × 10−40.96781/326.7308 × 10−41.99161.0781 × 10−31.70324.2689 × 10−40.96011/641.6666 × 10−42.01383.2889 × 10−41.71282.1515 × 10−40.9885Table 12
The temporal convergence rate can reach toOτ3−α; α=1.8, C′=0.6,andT=1.
τβ=1.15Rateβ=1.49Rateβ=1.999Rate1/42.2799 × 10−3—3.6088 × 10−3—6.5082 × 10−3—1/81.0716 × 10−31.08911.6547 × 10−31.12493.0294 × 10−31.10321/164.9550 × 10−41.11297.2342 × 10−41.19371.3576 × 10−31.15791/322.1798 × 10−41.18473.1736 × 10−41.18875.9499 × 10−41.19021/649.4494 × 10−51.20591.3783 × 10−41.20322.5780 × 10−41.2066


## 6.1. Numerical Examples

As examples, we consider the following nonlinear time-space fractional subdiffusion/superdiffusion equations:(80)∂αux,t∂tα+ux,t∂ux,t∂x−axε−−Δβ/2ux,t=hx,t,x∈Ω,t∈0,T,ux,0=u0x,x∈Ω,ux,t=0,x∈∂Ω,t∈0,T,where Ω=a,b, fu=1/2u2, and T>0; when 0<α<1, it is the nonlinear subdiffusion case, and when the derivative parameter is 1<α<2, it is the nonlinear superdiffusion case, and 1<β<2, ϵ>0.For brevity, we only derive the linear system obtained from the fully discrete LDG scheme (24). Similarly, we can get the corresponding linear system of scheme (28).We consider the case of discontinuous piecewise polynomial basis function sequencesϕixi=1k in the space Vkh. By the definition of the space Vkh, we know that, for all v∈Vkh, there holds vx=∑i=1Nviϕix, where vi=vxi.From the LDG scheme (24), we obtain the fully discrete scheme of equation (80) as follows:(81)∫Ωuhnvdx−α0∫Ωfuhnvxdx+α0∑j=1Nf^uhnvj+1/2−−f^uhnvj−1/2++α0εax∫Ωqhnvxdx−∑j=1Nqhn^vj+1/2−−qhn^vj−1/2+=∑i=1n−1bi−1−bi∫Ωuhn−ivdx+bn−1∫Ωuh0vdx+∫Ωhx,tvdx,∫Ωqhnwdx−∫ΩΔβ−2/2phnwdx=0,∫Ωphnzdx+εax∫Ωuhnzxdx−∑j=1Nuhn^zj+1/2−−uhn^zj−1/2+=0,∫Ωuh0vdx−∫Ωu0xvdx=0,where α0=ταΓ2−α.In terms of the basisϕixi=1N, we denote approximate solutions by uhn=∑i=1Nujtϕjx, phn=∑i=1Npjtϕjx, and qhn=∑i=1Nqjtϕjx, and take test functions as v=∑i=1Nviϕix,vi=vxi, w=∑i=1Nwiϕix,wi=wxi, and z=∑i=1Nziϕix,zi=zxi. Plug these expressions in the LDG scheme (81) and take flux equation (25); we get the linear system of the scheme in the jth element at tn as follows:(82)M+α0B0α0εA−H11+H120RMεA−H13+H14M0UjnPjnQjn=Fnv+∑i=1n−1bi−1−biUn−iv+bn−1U0v,where B represents the combined part of nonlinear terms; here, we select the Lax–Friedrichs flux as a numerical flux function of nonlinear term f^uhn, M=ϕmjx,ϕkjxj=1,…,N is the mass matrix, m,k=0,1,…, denote the polynomials order, andA=ϕmjx,ϕkj′xj=1,…,N,(83)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 a quadrature for fractional integral operator △β−2/2phn and test function w in the jth element. Recalling the works in [34], we know that(84)Δβ−2/2phn=−Iax2−βphn+Ixb2−βphn2cosβπ/2,where(85)Iax2−βphnx,t=1Γ2−β∫axx−ξ1−βphnξ,tdξ,Ixb2−βphnx,t=1Γ2−β∫xbξ−x1−βphnξ,tdξ,and Fn=h1,h2,…,hnT is a vector-valued function, hjtn=hxj,tn. Then, by solving linear system (82), we can solve the solution Un. This study uses MATLAB programs to execute the mathematical calculations involved in the proposed numerical method.


## 6.2. Numerical Results

Example 1.
Consider the nonlinear fractional time-space subdiffusion0<α<1 equation (80), and satisfy initial boundary conditions:(86)ux,0=x21−x,x∈0,1,u0,t=u1,t=0,t>0.
In order to validate the optimal convergence of the proposed method, we selectax=Γ2.8x/2,ε=1 and the source\term as(87)hx,t=2x21−xt1.2Γ2.2+0.2x1.81+t2.
Thus, equation (80) has an exact solution ux,t=x21−x1+t2.
By choosing the fluxes in equation (25) for the linear term and Lax–Friedrichs flux for the nonlinear term f^ and adopting the linear piecewise polynomials, then we numerically solve equation (80). We investigate the L2 errors and the corresponding convergence rate of the solutions with respect to different fractional-order parameters α and β. First, we select 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 1 and 2 list out the L2 errors and the corresponding convergence order, respectively, which is fixed the value of one parameter β (or α), with the values change of the other parameter α (or β). As can be seen from the data in the two tables, our schemes can achieve the optimal convergence Ohk+1 in spatial direction; this matches the theoretical results of Theorem 3.
The numerical data 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, fix the values of one derivative parameter β (or α), and then change the values of the other derivative constant α (or β), when 0<α<1, the accuracy in time direction of the nonlinear time-space fractional subdiffusion equation can achieve to Oτ2−α under the L2 norm and can reach to 1 when α⟶1; this experiment result agrees with the theoretical computation results in Theorem 3.
Although an order of convergenceOhk+1 is predicted by Theorem 3 in the spatial direction for any monotone flux of the nonlinear term, however, optimal order of superconvergence is observed when we choose a Lax–Friedrichs flux for the nonlinear term f^ and an alternating direction flux for the linear term. We do not understand the reason for this improved rate at present, but it may be associated with the “global dependency” of the fractional operator which exhibits optimal convergence.Table 1
The spatial convergence rate can achieve toOhk+1; here, β=1.6 and T=1.
hα=0.1Rateα=0.6Rateα=0.99Rate1/47.1268 × 10−2—6.3749 × 10−2—6.7888 × 10−2—1/81.9575 × 10−21.86421.6376 × 10−21.96081.8660 × 10−21.86321/165.2918 × 10−31.88724.1156 × 10−31.99244.9554 × 10−31.91291/321.3757 × 10−31.94361.0279 × 10−32.00131.2498 × 10−31.98731/643.4718 × 10−41.98642.5439 × 10−42.01473.1373 × 10−41.9941Table 2
The spatial convergence rate can achieve toOhk+1; here, α=0.3 and T=1.
hβ=1.01Rateβ=1.48Rateβ=1.97Rate1/45.6314 × 10−2—6.3326 × 10−2—4.7829 × 10−2—1/81.4588 × 10−21.94871.6475 × 10−21.94251.3025 × 10−21.87661/163.7356 × 10−31.96154.2249 × 10−31.96333.2786 × 10−31.99011/329.4929 × 10−41.98031.0590 × 10−31.99628.2262 × 10−41.99481/642.3761 × 10−41.99822.6413 × 10−42.00342.0523 × 10−42.0003Table 3
The temporal convergence rate can arrive atOτ2−α; β=1.1 and T=1.
τα=0.03Rateα=0.7Rateα=0.99Rate1/46.0622 × 10−3—6.3740 × 10−3—4.0562 × 10−3—1/81.7016 × 10−31.83292.7488 × 10−31.21342.2067 × 10−30.87821/164.5762 × 10−41.89471.1518 × 10−31.25491.1845 × 10−30.89761/321.1936 × 10−41.93884.6859 × 10−41.29756.1069 × 10−40.95581/643.0484 × 10−51.96921.9025 × 10−41.30043.0668 × 10−40.9973Table 4
The temporal convergence rate can arrive atOτ2−α; α=0.9 and T=1.
τβ=1.1Rateβ=1.6Rateβ=1.98Rate1/45.3762 × 10−2—4.1429 × 10−2—3.2679 × 10−2—1/82.6873 × 10−21.00042.0677 × 10−21.00261.6281 × 10−21.00511/161.3400 × 10−21.00291.0239 × 10−21.01398.0914 × 10−31.00881/326.3922 × 10−31.06794.8191 × 10−31.08733.9726 × 10−31.02631/642.9818 × 10−31.10012.2427 × 10−31.10351.8549 × 10−31.0987Example 2.
Consider the following nonlinear time-space fractional subdiffusion0<α<1 equation with Dirichlet boundary conditions:(88)∂αux,t∂tα+ux,t∂ux,t∂x+cxux,tε−−△β2ux,t,inx∈Ω,t∈0,T,ux,0=0,inx∈Ω,u−1,t=u1,t=0,inx∈∂Ω,t∈0,T.
Here,Ω=−1,1,T=1,cx=x2+1, ϵ=0.5; the equation has an exact solution ux,t=sintCx+12,2π. Here,(89)Cx+nγ,ω≜x+nγ∑j=0∞−1jωx+n2jΓγ+2j+1.
In order to complete the scheme, we choose an alternating direction flux for the linear term and Lax–Friedrichs flux for the nonlinear term. Tables5 and 6 reveal the L2 norm and the convergence rate, respectively, in which one of the parameters β (or α) is fixed and the other parameter α (or β) is changed, with an appropriate time step size τ. Analogously, if we choose an appropriate space step size h and fix one parameter β (or α) and then change the other constant α (or β), the results in Tables 7 and 8 show that when 0<α<1, the convergence rate in time direction of the nonlinear time-space fractional subdiffusion equation can converge to Oτ2−α under the L2 norm, and the rate can be close to 1 when α⟶1; this result agrees with the theoretical analysis conclusion in Theorem 3.Table 5
The spatial convergence rate can achieve toOhk+1; β=1.6 and T=1.
hα=0.1Rateα=0.5Rateα=0.9Rate1/41.7936 × 10−2—3.4408 × 10−2—5.3221 × 10−2—1/84.7607 × 10−31.91368.6408 × 10−31.99351.4560 × 10−21.86991/161.2281 × 10−31.95472.1630 × 10−31.99813.7084 × 10−31.97321/323.1192 × 10−41.97725.4001 × 10−42.00209.2814 × 10−41.99841/647.8041 × 10−51.99891.2619 × 10−42.09732.2497 × 10−42.0446Table 6
The spatial convergence rate can achieve toOhk+1; α=0.5 and T=1.
hβ=1.02Rateβ=1.4Rateβ=1.9Rate1/44.5016 × 10−2—5.1678 × 10−3—8.4526 × 10−3—1/81.1740 × 10−21.93891.3028 × 10−31.98792.1608 × 10−31.96781/163.0014 × 10−31.96783.2606 × 10−41.99845.4081 × 10−41.99841/327.2299 × 10−42.05367.7774 × 10−52.06781.3516 × 10−42.00041/641.6926 × 10−42.09471.8152 × 10−52.09913.9256 × 10−52.0361Table 7
The temporal convergence rate can arrive atOτ2−α; β=1.6 and T=1.
τα=0.1Rateα=0.7Rateα=0.97Rate1/46.0398 × 10−3—5.5331 × 10−3—8.9365 × 10−3—1/81.6957 × 10−31.83262.3184 × 10−31.25494.8416 × 10−30.88421/164.6458 × 10−41.86799.5726 × 10−41.27622.5280 × 10−30.93751/321.2436 × 10−41.90143.8866 × 10−41.30041.2411 × 10−31.02631/643.3314 × 10−51.90031.5764 × 10−41.30186.1737 × 10−41.0075Table 8
The temporal convergence rate can arrive atOτ2−α; α=0.1 and T=1.
τβ=1.2Rateβ=1.6Rateβ=1.99Rate1/43.2985 × 10−3—7.0233 × 10−3—6.4021 × 10−3—1/88.9844 × 10−41.87631.8923 × 10−31.89201.7360 × 10−31.88271/162.4326 × 10−41.88495.1141 × 10−41.88764.6633 × 10−41.89641/326.5050×10−51.90291.3624 × 10−41.90831.2490 × 10−41.90051/641.7311 × 10−51.90983.5915 × 10−51.92353.3368 × 10−51.9043Example 3.
In this example, we consider the nonlinear time-space fractional superdiffusion case1<α<2 with the nonhomogeneous boundary condition:(90)∂αux,t∂tα+ux,t∂ux,t∂x−axε−−Δβ/2ux,t=hx,t,x∈Ω,t∈0,T,ux,0=x3,utx,0=−4x3,x∈Ω,u0,t=0,u1,t=2t−12,x∈∂Ω,where Ω=0,1,T=1,1<α<2,1<β<2, ax=1,ϵ=1; we choose the exact solution of the equation as ux,t=2t−12x3, and hx,t is from equation (90). Similar to the nonlinear subdiffusion model, we also choose alternating direction flux for the linear term and Lax–Friedrichs flux for the nonlinear term and select the appropriate time step size τ such that the time discretization errors are negligible as compared to the spatial errors. The L2 error and convergence rate in spatial direction are indicated in Tables 9 and 10, in which one of the parameters β (or α) is fixed and the other parameter α (or β) is changed. As can be seen from the data in tables, both schemes have optimal convergence rate Ohk+1 in the spatial direction; this conclusion is consistent with the theoretical results predicted in Theorem 4. When h=C′τ (C′), we choose the time step size sequence as τ=1/2i,i=2,⋯,6 and fix β (or α), and then change α (or β); the results in Tables 11 and 12 show that when 1<α<2, the convergence rate in time direction for the time-space fractional superdiffusion can reach to Oτ3−α and tends to 1 when α⟶2. These results are consistent with the theoretical consequence proved in Theorem 4.Table 9
The spatial convergence rate can achieve toOhk+1; β=1.6 and T=1.
hα=1.05Rateα=1.6Rateα=1.999Rate1/43.8952 × 10−2—7.9493 × 10−2—5.0248 × 10−2—1/89.9577 × 10−31.96782.1840 × 10−21.86381.3661 × 10−21.87891/162.5203 × 10−31.98225.7575 × 10−31.92353.5435 × 10−31.94691/326.2991 × 10−42.00041.4499 × 10−31.98958.9173 × 10−41.99051/641.5363 × 10−42.03563.5356 × 10−42.03592.2245 × 10−42.0031Table 10
The spatial convergence rate can achieve toOhk+1; α=1.3 and T=1.
hβ=1.05Rateβ=1.5Rateβ=1.9Rate1/46.0285 × 10−2—4.9039 × 10−2—6.7648 × 10−2—1/81.6179 × 10−21.89761.3190 × 10−21.89441.7680 × 10−21.93591/164.2071 × 10−31.94333.4816 × 10−31.92174.5017 × 10−31.97361/321.0652 × 10−31.98168.8475 × 10−41.97641.1318 × 10−31.99181/642.6718 × 10−41.99532.2166 × 10−41.99692.6435 × 10−42.0981Table 11
The temporal convergence rate can reach toOτ3−α; β=1.5, C′=0.5,andT=1.
τα=1.1Rateα=1.3Rateα=1.99Rate1/43.7566 × 10−2—3.2989 × 10−2—3.0408 × 10−3—1/81.0324 × 10−21.86341.0920 × 10−21.59491.6243 × 10−30.90461/162.6767 × 10−31.94753.5105 × 10−31.63738.3049 × 10−40.96781/326.7308 × 10−41.99161.0781 × 10−31.70324.2689 × 10−40.96011/641.6666 × 10−42.01383.2889 × 10−41.71282.1515 × 10−40.9885Table 12
The temporal convergence rate can reach toOτ3−α; α=1.8, C′=0.6,andT=1.
τβ=1.15Rateβ=1.49Rateβ=1.999Rate1/42.2799 × 10−3—3.6088 × 10−3—6.5082 × 10−3—1/81.0716 × 10−31.08911.6547 × 10−31.12493.0294 × 10−31.10321/164.9550 × 10−41.11297.2342 × 10−41.19371.3576 × 10−31.15791/322.1798 × 10−41.18473.1736 × 10−41.18875.9499 × 10−41.19021/649.4494 × 10−51.20591.3783 × 10−41.20322.5780 × 10−41.2066


## 7. Conclusions

In this paper, we proposed two fully discrete local discontinuous Galerkin schemes for solving nonlinear time-space fractional subdiffusion/superdiffusion equations, respectively. Through carefully choosing the numerical flux on boundary terms and making a concrete theoretical analysis, we proved that our numerical schemes are unconditional stable and convergent under theL2 norm. According to the numerical results, we observe that when 0<α<1, a convergence rate Oτ2−α+hk+1 can be achieved. And when 1<α<2 and if we choose the space step size h and the time step size τ satisfies the relationship h=C′τ (C′ is a constant), the convergence order of the proposed schemes can approach to Oτ3−α+hk+1, which is consistent with the theoretical consequence proved in our theorems. The numerical results show that the fully discrete local discontinuous Galerkin (LDG) approximation is an effective and powerful method for solving the time-space fractional subdiffusion/superdiffusion equations.



---

*Source: 6954239-2020-06-22.xml*