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 () or second-order central difference method for (), 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 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.
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 in theoretical; for the superdiffusion linear problem, they got the approximation of order . 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 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 for the space fractional diffusion problem and suboptimal order of convergence 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: and the time-space fractional superdiffusion equation as follows:where , , denotes the order number of space fractional derivatives, describes the fractional subdiffusion equation (1), describes the fractional superdiffusion equation (2), is a real positive parameter, and 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] and the time fractional Caputo derivative in (2) defined as [24]where is the Gamma function.
The space fractional terms in (1) and (2) are defined through the fractional Laplacian operator , , which is also known as the Riesz fractional derivative under the assumption that vanishes at the end point of an infinite domain [32].
The spatial fractional differential operator is also an important tool to describe the anomalous diffusion phenomenon; when , it represents a Lévy -stable flight [33]; when , 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 , 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 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 function given in the domain ( can be finite or infinite), the expressionsare called left- and right-side Riemann-Liouville fractional integrals of order (, being an integer), respectively. is a Gamma function.
Definition 2. The left and right of Riemann-Liouville fractional derivatives of order 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 in the domain ( can be finite or infinite) are defined as follows:where , , denote the left and the right Riemann-Liouville fractional derivatives of order , respectively.
Definition 4. For a function defined on the infinite domain , 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 for , , , 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].
Lemma 7 (linearity). One has
Lemma 8 (semigroup property).
By equalities (10), for any continuous function with and , we can obtain the following result [29]: In the above equation, if , the fractional Laplacian operator becomes the fractional integral operator. Similar to [29], since , then , and we define
Lemma 9. When , we can rewrite the space fractional Laplacian operator in a composite form of first-order derivative and a fractional integral of order as follows [29]:
Definition 10. We define the left- and right-fractional space norms as [18]where is the seminorm of fractional space and , refer to the closure of with respect to and , respectively.
Lemma 11 (adjoint property). One has
Lemma 12 (see [18]). One has
Lemma 13 (see [19]). Let . Then, for any , then
From the above definitions and Lemma 8, together with the properties of Riemann-Liouville fractional integrals, the following lemmas hold.
Lemma 14. For any , the fractional integral satisfies the following property:
Proof. Since , by Lemmas 9, 11, 12, and 13, we obtain
Theorem 15 (see [26]). If , then (or or ) is embedded into (or or ), and is embedded into both of them.
Lemma 16 (fractional Poincaré-Friedrichs inequality [18]). One has
Lemma 17. Suppose the fractional integral is defined on and let . Then
Proof. One has
Lemma 18 (see [29]). Suppose is a smooth function defined on . is discretization of the domain with interval width , and is an approximation of u in . For all , is a polynomial of degree up to order , and . is the degree of the polynomial. Then for , we obtain where is a constant independent of .
Assume the following mesh to cover the computational domain , consisting of cells for , where , denoting the cell lengths , , . Define 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 , that is, for each ,and special projection , that is, for each ,The two projections satisfy the following approximation inequality [36]:where , or . 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 be denoted by and the associated norm by . If , we drop .
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 be the time mesh size, is a positive integer, let be mesh point, and let , be mid-mesh points. Here and below an unmarked norm refers to the usual 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 of subdiffusion equation.
For the case of subdiffusion equation (1), we adopt the backward Euler method to numerically approximate the time fractional derivative at as follows [21, 27]:Here satisfy the following properties: is the truncation error and satisfies the following inequality [21]:where is a constant only depending on . We rewrite (1) as a first-order systemAnd with (30), we give the weak form of system (33) at as follows:where
Let be the approximation of , respectively. We define a fully discrete local discontinuous Galerkin (LDG) scheme as follows: Find , such that for all 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 thatWe 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).
3.2. Superdiffusion Case
Similar to the subdiffusion case, we adopt the center difference scheme to estimate the time fractional derivative of superdiffusion case (2) at as follows [24]:where , satisfy the following properties: is the truncation error; that is, which satisfies the following estimate:where is the constant which only depends on .
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 :where
Let be the approximation of , respectively. We define a fully discrete local discontinuous Galerkin (LDG) scheme as follows: find , such that for all
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 , 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 as follows:
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 satisfies
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.
Proof. First, consider ; the LDG scheme (35) isTaking the test functions , we getConsidering the integration by parts formulation: , we obtain the interface condition, Thus, we obtainFrom the fact that we have Combining the boundary condition and Lemma 14, we getNow we assume that the following inequality holds:We need to prove . Let and take the test functions in scheme (35); we obtain Calling Lemma 14 again and considering the boundary condition, we have
Theorem 20. For a sufficiently small step size , the fully discrete LDG scheme (42) is unconditional stable and the numerical solution satisfies
Proof. Similar to the subdiffusion case, using the inequality 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. 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:
When ,When ,where is a constant only depending on .
Proof. DenoteSubtracting 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 asTaking the test functions , , in (61), using the properties (27)-(28), then the following equality holds:That is,We prove Theorem 21 by mathematical induction.
First, we consider the case when , (63) becomesNotice the fact that
Together with property (29) and Lemma 18, we can get Combining this with Definition 4 and (59), using Young’s inequality we obtainUsing Young’s inequality again in (66), we getIf we choose and recall the fractional Poincaré-Friedrichs Lemma 16, we havewhere is a constant depending on .
Next, we suppose that the following inequality holds:When , from (63) and Young’s inequality, we can obtain SincethenUsing Definition (7) into above inequalities and combining Young’s inequalities, we haveAnalogously, if we choose small enough and recalling the fractional Poincaré-Friedrichs Lemma 16 again, we obtainwhere is a constant related to .
According to the principle of mathematical induction, we getBy [21, 27], we know that is monotonly increasing and tends to ; thus, where is a constant only depending on .
When , , it is meaningless for (76), so as for the case of , we should give out the estimate again.
We prove the following estimate by using the mathematical induction:Obviously, when , (77) is workable clearly.
When , by (63), we obtain Using Young’s inequality into the last inequality above, and choosing , then combining the property of and Lemma 16 (Poincaré-Friedrichs inequality), we obtain By mathematical induction, we get inequality (77).
Since , , that is, , hence, when , inequality (77) becomeswhere is a constant only depending on .
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 ,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 haveWhen , the denominator of equality (88) is strictly greater than zero, so we only need to consider the condition of the molecular, because When and , always , so we only need to estimate the positive or negative property of the molecular of the last equality above. For , always holds.
For any given , we define the function as follows:Differentiating the function with respect to , we obtainLet ; there is only one zero root of the derived function on interval : Since , for any given , the function with respect to can only increase first and then decrease or decrease first and then increase on interval .
Now if we can prove that the value of any point on interval is positive, then the function must be increasing first and then decreasing. For example, we choose the parameter and estimate the positive or negative property of the following equality: . Since , we get ; hence, the function with respect to can only increase first and then decrease on interval . That is, for any given , when , we have . When , .
From the above discussion we obtain ; that is, the function is a monotone increasing function on interval .
Considering the function limit , by Hospital’s rule, we have Hence, Therefore, the sequence monotone increasing tends to . We obtainwhere is a constant only depending on .
When , , estimate (96) is meaningless, so for the case of , we need to estimate again.
Analogous to the subdiffusion case, we can prove the following estimate still using the mathematical induction:Because , , that is, , as a result, when , the inequality (97) becomeswhere is a constant only depending on .
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 . 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:where , , for the subdiffusion case (1), or for the superdiffusion case (2), and , .
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 in space . By the definition of the space , we know that, for all , holds, where
By the LDG scheme (35), we obtain the fully discrete scheme of (99) as follows:where
In terms of the basis , we denote approximate solution functions by , , and and let test functions , , and . Putting these expressions into the LDG scheme (100) and taking the flux in (36), we get the linear system of the scheme in th element at as follows:where is the mass matrix, denote the polynomials order, and is the quadrature of the fractional integral operator and test function in th element. From [29], we know that whereand is a vector valued function, . Then, by solving the linear system (101), we can obtain the solution .
6.2. Numerical Results
Example 1. Consider the fractional time-space subdiffusion equation (99) in domain which satisfies the following initial boundary conditions: In order to obtain the exact solution , 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 ; Tables 1 and 2 list out the 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 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 , fixing (or ), and changing (or ), the temporal accuracy of the time-space fractional subdiffusion equation under the norm can converge to for , and close to 1 when , the results agree with the theoretical computation results in Theorem 21.
Example 2. Consider the time-space fractional subdiffusion equation (99); for , we choose the exact solution as and .
Similar to Example 1, we give out the 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 , 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 . The data in Table 7 display that, for , the convergence order in time direction can be up to and close to 1 order for , which are consistent with Theorem 21.
Example 3. In this example, we consider the case for ; that is, (99) is a time-space fractional superdiffusion model. We assume that the equation satisfies the following initial boundary conditions in the domain :Here we choose , and the source term is Then the time-space fractional superdiffusion obtains the same exact solution as Example 1: .
Analogous to the case of fractional subdiffusion equation, we consider the errors under the 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 , 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 , 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 , fix (or ), and then change (or ); the temporal accuracy of the time-space fractional superdiffusion can reach for and be close to 1 when . The numerical results are consistent with the theoretical results proved in Theorem 22.
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 the norm. From the numerical experiment results, we can observe that a convergence rate is achieved for , and for , when the space step size and the time step size satisfy the relationship ( is a constant), the convergence rate of our scheme can approach . 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.
Competing Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The project is supported by NSF of China (11371289, 11426068, and 11501441) and the talent project of Huizhou University of China (2015JB018).