Discontinuous Galerkin methods for nonlinear scalar hyperbolic conservation laws: divided difference estimates and accuracy enhancement.

In this paper, an analysis of the accuracy-enhancement for the discontinuous Galerkin (DG) method applied to one-dimensional scalar nonlinear hyperbolic conservation laws is carried out. This requires analyzing the divided difference of the errors for the DG solution. We therefore first prove that the [Formula: see text]-th order [Formula: see text] divided difference of the DG error in the [Formula: see text] norm is of order [Formula: see text] when upwind fluxes are used, under the condition that [Formula: see text] possesses a uniform positive lower bound. By the duality argument, we then derive superconvergence results of order [Formula: see text] in the negative-order norm, demonstrating that it is possible to extend the Smoothness-Increasing Accuracy-Conserving filter to nonlinear conservation laws to obtain at least [Formula: see text]th order superconvergence for post-processed solutions. As a by-product, for variable coefficient hyperbolic equations, we provide an explicit proof for optimal convergence results of order [Formula: see text] in the [Formula: see text] norm for the divided differences of DG errors and thus [Formula: see text]th order superconvergence in negative-order norm holds. Numerical experiments are given that confirm the theoretical results.

Introduction

In this paper, we study the accuracy-enhancement of semi-discrete discontinuous Galerkin (DG) methods for solving one-dimensional scalar conservation laws where is a given smooth function. We assume that the nonlinear flux function f(u) is sufficiently smooth with respect to the variable u and the exact solution is also smooth. For the sake of simplicity and ease in presentation, we only consider periodic boundary conditions. We show that the -th order divided difference of the DG error in the norm achieves th order when upwind fluxes are used, under the condition that possesses a uniform positive lower bound. By using a duality argument, we then derive superconvergence results of order in the negative-order norm. This allows us to demonstrate that it is possible to extend the post-processing technique to nonlinear conservation laws to obtain at least th order of accuracy. In addition, for variable coefficient hyperbolic equations that have been discussed in [19], we provide an explicit proof for optimal error estimates of order in the norm for the divided differences of the DG errors and thus in the negative-order norm.

Various superconvergence properties of DG methods have been studied in the past two decades, which not only provide a deeper understanding about DG solutions but are useful for practitioners. According to different measurements of the error, the superconvergence of DG methods is mainly divided into three categories. The first one is superconvergence of the DG error at Radau points, which is typically measured in the discrete norm and is useful to resolve waves. The second one is superconvergence of the DG solution towards a particular projection of the exact solution (supercloseness) measured in the standard norm, which lays a firm theoretical foundation for the excellent behaviour of DG methods for long-time simulations as well as adaptive computations. The last one is the superconvergence of post-processed solution by establishing negative-order norm error estimates, which enables us to obtain a higher order approximation by simply post-processing the DG solution with a specially designed kernel at the very end of the computation. In what follows, we shall review some superconvergence results for the aforementioned three properties and restrict ourselves to hyperbolic equations to save space. For superconvergence of DG methods for other types of PDEs, we refer to [21].

Let us briefly mention some superconvergence results related to the Radau points and supercloseness of DG methods for hyperbolic equations. Adjerid and Baccouch [1-3] studied the superconvergence property as well as the a posteriori error estimates of the DG methods for one- and two-dimensional linear steady-state hyperbolic equations, in which superconvergence of order and are obtained at downwind-biased Radau points and downwind points, respectively. Here and below, k is the highest polynomial degree of the discontinuous finite element space. For time-dependent linear hyperbolic equations, Cheng and Shu [9] investigated supercloseness for linear hyperbolic equations, and they obtained superconvergence of order towards a particular projection of the exact solution, by virtue of construction and analysis of the so-called generalized slopes. Later, by using a duality argument, Yang and Shu [24] proved superconvergence results of order of the DG error at downwind-biased points as well as cell averages, which imply a sharp th order supercloseness result. By constructing a special correction function and choosing a suitable initial discretization, Cao et al. [7] established a supercloseness result towards a newly designed interpolation function. Further, based on this supercloseness result, for the DG error they proved the th order superconvergence at the downwind points as well as domain average, -th order superconvergence at the downwind-biased Radau points, and superconvergent rate for the derivative at interior Radau points. We would like to remark that the work of [7, 24] somewhat indicates the possible link between supercloseness and superconvergence at Radau points. For time-dependent nonlinear hyperbolic equations, Meng et al. [18] proved a supercloseness result of order . For superconvergent posteriori error estimates of spatial derivative of DG error for nonlinear hyperbolic equations, see [4].

Let us now mention in particular some superconvergence results of DG methods regarding negative-order norm estimates and post-processing for hyperbolic equations. The basic idea of post-processing is to convolve the numerical solution by a local averaging operator with the goal of obtaining a better approximation and typically of a higher order. Motivated by the work of Bramble and Schatz in the framework of continuous Galerkin methods for elliptic equations [5], Cockburn et al. [11] established the theory of post-processing techniques for DG methods for hyperbolic equations by the aid of negative-order norm estimates. The extension of this post-processing technique was later fully studied by Ryan et al. in different aspects of problems, e.g. for general boundary condition [20], for nonuniform meshes [13] and for applications in improving the visualization of streamlines [22] in which it is labeled as a Smoothness-Increasing Accuracy-Conserving (SIAC) filter. For the extension of the SIAC filter to linear convection-diffusion equations, see [15].

By the post-processing theory [5, 11], it is well known that negative-order norm estimates of divided differences of the DG error are important tools to derive superconvergent error estimates of the post-processed solution in the norm. Note that for purely linear equations [11, 15], once negative-order norm estimates of the DG error itself are obtained, they trivially imply negative-order norm estimates for the divided differences of the DG error. However, the above framework is no longer valid for variable coefficient or nonlinear equations. In this case, in order to derive superconvergent estimates about the post-processed solution, both the norm and negative-order norm error estimates of divided differences should be established. In particular, for variable coefficient hyperbolic equations, although negative-order norm error estimates of divided differences are given in [19], the corresponding norm estimates are not provided. For nonlinear hyperbolic conservation laws, Ji et al. [16] showed negative-order norm estimates for the DG error itself, leaving the estimates of divided differences for future work.

For nonlinear hyperbolic equations under consideration in this paper, it is therefore important and interesting to address the above issues by establishing both the norm and negative-order norm error estimates for the divided differences. The major part of this paper is to show norm error estimates for divided differences, which are helpful for us to obtain a higher order of accuracy in the negative-order norm and thus the superconvergence of the post-processed solutions. We remark that it requires having a uniform positive lower bound due to the technicality of the proof. The generalization from purely linear problems [11, 15] to nonlinear hyperbolic equations in this paper involves several technical difficulties. One of these is how to establish important relations between the spatial derivatives and time derivatives of a particular projection of divided differences of DG errors. Even if the spatial derivatives of the error are switched to their time derivatives, it is still difficult to analyze the time derivatives of the error; for more details, see Sect. 3.2 and also the appendix. Another important technicality is how to construct a suitable dual problem for the divided difference of the nonlinear hyperbolic equations. However, it seems that it is not trivial for the two-dimensional extension, especially for establishing the relations between spatial derivatives and time derivatives of the errors. The main tool employed in deriving norm error estimates for the divided differences is an energy analysis. To deal with the nonlinearity of the flux functions, Taylor expansion is used following a standard technique in error estimates for nonlinear problems [25]. We would like to remark that the superconvergence analysis in this paper indicates a possible link between supercloseness and negative-order norm estimates.

This paper is organized as follows. In Sect. 2, we give the DG scheme for divided differences of nonlinear hyperbolic equations, and present some preliminaries about the discontinuous finite element space. In Sect. 3, we state and discuss the norm error estimates for divided differences of nonlinear hyperbolic equations, and then display the main proofs followed by discussion of variable coefficient hyperbolic equations. Section 4 is devoted to the accuracy-enhancement superconvergence analysis based on negative-order norm error estimates of divided differences. In Sect. 5, numerical experiments are shown to demonstrate the theoretical results. Concluding remarks and comments on future work are given in Sect. 6. Finally, in the appendix we provide the proofs for some of the more technical lemmas.

The DG scheme and some preliminaries

The DG scheme

In this section, we follow [11, 12] and present the DG scheme for divided differences of the problem (1.1).

Let be a partition of , and set . Since we are focused on error analysis of both the norm and the negative-order norm for divided differences of the DG solution and the problem under consideration is assumed to be periodic, we shall introduce two overlapping uniform (translation invariant) meshes for , namely and with mesh size . Associated with these meshes, we define the discontinuous finite element spacewhere denotes the set of polynomials of degree up to k defined on the cell . Here and below, represents the -th order divided difference of a given function, whose definition is given in (2.6a). Clearly, is a piecewise polynomial space on mesh for even (including ) and a piecewise polynomial space on mesh for odd of the DG scheme. For simplicity, for even , we would like to use to denote the standard finite element space of degree k defined on the cell , if there is no confusion. Since functions in may have discontinuities across element interfaces, we denote by and the values of w(x) at the discontinuity point from the left cell and the right cell, respectively. Moreover, we use and to represent the jump and the mean of w(x) at each element boundary point.

The -th order divided difference of the nonlinear hyperbolic conservation law is where with . Clearly, (2.1) reduces to (1.1) when . Then the approximation of the semi-discrete DG method for solving (2.1) becomes: find the unique function (and thus ) such that the weak formulationholds for all and all Note that, by (2.6a), for any and t, is a linear combination of the values of at equally spaced points of length h, namely . Here and in what follows, denotes the usual inner product in , and is the DG spatial discretization operator defined on each cell , namelyWe point out that in order to obtain a useful bound for the norm error estimates of divided differences, the numerical flux is chosen to be an upwind flux, for example, the well-known Godunov flux. Moreover, the analysis requires a condition that has a uniform positive lower bound. Without loss of generality, throughout the paper, we only consider , and thus . Therefore,

For periodic boundary conditions under consideration in this paper, we use to denote the summation of with respect to cell , that is where represents the inner product in . Note that we have used the summation with respect to j instead of to distinguish two overlapping meshes, since for even and for odd .

Preliminaries

We will adopt the following convention for different constants. We denote by C a positive constant independent of h but may depend on the exact solution of the Eq. (2.1), which could have a different value in each occurrence. To emphasize the nonlinearity of the flux f(u), we employ to denote a nonnegative constant depending solely on the maximum of a high order derivative (). We remark that for a linear flux function or a variable coefficient flux function , where c is a constant and a(x) is a given smooth function.

Prior to analyzing the norm and the negative-order norm error estimates of divided differences, let us present some notation, definitions, properties of DG discretization operator, and basic properties about SIAC filters. These preliminary results will be used later in the proof of superconvergence property.

Sobolev spaces and norms

We adopt standard notation for Sobolev spaces. For any integer , we denote by the Sobolev space on subdomain equipped with the norm . In particular, if , we set , and , and further if , we set . Throughout the paper, when , we will omit the index D for convenience. Furthermore, the norms of the broken Sobolev spaces with being the union of all cells can be defined analogously. The Bochner space can also be easily defined. For example, the space is equipped with the norm .

Additionally, we denote by the standard norm on the cell interfaces of the mesh . Specifically, for the one-dimensional case under consideration in this paper, with . To simplify notation in our later analysis, following [23], we would like to introduce the so-called jump seminorm for .

In the post-processing framework, it is useful to consider the negative-order norm, defined as: Given and domain ,

Properties for divided differences

For any function w and integer , the following standard notation of central divided difference is used Note that the above notation is still valid even if w is a piecewise function with possible discontinuities at cell interfaces. In later analysis, we will repeatedly use the properties about divided differences, which are given as follows. For any functions w and v which is the so-called Leibniz rule for the divided difference. Moreover, for sufficiently smooth functions w(x), by using Taylor expansion with integral form of the remainder, we can easily verify that is a second order approximation to , namelywhere is a positive constant and is a smooth function; for example, for , andHere and below, denotes the -th order partial derivative of a function with respect to the variable x; likewise for . The last identity is the so-called summation by parts, namelyIn addition to the properties of divided differences for a single function w(x), the properties of divided differences for a composition of two or more functions are also needed. However, we would like to emphasize that properties (2.6a), (2.6b), (2.6d) are always valid whether w is a single function or w is a composition of two or more functions. As an extension from the single function case in (2.6c) to the composite function case, the following property (2.6e) is subtle, it requires a more delicate argument for composite functions. Without loss of generality, if w is the composition of two smooth functions r and u, i.e., , we can prove the following identitywhere is a positive constant and is a smooth function. We can see that, unlike (2.6c), the divided difference of a composite function is a first order approximation to its derivative of the same order. This finding, however, is sufficient in our analysis; see Corollary 1.

It is worth pointing out that in (2.6e), and should be understood in the sense of the chain rule for high order derivatives and divided differences of composite functions, respectively. In what follows, we use to denote the standard -th order Newton divided difference, that isIt is easy to verify thatwhere .

For completeness, we shall list the chain rule for the derivatives (the well-known Faà di Bruno’s Formula) and also for the divided differences [14]; it readswhere , and the sum is over all and non-negative integer solutions toandwith the sum being over integers such that .

It follows from the mean value theorem for divided differences thatConsequently, by (2.7),We are now ready to prove (2.6e) for the relation between the divided difference and the derivative of composite functions. Using a similar argument as that in the proof of (2.6c) for and the relation thatdue to the smoothness of and the fact that may not necessarily be equally spaced, with and being smooth functions, we can obtain the relation (2.6e).

The inverse and projection properties

Now we list some inverse properties of the finite element space . For any , there exists a positive constant C independent of p and h, such thatNext, we introduce the standard projection of a function into the finite element space , denoted by , which is a unique function in satisfyingNote that the proof of accuracy-enhancement of DG solutions for linear equations relies only on an projection of the initial condition [11, 15]. However, for variable coefficient and nonlinear hyperbolic equations, a suitable choice of the initial condition and a superconvergence relation between the spatial derivative and time derivative of a particular projection of the error should be established, since both the norm and negative-order norm error estimates of divided differences need to be analyzed. In what follows, we recall two kinds of Gauss–Radau projections into following a standard technique in DG analysis [8, 25]. For any given function and an arbitrary element the special Gauss–Radau projection of q, denoted by , is the unique function in satisfying, for each , We would like to remark that the exact collocation at one of the end points on each cell plus the orthogonality property for polynomials of degree up to of the Gauss–Radau projections play an important role and are used repeatedly in the proof. We denote by or the projection error, then by a standard scaling argument [6, 10], it is easy to obtain, for smooth enough q(x), that, Moreover,

The properties of the DG spatial discretization

To perform the error estimates of divided differences, several properties of the DG operator are helpful, which are used repeatedly in our proof; see Sect. 3.

Lemma 1

Suppose that r(u(x, t)) etc) is smooth with respect to each variable. Then, for any , there holds the following inequality and in particular, if , there holds

Proof

Let us first prove (2.11b), which is a straightforward consequence of the definition of , since, after a simple integration by partsWe would like to emphasize that (2.11b) is still valid even if the smooth function r and depend on different x, e.g. etc, as only integration by parts as well as the boundedness of r is used here.

To prove (2.11a), we consider the equivalent strong form of (2.4b). An application of Cauchy–Schwarz inequality and inverse inequality (ii) leads toThis completes the proof of Lemma 1.

Corollary 1

Under the same conditions as in Lemma 1, we have, for small enough h,

The case has been proved in Lemma 1. For general , let us start by using the relation (2.6e) for to obtainwith C a positive constant and a smooth function. Next, applying (2.11a) in Lemma 1 to and , we have for small enough h This finishes the proof of Corollary 1.

Lemma 2

Suppose that r(u(x, t)) is smooth with respect to each variable. Then, for any and , there holds

Using the definition of the projection (2.9a), we have that , and thusNext, on each cell , we rewrite r(u(x, t)) as with . Clearly, on each element , due to the smoothness of r and u. Using the orthogonality property of again (2.9b), we arrive atwhere we have used Cauchy–Schwarz inequality, inverse inequality (i) and the approximation property (2.10a) consecutively.

Corollary 2

Suppose that r(u(x, t)) is smooth with respect to each variable. Then, for any , , there holds

The case has been proved in Lemma 2. For , by the Leibniz rule (2.6b) and taking into account the fact that both the divided difference operator and the projection operator are linear, we rewrite aswithThus,Clearly, by (2.6e), is also a smooth function with respect to each variable with leading term . To complete the proof, we need only apply the same procedure as that in the proof of Lemma 2 to each term on the right side of (2.15).

Regularity for the variable coefficient hyperbolic equations

Since the dual problem for the nonlinear hyperbolic equation is a variable coefficient equation, we need to recall a regularity result.

Lemma 3

[16] Consider the variable coefficient hyperbolic equation with a periodic boundary condition for all where a(x, t) is a given smooth periodic function. For any , fix time t and , then the solution of (2.16) satisfies the following regularity propertywhere C is a constant depending on .

SIAC filters

The SIAC filters are used to extract the hidden accuracy of DG methods, by means of a post-processing technique, which enhances the accuracy and reduces oscillations of the DG errors. The post-processing is a convolution with a kernel function that is of compact support and is a linear combination of B-splines, scaled by the uniform mesh size,where is the B-spline of order obtained by convolving the characteristic function of the interval with itself k times. Additionally, the kernel function should reproduce polynomials of degree by convolution, which is used to determine the weights . For more details, see [11].

The post-processing theory of SIAC filters is given in the following theorem.

Theorem 1

(Bramble and Schatz [5]) For , where is the maximal time of existence of the smooth solution, let be the exact solution of (1.1). Let and U be any approximation to u, thenwhere and depend on , but is independent of h.

norm error estimates for divided differences

By the post-processing theory [5, 11] (also see Theorem 1), it is essential to derive negative-order norm error estimates for divided differences, which depend heavily on their norm estimates. However, for both variable coefficient equations and nonlinear equations, it is highly nontrivial to derive norm error estimates for divided differences, and the technique used to prove convergence results for the DG error itself needs to be significantly changed.

The main results in norm

Let us begin by denoting to be the error between the exact solution and numerical solution. Next, we split it into two parts; one is the projection error, denoted by , and the other is the projection of the error, denoted by . Here the projection is defined at each time level t corresponding to the sign variation of ; specifically, for any and , if we choose , and if , we take .

We are now ready to state the main theorem for the norm error estimates.

Theorem 2

For any , let be the exact solution of Eq. (2.1), which is assumed to be sufficiently smooth with bounded derivatives, and assume that is uniformly lower bounded by a positive constant. Let be the numerical solution of scheme (2.2) with initial condition when the upwind flux is used. For a uniform mesh of , if the finite element space of piecewise polynomials with arbitrary degree is used, then for small enough h and any there holds the following error estimatewhere the positive constant depends on the u, , T and f, but is independent of h.

Corollary 3

Under the same conditions as in Theorem 2, if in addition we have the following error estimates:

Proof

As shown in Corollary 2, we have that , and thusby the approximation error estimate (2.10a). Now, the error estimate (3.2) follows by combining the triangle inequality and (3.1).

Remark 1

Clearly, the error estimates for the divided differences in Theorem 2 and Corollary 3 also hold for the variable coefficient equation (2.1) with and . In fact, for variable coefficient equations, we can obtain optimal th order in the norm and thus th order in the negative-order norm; see Sect. 3.3.

Remark 2

The result with in Theorem 2 is indeed a superconvergence result towards a particular projection of the exact solution (supercloseness) that has been established in [18], which is a starting point for proving with . For completeness, we list the superconvergence result for (zeroth order divided difference) as follows where, on each element , we have used with being a constant and . Note that the proof of such superconvergence results requires that is uniformly lower bounded by a positive constant; for more details, see [18].

In the proof of Theorem 2, we have also obtained a generalized version about the norm estimates of in terms of the divided differences, their time derivatives, and spatial derivatives. To simplify notation, for an arbitrary multi-index , we denote by the mixed operator containing divided differences and time derivatives of a given function, namely

Corollary 4

Under the same conditions as in Theorem 2, for and a multi-index with , we have the following unified error estimate:where .

Proof of the main results in the norm

Similar to the discussion of the DG discretization operator properties in Sect. 2.2.4, without loss of generality, we will only consider the case for all ; the case of is analogous. Therefore, we take the upwind numerical flux as on each cell interface and choose the projection as on each cell, and the initial condition is chosen as . Since the case has already been proven in [18] (see (3.4a)), we need only to consider . In order to clearly display the main ideas of how to perform the norm error estimates for divided differences, in the following two sections we present the detailed proof for Theorem 2 with and , respectively; the general cases with can be proven by induction, which are omitted to save space.

Analysis for the first order divided difference

For , the DG scheme (2.2) becomeswith , which holds for any and . By Galerkin orthogonality and summing over all , we have the error equationfor all To simplify notation, we would like to denote with . If we now take , we get the following identityThe estimate for the right side of (3.7) is complicated, since it contains some integral terms involving mixed order divided differences of , namely and , which is given in the following lemma.

Lemma 4

Suppose that the conditions in Theorem 2 hold. Then we havewhere the positive constants C and are independent of h and

Let us start by using the second order Taylor expansion with respect to the variable u to write out the nonlinear terms, namely and , as where and are the integral form of the remainders of the second order Taylor expansion. We would like to emphasize that the various order spatial derivatives, time derivatives and divided differences of are all bounded uniformly due to the smoothness of f and u. Thus,which will be estimated separately below.

To estimate , we employ the Leibniz rule (2.6b), and rewrite asand thus,where we have omitted the dependence of x for convenience if there is no confusion, since the proof of (2.11b) is still valid even if and are evaluated at different x; see proof of (2.11b) in Sect. 2.2.4. A direct application of Lemma 1 together with the assumption that , (2.11b), leads to the estimate for : By Corollary 1, we arrive at the estimate for : Substituting (3.4a)–(3.4c) into (3.10b), and combining with (3.10a), we have, after a straightforward application of Young’s inequality, thatLet us now move on to the estimate of . By Corollary 2, we haveTo estimate , let us first employ the identity (2.6b) and rewrite asConsequently,It is easy to show, for the high order nonlinear term , thatwhere in the first step we have used the Cauchy–Schwarz inequality, in the second step we have used the inverse properties (i) and (ii), and in the last step we have employed the interpolation properties (3.3). We see that in order to deal with the nonlinearity of f we still need to have a bound for . Due to the superconvergence result (3.4a), we conclude, by combining inverse inequality (iii) and the approximation property (2.10b), thatTherefore, for small enough h, we have By using analysis similar to that in the proof of (3.13), we have, for , thatAs a consequence, by (3.14) and (3.4a) A combination of (3.15a) and (3.15b) produces a bound for :To complete the proof of Lemma 4, we need only combine (3.11), (3.12), (3.16) and use Young’s inequality.

We are now ready to derive the norm estimate for . To do this, let us begin by inserting the estimate (3.8) into (3.7) and taking into account the bound for in (3.3) and thus to get, after an application of Cauchy–Schwarz inequality and Young’s inequality, thatNext, we integrate the above inequality with respect to time between 0 and T and note the fact that due to to obtainwhere we have used the superconvergence result (3.4a). An application of Gronwall’s inequality leads to the desired resultThis finishes the proof of Theorem 2 for .

Remark 3

We can see that the estimates (3.17) for the norm and the jump seminorm of are based on the corresponding results for in Remark 2, which are half an order lower than that of . This is mainly due to the hybrid of different order divided differences of , namely and , and thus the application of inverse property (ii). It is natural that the proof for the high order divided difference of , say , should be based on the corresponding lower order divided difference results of ( and ) that have already been established; see Sect. 3.2.2 below.

Analysis for the second order divided difference

For , the DG scheme (2.2) becomeswith , which holds for any and . By Galerkin orthogonality and summing over all j, we have the error equationfor all To simplify notation, we would like to denote with . If we now take , we get the following identityThe estimate for right side of (3.19) is rather complicated, since it contains some integral terms involving mixed order divided differences of , namely , and , which is given in the following Proposition.

Proposition 1

Suppose that the conditions in Theorem 2 hold. Then we havewhere the positive constants C and are independent of h and

By the second order Taylor expansion (3.9), we havewhich will be estimated one by one below.

To estimate , we use the Leibniz rule (2.6b), to rewrite asand thus,where we have omitted the dependence of x for convenience if there is no confusion. A direct application of Lemma 1 together with the assumption that , (2.11b), produces the estimate for : By Corollary 1, we arrive at the estimates for and : Substituting (3.4a)–(3.4c), (3.17) into (3.22b), (3.22c), and combining with (3.22a), we have, after a straightforward application of Young’s inequality, thatFor terms on the right side of (3.23), although we have information about and as shown in (3.4a) and (3.17), we still need a suitable bound for , which is given in the following lemma.

Lemma 5

Suppose that the conditions in Theorem 2 hold. Then we havewhere depends on u and but is independent of h and .

The proof of this lemma is given in the appendix. Up to now, we see that we still need to have a bound for . In fact, the proof for would require additional bounds for and , whose results are shown in Lemmas 6 and 7.

Lemma 6

Suppose that the conditions in Theorem 2 hold. Then we have

The proof of Lemma 6 follows along a similar argument as that in the proof of Lemma 5, so we omit the details here.

Lemma 7

Suppose that the conditions in Theorem 2 hold. Then we have

The proof of this lemma is deferred to the appendix. Based on the above two lemmas, we are able to prove the bound for in Lemma 8, whose proof is deferred to the appendix.

Lemma 8

Suppose that the conditions in Theorem 2 hold. Then we havewhere depends on u and but is independent of h and .

We now collect the estimates in Lemmas 5 and 8 into (3.23) to getLet us now move on to the estimate of . By Corollary 2, we haveTo estimate , let us first employ the identity (2.6b) and rewrite aswhereThus,By using analysis similar to that in the proof of (3.13), we getwhere we have used the fact that for and small enough h, ; for more details, see the appendix. ConsequentlyCollecting the estimates (3.28)–(3.30) into (3.21) and taking into account (3.4a) and (3.17), we getThis finishes the proof of Proposition 1.

We are now ready to derive the norm estimate for . To do this, we begin by combining (3.19) and (3.20) to getNext, integrate the above inequality with respect to time between 0 and T and use (and thus ) to obtainby the estimates (3.4a) and (3.17). An application of Gronwall’s inequality leads to the desired resultThis completes the proof of Theorem 2 with .

Remark 4

Through the proof of Theorem 2 with , , we can see that apart from the bounds for that have already been obtained for proving , we require additional bounds for and , as shown in Lemmas 5–8. The proof for the norm estimates for higher order divided differences are more technical and complicated, and it would require bounds regarding lower order divided differences as well as its corresponding spatial and time derivatives. For example, when , in addition to the abounds aforementioned, we need to establish the bounds for and . Thus, Theorem 2 can be proven along the same lines for general . Finally, we would like to point out that the corresponding results on the jump seminorm for various order divided differences and time derivatives of are useful, which play an important role in deriving Theorem 2.

Variable coefficient case

The main results

In this section we consider the error estimates for divided differences for the variable coefficient equation (1.1) with . Similar to the nonlinear hyperbolic case, to obtain a suitable bound for the norm the numerical flux should be chosen as an upwind flux. Moreover, the analysis requires a condition that |a(x)| is uniformly lower bounded by a positive constant. Without loss of generality, we only consider , and thus the DG scheme isfor . We will use the same notation as before.

For nonlinear hyperbolic equations, the loss of order in Theorem 2 is mainly due to the lack of control for the interface jump terms arising from (2.11a) in the superconvergence relation, for example, (3.4b), (3.24) and (3.25). Fortunately, for variable coefficient hyperbolic equations, we can establish a stronger superconvergence relation between the spatial derivative as well as interface jumps of the various order divided difference of and its time derivatives; see (3.37b) below. Thus, optimal error estimates of order are obtained.

Prior to stating our main theorem, we would like to present convergence results for time derivatives of , which is slightly different to those for nonlinear hyperbolic equations.

Lemma 9

Let u be the exact solution of the variable coefficient hyperbolic Eq. (1.1) with , which is assumed to be sufficiently smooth with bounded derivatives. Let be the numerical solution of scheme (3.32) with initial condition when the upwind flux is used. For regular triangulations of , if the finite element space of piecewise polynomials with arbitrary degree is used, then for any and any there holds the following error estimatewhere the positive constant C depends on u, T and a, but is independent of h.

The proof of this lemma is postponed to the appendix.

We are now ready to state our main theorem.

Theorem 3

For any , let be the exact solution of the problem (2.1) with , which is assumed to be sufficiently smooth with bounded derivatives, and assume that |a(x)| is uniformly lower bounded by a positive constant. Let be the numerical solution of scheme (3.32) with initial condition when the upwind flux is used. For a uniform mesh of , if the finite element space of piecewise polynomials with arbitrary degree is used, then for any there holds the following error estimatewhere the positive constant C depends on u, , T and a, but is independent of h.

Remark 5

Based on the optimal error estimates for together with approximation error estimates (3.3) and using the duality argument in [19], we can obtain the negative-order norm estimatesand thusFor more details, see [5, 19] and also Sect. 4 below.

Proof of main results

We shall prove Theorem 3 for general . First we claim that if we can prove the following three inequalities where represents the mixed operator containing divided differences and time derivatives of that has already been defined in (3.5), then . In what follows, we sketch the verification of this claim. To do that, we start by taking in the following error equationwhich isNext, consider the term . Use Leibniz rule (2.6b) to rewrite and employ (2.11a), (2.11b) in Lemma 1 to get the boundwhere we have also used the relations (3.37a)–(3.37c). For the estimate of , we need only use Corollary 2 to getCollecting above two estimates into (3.38) and using Cauchy–Schwarz inequality as well as Gronwall’s inequality, we finally getThe claim is thus verified.

In what follows, we will prove (3.37) by induction.

Step 1  For , is well known, and thus (3.37a) is valid for . Moreover, (3.37c), namely has been given in (3.4c); see [18]. To complete the proof for , we need only to establish the following relation

Noting the relation (3.4b), we need only to proveTo do that, we consider the cell error equationwhich holds for any and . If we now take in the above identity and use the strong form (2.3b) for , we getIt follows from the assumption thatBy Cauchy–Schwarz inequality, we haveBy the definition of the projection , (2.9b)Inserting the above two estimates into (3.41), we arrive atwhich iswhere we have used the bound for , the relation (3.4b) and approximation error estimates (2.10a), and thus (3.40) follows. Therefore, (3.37) is valid for .

Step 2  Suppose that (3.37) is true for . That is let us prove that it also holds for .

First, as shown in our claim, (3.42) implies thatThe above estimate together with (3.42a) producesTherefore, (3.37a) is valid for .

Next, by assumption (3.42b), we can see that to show (3.37b) for , we need only to showWithout loss of generality, let us take for example. To this end, we consider the following error equationwhich holds for any . We use Leibniz rule (2.6b) to write out asTherefore, the error equation becomeswhere for . Let us now work on . By the strong form of , (2.4b), we haveDenote the standard Legendre polynomials of degree k in . If we now let with being a constant and , we getsince . Substituting above expression into (3.44) and taking into account the assumption that , we have It is easy to show by Corollary 1 thatwhere we have used (3.42a)–(3.42c), since for . By Corollary 2, we haveBy (3.43) and inverse property (i), we arrive at a bound for and The triangle inequality and the approximation error estimate (3.3) yield Collecting the estimates (3.46a)–(3.46d) into (3.45) and using the fact that , we arrive atIf we take in the cell error equation and use an analysis similar to that in the proof of (3.40), we will get the following relationA combination of (3.47) and (3.48) gives usTherefore, (3.37b) still holds for .

Finally, let us verify that (3.37c) is valid for . Noting the assumption (3.42c), we need only consider and . To do that, we start from the estimate for with that has already been established in (3.33). By an analysis similar to that in the proof of Lemma 8 and taking into account relations (3.37a) and (3.37b) for , we conclude that (3.37c) is valid for . Repeating the above procedure, we can easily verify that (3.37c) is also valid for . Therefore, (3.37c) holds true for , and thus (3.34) in Theorem 3 is valid for general .

Superconvergent error estimates

For nonlinear hyperbolic equations, the negative-order norm estimate of the DG error itself has been established in [16]. However, by post-processing theory [5, 11], negative-order norm estimates of divided differences of the DG error are also needed to obtain superconvergent error estimates for the post-processed solution in the norm. Using a duality argument together with norm estimates established in Sect. 3, we show that for a given time T, the -th order divided difference of the DG error in the negative-order norm achieves th order superconvergence. As a consequence, the DG solution , converges with at least th order in the norm when convolved with a particularly designed kernel.

We are now ready to state our main theorem about the negative-order norm estimates of divided differences of the DG error.

Theorem 4

For any , let be the exact solution of the problem (2.1), which is assumed to be sufficiently smooth with bounded derivatives, and assume that is uniformly lower bounded by a positive constant. Let be the numerical solution of scheme (2.2) with initial condition when the upwind flux is used. For a uniform mesh of , if the finite element space of piecewise polynomials with arbitrary degree is used, then for small enough h and any there holds the following error estimatewhere the positive constant C depends on u, , T and f, but is independent of h.

Combining Theorems 4 and 1, we have

Corollary 5

Under the same conditions as in Theorem 4, if in addition is a convolution kernel consisting of B-splines of order such that it reproduces polynomials of degree , then we havewhere .

Remark 6

The th order superconvergence is shown for the negative norm, and thus is valid for B-splines of order (by Theorem 1). For general order of B-splines and , using similar argument for the proof of the negative norm estimates (see Sect. 4.1), we can prove the following superconvergent error estimateTherefore, from the theoretical point of view, a higher order of B-splines may lead to a superconvergence result of higher order, for example and thus th order in Corollary 5. However, from the practical point of view, changing the order of B-splines does not affect the order of superconvergence; see Sect. 5 below and also [17].

Proof of the main results in the negative-order norm

Similar to the proof for the norm estimates of the divided differences in Sect. 3.2, we will only consider the case for all . To perform the analysis for the negative-order norm, by (2.5), we need to concentrate on the estimate offor . To do that, we use the duality argument, following [16, 19]. For the nonlinear hyperbolic Eq. (2.1), we choose the dual equation as: Find a function such that is periodic for all and

Unlike the purely linear case [11, 15] or the variable coefficient case [19], the dual equations for nonlinear problems will no longer preserve the inner product of original solution and its dual solution , namely . In fact, if we multiply (2.1a) by and (4.4a) by and integrate over , we get, after using integration by parts and summation by parts (2.6d), thatwhereNote that is the same as that in [16] when . We now integrate (4.5) with respect to time between 0 and T to obtain a relation in different time levelIn what follows, we work on the estimate of (4.3). To do that, let us begin by using the relation (4.6) to get an equivalent form of (4.3). It reads, for any wherewill be estimated one by one below.

Note that in our analysis for in Theorem 2, we need to choose a particular initial condition, namely instead of for purely linear equations [11, 15]. Thus, we arrive at a slightly different bound for , as shown in the following lemma. We note that using the projection in the numerical examples is still sufficient to obtain superconvergence.

Lemma 10

(Projection estimate) There exists a positive constant C, independent of h, such that

Since , we have the following identitywhere is the projection into . A combination of Cauchy–Schwarz inequality and approximation error estimates (2.10a) leads to the desired result (4.7).

The bound for is given in the following lemma.

Lemma 11

(Residual) There exists a positive constant C, independent of h, such that

Denoting by G the term inside the time integral of , we get, by taking , the following expression for G,which is equivalent towhere we have added and subtracted the term and used integration by parts.

Let us now consider the estimates of . For , by using the second order Taylor expansion for , (3.9), we getwhere and , respectively, represent the linear part and the nonlinear part of . It is easy to show, by using the Leibniz rule (2.6b) and Cauchy–Schwarz inequality, that where we have used the estimate for in Corollary 3 and the approximation error estimate (2.10a). Analogously, for high order nonlinear term , we havewhere we have used the (2.6b) twice, the inverse property (iii), the norm estimate (3.2), and the approximation error estimate (2.10a). A combination of above two estimates yieldsTo estimate , we use an analysis similar to that in the proof of in Lemma 10 and make use of the orthogonal property of the projection to getwhere we have used the approximation error estimate (2.10a).

We proceed to estimate . It follows from the Taylor expansion (3.9), the Leibniz rule (2.6b), the Cauchy–Schwarz inequality and the inverse properties (ii), (iii) thatwhere we have also used (3.2) and (2.10a). Collecting the estimates (4.10)–(4.12), we getConsequently, the estimate for follows by integrating the above inequality with respect to time.

We move on to the estimate of , which is given in the following lemma.

Lemma 12

(Consistency) There exists a positive constant C, independent of h, such that

To do that, let us denote by the term inside the integral and take into account (2.6d) to obtain an equivalent form of where we have used the dual problem (4.4) and the fact that due to the smoothness of . Next, by using the second order the Taylor expansion (3.9) and (2.6d) again, we arrive atIf we now use (2.6b) twice for and the Cauchy–Schwarz inequality together with the error estimate (3.2), we getwhere we have also used the Sobolev inequality , under the condition that . The bound for follows immediately by integrating the above inequality with respect to time.

We are now ready to obtain the final negative-order norm error estimates for the divided differences. By collecting the results in Lemmas 10–12 and taking into account the regularity result in Lemma 3, namely , we get a bound for Thus, by (2.5), we have the bound for the negative-order normThis finishes the proof of Theorem 4.

Numerical examples

For nonlinear hyperbolic equations, we proved norm superconvergence results of order for post-processed errors, as shown in Corollary 5. The superconvergence results together with the post-processing theory by Bramble and Schatz in Theorem 1 entail us to design a more compact kernel to achieve the desired superconvergence order. We note that superconvergence of post-processed errors using the standard kernel (a kernel function composed of a linear combination of B-splines of order ) for nonlinear hyperbolic equations has been numerically studied in [11, 16]. Note that the order of B-splines does not have significant effect on the rate of convergence numerically and that it is the number of B-splines that has greater effect to the convergence order theoretically [11], we will only focus on the effect of different total numbers (denoted by with ) of B-splines of the kernel in our numerical experiments. For more numerical results using different orders of B-splines, we refer the readers to [17].

We consider the DG method combined with the third-order Runge–Kutta method in time. We take a small enough time step such that the spatial errors dominate. We present the results for and polynomials only to save space, in which a specific value of is chosen to match the orders given in Corollary 5. For the numerical initial condition, we take the standard projection of the initial condition and we have observed little difference if the projection is used instead. Uniform meshes are used in all experiments. Only one-dimensional scalar equations are tested, whose theoretical results are covered in our main theorems.

Example 1

We consider the Burgers quation on the domain with periodic boundary conditions.

Noting that changes its sign in the computational domain, we use the Godunov flux, which is an upwind flux. The errors at , when the solution is still smooth, are given in Table 1. From the table, we can see that one can improve the order of convergence from to at least , which is similar to the results for Burgers equations in [11]. Moreover, superconvergence of order 2k can be observed for the compact kernel with , as, in general, a symmetric kernel could yield one additional order. This is why instead of , is chosen in our kernel. The pointwise errors are plotted in Fig. 1, which show that the post-processed errors are less oscillatory and much smaller in magnitude for a large number of elements as observed in [11], and that the errors of our more compact kernel with are less oscillatory than that for the standard kernel with , although the magnitude of the errors increase. This example demonstrates that the superconvergence result also holds for conservation laws with a general flux function.

Table 1

Before post-processing (left), after post-processing (middle) and post-processing with the more compact kernel (right). T = 0.3 - and errors for Example 1

Mesh	Before post-processing				Post-processed (ω=0)				Post-processed (ω=-2)
	L^2 error	Order	L^∞ error	Order	L^2 error	Order	L^∞ error	Order	L^2 error	Order	L^∞ error	Order
P^2
20	1.54E−04	–	5.70E−04	–	1.04E−04	–	3.16E−04	–	5.36E−04	–	1.40E−03	–
40	2.06E−05	2.90	1.03E−04	2.47	2.28E−06	5.52	7.53E−06	5.39	3.69E−05	3.86	9.93E−05	3.81
80	2.73E−06	2.92	1.55E−05	2.73	3.97E−08	5.84	1.38E−07	5.77	2.37E−06	3.96	6.43E−06	3.95
160	3.56E−07	2.93	2.25E−06	2.78	1.13E−09	5.13	9.86E−09	3.81	1.49E−07	3.99	4.06E−07	3.99
P^3
20	7.68E−06	–	2.91E−05	–	5.88E−05	–	1.88E−04	–	1.59E−04	–	4.80E−04	–
40	5.21E−07	3.88	2.36E−06	3.62	5.47E−07	6.75	1.97E−06	6.58	3.71E−06	5.42	1.21E−05	5.31
80	3.45E−08	3.92	1.74E−07	3.76	2.87E−09	7.57	1.09E−08	7.50	6.56E−08	5.82	2.20E−07	5.78
160	2.23E−09	3.95	1.19E−08	3.87	1.22E−11	7.88	4.70E−11	7.86	1.06E−09	5.95	3.58E−09	5.94

Fig. 1

The errors in absolute value and in logarithmic scale for (top) and (bottom) polynomials with and 160 elements for Example 1 where . Before post-processing (left), after post-processing (middle) and post-processing with the more compact kernel (right).

Example 2

In this example we consider the conservation laws with more general flux functions on the domain with periodic boundary conditions.

We test the Example 2 at before the shock is developed. The orders of convergence with different kernels are listed in Table 2 and pointwise errors are plotted in Fig. 2. We can see that the post-processed errors are less oscillatory and much smaller in magnitude for most of elements as observed in [16], and that the errors of our more compact kernel with are slightly less oscillatory than that for the standard kernel with . This example demonstrates that the accuracy-enhancement technique also holds true for conservation laws with a strong nonlinearity that is not a polynomial of u.

Table 2

Before post-processing (left), after post-processing (middle) and post-processing with the more compact kernel (right). T = 0.1 - and errors for Example 2

Mesh	Before post-processing				Post-processed (ω=0)				Post-processed (ω=-2)
	L^2 error	Order	L^∞ error	Order	L^2 error	Order	L^∞ error	Order	L^2 error	Order	L^∞ error	Order
P^2
20	1.25E−04	–	5.76E−04	–	4.45E−05	–	1.61E−04	–	2.49E−04	–	7.98E−04	–
40	1.61E−05	2.95	7.64E−05	2.91	1.01E−06	5.46	4.03E−06	5.32	1.68E−05	3.88	5.73E−05	3.80
80	1.96E−06	3.04	1.02E−05	2.91	1.80E−08	5.81	7.35E−08	5.78	1.08E−06	3.97	3.72E−06	3.95
160	2.45E−07	3.00	1.32E−06	2.95	3.02E−10	5.90	1.25E−09	5.88	6.77E−08	3.99	2.35E−07	3.99
P^3
20	3.99E−06	–	2.52E−05	–	2.50E−05	–	9.12E−05	–	6.64E−05	–	2.38E−04	–
40	2.62E−07	3.93	1.67E−06	3.91	2.41E−07	6.70	1.00E−06	6.51	1.57E−06	5.40	6.17E−06	5.27
80	1.68E−08	3.96	1.13E−07	3.89	1.29E−09	7.55	5.66E−09	7.47	2.79E−08	5.81	1.14E−07	5.76
160	1.04E−09	4.01	7.38E−09	3.93	5.45E−12	7.88	2.45E−11	7.85	4.51E−10	5.95	1.86E−09	5.94

Fig. 2

The errors in absolute value and in logarithmic scale for (top) and (bottom) polynomials with and 160 elements for Example 2 where . Before post-processing (left), after post-processing (middle) and post-processing with the more compact kernel (right).

Concluding remarks

In this paper, the accuracy-enhancement of the DG method for nonlinear hyperbolic conservation laws is studied. We first prove that the -th order divided difference of the DG error in the norm is of order when piecewise polynomials of degree k and upwind fluxes are used, provided that is uniformly lower bounded by a positive constant. Then, by a duality argument, the corresponding negative-order norm estimates of order are obtained, ensuring that the SIAC filter will achieve at least th order superconvergence. As a by-product, we show, for variable coefficient hyperbolic equations with , the optimal error estimates of order for the norm of divided differences of the DG error, provided that |a(x)| is uniformly lower bounded by a positive constant. Consequently, the superconvergence result of order is obtained for the negative-order norm. Numerical experiments are given which show that using more compact kernels are less oscillatory and that the superconvergence property holds true for nonlinear conservation laws with general flux functions, indicating that the restriction on f(u) is artificial. Based on our numerical results we can see that these estimates are not sharp. However, they indicate that a more compact kernel can be used in obtaining superconvergence results.

Future work includes the study of accuracy-enhancement of the DG method for one-dimensional nonlinear symmetric/symmetrizable systems and scalar nonlinear conservation laws in multi-dimensional cases on structured as well as unstructured meshes. Analysis of the superconvergence property of the local DG (LDG) method for nonlinear diffusion equations is also on-going work.