OPTIMAL ERROR ESTIMATES TO SMOOTH SOLUTIONS OF THE CENTRAL DISCONTINUOUS GALERKIN METHODS FOR NONLINEAR SCALAR CONSERVATION LAWS ⋆

. In this paper, we study the error estimates to sufficiently smooth solutions of the nonlinear scalar conservation laws for the semi-discrete central discontinuous Galerkin (DG) finite element methods on uniform Cartesian meshes. A general approach with an explicitly checkable condition is established for the proof of optimal 𝐿 2 error estimates of the semi-discrete CDG schemes, and this condition is checked to be valid in one and two dimensions for polynomials of degree up to 𝑘 = 8. Numerical experiments are given to verify the theoretical results


Introduction
In this paper, we study the central discontinuous Galerkin (DG) finite element method for solving scalar conservation laws [10].The optimal error estimates of the central DG methods have been proved for linear conservation laws in Liu et al. [12].In this paper, we present the optimal error estimates of central DG approximation based on tensor-product polynomials under suitable assumptions for the general nonlinear scalar conservation laws ⎧ ⎪ ⎨ ⎪ ⎩ (  ())  = 0, (x, ) ∈ Ω × (0,  ], where x = ( 1 ,  2 , . ..,   ) and Ω is a bounded rectangular domain in R  .Here  0 (x) is a given smooth function.
We do not pay attention to boundary conditions in this paper; hence the exact solution is considered to be either periodic or compactly supported.We also assume the flux  () is smooth in the variable ; for example,  ∈  2 is enough for our proof.The analysis in this paper is for the smooth solutions of (1.1).Discontinuous solutions with shocks are not considered here.We study the cases with  = 1 and 2, but the approach is applicable to any .
The central scheme of Nessyahu and Tadmor [14] computes hyperbolic conservation laws on a staggered mesh and avoids the Riemann solver.In [3], Kurganov and Tadmor introduced a new type of central scheme without the large dissipation error related to the small time step size by using a variable control volume whose size depends on the time step size.To avoid the excessive numerical dissipation for small time steps, Liu [8] uses another coupling technique.The overlapping cell approach evolves two independent cell averages on overlapping cells, which opens up many new possibilities.The advantages of overlapping cells motivate the combination of the central scheme and the DG method, which results in the central DG methods [7,9,10].The central DG method evolves two copies of approximating solutions defined on staggered meshes and avoids using numerical fluxes which can be complicated and costly [4].Like some previous central schemes, the central DG method also avoids the excessive numerical dissipation for small time steps by a suitable choice of the numerical dissipation term.Besides, the central method carries many features of standard DG methods, such as compact stencil, easy parallel implementation, etc.The central DG method with Runge-Kutta time stepping has a larger CFL number for stability than the standard Runge-Kutta DG method with the same polynomial order .Also the central DG method has a smaller error than the standard DG method on the same mesh.See Liu et al. [10], Reyna and Li [15] for more details.Later in Liu et al. [11], the central local discontinuous Galerkin method was introduced to solve diffusion equations, which is formulated based on the local discontinuous Galerkin scheme on overlapping cells.Recently, the central DG method has been used to solve systems of conservation laws in many applications [5,6,16,18,22].
In Liu et al. [12], suitable special projections for central DG methods were proposed to yield optimal error estimates for scalar linear conservation laws.The proper local projections were constructed according to the superconvergence property and the duality of overlapping cells, which also required uniform Cartesian meshes.Zhang and Shu firstly presented a priori error estimates for the fully discrete second order Runge-Kutta DG methods with smooth solutions for scalar nonlinear conservation laws [19] and symmetrizable systems [20].The main techniques they used are Taylor expansion and energy estimates.Later these techniques are widely used in error estimates for DG-type methods of nonlinear equations, like the local DG methods for convection-diffusion and KdV equations [17], the ultra weak DG methods for equations with higher order derivatives [1], the third order Runge-Kutta DG methods for scalar conservation laws [21] and for symmetrizable systems [13].
In this paper, we combine the special projections in Liu et al. [12] and the techniques used in Zhang and Shu [19] to construct new projections to provide the optimal error estimates of the central DG methods on uniform Cartesian meshes for nonlinear scalar conservation laws with smooth solutions.In one dimension, we construct a proper local projection P * ℎ similar to Liu et al. [12].The existence and optimal approximation properties of this projection are proved by standard finite element techniques.Moreover, this projection has similar superconvergence property as the projections in Liu et al. [12].By using this property we develop a general approach with an explicitly checkable condition, and this condition is checked to be valid in one dimension for polynomials of degree up to  = 8.This condition could also be checked for larger  with increased algebraic complexity, but it is not carried out in this paper.The optimal convergence results is valid for uniform meshes and for polynomials of degree  ≥ 1, while for  = 0 we need the convection flux to be linear to get the optimal results.For two-dimensional conservation laws, we follow the same arguments as in the one-dimensional case to construct a suitable projection P * ℎ and to analyze its existence and approximation properties.This new projection utilizes   , the space of tensor-product polynomials of degree at most  in each variable.Similarly, the optimal convergence result is valid for uniform meshes and for polynomials of degree  ≥ 2 in the twodimensional case, while for  = 0, 1 we need the convection flux to be linear to get the optimal results.The superconvergence result of P * ℎ on uniform Cartesian meshes will help to yield optimal convergence results under some suitable assumptions.Similar approach with an explicitly checkable condition is established, and here we also check this condition for polynomials of degree up to  = 8.Likewise, this condition could also be checked for larger  with increased algebraic complexity, but we will not carried it out.The approach is applicable to higher dimension , but it will not be discussed in this paper.
The rest of the paper is organized as follows.In Section 2, we recall the central DG method for one-dimensional conservation laws.Then we construct a special projection and study its existence, uniqueness and optimal approximation properties.With the help of this projection, we will prove the optimal error estimate for the semi-discrete central DG methods on uniform meshes for the nonlinear conservation laws in one dimension.In Section 3, we extend the analysis to two-dimensions.Optimal error estimates are proved by following the same lines of the one dimensional case.We provide numerical examples to show our theoretical results in Section 4. In Section 5, we give a few concluding remarks and perspectives for future work.Finally, in the appendix we provide proofs for some of the more technical results of the error estimates.

Basic notations
For a given interval  = [, ], we divide it into  cells as follows: We denote and similarly for the dual mesh We let ℎ = max  ℎ + 1 2 and assume the mesh is regular.Define the approximation space as Here   (  ) denotes the set of all polynomials of degree at most  on   .For a function  ℎ ∈   ℎ , we use ( ℎ ) to refer to the value of  ℎ at  + 1 2 from the left cell   and the right cell  +1 , respectively.For . the jump of  ℎ or  ℎ at cell interfaces.We denote by  a positive constant independent of ℎ, which may depend on the solution of the problem and other parameters.For our analysis, we need the uniform boundedness of  ′ and  ′′ .We shall take this as an assumption for simplicity, although such boundedness can be shown a posteriori by the eventual boundedness of the numerical solution through the verification of the a priori assumptions at the end of Sections 2 and 3. Similar to Xu and Shu [17], Zhang and Shu [19], to emphasize the nonlinearity of the flux  (), we use  * to denote a non-negative constant depending on the maximum of | ′′ |.We remark  * = 0 for linear fluxes  () =  with a constant .

The central DG scheme
We propose the following semi-discrete central DG scheme for periodic boundary condition: find  ℎ ∈   ℎ and  ℎ ∈   ℎ , such that for any  ℎ ∈   ℎ and where  max is an upper bound for the time step size due to the CFL restriction, that is,  max =  ℎ with a given constant CFL number  dictated by stability.For the initial condition, we simply take  ℎ (•, 0 , where P ℎ and Q ℎ are the  2 projections into   ℎ and   ℎ , respectively, and we have (2.7)

𝐿 2 stability for the linear equation
In Liu et al. [10], the following stability result is proved for this scheme if  () is linear.Without loss of generality, we take  () = .Hence, we have with periodic boundary condition.

Optimal 𝐿 2 error estimate
It is worth noting that the  2 stability for CDG scheme for nonlinear problem is generally not available [10].But under the assumption of the smoothness of the exact solution, we can still get the error estimate of the nonlinear case.In this subsection, we show a priori  2 error estimate of the scheme (2.6) for the equation (2.1).
Here and below, we use ‖ • ‖ to denote the standard  2 norm.For the proof, we recall the classical inverse and trace inequalities [2].For any  ℎ ∈   ℎ or  ℎ ∈   ℎ , there exists a positive constant  independent of  ℎ and ℎ, such that where Γ is the set of boundary points of all elements   or  + 1 2 .First we introduce some notations.For the numerical solutions  ℎ and  ℎ of the CDG scheme (2.6) for equation (2.1), we define and Obviously, we have (2.14) It is also clear that the exact solution  of (2.1) satisfies (2.15) Subtracting (2.14) from (2.15), we obtain the error equation (2.16) Summing over all , the error equation becomes (2.17) Next, we will discuss the properties of the projections P * ℎ and Q * ℎ .Without loss of generality we will only consider P * ℎ .The equation (2.18a) is required by conservation.Note that Pℎ (;  ℎ ; , )  = 0 for ∀ when  ℎ is a constant, so (2.18b) alone misses one condition which is provided by (2.18a).The following lemma gives the existence and uniqueness of the special projection P * ℎ .
Lemma 2.2.The projection P * ℎ defined by (2.18) exists and is unique for any smooth function (), and the following inequality holds for all .The positive constant  depends on , the bound of  ′ (), the constant  in the scheme (2.6) and is independent of ℎ and .
Proof.The proof of this lemma is given in Appendix A.1.
Since P * ℎ and Q * ℎ are -th degree polynomial preserving local projections, standard approximation theory [2] implies, for smooth function , where  is the polynomial degree in the finite element spaces   ℎ and   ℎ , and the constant  depends on , the final time  , ‖‖  +2 and the bounds on the derivatives |  |,  = 1, 2, but is independent of the mesh size ℎ.
Proof.Let   =  −  ℎ ,   =  −  ℎ be the error between the numerical and exact solutions.To deal with the nonlinearity of  (), we would like to first make the a priori assumption that, for small enough ℎ, we have which also establishes the Lipschitz continuity of the right-hand side of the method of lines semi-discrete ordinary differential equation system, hence the very existence of  ℎ and  ℎ .By the interpolation property, we then have This assumption is not necessary for linear  .We will verify this assumption for  ≥ 1 later.By taking (2.30) For the left-hand side of (2.30), we follow the  2 stability proof in Theorem 2.1 for linear case to conclude Similar to Zhang and Shu [19] and Xu and Shu [17], to deal with the nonlinear part of (2.30) we would like to use the following Taylor expansions: By Young's inequality and (2.23), we have (2.34) Next we estimate the nonlinear part.First for the  1 term, we can rewrite it in the form .
By the inequality in (2.10), (2.23) and ‖ ′ ((  )) −  ′ ()‖  ∞ ( ) = (ℎ), we have For B (  ,   ;  ℎ ; , ), let ̂︁   be the Taylor polynomial of order  + 1 of  near Recalling the Bramble-Hilbert lemma [2], we have (2.36) Then we rewrite   and (2.37) Hence, using Proposition 2.3, we have Therefore, by using Young's inequality, (2.23), the inequality in (2.10) and (2.36), we have (2.39) Hence, for  1 we have Similarly, for  2 we have (2.41) The  3 term can be rewritten as the following form 4 is the high order term in Taylor expansion, it is easy to show that ( (2.44) When  ≥ 1, by using a priori assumption (2.26) we have Finally, by Gronwall's inequality and the fact that This, together with the approximation result (2.23), implies the desired error estimate.
For the case of  = 0, we assume that the convection term is linear, namely  () = .This is to avoid the need of the a priori assumption (2.26) which is no longer justifiable since our  2 error estimate is only of order (ℎ) in this case.The proof is similar to that for  ≥ 1 case given above, and the only difference is  * = 0 in this case.By similar lines of proof, we have (2.47) An application of Gronwall's inequality give us that This, together with the approximation result (2.23), implies the desired error estimate.Finally, let us justify the a priori assumption (2.26) for  ≥ 1.Similar to Zhang and Shu [19] and Cheng and Shu [1], we can verify this by a proof by contradiction.By (2.25), we can consider ℎ small enough so that ℎ +1 < 1 2 ℎ 3 2 , where  is the constant in (2.25) determined by the final time  .Define . This is a contradiction if  * <  .Hence,  * ≥  and our a priori assumption is justified.

The central DG method in multi-dimensions
In this section, we consider the semi-discrete central DG method for multidimensional nonlinear conservation laws.Without loss of generality, we will show our central DG scheme and prove the optimal a priori error estimates in two dimensions ( = 2); all the arguments we present in our analysis depend on the tensor product structure of the mesh and finite element space and can be easily extended to the more general cases  > 2. Now we consider the following two-dimensional problem, with periodic boundary condition or compactly supported boundary condition.

Basic notations
]︁}︁ be a partition of Ω into uniform square cells, depicted by the solid lines in Figure 1, and tagged by their cell centroid at (  ,   ).Define , where   ( , ) is the tensor-product polynomials of degrees at most  in each variable defined on  , and no continuity is assumed across cell boundaries.Let  + 1 2 ,+ 1 2 be the dual mesh which consists of a ℎ 2 shift of the  , , depicted by the dashed lines in Figure 1.Let be the cell centroid of the cell  +  1  2 and [ ℎ ] , .

The central DG scheme
We propose the following semi-discrete CDG scheme for periodic boundary condition: find  ℎ ∈   ℎ and  ℎ ∈   ℎ , such that for any ,+ 1 2 where  max is a max step size, determined by  max = (CFL factor) × ℎ/(maximum characteristic speed), in which the CFL constant should be less than 1/2.Similarly, for the initial condition we simply take , where P ℎ and Q ℎ are the  2 projections into   ℎ and   ℎ , respectively, and we have )︂ . (3.3)

𝐿 2 Stability for linear equation
The  2 -stability is proved for the CDG scheme (3.2) in Liu et al. [10] if  () and () are linear.Without loss of generality, we take  () = () = .Hence, we have with periodic boundary condition.
Theorem 3.1.The numerical solutions  ℎ and  ℎ of the semi-discrete CDG scheme (3.2) for the equation (3.4) have the following  2 stability property

Optimal 𝐿 2 error estimate
In this subsection, we show the a priori  2 error estimate of the scheme (3.2) for the equation (3.1).
Here and below, we again use ‖ • ‖ to denote the standard  2 norm.Similar to the one-dimensional case, we recall the classical inverse and trace inequalities [2].For any  ℎ ∈   ℎ or  ℎ ∈   ℎ , there exists a positive constant  independent of  ℎ and ℎ, such that where Γ is the set of boundaries of all elements  , or  + 1 2 ,+ 1 2 .
Let  be the exact solution of equation (3.1), clearly we have 2 )︁ .

Projection operators
To prove the error estimates for two-dimensional problems in uniform Cartesian meshes, we need two suitable projections P * ℎ and Q * ℎ similar to the one-dimensional case.By applying the shifting technique in the twodimensional case, for  and  variables respectively, for a given function () we define P * ℎ  ∈   ( Similarly, we can define the projection 2 2 2 .Next we will discuss the properties of these two special projections.Without loss of generality we will only consider P * ℎ .The equation (3.13a) is required by conservation.Note that Pℎ (;  ℎ ) , = 0 for ∀ when  ℎ is a constant, so (3.13b) alone misses one condition which is provided by (3.13a), just like the one-dimensional case.Existence and optimal approximate property of the projection P * ℎ are established in the following lemma.
Lemma 3.2.The projection P * ℎ defined by (3.13) exists and is unique for any smooth function (), and the following inequality holds for all .The positive constant  depends on , the bound of  ′ (),  ′ (), the constant  and is independent of ℎ and .
Proof.The proof of this lemma is given in Appendix A.3.
Similarly, for Q * ℎ we have if  is a smooth function.

A priori 𝐿 2 error estimates
Now let us give the a priori error estimate for the two-dimensional case.
Proof.Let   =  −  ℎ ,   =  −  ℎ be the error between the numerical and exact solutions.Similar to the onedimensional case, to deal with the nonlinearity of  () and (), we would like first make a priori assumption that, for small enough ℎ, we have which also establishes the Lipschitz continuity of the right-hand side of the method of lines semi-discrete ordinary differential equation system, hence the very existence of  ℎ and  ℎ .By the interpolation property, we then have This assumption is not necessary for linear  and .We will verify this assumption for  ≥ 2 later.By taking For the left-hand side of (3.24), we follow the  2 stability proof in Theorem 3.1 for linear case to conclude Similar to the proof in Zhang and Shu [19] and Xu and Shu [17], to deal with the nonlinear part of (3.24) we would like to use the following Taylor expansions: where (  )   ℎ d d, ,+ 1 2 By Young's inequality and (3.15), (3.16) we have (3.28) Next we estimate the nonlinear part.First for the  1 term, we can rewrite it as  where  Hence, for  1 we have Similarly, for  2 we have Similar to the one-dimensional case, the  3 term can be rewritten as 2 4 is the high order term in Taylor expansion, it's easy to show that (3.41) When  ≥ 2, by using a priori assumption (3.21) we have Finally, by Gronwall's inequality and the fact that This, together with the approximation result (3.15), (3.16) implies the desired error estimate.
For the case of  = 0 or 1, we assume that  () and () are linear fluxes, namely  () =  1 , () =  2  with constants  1 ,  2 .This is to avoid the need of the a priori assumption (3.20) which is no longer justifiable in this case.By similar lines of proof and noting that  * = 0 in this case, we can obtain By using the Gronwall's inequality we have This, together with the approximation result (3.15), (3.16), implies the desired error estimate for  = 0, 1 with linear fluxes.Just like the one-dimensional case, let us justify the a priori assumption (3.20) with  ≥ 2. Similar to Zhang and Shu [19] and Cheng and Shu [1], we can verify this by a proof by contradiction.By (3.19), we can consider ℎ small enough so that ℎ +1 < 1 2 ℎ 2 , where  is the constant in (3.19) determined by the final time  .
. This is a contradiction if  * <  .Hence,  * ≥  and our a priori assumption is justified.

Numerical examples
In this section, we present numerical examples to verify our theoretical findings.Uniform meshes are used in all examples.The schemes are integrated in time with the third order SSP Runge-Kutta method.We would like to compute on elements of degree  = 0, 1, 2, 3. We set the CFL number to be 0.05.For  = 0, 1, 2 we let ∆ = CFL • ℎ and ∆ = CFL • ℎ 4 3 for  = 3 where ℎ is the characteristic length of the mesh, so that the time error will be dominated by the spatial error.
Example 4.1.We solve the one-dimensional Burgers equation given by The exact solution is obtained by Newton iteration.In this example, we use  max = ℎ 2+1 , ℎ = 2  to test the numerical schemes.The errors and numerical order of accuracy at  = 0.5 with 0 ≤  ≤ 3 are listed in Tables 1.
Table 1 shows that the order of convergence of the error achieves the expected ( + 1)-th order of accuracy.
Example 4.2.We solve the two-dimensional Burgers equation given by with periodic boundary condition.The exact solution follows from the solution of one-dimensional Burgers equation with  =  + .In this example, we use  max = ℎ 2+1 , ℎ = 2  to test the numerical schemes.The central DG scheme is evolved up to  = 0.2 when the solution is still smooth.The errors and numerical order of accuracy with 0 ≤  ≤ 3 are listed in Table 2.
Table 2 shows that the order of convergence of the error achieves the expected ( + 1)-th order of accuracy.

Concluding remarks
In this paper, a priori optimal  2 error estimates to central DG methods on uniform meshes applied to nonlinear conservation laws with smooth solutions are proved with polynomial degrees of  ≤ 8.The main techniques used in this paper are special projections and Taylor expansions.Our analysis is carried out both in one dimension and in two-dimensions for uniform Cartesian meshes and tensor-product polynomial spaces.We also give some numerical examples to verify the results of our theoretical analysis.The error estimates for nonlinear conservation laws in this paper were obtained using stability for the linear case and the smoothness of the exact solution.It is not clear whether stability holds for the scalar nonlinear conservation laws with

Appendix A. Collection of technical proofs
In this appendix, we collect the proofs of some technical lemmas and propositions.
A.1.Proof of Lemma 2.2 Proof.We only consider P * ℎ , while the proof for Q * ℎ follows similar lines.For ∀, we let  = 2(− ) ℎ on   , for a smooth function () and a -th order polynomial  ℎ () on   , and define Note that the procedure to find the P * ℎ ω ∈ P  ([−1, 1]) is to solve for a linear system, so existence of the projection can be proved by proving its uniqueness.Thus, we only need to prove the uniqueness of the projection P * ℎ .We set   () = P * ℎ ω() = P * ℎ () with ω() = () = 0, and would like to prove   () = 0. Then by the definition of the projection P * ℎ , we have: We rewrite Pℎ (  ;   ; , )  by a change of variable  →  + 1 for the integrations over [−1, 0] to get For  = 0 it clearly holds.For  ≥ 1, now from (A.5) we have Assume   , 1 ≤  ≤  are not all zeros, then () is a non-zero polynomial of degree at most  − 1, thus it has at most  − 1 roots, which is a contradiction to (A.6).Hence, we have . Hence, by (A.2b), we have We have now finished the proof of uniqueness.Next we move to prove the boundedness.Let   () = P * ℎ () = ∑︀  =0     and set the test functions  ℎ = ,  2 , . ..,   .Then we have By calculation, for 1 ≤  ≤  we have where From the formulation of the scheme (2.6) we have  max =  ℎ, here  is a constant dictated by stability.Then we have From (A.11) we know that and from (A.9) we have Hence, we can estimate the infinity norm of   , Since   > 0 for ( + ) is even and  ′ ((  )) is bounded, then we have where ℰ is a constant which depends on polynomial degree , the bound of  ′ ((  )) and constant .Since the first row of the matrix   are constants  0 which only depends on degree  and the other elements of   either only contain )︀  .From the previous proof of the existence and uniqueness of the projection, we know that   is always invertible which means det (︀   )︀ ̸ = 0 holds for any value of  ′ ((  )).Hence, here we have   () ̸ = 0. Therefore, we can take  small enough so that )︂  ( ′ ((  ))) We emphasize that this choice of  is only a sufficient condition for our proof, in numerical computation  should be chosen as the largest CFL number for linear stability to avoid excessive numerical dissipation.We now have (A.26) By the equivalence of norms we have It is obvious that ‖‖ ∞ ≤ C‖‖ ∞ due to the boundedness of  ′ ((  )).Here C is a constant which depends on degree  and the bound of  ′ ((  )).Hence, for the coefficients  we have which immediately results in the boundedness of P * ℎ .
For  = 0, 1. .., 8, by using the definition of the projection and the property that = (ℎ) we have the following results.For  =  +1 , by the definition (for  = 0 we only have the first equation in the definition), (A.34) Here   and   are the coefficients obtained by solving the local linear system (A.33).We leave the detailed calculations and formulas for  up to 8 in a separate file, as a supplement to this paper, since they are too lengthy.We then have,  = 0, 1, . . ., 8, that and therefore we can prove that where we have again used change of variable to shift all the integration regions to the same subcell Thus   (, ) ≡  0 on  , ,  0 is a constant.By (3.13) we immediately get   ≡ 0, and we have finished the proof of uniqueness, hence also existence.We note that this projection is a local projection, hence we can make a change of variables to the reference element . Taking a similar derivation as in the proof (Sect.A.1), we obtain (A.40) Again standard approximation theory [2] implies the optimal approximating estimates.
A.4. Proof of Lemma 3.3 , then by the definition of B, and B+  = (ℎ) we have the following results:

Figure 1 .
Figure 1.2D overlapping cells formed by collapsing the staggered dual cells on two adjacent time levels to one time level.

Table 2 .
Errors and numerical orders of accuracy for Example 4.2 on a uniform mesh of  × cells.Here  max = ℎ 2+1 and final time  = 0.2.general non-smooth solutions.Such a stability proof for the central DG schemes and the extension of this work to non-uniform meshes and unstructured triangular meshes are interesting and challenging, and constitutes our ongoing work.