Optimal a priori error estimates for the -version of the local discontinuous Galerkin method for convection–diffusion problems
Abstract
We study the convergence properties of the -version of the local discontinuous Galerkin finite element method for convection-diffusion problems; we consider a model problem in a one-dimensional space domain. We allow arbitrary meshes and polynomial degree distributions and obtain upper bounds for the energy norm of the error which are explicit in the mesh-width , in the polynomial degree , and in the regularity of the exact solution. We identify a special numerical flux for which the estimates are optimal in both and . The theoretical results are confirmed in a series of numerical examples.
1. Introduction
This paper contains the first a priori error estimate of the -version of the so-called local discontinuous Galerkin (LDG) finite element method for convection-diffusion problems. Such an error analysis, which takes into account both the mesh-size of the element, , and the degree of the approximating polynomial in it, , is quite relevant for the LDG method since, being a locally conservative method that does not require any inter-element continuity, it is ideally suited for -adaptivity. In this paper, we consider a model convection-diffusion equation in one space dimension with Dirichlet boundary conditions and obtain, for a special choice of the numerical fluxes defining the LDG method, a priori error estimates that are optimal both in and , even for ; all other error estimates available in the current literature are suboptimal in both and and do not give a rate of convergence for .
The LDG method was introduced by Cockburn and Shu in Reference 11 as an extension to general convection-diffusion problems of the numerical scheme for the compressible Navier-Stokes proposed by Bassi and Rebay in Reference 1. This scheme was in turn an extension of the Runge-Kutta discontinuous Galerkin (RKDG) method developed by Cockburn and Shu Reference 10Reference 9Reference 8Reference 6Reference 12 for nonlinear hyperbolic systems. For a fairly complete set of references on RKDG and LDG methods, see the short monograph by Cockburn Reference 4; see also the review of the development of discontinuous Galerkin methods by Cockburn, Karniadakis, and Shu Reference 7.
To put our result in proper perspective, let us briefly describe the relevant results available in the current literature. There are only a few a priori error estimates for the LDG method and they are all for the -version of the method. The first a priori error estimate for the LDG method was obtained in 1998 by Cockburn and Shu Reference 11 who proved that, when polynomials of degree are used, the LDG method converges in the energy norm at a rate of order . This rate of convergence was obtained for the general form of the so-called numerical fluxes that appear in the definition of the LDG method and is sharp since for the numerical flux proposed by Bassi and Rebay Reference 1 this rate is actually achieved. Later, this analysis was extended by Cockburn and Dawson Reference 5 to the case in which the convective velocity and the diffusion tensor depend on and the domain is bounded; the rate of convergence of order was once again obtained.
Although the rate of convergence of order is sharp for general fluxes, Cockburn and Shu Reference 11 reported numerical experiments in the one-dimensional case indicating that, for a special numerical flux, a rate of convergence of order is achieved for very smooth solutions. This indication was later put on firm mathematical grounds by Castillo Reference 3 who showed, for the model problem of constant-coefficient, linear convection-diffusion in one space dimension, that the LDG method with a particular numerical flux converges with the optimal rate of convergence of order . Castillo’s result can be viewed as an extension to the convection-diffusion setting of the a priori error estimate for the discontinuous Galerkin (DG) method for the purely convective case obtained in 1974 by LeSaint and Raviart Reference 15 who prove that the rate of convergence is of order .
In this paper, we obtain an a priori error estimate for the -version of the LDG method for general numerical fluxes which is explicit in the mesh-width and the polynomial degree . Assuming that the -th derivative of the exact solution in the energy norm plus the -norm of the time derivative is finite, we show that, for general numerical fluxes, the energy norm of the error has a rate of convergence of order in the purely convective case and of in the convection-diffusion case. Moreover, by using the special numerical flux studied by Castillo Reference 3, we obtain the optimal rate of convergence of order for totally arbitrary meshes and polynomials of degree in all elements. This result holds in the purely convective case as well as in the purely parabolic case.
Let us give an idea of how the error estimate is obtained. First, using the technique employed by Cockburn and Shu Reference 11, we find an upper bound for the energy norm of a projection into the finite element space of the error. Then, following Castillo Reference 3, we eliminate as many as possible terms in the upper bound of the error by carefully defining the numerical flux of the LDG method and by suitably choosing such a projection. Indeed, instead of using the L-projection operator used by Cockburn and Shu Reference 11, the projections used by Houston, Schwab and Süli Reference 14Reference 30Reference 29, or the Lagrange interpolation of Gauss-Radau points used by Castillo Reference 3, we pick the more advantageous projection used in 1985 by Thomeé Reference 31 in his study of discontinuous Galerkin time-discretizations for parabolic problems and recently by Schötzau and Schwab Reference 23Reference 22 in their study of the -version of this method. Indeed, with this projection, many terms in the upper bound of the error become identically zero which allows us to obtain an optimal rate of convergence after a simple application of the sharp -approximation results for this projection.
Recent work on other discontinuous finite element methods for convection-diffusion (and for pure diffusion) problems has been reviewed by Cockburn, Karniadakis and Shu Reference 7. See, in particular, the numerical method of Baumann and Oden Reference 2, the optimal error estimates for the method as applied to nonlinear convection-diffusion equations by Rivière and Wheeler Reference 20, the analysis of several versions of the Baumann and Oden method for elliptic problems by Rivière, Wheeler and Girault Reference 21, and the -version analyses by Houston, Süli and Schwab Reference 14Reference 30Reference 29 of discontinuous Galerkin methods for second-order problems with nonnegative characteristic form. We also mention the recent work of Wihler and Schwab Reference 32 in which robust exponential rates of convergence of DG methods for (stationary) convection-diffusion problems in one space dimension are proven.
The organization of the paper is as follows. In section 2, we describe the LDG method. In section 3, we state and prove the a priori error estimate for the constant-coefficient convection-diffusion problem and the special numerical flux for which the estimates are optimal in both and . In section 4, we discuss several extensions and, in section 5, we perform numerical experiments that verify the theoretical results. We end our presentation with some concluding remarks in section 6.
2. The LDG method
In this section, we introduce and briefly discuss the various key elements of the LDG method for a simple model problem.
2.1. The model problem and its weak formulation
In this paper we consider the following model convection-diffusion equation in one space dimension:
with the initial condition
and the Dirichlet boundary conditions
The unknown function is a scalar, and we assume the velocity to be a positive and the diffusion coefficient to be a nonnegative number; we choose to work with a positive velocity simply to fix the location of the possible boundary layer at . Note that in the purely convective case , only the Dirichlet boundary condition at is taken into account.
The weak formulation we are going to use is obtained as follows. First, we introduce the new variable and the “flux” function
and rewrite Equation 2.1–Equation 2.3 in the form
Next, given the nodes , we define the mesh and set ; furthermore, we define . To the mesh , we associate the so-called broken Sobolev space
For a function the one-sided limits at the nodes are denoted as follows:
Throughout, we assume the exact solution of (Equation 2.1)–(Equation 2.3) belongs to . Then, it satisfies the following equations:
for all test functions and for . Here, the time derivative is to be understood in the weak sense and .
2.2. The method
A discrete version of the above mixed formulation is obtained by restricting the trial and test functions to finite-dimensional subspaces and by replacing the flux function at the nodes by a numerical flux : Find such that for all and for the following equations hold:
Upon a choice of basis for the subspaces , and, more importantly, of the numerical fluxes, the semidiscrete problem Equation 2.5 becomes an ODE system of dimension on , where . In what follows we do not consider the impact of the time discretization and refer to Shu Reference 28, Shu and Osher Reference 26Reference 27 and Gottlieb and Shu Reference 13 for the analysis of certain TVD Runge-Kutta methods for the solution of ODE systems.
Our choice of the space is the space of discontinuous, piecewise polynomial functions
where denotes the set of all polynomials of degree less or equal than on . Notice that the polynomial degrees can vary from element to element.
To complete the definition of the LDG method, it remains to define the numerical flux .
2.3. The numerical flux
Crucial for the stability as well as for the accuracy of the LDG method is the choice of the numerical flux . To define it, we introduce with the notation in (Equation 2.4) the following quantities:
The numerical flux has the following general form:
which also holds at the boundary if we define
Let us stress several important points concerning this numerical flux:
- •
In the purely hyperbolic case, i.e., in the case , if we take , we obtain the well-known “upwinding” flux of the original DG method; see, e.g., Reference 4.
- •
Note that . This is so because we want to be able to solve for in terms of element by element. This local solvability, which gives the name to the LDG method, is not shared by most mixed methods and allows us to easily eliminate the unknown from the equations.
- •
The main purpose of the coefficient is to enhance the stability of the method; that is why it must be a positive number. This results in an improvement of the accuracy of the method too.
- •
The choice ensures the stability of the LDG method.
- •
The main purpose of the coefficients is to enhance the accuracy of the method. Thus, if we take, following Bassi and Rebay Reference 1, , the rate of convergence of the energy norm is of order for smooth functions. If instead, following Castillo Reference 3, we take , we obtain the optimal rate of order .
Let us point out that, if we consider the general form of the numerical fluxes, it is possible to obtain exponential convergence for piecewise analytic exact solutions, but not optimality in both and . As we shall see, this optimality is guaranteed for completely arbitrary meshes if we take (an extension of) Castillo’s Reference 3 choice of the numerical flux , namely,
where .
This flux is obtained by setting and
Note again that in the purely convective case, , this numerical flux is nothing but the standard upwinding flux used by the original DG method. Note also that the fact that the coefficient has a special form is a reflection of the fact that at there might be a boundary layer which requires special treatment. The factor ensures the optimality in both and of the energy norm of the error.
The definition of the LDG method is now complete.
3. Error analysis
This section is devoted to our main a priori error estimates. First, we state and briefly discuss the results; the remainder of the section is devoted to their proof.
3.1. A priori estimates
Our main result follows naturally from an estimate of a suitably defined projection of the error and from the -approximation properties of the projection . Each of these results is contained in several lemmas that we state next. To do that, we need to introduce the projection and the norm in which we measure .
The projection is the operator from to that associates to where, for each interval , is defined by the following conditions:
Let us now introduce a norm that appears naturally in the analysis of the LDG method. In what follows, we denote by the L-norm in the subdomain and omit when . For the norm is defined as follows:
where the energy norm is given by
and
Note that information about the numerical flux is contained in the norm only through . We are now ready to state our results.
The error between the exact solution and the approximation given by the LDG method with numerical flux (Equation 2.6) satisfies the inequality
where
and
If we combine the above result with the estimates of , for and , we immediately obtain our desired error bound. To state the corresponding approximation results, we introduce on the reference interval and for the following weighted seminorm:
We can now state our approximation result.
Let be the reference interval and a generic polynomial degree on . Assume for and . Then we have the following estimates:
where depends on but is independent of and .
For standard finite element methods, the introduction of the weighted norms enables one to show that for singular solutions the -version of the method, i.e., when the mesh is fixed and increases unboundedly, yields twice the convergence rate than the -version provided that the singularity lies at a mesh point ; see, e.g., Schwab Reference 24. The results in Lemma 3.2 are slightly suboptimal with regard to these aspects as will be shown in the numerical experiments in section 5 below. However, for smooth solutions we obtain optimal approximation properties for in and . This can be inferred immediately from Lemma 3.2 and standard scaling and interpolation arguments.
Let for and . Then we have the following estimates:
where depends on but is independent of , and .
Now, assume that , where
Our main result is a simple consequence of the above lemmas. Indeed, since
from Lemmas 3.1 and 3.3, we obtain the following result.
Let be the error between the exact solution and the approximation given by the LDG method with numerical flux (Equation 2.6) and polynomial degree on each interval. Then, for totally arbitrary meshes, the energy norm of the error satisfies the inequality
The error estimates in Theorem 3.4 are optimal in and for smooth solutions, even for the case in which piecewise constant approximations () are used. Note also that in the purely convective case, , the above error estimate is nothing but the extension of the super-convergence error estimate of LeSaint and Raviart Reference 15 for the -version of the DG method for purely convective problems.
The proof of Lemma 3.3 actually gives us estimates which are completely explicit in the mesh-width , in the polynomial degree , and in the regularity of the exact solution (see Proposition 3.12 below). Hence, completely explicit error estimates in the energy norm can be obtained. In conjunction with geometric meshes and linearly increasing polynomial degrees, such estimates can be used in the -version to prove exponential rates of convergence in the presence of solution singularities; see, e.g., the recent monograph by Schwab Reference 24 and the references therein. However, since the corresponding analytic regularity in space-time still remains to be found, we do not further pursue these issues here.
From Theorem 3.4, we conclude that in the -version of the LDG method, where the mesh is kept fixed and the polynomial degree is increased, we have . Hence, for smooth solutions convergence rates of arbitrarily high algebraic order in are possible. This is sometimes referred to as spectral convergence. Furthermore, for solutions which are analytic in , even exponential rates of convergence are obtained in the -version, i.e.,
with constants independent of . This result can immediately be derived from Lemma 3.1, properties of (see, e.g., (Equation 3.15) and (Equation 3.16) below) and standard approximation theory for analytic functions. We note that (Equation 3.4) holds true for general fluxes as well.
For small diffusivities , i.e., for in (Equation 2.1), the solutions typically exhibit viscous boundary layers (or shock profiles) of length scale or . In principle, layer components in the solutions are still analytic (see the work of Melenk Reference 16 and Melenk and Schwab Reference 19Reference 18 for a complete characterization of boundary layers in stationary problems with analytic input data) and can thus be approximated at exponential rates of convergence, in agreement with (Equation 3.4). However, the estimate (Equation 3.4) is not robust with respect to the diffusivity and deteriorates as . A remedy is to employ a needle-element of the appropriate width or geometric mesh refinement near the boundary. It has been shown recently by Melenk Reference 16, Melenk and Schwab Reference 17Reference 19, Schwab and Suri Reference 25 and Wihler and Schwab Reference 32 that the use of these mesh-design principles yields exponential rates of convergence that are robust with respect to the diffusivity parameter . We demonstrate this robustness in our numerical examples in section 5.5 below.
3.2. Proof of the basic estimate
This section is devoted to the proof of Lemma 3.1. To do so, we follow the technique used by Cockburn and Shu Reference 11; see also Cockburn Reference 4, Castillo Reference 3 and Cockburn and Dawson Reference 5.
We start by rewriting the definition of the LDG method in compact form for general numerical fluxes. Integrating Equation 2.5 with respect to from to and summing over all elements, it turns out that the LDG solution is defined as follows:
Find such that
where the discrete bilinear form is given by
and the discrete linear form is given by
The basic error estimate now follows by standard manipulations. Indeed, since we have that
we obtain that
where , and this implies that
It only remains to obtain a suitable expression for and an upper bound for .
For any , it holds that
where
and is defined by (Equation 3.3).
This is a direct consequence of the definition of the form .
■For any and the numerical flux (Equation 2.6), we have
Taking into account the definition of the form , (Equation 3.6), and the definition of the projection , (Equation 3.1) and (Equation 3.2), we easily get that
The result follows from simple applications of Cauchy-Schwarz’s and Young’s inequalities and from the definition of the functional defined in Lemma 3.9.
■Now, inserting the results of Lemmas 3.9 and 3.10 into (Equation 3.8), we get the inequality
which is of the form
with
Since inequality (Equation 3.9) holds true for all , Lemma 3.1 now follows after a simple application of the following result.
Suppose that for all we have
where , , and are nonnegative functions. Then, for any ,
Define and fix . Setting , the hypothesis implies that for
Integrating over yields
Hence,
Since , the assertion in Lemma 3.11 follows.
■This completes the proof of Lemma 3.1.
3.3. Proof of the -approximation results
This section is devoted to the proof of Lemma 3.2 which follows from the subsequent finer approximation results after an application of Stirling’s formula.
Let and recall that
We have:
Let for . Let be defined by
Then we have
and
for any .
The first estimate was obtained by Schötzau Reference 22 and Schötzau and Schwab Reference 23. Nevertheless, we present a detailed proof for the sake of completeness. We consider only since the proof for is similar. We proceed in several steps.
Step 1: First, we derive bounds on the difference in terms of the Legendre coefficients of . To do so, denote by , , the Legendre polynomial of degree on and expand the function into the series with . Since , it can be seen from (Equation 3.10) that is uniquely given by the series
Hence, the difference can be written as
Let be the -projection from onto . Since , and , we obtain from the definition of the following bounds:
Step 2: To estimate the sums in the above equalities, we start by expanding into the series . Integrating this expression yields
and employing for the identity
and rearranging terms, we obtain
Comparing coefficients in the Legendre expansions, we conclude that
Hence, after some simple algebraic manipulations,
and
As a consequence,
and since , we get
Step 3: Now, note that after inserting the estimates (Equation 3.13) and (Equation 3.14) into (Equation 3.11) and (Equation 3.12), respectively, we get
Replacing in these inequalities by , where is an arbitrary polynomial of degree , and taking into account that and give
Schwab Reference 24 proved that
for any , and the existence of a polynomial such that
for any . We now simply insert these estimates into (Equation 3.15) and (Equation 3.16) to conclude that
for any .
The corresponding estimates for are obtained by symmetry. Since
this proves Proposition 3.12.
■From Proposition 3.12 we obtain by standard scaling and interpolation arguments the following -approximation properties of :
For each interval , , we have for , real, the estimates
and
for any . The constant is independent of , and .
4. Extensions
The a priori error estimate of Theorem 3.4 can be easily extended to the case of general boundary conditions and to general numerical fluxes.
4.1. Other boundary conditions
Theorem 3.4 holds unchanged for Neumann, Robin or mixed boundary conditions. To see this, let us consider, for example, the following Neumann boundary conditions:
First, we take
Then we redefine the numerical flux as follows:
Note that this flux is obtained by setting and
Now, we simply have to go through the proof of Lemma 3.1 in section 3.2 to verify that the basic error estimate of Lemma 3.1 still holds with as before and
This and the approximation result of Lemma 3.3 imply that the estimate of Theorem 3.4 holds in this case.
4.2. General numerical fluxes
In the case of general numerical fluxes, the optimality of the estimate of Theorem 3.4 is lost in both and . Theoretically, the main reason is that now there are new terms in which were equal to zero for the flux (Equation 2.6).
Indeed, in the purely convective case, , the new term is
which is estimated as follows:
The last two terms are absorbed by the term and the first two are bounded, using Lemma 3.3, by
where
Hence, the error estimate is
Note the loss of half a power in both and .
In the case in which , the following additional term appears:
By using the approximation results of Lemma 3.3, we see that we can bound as follows:
where
and
Next, we use the following inverse inequality:
for in order to get
Note that this produces an additional loss of half power in and a full power in . Thus, after a few simple manipulations, we obtain the following estimate for general numerical fluxes:
5. Numerical results
The purpose of this section is to numerically validate the a priori error estimates given in section 3. In all our experiments, we use a TVD Runge–Kutta time stepping method (see Shu and Osher Reference 26Reference 27), with sufficiently small time steps, such that the overall error is governed by the spatial error.
5.1. Exponential convergence
Our first example illustrates the exponential convergence in for analytic solutions. We solve (Equation 2.1) on the space-time domain , with exact solution . We use a fixed grid consisting of a uniform mesh with 4 elements and increase the polynomial degree . The corresponding errors in the energy norm at time are shown in Figure 1. The diffusion coefficient is and the convection coefficient is chosen as (left) and (right). The curves clearly show exponential rates of convergence as predicted in (Equation 3.4) of section 3. Since the quadrature points and weights used to determine the LDG solution are computed only with an accuracy of , the curves bottom out for .
5.2. Optimal order of convergence in
In these examples, we show that an optimal order of convergence of is achieved when using the numerical flux in (Equation 2.6). For this set of tests, we solve (Equation 2.1) on , again with exact solution . To determine numerically the convergence order we consider the two sequences of successively refined meshes shown in Figure 2. Since our analysis is valid for arbitrary meshes, we choose the second sequence to consist of nonuniform meshes whereas the first one contains uniform meshes. Note that in both cases the mesh-size parameter of is half of the one of .
If denotes the error on the th mesh in the energy norm, then the numerical rate of convergence is defined as
In Tables 1, 2 and 3, we present these numerical orders in the energy norm at for polynomials of degree to on the above two mesh-sequences. In all the experiments we use the same convection coefficient and increase the diffusion coefficient from to . The results show that our estimates are optimal in not only for convection dominated problems, but also for diffusion dominated problems. In all the cases the numerical orders agree with the theoretical orders of our error estimates in Theorem 3.4.
5.3. Nonsmooth solutions
In this subsection, we present some numerical results to illustrate the performance of the LDG method for a solution that is nonsmooth in space.
We consider first the -version and start by solving the purely convective () problem (Equation 2.1), on with and with data chosen in such a way that the exact solution is . The corresponding uniform and nonuniform spatial discretizations are similar to those used in section 5.2 (cf., Figure 2). In this purely convective problem an order of convergence of is expected from our error estimate in section 3. These orders can clearly be seen in Table 4. Again, they are calculated for the energy norm at .
Now, we consider the linear problem with diffusion, i.e., with and . Again, we choose the data in such a way that the exact solution is . The results at are shown in Table 5. From our a priori error estimate, an order of convergence of is expected. This is what we actually see for all values of except for . In this case, we observe an order of convergence of instead of the predicted . Since the order of convergence for is smaller than , we believe that an error cancellation might be taking place which occurs only for .
Since the solution is singular at the mesh point , we expect a doubling of the convergence rate in the -version where is increased on a fixed mesh ; see, e.g., Schwab Reference 24. This is shown in Figure 3 for the same model problems as above. We can see a convergence rate of in the purely hyperbolic case () and of in the convection-diffusion case, respectively, which corresponds to an exact doubling of the rates. However, in our theoretical results, if we insert the weighted bounds from Lemma 3.2 in the proof of Theorem 3.4, we obtain rates of in the hyperbolic case and of in the convection-diffusion case, resulting in a loss of one power of and indicating the suboptimality of Lemma 3.2 with respect to the weighted spaces .
5.4. Testing the optimality of the smoothness requirement
To test if the smoothness on the exact solution required by our Theorem 3.4 when is optimal, we consider problem (Equation 2.1) with , , homogeneous Dirichlet boundary conditions and initial data . The results at are given in Table 6. Theorem 3.4 predicts an order of convergence of but we actually see an order of convergence of . This gives a strong indication that, to obtain optimal orders of convergence at least in , less smoothness of the exact solution than required in Theorem 3.4 for is sufficient. However, obtaining optimality in the smoothness of the exact solution seems to ask for more sophisticated theoretical techniques than the ones available in the current literature and has to be addressed in future work.
5.5. Robust exponential convergence
Our last example shows that robust exponential convergence can be obtained in the presence of a boundary layer, when suitable meshes are used; see Remark 3.8. We solve (Equation 2.1) on with , and right-hand side such that the exact solution is . For small , this solution has an exponential boundary layer of strength at the outflow boundary . In Figure 4, we compare the -version of the LDG method when using uniform and geometric meshes. Both meshes are chosen in such a way that they have the same number of elements; however, the distribution of the grid points is different: In the geometric mesh the size of the first element near is on the order of the length of the boundary layer, , the size of the next element is twice the size of the previous and so forth. The errors in the energy norm at for and are depicted in Figure 4. All the curves show exponential rates of convergence. However, the robustness of the rates for the geometric boundary layer meshes can clearly be observed, whereas the uniform mesh performs orders of magnitude worse for .
6. Concluding remarks
In this paper, we have obtained optimal error estimates for the -version of the LDG method for the model problem of the initial boundary value problem for a one-dimensional convection-diffusion equation. We have shown that this is possible by a careful choice of the numerical fluxes and the associated projections and ; we have also shown how this optimality in and is lost, at least theoretically, when general numerical fluxes are used.
Our numerical results confirm the optimality in of our main result and the exponential convergence that follows when the solution is analytic. These results also indicate that the smoothness requirement on the exact solution is too stringent. The problem of obtaining optimality in the smoothness of the exact solution seems to ask for more sophisticated theoretical techniques than the ones available in the current literature and constitute the subject of ongoing work. Also, extensions of our main result to the more challenging cases of nonconstant coefficients and , and to the multidimensional case will be considered elsewhere.
