Research Article | Open Access
A Computational Study of an Implicit Local Discontinuous Galerkin Method for Time-Fractional Diffusion Equations
We propose, analyze, and test a fully discrete local discontinuous Galerkin (LDG) finite element method for a time-fractional diffusion equation. The proposed method is based on a finite difference scheme in time and local discontinuous Galerkin methods in space. By choosing the numerical fluxes carefully, we prove that our scheme is unconditionally stable and convergent. Finally, numerical examples are performed to illustrate the effectiveness and the accuracy of the method.
Fractional calculus which is considered as the generalization of the integer order calculus attracts much attention recently its of their numerous applications in physics and engineering. They provide an excellent instrument for the description of memory and hereditary properties of various materials and processes. This is the main advantage of fractional derivatives in comparison with classic integral-order models, in which such effects are, in fact, neglected. Interest of some scholars has been shown in research on the problems involving the fractional order partial differential equations (PDEs) [1–19]. Machado et al.  introduced the recent history of fractional calculus; as for the detailed theory and applications of fractional integrals and derivatives, we can refer to [21, 22] and the references therein. Due to their numerous applications in the areas of physics and engineering, solving such equations and numerical schemes for fractional differential equations has been stimulated.
Fractional equations arise in continuous-time random walks, modeling of anomalous diffusive and subdiffusive systems, unification of diffusion and wave propagation phenomenon, and simplification of the results. There are only a few numerical works in the literature to solve fractional diffusion equations. Liu et al.  use a first-order finite difference scheme in both time and space directions for this equation, where some stability conditions are derived. In  Lin and Xu examine a practical finite difference/Legendre spectral method to solve the initial-boundary value time-fractional diffusion problem on a finite domain. In , Jiang and Ma use high-order finite element methods to solve the equation and prove an optimal convergence rate.
In this paper, we consider the following time-fractional diffusion equation: where is the order of the time-fractional derivatives. and are given smooth functions. We do not pay attention to boundary condition in this paper; hence, the solution is considered to be either periodic or compactly supported.
In the present paper, we propose a fully discrete local discontinuous Galerkin (LDG) finite element method for solving the time-fractional diffusion equation. Our fully discrete scheme is based on a finite difference scheme in time and local discontinuous Galerkin methods in space. By choosing the numerical fluxes carefully, we prove that our scheme is unconditionally stable and gives an error estimate.
What remains of this paper is organized as follows. We begin by introducing some basic notations and mathematical preliminaries which are required for establishing our results. In Section 3, we discuss the LDG scheme for the fractional equation (1), and we prove that the fully discrete scheme is unconditionally stable and convergent. Numerical experiments to illustrate the accuracy and capability of the method are given in Section 4. Finally, in Section 5, concluding remarks are provided.
2. Notations and Auxiliary Results
In this section, we introduce notations and definitions to be used later in the paper and also present some auxiliary results.
Given a spatial grid , define the mesh , for and the cell lengths , and .
We denote by and the values of at , from the right cell and from the left cell . is used to denote , that is, the jump of at cell interfaces.
We define the piecewise-polynomial space as the space of polynomials of the degree up to in each cell ; that is,
For error estimates, we will be using two projections in one dimension , denoted by ; that is, for each , and special projection ; that is, for each ,
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. Fully Discrete LDG Scheme
Let be the approximation of ; respectively, . We define a fully discrete local discontinuous Galerkin scheme as follows: find , such that, for all test functions ,
The “hat” terms in (11) in the cell boundary terms from integration by parts are the so-called “numerical fluxes,” which are single valued functions defined on the edges and should be designed based on different guiding principles for different PDEs to ensure stability. It turns out that we can take the simple choices such that
In order to simplify the notations and without loss of generality, we consider the case in its numerical analysis. Now, we consider the stability for the scheme (11), we have the following result.
Theorem 1. For periodic or compactly supported boundary conditions, the fully-discrete LDG scheme (11) is unconditionally stable, and the numerical solution satisfies
Proof. We will prove Theorem 1 by mathematical induction. When , scheme (11) is
Taking the test functions , we obtain Here,
If we take fluxes (12) and after some manual calculation, we can easily obtain .
From the fact that we can get
Now, suppose the following inequality holds We need to prove . Let and take the test functions in scheme (11); we can obtain
Taking fluxes (12), we can easily obtain . Then, the last inequality gives
This finishes the proof of the stability result.
Theorem 2. Let be the exact solution of the problem (1), which is sufficiently smooth with bounded derivatives. Let be the numerical solution of the fully discrete LDG scheme (11); then, there hold the following error estimates when : and when :
Proof. We denote
Subtracting (11) from (10) and with fluxes (12), we can obtain the error equation
Using (24), the error equation (25) can be written as
Taking the test functions in (26), using the properties (4) and (6), then the following equality holds:
Based on the fact that , we can obtain
For the sake of convenience, we denote
(1) We start with the following estimate:
When , (35) becomes
Denoting , then we can obtain
Next, we suppose the following inequality holds:
Let in the inequality (35); we deduce
Notice the fact that we know that is,
Inequality (31) follows.
By some calculations and analyses, we know that increasingly tends to . For more details of the proof, we refer to . So we can obtain
(2) The above estimate has no meaning when due to . So we must reconsider it for the case .
We suppose the following estimate holds:
By the similar techniques used in (1) and that in , we can obtain (40) easily. Here, we omitted the proof to save space. Then, we know that when ,
Thus, Theorem 2 follows by the triangle inequality and the interpolating property (7).
4. Numerical Examples
In this section, we offer some numerical examples to illustrate the accuracy and capability of the method. For this purpose, we calculate the numerical results of the exact solutions (for the cases where exact solutions are available). We mainly focus on the spatial accuracy, so a small time step is used such that errors stemming from the temporal approximation are negligible. With the aid of successive mesh refinements, we have verified that the results shown are numerically convergent.
Example 1. We consider time-fractional equation (1) in ; the corresponding forcing term is of the form Then, the exact solution is . The space and time step is , respectively. We check the spatial accuracy by fixing the time step sufficiently small to avoid contamination of the temporal error. From Tables 1, 2, 3, 4, 5, 6, 7, 8, and 9, we can see that the errors in -norm and -norm attain optimal order of accuracy for piecewise polynomials for . In Table 10, we show the errors in -norm and -norm attains order of accuracy for two values of and .