Superconvergence of the semi-discrete local discontinuous Galerkin method for nonlinear KdV-type problems

In this paper, we present and analyze a superconvergent local discontinuous Galerkin (LDG) scheme for the numerical solution of nonlinear KdV-type partial differential equations. Optimal a priori error estimates for the LDG solution and for the two auxiliary variables that approximate the first-and second-order derivative are derived in the L2-norm for the semi-discrete formulation. The order of convergence is proved to be p+1, when piecewise polynomials of degree at most p are used. We further prove that the derivative of the LDG solution is superconvergent with order p+1 towards the derivative of a special projection of the exact solution. We use this results to prove that the LDG solution is superconvergent with order p+3/2 toward a special Gauss-Radau projection of the exact solution. Finally, several numerical examples are given to validate the theoretical results. Our proofs are valid for arbitrary regular meshes using Pp polynomials with p≥1 and under the condition that |f'(u)| possesses a uniform positive lower bound, where f(u) is the nonlinear flux function. Our experiments demonstrate that our results hold true for KdV equations with general flux functions.

1. Introduction. In this paper, we present and analyze a superconvergent local discontinuous Galerkin (LDG) method for the following Korteweg-de Vries (KdV) equation with periodic boundary conditions and initial condition where f (u), g(x, t), and u 0 (x) are arbitrary smooth and given functions. For the sake of simplicity, we only consider the case of periodic boundary conditions. This assumption is not essential and the method can be easily designed for non-periodic boundary conditions (e.g., Dirichlet or Neumann or mixed boundary conditions); see [10,7,6,18,43] for some discussion. In our analysis, the exact solution is assumed to be smooth function on [a, b] × [0, T ] for a fixed T . The nonlinear KdV equation was derived by Diederik Korteweg and Gustav de Vries in 1895. It describes the propagation of waves in a variety of nonlinear dispersive media. This equation is a generic equation for the study of weakly nonlinear long waves and arises in many physical situations, such as surface water waves and plasma waves; see e.g., [11]. It has been shown that the KdV equation describes a large class of solitons observed in various situations including acoustic waves on a crystal lattice, plasma waves, hydrodynamics internal or surface waves, elastic surface waves, and waves in optical fibers [25].
The LDG method is an extension of the discontinuous Galerkin (DG) method aimed at solving PDEs containing higher than first order spatial derivatives. The LDG method was first introduced by Cockburn and Shu in [24] for solving convection-diffusion problems. The main idea of the LDG method is to rewrite a higher order equation into a system of first order PDE and then discretize it by the standard DG method. The LDG method has several advantages over the classical numerical methods available in the literature. One of the advantages of the LDG method over the existing numerical schemes in the literature is that the LDG scheme achieves optimal (p + 1)th order convergence for the solution and its spatial derivative in the L 2 -norm. Furthermore, as we shall see below, it achieves superconvergence results, which can be used to construct asymptotically exact a posteriori error estimates by solving a local residual problem on each element. Moreover, LDG methods are extremely flexible in the mesh-design; they can easily handle meshes with hanging nodes, elements of various types and shapes, and local spaces of different orders. LDG schemes have been successfully applied to hyperbolic, elliptic, and parabolic PDEs [2,3,4,5,23,32,26,27,18,42,30,24,14,17,6,35,16], to mention a few. Furthermore, LDG methods can be easily designed for many other higher order nonlinear wave and diffusion equations such as the KdV equations, the Kadomtsev-Petviashvili equations, the Zakharov-Kuznetsov equations, the Kuramoto-Sivashinsky-type equations, the Cahn-Hilliard equation, and the equations for surface diffusion and Willmore flow of graphs. A review of the LDG methods is given in [6,8,16,22,20,15,21,18,37,43,13,36].
The LDG method for the third-order KdV equations has been studied in [39,41,38,40,28,9]. In [41], the authors developed an LDG method for solving KdVtype equations containing third derivative terms in one and two space dimensions. They proved the L 2 stability and a cell entropy inequality for the square entropy for a class of nonlinear KdV equations in both one and multiple space dimensions. They further obtained, for the linear KdV equation, the sub-optimal error estimate of O(h p+1/2 ) in the L 2 norm. Numerical examples in [41] suggest that the LDG solution converges at O(h p+1 ). In [38], Xu and Shu provided L 2 error estimates for the semi-discrete LDG method for nonlinear KdV equations with smooth solutions. The proved sub-optimal error estimate of O(h p+1/2 ) in the L 2 norm. However, they could not obtain the optimal error estimates O(h p+1 ) due to some extra boundary terms arising from high order derivatives. Later, Xu and Shu [40] proved optimal L 2 error estimates of the semi-discrete LDG methods for solving linear high order wave equations including the linearized KdV equation. More recently, Hufford and Xing [28] investigated the superconvergence property of the LDG method for solving the linearized KdV equation. They selected a special projection of the initial condition and proved that the LDG solution is O(h p+3/2 ) super close to a particular projection of the exact solution, when the upwind flux is used for the convection term and the alternating flux is used for the dispersive term. In [9], we presented and analyzed a posteriori error estimates for the LDG method for the linearized KdV equation in one space dimension. These error estimates are computationally simple and are obtained by solving a local steady problem with no boundary condition on each element. We extended the work of Hufford and Xing [28] to prove new superconvergence results towards particular projections of the exact solutions for the two auxiliary variables in the LDG method that approximate the first and second derivatives of the solution. The superconvergence rate is proved to be of order p + 3/2. These results allow us to prove that the significant parts of the spatial discretization errors for the LDG solution and its spatial derivatives (up to second order) are proportional to (p + 1)-degree Radau polynomials. We used these results to construct asymptotically exact a posteriori error estimates. Finally, we proved that, for smooth solutions, our a posteriori LDG error estimates for the solution and its spatial derivatives at a fixed time t converge to the true errors at O(h p+3/2 ) rate. The purpose of this paper is to extend the analysis to the nonlinear KdV equation. To the best knowledge of the author, this is the first study investigating the superconvergence properties of the LDG method for the nonlinear KdV equation.
In this paper, we study the superconvergence properties for the LDG method for the nonlinear KdV equation (1), extending the convergence and superconvergence results in [28,9] for the linearized KdV equation. We prove that the LDG solutions are (p + 1)th order convergent in the L 2 -norm, when the space of piecewise polynomials of degree p ≥ 1 is used. Computational results indicate that the theoretical order of convergence is sharp. Moreover, we show that the derivative of the LDG solution is superconvergent with order p + 1 towards the derivative of a Gauss-Radau projection of the exact solution. Finally, we prove that the LDG solution is superconvergent with order p + 3/2 towards a Gauss-Radau projection of the exact solution, while computational results show higher O(h p+2 ) convergence rate. In our analysis we proved these convergence results under mesh refinement and at a fixed time t and time discretization is assumed to be exact. Our proofs are valid for any regular meshes and using piecewise polynomials of degree p ≥ 1.
An outline of the paper is given as follows: In section 2 we present the semidiscrete LDG method for solving the nonlinear KdV equation. We also introduce some notation and definitions. In section 3 we present the LDG error analysis and prove several a priori error estimates in the L 2 -norm. In section 4 we state and prove the main superconvergence results. In section 5 we present numerical results to validate the theoretical results. We conclude and discuss our results in section 6.
2. The LDG method for the KdV equations. To define the LDG method, we introduce two auxiliary variables q = u x , r = q x , and convert equation (1a) into the following first-order system We divide the computational domain Ω = [a, b] into N intervals h i to be the length of the largest interval. In our analysis, we assume that the mesh is regular in the sense that the ratio between the maximum and the minimum mesh sizes stays bounded during mesh refinements. For simplicity, we use v i to denote the value of the continuous function v = v(x, t) at x = x i . We also use v − i and v + i to denote the left limit and the right limit of v at the discontinuity point Multiplying the three equations in (2) by test functions v, w, and z, respectively, integrating over an arbitrary element I i , and using integration by parts, we obtain We introduce the following discontinuous finite element space . . , N }, where P p (I i ) denotes the space of polynomials of degree at most p on I i with coefficients as functions of t. We note that polynomials in the space V p h are allowed to be completely discontinuous at the mesh points.
Next, we approximate the exact solutions u, q, and r, at any fixed time t, by piecewise polynomials of degree at most p. For simplicity, we denote the approximations of u, q, r by u h , q h , r h ∈ V p h , respectively. The LDG scheme can now be defined as: where the termsf ,û h , andq h are the so-called numerical fluxes, which are, respectively, the discrete approximations to f (u), u, and q at the nodes. The numerical fluxf is a single-valued function defined at the nodes and in general depends on the values of u h from both sides i.e.,f =f (u − h , u + h ). We would like to mention that the numerical fluxes have to be suitably chosen in order to ensure the stability of the method and also to improve the order of convergence. In this paper,f is a monotone numerical flux for f (u), i.e., it satisfies the following three conditions: is a nondecreasing function of its first argument u − h and a nonincreasing function of its second argument u + h . The popular monotone numerical fluxes are the Godunov flux, the Engquist-Osher flux, the Lax-Friedrichs flux, etc (see [31]). It turns out that we can take the following simple choices: • The numerical fluxf associated with the convection is taken as the upwind flux which depends on the sign of f ′ (u) i.e., (4d) • The numerical fluxesû h ,q h , andr h can be taken as (e.g., see [39,41,38]) We remark that the choice for the fluxes (4e) is not unique. In fact, the crucial part is to takeû h andr h from opposite sides. Thus, the same results can be proved using the following numerical fluxes with only minor modificationŝ Remark 1. In our analysis, we will assume that f ′ (u) ≥ 0 so thatf i = f (u − h ) i . We note that the other case f ′ (u) ≤ 0 can be handled in a very similar way and hence we omit the details. Our numerical experiments (see section 5) demonstrate that the superconvergence property holds true for nonlinear problem with general flux functions, indicating that the assumption f ′ (u) ≥ 0 is artificial.
To complete the definition of the LDG scheme, we still need to define the discrete initial condition u h (x, 0) ∈ V p h . In this paper, we use a special projection of the exact initial condition u 0 (x). This particular projection will be defined later and is needed to achieve global superconvergence result towards the Gauss-Radau projection, which will be defined later.
Implementation: We note that the approximate solutions u h , q h and r h can be efficiently computed in an element-by-element fashion. Specifically, we can obtain the LDG solutions in the following ordering: first we find u h (x, 0) ∈ V p h using a special projection of the exact initial condition u 0 (x). Then, we use (4c) to obtain q h (x, 0) ∈ V p h since u h (x, 0) is given. Once we have q h (x, 0), we use (4b) to figure out r h (x, 0). Finally, we discretize the system of ODEs (4a) using an ODE solver to obtain u h at the next time step ∆t. This process can be repeated to obtain u h , q h and r h for t ∈ [0, T ]. L 2 stability: We have the following L 2 stability for the LDG scheme (4).
Corollary 1 (L 2 stability). The solution to the LDG scheme (4) satisfies the L 2 stability 1 2 Proof. This proof can be found in [41]. More precisely, in its Corollary 2.2 page 776.
Norms: The standard L 2 -norm of u = u(x, t) over I i is denoted by u(·, t) 0,Ii = Ii u 2 (x, t)dx 1/2 . Moreover, the L ∞ -norm of u(·, t) on I i at time t is defined by u(·, t) ∞,Ii = sup . The H s (I i )-seminorm of u on I i is given by |u(·, t)| s,Ii = ∂ s x u(·, t) 0,Ii . We also define the norms on the whole computational domain Ω as follows: The seminorm on the whole computational domain Ω is defined as |u(·, t)| s,Ω = N i=1 |u| 2 s,Ii . We note that if u(·, t) ∈ H s (Ω), the norm u(·, t) s,Ω on the . For simplicity, if we consider the norm on the whole computational domain Ω, then the corresponding index will be omitted. Thus, we use u , u s , and u ∞ to denote u 0,Ω , u s,Ω , and u ∞,Ω , respectively. We also use u(t) to denote the value of u(·, t) at time t. In particular, we use u(0) to denote u(·, 0) . Finally, we omit the argument t and we use u to denote u(t) whenever confusion is unlikely.
Projections: For p ≥ 1, we define P ± h u as two special Gauss-Radau projections of u onto V p h as follows [18]: The restrictions of P + h u and P − h u to I i are polynomials in P p (I i ) satisfying These special projections are used in the error estimates of the DG methods to derive optimal L 2 error bounds in the literature, e.g., in [18]. They are mainly used to eliminate the jump terms at the element boundaries in the error estimates in order to prove the optimal L 2 error estimates. In our analysis, we need the following projection results [19]: If u ∈ H p+1 (I i ), then there exists a positive constant C independent of the mesh size h, such that In the rest of the paper, we will not differentiate between various constants, and instead will use a generic constant C (or accompanied by lower indices) to represent a positive constant independent of the mesh size h, but which may depend upon the exact smooth solution of the PDE (1a) and its derivatives.
Properties of the finite element space: We recall some inverse properties of the space V p h that will be used in our error analysis [30]: For any v ∈ V p h , there exists a positive constant C independent of v and h such that (8) Finally, the Sobolev's inequality implies that there exists a positive constant C is independent of h such that The initial condition of the LDG method: In 2014, Hufford and Xing [28] studied the superconvergence property for the LDG method for solving the linearized KdV equation. They selected a special projection of the initial condition u h (x, 0) = P h u(x, 0) and proved that the LDG solution is O(h p+3/2 ) super close to P − h u. The projection P h u ∈ V p h is designed to better control the error of the initial condition. It is defined as follows: let q = u x and r = q x . Suppose q h , r h ∈ V p h are the unique solutions (with given P h u) to then we require Proof for the existence and uniqueness of P h u is provided in [28].
Lemma 2.1. The operator P h exists and is unique. Moreover, there exists a constant C independent of h such that Proof. The proof of this lemma is given in [28], more precisely in its Lemma 2.1.
In our mathematical error analysis and numerical examples we will approximate the initial conditions of our numerical scheme on each interval as As discussed in [28], P h is only needed for technical purposes in the proof of the error estimates. In our numerical experiments we used the special projection P h and the projection P − h and observed similar results.

3.
A priori error estimates of the LDG method for the KdV equations.
In this section, we will derive optimal L 2 error estimates for the LDG method. In our analysis, we always assume that the exact solution u is smooth enough. We also assume that the flux function f in (1) is smooth enough. More precisely, we always assume that f (u) satisfies the following conditions Assumption 1. f and its derivatives with respect to u are smooth functions and u and its derivatives are bounded. Thus, f (u) and its derivatives are bounded functions.
To be more precise, we always assume that there exists a constant C such that f (k) (u) ≤ C for k = 0, 1, 2, 3. Throughout this paper, e u , e q and e r , respectively, denote the errors between the exact solutions of (2) and the LDG solutions defined in (4), i.e., e u = u − u h , e q = q − q h and e r = r − r h . Let the projection errors be defined as ǫ u = u − P − h u, ǫ q = q − P + h q and ǫ r = r − P + h r and the errors between the numerical solutions and the projection of the exact solutions be defined asē We note that the true errors can be split as e u = ǫ u +ē u , e q = ǫ q +ē q , e r = ǫ r +ē r .
We also note that, by the definitions of the projections P ± h (6), the following hold: since v is a polynomial of degree at most p and thus v x is a polynomial of degree at most p − 1.
Next, we derive some error equations which will be used repeatedly throughout this and the next sections. Subtracting (4) from (3) with v, w ∈ V p h and using the numerical fluxes (4d) and (4e), we obtain the following error equations Using the classical Taylor's series with integral remainder in the variable u and using the relation Substituting (17) into (16), using (14), and applying (15), we get which, after integrating by parts, are equivalent to To deal with the nonlinearity of the flux f (u), we would like to make an a priori assumption that, for small enough h and p ≥ 1, there holds where C is a constant independent of h. This is obviously satisfied at time t = 0 since, initially by (12), ||ē u (0)|| ≤ Ch p+3/2 ≤ Ch 3/2 for p ≥ 0. We will justify this a priori assumption for piecewise polynomials of degree p ≥ 1 at the end of the proof of Theorem 3.1. In remark 2, we will explain that this a priori assumption is unnecessary for the linear flux f (u) = cu, where c is a constant.
As a consequence of the a priori assumption (20), we have the following results.
Proof. Using the inverse property (8) and the a priori assumption (20), we obtain which completes the proof of (21). Finally, using (14), the triangle inequality, the estimate (9), and the estimate (21), we get, which completes the proof of (22). Now we are ready to derive several optimal a priori error estimates in the L 2norm.
Theorem 3.1. Let p ≥ 1. Let (u, q, r) and (u h , q h , r h ) be the exact and numerical solutions of (2) and (4), where u h (x, 0) is defined in (13). Also, we assume that is the set of real m-times continuously differentiable functions which are bounded together with their derivatives up to the mth order. Then, for sufficiently small h, there exists a positive constant C independent of h such that, ∀ t ∈ [0, T ], Proof. We follow the proof in [28], where the optimal error estimate of LDG method for the third order linear equation is provided. We first use the error equations (18) and (19) to derive four important estimates.

MAHBOUB BACCOUCH
Adding the above equations, we get Summing over all elements, we obtain Next, we will estimate T k , k = 1, . . . , 5 one by one.
Estimate of T 1 . Using the Fundamental Theorem of Calculus and the fact θ i = f ′ (u i ) ≥ 0, we obtain Summing over all elements and using the periodic boundary conditions yields Next, we will estimate θ i+1 − θ i . Using Taylor's Theorem and the assumption |f ′′ | ≤ C 1 on R, we get where we used the smoothness of the exact solution u (i.e., |u x | ≤ C 2 ). Thus, Combining (27) and (28) and applying the inverse property (8), we conclude that Estimate of T 2 . We first rewrite θ − θ i on each element I i as follows Adding and subtracting f ′ (u), we write Using Taylor's Theorem, we can bound the interpolation error Using the smoothness of the exact solution u and f , we get Similarly, we can bound the error Using the estimate (22) and the assumption |f ′′ | ≤ C 2 , we get Combining (30), (31), and (32), we conclude that We can use the same argument the show the following estimates Now, applying the Cauchy-Schwarz inequality, (33), (34), the inverse inequalities in (8), we obtain Estimate of T 3 . Since ǫ u is orthogonal to (ē u ) x ∈ P p−1 (I i ) (due to the properties in (6a)), T 3,i simplifies to Summing over all element, using the estimate (33), the Cauchy-Schwarz inequality, the inverse inequality (8), the projection result (7), and Young's inequality, we get Estimate of T 4 . Applying the Cauchy-Schwarz inequality, the projection result (7), and Young's inequality, we get Estimate of T 5 . Using the periodic boundary conditions, we get Now, combining (26) with (29), (35), (37), (38), and (39), we conclude that Second estimate. Taking w = (ē u ) t in (18b), and then taking the first time derivative of (19c) and letting z =ē q , we obtain Adding these two equations together, we obtain the error equation Summing over all elements and using the periodic boundary conditions, we get Applying the Cauchy-Schwarz inequality, the projection result (7), and Young's inequality, we get Third estimate. Taking the first time derivative of (19a), (18b), (18c), and then choosing the test functions as Adding the above three equations, we get Adding and subtracting Summing over all elements, we obtain where Next, we will estimate A k , k = 1, . . . , 9 one by one.
Consequently, we have Using (43), we rewrite A 1,i as Summing over all elements, using the assumption |f ′′ | ≤ C on R, applying the Cauchy-Schwarz inequality, the projection result (7), and Young's inequality, we get t)) and u, f are assumed to be smooth functions.
Estimate of A 2 . Using the Fundamental Theorem of Calculus and the assumption Summing over all elements and using the periodic boundary conditions yields Combining (45) and (28) and applying the inverse property (8), we conclude that Estimate of A 3 . Adding and subtracting u t f ′′ (u i ), we write (θ−θ i ) t on each element I i as In order to estimate Θ 2 , we add and subtract f ′′ (u) to write Using Taylor's Theorem, we can bound the interpolation error f ′′ (u) − f ′′ (u i ) on each element I i as , and λ k ∈ [0, 1], k = 1, 2. Using the smoothness of the exact solution u and f , we get Similarly, we can bound the error f ′′ (u − se u ) − f ′′ (u) on each element I i as Using the estimate (22) and the assumption |f ′′′ | ≤ C 3 , we get Using the smoothness of the exact solution u, we get Finally, using Taylor's Theorem, we estimate the interpolation error (u − u i ) t on each element I i as Using the smoothness of the exact solution u and f , we get Using the same argument one can show the following estimates Now, following the proof of the estimate (35), we find that Estimate of A 4 . The estimate of A 4 is similar to the estimate of T 2 (35) exceptē u is replaced by (ē u ) t . Thus, Estimate of A 5 . Using the periodic boundary conditions, we get Estimate of A 6 . We first write A 6,i as Summing over all elements and using the periodic boundary conditions, we get Estimate of A 7 . The estimate of A 7 is similar to the estimate of T 3 derived in (37) except thatē u is replaced by (ē u ) t . Hence, following the same line of reasoning used to prove (37), we can show the estimate Estimate of A 8 . The estimate of A 8 is also similar to the estimate of T 3 derived in (37). Thus, we can show Estimate of A 9 . Applying the Cauchy-Schwarz inequality, the projection result (7), and Young's inequality, we get Now, combining (42) with (44), (46), and (49)-(55), we arrive at Fourth estimate. Taking v = −ē q in (19a) and w =ē r in (18b), we get Adding these two equations, we get Adding and subtracting θ i = f ′ (u i ) to θ, we get ē r Summing over all elements, we obtain where B k = N i=1 B k,i , k = 1, . . . , 5. Next, we will estimate these terms separately.
Estimate of B 1 . Taking z =ē q in (18c), we obtain Substituting (59) into (57a) yields Summing over all elements, using the assumption |f ′ | ≤ C on R, the Cauchy-Schwarz inequality, the projection result (7), and Young's inequality yields Estimate of B 2 . Applying the Cauchy-Schwarz inequality, (33), (34), the inverse inequalities in (8), we obtain Estimate of B 3 . Since ǫ u is orthogonal to (ē q ) x ∈ P p−1 (I i ) (due to the properties in (6a)), B 3,i simplifies to Summing over all element, using the estimate (33), the Cauchy-Schwarz inequality, the inverse inequality (8), the projection result (7), and Young's inequality, we obtain Estimate of B 4 . Applying the Cauchy-Schwarz inequality, the projection result (7), and Hölder's inequality, we get Estimate of B 5 . We first write B 5,i as i . Summing over all elements and using the periodic boundary conditions, we get Combining (58) with (60), (61), (62), (63), and (64), we get Next, we will estimate 1 Summing this error equation over all elements and using the periodic boundary conditions, we get Applying the Cauchy-Schwarz inequality, the projection result (7), and Hölder's inequality yields Adding and subtracting Summing the error equation (67) over all elements, we obtain where D k = N i=1 D k,i , k = 1, . . . , 5. Next, we will estimate D k , k = 1, . . . , 5 separately.
Estimate of D 1 . Taking z =ē r in (18c), we obtain Substituting (69) into (67a), we get Summing over all elements, using the assumption |θ i | = |f ′ (u i )| ≤ C on R, the Cauchy-Schwarz inequality, the projection result (7), and Hölder's inequality yields Estimate of D 2 . The estimate of D 2 is similar to the estimate of B 2 (61) exceptē q is replaced byē r . Thus, one can prove the following estimate Estimate of D 3 . The estimate of D 3 is similar to the estimate of B 3 (62). Thus, we can easily show Estimate of D 4 . Applying the Cauchy-Schwarz inequality, the projection result (7), and Hölder's inequality, we get Estimate of D 5 . Using the periodic boundary conditions, we get Combining (68) Now, we are ready to estimate ē r 2 . Combining (65), (66), and (75), we arrive at which yields Integrating over the interval [0, t] and letting Λ(t) = ē u 2 + ē q 2 + (ē u ) t 2 , gives where Using the estimates in (12), we see that Λ(0) = ē u (0) 2 + ē q (0) 2 + (ē u ) t (0) 2 ≤ Ch 2p+2 . Thus, the following holds Using integration by parts with respect to time, we write K(t) as
Using (14), Young's inequality, and the projection result (7), we obtain which completes the proof of the theorem.
Justification of the a priori assumption: Now, to complete the proof of Theorem 3.1, let us justify the a priori assumption (20). We will follow the same arguments used in [30,38,29]. First of all, the a prior assumption is satisfied at t = 0 since initially (due to (12)) On the other hand, our proof implies that (23a) holds for t ≤ t * , in particular ē u (t * ) ≤ Ch p+1 < 1 2 h 3/2 . This is a contradiction if t * < T . Hence t * ≥ T and our a priori assumption (20) is justified.

Remark 2.
We remark that when f (u) is linear, we have f ′ (u) is a constant. Therefore, θ − θ i = 0. Consequently, the a priori assumption (20) is unnecessary. 4. Superconvergence error analysis. In this section, we will prove that the derivative of the LDG solution (u h ) x is O(h p+1 ) superconvergent towards (P − h u) x in the L 2 -norm. This result allows us to prove that the LDG solution is O(h p+3/2 ) superconvergent to P − h u, when p-degree piecewise polynomials with p ≥ 1 are used. In our analysis, we need some properties of Legendre polynomials. We denote bỹ L p the Legendre polynomial of degree p on [−1, 1], which can be defined by the Rodrigues formula [1] The Legendre polynomial satisfies the propertiesL p (1) = 1,L p (−1) = (−1) p , and the orthogonality relation Mapping the physical element I i into the reference element [−1, 1] by the standard affine mapping we obtain the p-degree shifted Legendre polynomial on Using the mapping (82) and the orthogonality relation (81), we obtain Proof. Taking the test function z = (ē u ) x −(−1) p (ē u ) + x i−1 L p,i (x) ∈ P p (I i ) in (19c), using the property L p,i (x i−1 ) = (−1) p , and applying (81), we get since for this choice z + i−1 = 0 and Ii (ē u ) x L p,i dx = 0. Applying the Cauchy-Schwarz inequality, the inverse inequality, and the estimate (83) yields Consequently, we have (ē u ) x 0,Ii ≤ C 4 e q 0,Ii . Squaring both sides, summing over all elements, and using the estimate (23b), we obtain which completes the proof of (84).

4.2.
Superconvergence for the LDG solution towards P − h u. We will prove that the LDG solution is O(h p+3/2 ) superconvergent to P − h u. We first recall the following lemma which will be needed in our analysis.
Proof. The proof of this lemma can be found in [6], more precisely in its Lemma 2.3.
Sinceē u ,ē q ,ē r ∈ V p h are piecewise polynomials, they can be written on each element I i as where a i =ē − u i , b i =ē + q i−1 , c i =ē + r i−1 , and P u (·, t), P q (·, t), P r (·, t) ∈ P p−1 (I i ). Next, we prove the following theorem which will be needed to prove our main superconvergence result.
into (19b), then choosing z = x−xi−1 hi P u and w = x−xi hi P q , we get , and e q =ē q + ǫ q . Using (85a) with f = P u and (85b) with f = P q , we obtain 1 Summing over all elements, applying the Cauchy-Schwarz inequality, and invoking the estimate in (23b), we obtain P u 2 ≤ 2h e q P u ≤ Ch p+2 P u , P q 2 ≤ 2h e r P q ≤ Ch p+2 P q , which completes the proof of the first two estimates in (87). Next, we show the third estimate in (87). Substitutingē r = c i (t)+ x−xi−1 hi P r (x, t) into (19a) and taking v = x−xi hi P r , we obtain where we used the fact that (c i (t))

Consequently, we have
Adding and subtracting the constant θ i = f ′ (u(x i−1 , t)) to θ then summing over all elements, we obtain where δ k = N i=1 δ k,i , k = 1, . . . , 4 and Next, we will estimate each of these terms separately.
Estimate of δ 1 . Using the fact that |x − x i | ≤ h, x ∈ I i , applying the Cauchy-Schwarz inequality, and using the estimate (23b), we get Estimate of δ 2 . Applying (33), using the Cauchy-Schwarz inequality, the inverse inequality, the estimate |x − x i | ≤ h i , x ∈ I i , and invoking the estimate (23b), we obtain Estimate of δ 3 . Since e u =ē u + ǫ u and ǫ u is orthogonal to ((x − x i )P r ) x ∈ P p−1 (I i ) (due to the properties in (6a)), δ 3 simplifies to Multiplying by θ i and summing over all elements, we get Combining (91) and (92), we obtain Using the assumption |f ′ | ≤ C 1 on R, the estimate |x − x i | ≤ h, x ∈ I i , the Cauchy-Schwarz inequality, and the estimate (23b) gives Estimate of δ 4 . Using (34), we have On the other hand, we take z = (x − x i )P r in (18c) to get Combining (94) and (95), we obtain Applying the Cauchy-Schwarz inequality, the inverse inequality, using the fact that |x − x i | ≤ h i , x ∈ I i , and invoking the estimates in (23) yields Now, combining (88a) with the estimates (89), (90), (93), and (96), we deduce that which completes the proof of (87). Now, we are ready to prove that u h is O(h p+3/2 ) superconvergent to P − h u.
Estimate of T 3 . We start from (36). Using the estimate (84) and the projection result (7), we obtain Estimate of T 4 . Taking w =ē q in (19b), we obtain Combining (99) and (25d), we get Summing over all elements and using the periodic boundary conditions yields Substituting the definitions ofē u ,ē q andē r , given in (86), into (101), and using the fact that (ǫ u ) t , ǫ q and ǫ r are orthogonal to any piecewise constant functions (due to the properties given in ( 6)), we get Using the fact that |x − x i | ≤ h i and |x − x i−1 | ≤ h i for all x ∈ I i , then applying the Cauchy-Schwarz inequality, the projection result (7), and (87), we get Now, combining (26) with (29), (35), (98), (102), and (39), we obtain Integrating with respect to time and using (12) By the continuous Gronwall's inequality (see, e.g., [33]), we conclude that, where C is independent of h. Thus, we completed the proof of the theorem.
Remark 3. We have also proved similar superconvergence rate for ē q and ē r . These results are not included here to save space.

Numerical examples.
In this section, we present some numerical experiments to validate our theoretical results. The initial condition is determined by u h (x, 0) = P h u(x, 0). Temporal integration is performed by the fourth-order classical implicit Runge-Kutta method. A time step ∆t is chosen so that temporal errors are small relative to spatial errors. We do not discuss the influence of the time discretization error in this paper. In the first two numerical examples, we examine the orders of convergence and superconvergence for the nonlinear KdV equation with periodic and non-periodic boundary conditions. In the third and fourth examples, we apply the LDG method to solve some famous dispersive wave problems. These tests are used by many authors; see e.g. [12,34].
Example 1. In this example, we consider the following nonlinear KdV problem subject to the periodic boundary condition We select g(x, t) such that the exact solution is u(x, t) = sin(x + t). We note that f ′ (u) = 6u changes sign in the computational domain. In this case, we use the Godunov flux which is an upwind flux. We solve this problem using the LDG method on uniform meshes obtained by partitioning the computational domain [0, 2π] into N subintervals with N = 10, 20, 30, 40, 50, 60, 70, 80 and using the spaces V p h with p = 1 − 4. Tables 1-3 show the L 2 errors ||e u ||, ||e q ||, and ||e r || at t = 1 and their orders of convergence. These results indicate that ||e u ||, ||e q ||, and ||e r || are O(h p+1 ). This is in full agreement with the theory. Thus, the error estimates proved in this paper are optimal in the exponent of the parameter h. In Table 4, we report the L 2 -norm of the errors ||ē u || as well as their orders of convergence. We observe that ||ē u || = O(h p+2 ). Thus, the LDG solution u h is superconvergent with order p + 2 towards the Gauss-Radau projection P − h u. Although the superconvergence rate is proved to be of order p + 3/2, our computational results indicate that the observed numerical convergence rate is higher than the theoretical rate. This example shows that the theoretical results actually holds true for nonlinear problems with general flux function f (u), indicating that the restriction on f ′ (u) is artificial.
Example 2. In this example, we consider the following nonlinear KdV equation  where g is chosen so that the exact solution is u(x, t) = sin(x + t) on the domain (x, t) ∈ [0, π] × [0, 1]. We extract the initial and boundary conditions from the exact solution. More precisely, the initial condition is u 0 (x) = sin(x) and the boundary conditions are u(0, t) = sin(t), u x (π, t) = − cos(t) and u xx (π, t) = sin(t). We solve (104) on uniform meshes having N = 10, 20, 30, 40, 50, 60, 70, 80 elements. We compute the LDG solutions in the spaces V p h , p = 1 − 4. The errors and their orders presented in Tables 5-7 suggest optimal O(h p+1 ) convergence rate for u, q, and r. In Table 8, we present the L 2 -norm of the errors ||ē u ||. We observe that ||ē u || = O(h p+2 ). Thus, the LDG solution u h is superconvergent with order p + 2 to the Gauss-Radau projection P − h u. Again this example indicates that the observed numerical superconvergence rate is higher than the theoretical rate which is proved to be of order p + 3/2.   The space-time graphs of the exact solution u and the LDG solution u h using N = 80 and p = 3 are displayed in Figure 1. We observe an excellent match between the exact and numerical solutions. These results are in full agreement with those published in [12,34].
The space-time graphs of the exact solution u and the LDG solution u h using N = 80 and p = 3 are shown in Figure 2. It can be seen from these graphs that the two solitary waves are moving toward the same direction. At time around t = 0.5, the two waves overlap. The faster soliton overtakes the slower one and then proceeds on its way. We also observe that both solitons continue to propagate, but the faster soliton is shifted forward while the slower soliton is shifted backward. These results in complete agreement with the results in [34]. 6. Concluding remarks. In this paper, we studied a superconvergent local discontinuous Galerkin (LDG) finite element method for nonlinear KdV-type problems. We proved several L 2 error estimates and superconvergence results toward special projections. In particular, we showed that the LDG solutions converge to the true solutions with order p + 1, when the space of piecewise polynomials of degree p is used. Our numerical experiments demonstrate that the rate of convergence is sharp. We further proved that the derivative of the LDG solution is superconvergent with order p + 1 towards the derivative of the Gauss-Radau projection of the exact solution. Finally, we proved that the LDG solution is O(h p+3/2 ) superconvergent towards the Gauss-Radau projection of the exact solution. All proofs are valid for arbitrary regular meshes using P p polynomials with p ≥ 1 and under the condition that nonlinear flux function, f ′ (u) possesses a uniform positive lower bound. As in [30], our numerical experiments demonstrate that the results in this paper hold true for nonlinear conservation laws with general flux functions, indicating that the restriction on f (u) is artificial. The generalization of our proofs to nonlinear equations with general flux functions involves several technical difficulties and will be investigated in the future. Superconvergence of the current paper will used to construct a posteriori error estimates by solving a local problem on each element. This will be reported in the second part of this work, where we will prove that the a posteriori LDG error estimates for the solution converges to the true error in the L 2 -norm under mesh refinement. We are currently investigating the superconvergence properties of the LDG method applied to 2D KdV-type problems on rectangular and triangular meshes. Also, because we observed superconvergence of order p + 2 in our numerical examples, future work will include investigating how to improve our proofs to obtain optimal superconvergence results. We expect that a new technique will be needed to obtain the optimal rate of superconvergence.