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

. 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 𝐿 2 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.


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][2][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][15][16][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 2 Advances in Mathematical Physics 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 Δ 2− + Δ − ℎ 2 in theoretical; for the superdiffusion linear problem, they got the approximation of order Δ 2 + Δ 1− ℎ 2 .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 convectiondiffusion 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 (ℎ +1 ) for the space fractional diffusion problem and suboptimal order of convergence (ℎ +1/2 ) for the general spatial fractional convection-diffusion problem are established, respectively.
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] and the time fractional Caputo derivative in (2) defined as [24]    (, ) where Γ(⋅) is the Gamma function.
The spatial fractional differential operator is also an important tool to describe the anomalous diffusion phenomenon; when 1 <  < 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  2 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.

Definitions and Auxiliary Results
In this section, we present a few definitions and recall some properties of fractional calculus [34,35].
Definition 2. The left and right of Riemann-Liouville fractional derivatives of order  ( − 1 <  < ) in the domain Ω = (, ) (,  can be finite or infinite) can be derived by Riemann-Liouville fractional integrals as follows: where Γ(⋅) is a Gamma function.When  is an integer,   =   =   .Definition 3 (see [34]).The fractional Riesz derivatives of order  ( − 1 <  ≤ ) in the domain Ω = (, ) (,  can be finite or infinite) are defined as follows: where   = 1/2 cos(/2),  ̸ = 1,     ,     denote the left and the right Riemann-Liouville fractional derivatives of order , respectively.Definition 4. For a function (, ) defined on the infinite domain (−∞ <  < ∞),  ∈ (0, ], the following equality holds [34]: 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 function (, ) is subject to homogeneous Dirichlet boundary conditions.In this paper, the developments and analysis are based on this definition.Lemma 6. Suppose that (  /)(, ) = 0 for  → ∓∞, ∀0 ≤  ≤ , ( − 1 <  < ), and then the relationship between Riemann-Liouville fractional integrals and Riemann-Liouville fractional derivatives satisfies [30] 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].
Assume the following mesh to cover the computational domain Ω = [, ], consisting of cells ℎ as the space of polynomials of the degree up to  in each cell   ; that is, To prepare for the analysis of the error estimates, we need two projections in one dimension [, ], denoted by Q, that is, for each , and special projection P ± , that is, for each , The two projections satisfy the following approximation inequality [36]: where   = P() − (), or   = P ± () − ().The positive constant , solely depending on , is independent of ℎ.  ℎ denotes the set of boundary points of all elements   .
In the present paper, we use  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  2 () be denoted by (⋅, ⋅)  and the associated norm by ‖ ⋅ ‖  .If  = Ω, we drop .

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.

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   (, )/  of subdiffusion equation.

Advances in Mathematical Physics
Let   ℎ ,   ℎ ,   ℎ ∈   ℎ be the approximation of (⋅,   ), (⋅,   ), (⋅,   ), respectively.We define a fully discrete local discontinuous Galerkin (LDG) scheme as follows: Find where 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 We remark that the choice for the fluxes (36) is not unique.In fact the crucial part is taking   ℎ and   ℎ from opposite sides.We know the truncation error is   () from (30).
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   : Advances in Mathematical Physics where Let   ℎ ,   ℎ ,   ℎ ∈   ℎ be the approximation of (⋅,   ), (⋅,   ), (⋅,   ), respectively.We define a fully discrete local discontinuous Galerkin (LDG) scheme as follows: find Remarks 1. Originally, we assume that  has compact support and restrict the problem to the bounded domain Ω.As a consequence, we impose homogeneous Dirichlet boundary conditions for  ∉ Ω.First, we define the fluxes at the boundary as ] for the right boundary and û1/2 = (, ), q1/2 =  − 1/2 + (/ℎ)[ 1/2 ] for the left boundary [29], where  is a positive constant and [⋅] is the jump term operator of function .
Next, we define the boundary term operator L as follows: (43)

Stability
In this subsection, we present the stability analysis for the fully discrete LDG schemes (35) 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 Now, we use the mathematical induction to proof Theorem 19.

Advances in Mathematical Physics
Proof.First, consider  = 1; the LDG scheme (35) is Taking the test functions Considering the integration by parts formulation: , we obtain the interface condition, Thus, we obtain From the fact that we have (51) Combining the boundary condition and Lemma 14, we get Now we assume that the following inequality holds: We need to prove ‖ +1 ℎ ‖ ≤ ‖ 0 ℎ ‖.Let  =  + 1 and take the test functions V =  +1 ℎ ,  = − 0  +1 ℎ ,  =  0  +1 ℎ in scheme (35); we obtain Proof.Similar to the subdiffusion case, using the inequality ‖ −1 ‖ ≤ (‖()‖+‖()‖) 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.

Error Estimates
In this section, we discuss the error estimates for the fully discrete LDG schemes (35), (42) respectively.
Theorem 21.Let (,   ) be the exact solution of the problem (1), which is sufficiently smooth with bounded derivatives, and   ℎ be the numerical solution of the fully discrete LDG scheme (35); then the following error estimates hold: (57) where   is a constant only depending on . ( Subtracting equation (35) from equation (34), with the fluxes on boundaries, we obtain the error equation: Using equations ( 59) and (43), after a simple rearrangement of terms, the error equation (60) can be rewritten as Advances in Mathematical Physics +  0 √ (P +  (,   ) −  (,   )) Taking the test functions V = P −    ,  = − 0 Q   ,  =  0 P +    in (61), using the properties ( 27)-( 28), then the following equality holds: We prove Theorem 21 by mathematical induction.First, we consider the case when  = 1, (63) becomes Notice the fact that ‖P −  0  ‖ ≤ ℎ +1 ,  ≤  (67) If we choose  ≤ min{1, 1/} and recall the fractional Poincaré-Friedrichs Lemma 16, we have where  is a constant depending on , , .
Next, we suppose that the following inequality holds: When  =  + 1, from (63) and Young's inequality, we can obtain Since (73) Analogously, if we choose  small enough and recalling the fractional Poincaré-Friedrichs Lemma 16 again, we obtain where  is a constant related to , , .
According to the principle of mathematical induction, we get By [21,27], we know that  −  −1 −1 is monotonly increasing and tends to 1/(1 − ),  → ∞; thus, where   is a constant only depending on .

Advances in Mathematical Physics 13
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: Using Young's inequality into the last inequality above, and choosing  ≤ min{1, 1/}, then combining the property of   and Lemma 16 (Poincaré-Friedrichs inequality), we obtain By mathematical induction, we get inequality (77).
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.Let (,   ) be the exact solution of problem (2), which is sufficiently smooth with bounded derivatives, and   ℎ be the numerical solution of the fully discrete LDG scheme (42); then the following error estimates hold. When where   is a constant only depending on .
Proof.Subtracting equation (42) from equation (41), and combining the fluxes on boundaries, we obtain the error equation as follows: Similar to the error estimate of subdiffusion equation, first we can prove the following error inequality: where  is a constant related to , , .
In fact, 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 Theorem 21 to obtain Now we consider the function, Differentiating the function () with respect to , we have When  > 1, the denominator of equality (88) is strictly greater than zero, so we only need to consider the condition of the molecular, because When 1 <  < 2 and  > 1, always  2 ( − 1)  > 0, so we only need to estimate the positive or negative property of the molecular of the last equality above.
For any given , we define the function as follows: Differentiating the function () with respect to , we obtain Let   () = 0; there is only one zero root of the derived function   () on interval (1, 2): Since (1) = (2) = 0, for any given , the function () 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 () 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: we get (3/2) > 0; hence, the function () with respect to  can only increase first and then decrease on interval (1, 2).That is, for any given , when 1 <  < 2, we have () > 0.
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.

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 ℎ.By choosing the appropriate time step size, we avoid the contamination of the temporal error compared with the spatial error.
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 {  ()}  =1 in space  ℎ  .By the definition of the space  ℎ  , we know that, for all V ∈  ℎ  , V() = ∑  =1 V    () holds, where V  = V(  ).

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: In order to obtain the exact solution (, ) = (2 − 1) 2  2 (1 − ) 2 , we choose the source term as 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 ℎ = 1/2  , ( = 2, . . ., 6); Tables 1 and 2 list out the  2 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(ℎ +1 ) in spatial direction; this matches the theoretical results of Theorem 21.
The numerical example results listed in Tables 3 and 4 show that when we choose the appropriate space step size ℎ and the time step size sequence  = 1/2  , ( = 2, . . ., 6), fixing  (or ), and changing  (or ), the temporal accuracy 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 ℎ = 1/2  , ( = 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 ℎ and the time step size sequence as  = 1/2  , ( = 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.Example 3. In this example, we consider the case for 1 <  < 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. Analogous to the case of fractional subdiffusion equation, we consider the errors under the  2 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 ℎ = 1/2  , ( = 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(ℎ +1 ), which agrees with the theoretical results proved in Theorem 22.The numerical results in Tables 11 and 12 show that when ℎ =    (  is a positive constant), we choose the time step size sequence  = 1/2  , ( = 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.

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 the  2 norm.From the numerical experiment results, we can observe that a convergence rate O( 2− + ℎ +1 ) is achieved for 0 <  < 1, and for 1 <  < 2, when the space step size ℎ and the time step size  satisfy the relationship ℎ =    (  is a constant), the convergence rate of our scheme can approach O( 3− + ℎ +1 ).The numerical results show that the fully discrete local discontinuous Galerkin (LDG) approximation is effective and powerful for solving timespace fractional subdiffusion/superdiffusion equations.