Analysis of optimal superconvergence of a local discontinuous Galerkin method for nonlinear second-order two-point boundary-value problems

Article history: Received 23 February 2019 Received in revised form 30 April 2019 Accepted 1 May 2019 Available online xxxx


Introduction
The purpose of this paper is to study the convergence and superconvergence properties of the local discontinuous Galerkin (LDG) method for the nonlinear two-point second-order boundary-value problems (BVPs) (1.1a) where u : [a, b] → R and f : D → R is a given smooth function on the set D = [a, b] × R 2 . In this paper, we consider one of the following set of boundary conditions, which are commonly encountered in practice: where α 1 , β 2 are given constants. In our analysis, we assume that the BVP (1.1) has one and only one solution. The conditions on f for the existence and uniqueness of the solution to the general BVP (1.1) are given in [22]. The nonlinear two-point BVP (1.1) arises in applied mathematics, theoretical physics, engineering, control and optimization theory; see e.g., [3,28]. Since the analytic solution to (1.1) is difficult to obtain for general f , numerical techniques are often needed to approximate its solution. Many authors have designed numerical schemes to solve second-order BVPs. We refer to [30,17,22,3,25,20,5,19,11,23,27,2,30,21] for some numerical methods including the shooting method, the finite difference method, the collocation method, the monotone iterative method, and the quasilinearization method.
Superconvergent numerical methods of BVPs are necessary in many important scientific and engineering applications such as boundary layer theory, the study of stellar interiors, control and optimization theory, and flow networks in biology. A knowledge of superconvergence properties can be used to (i) construct simple and asymptotically exact a posteriori estimates of discretization errors and (ii) help detect discontinuities to find elements needing limiting, stabilization and/or refinement. Typically, a posteriori error estimators employ the known numerical solution to derive estimates of the actual solution errors. They are also used to steer adaptive schemes where either the mesh is locally refined (h-refinement) or the polynomial degree is raised (p-refinement). In the past several decades, there also has been considerable interest in studying superconvergence properties of numerical methods. In this paper, we present new superconvergence results of the LDG method for solving (1.1). Discontinuous Galerkin (DG) methods form a class of high order numerical methods for solving ordinary differential equations (ODEs) and partial differential equations (PDEs). They combine many attractive features of the finite element and finite volume methods. DG schemes have been successfully applied to many problems arising from a wide range of applications. The DG method is a finite element method using a completely discontinuous piecewise polynomial space for the numerical solution and the test functions. DG methods are becoming important techniques for the computational solution of many real-world problems. They are known to have good stability properties when applied to hyperbolic PDEs. Furthermore, DG methods have been successfully applied to hyperbolic, elliptic, and parabolic problems arising from a wide range of applications. DG methods are highly accurate numerical methods with the advantage that they can handle problems having discontinuities such as those that arise in hyperbolic problems, can handle problems with complex geometries, simplify adaptive h−, p−, and r− refinement, and produce efficient parallel solution procedures. DG method was initially introduced by Reed and Hill in 1973 as a technique to solve neutron transport problems [29].
The local DG (LDG) methods are natural extension of the DG methods aimed at solving higher-order PDEs. The LDG method was first proposed by Cockburn and Shu in [16] for solving convection-diffusion problems. The LDG method consists of rewriting a higher order differential equation into a system of first-order equations and then discretizing it by the standard DG method. The success of LDG methods is due to the following properties: (i) they are robust and high order accurate, (ii) they can achieve stability without slope limiters, and (iii) they are element-wise conservative. This last feature is very important in the area of computational fluid dynamics, especially in situations where there are steep gradients or boundary layers or shocks. Moreover, LDG schemes are extremely flexible in the mesh-design. Thus, they can easily handle meshes with hanging nodes, elements of various types and shapes, and local spaces of different orders. Furthermore, they exhibit useful superconvergence properties that can be used to estimate the actual discretization errors. We refer the reader to e.g., [8,15,13,7,14] and references therein for a more complete survey of several LDG methods.
Several authors designed and analyzed the LDG method for BVPs of the form (1.1), see e.g., [9,24,33,36,32,31,37,4]. In [4], we proposed and analyzed a superconvergent and high order accurate LDG method for nonlinear two-point second-order BVPs of the form u = f (x, u) subject to some suitable boundary conditions. We proved optimal L 2 error estimates for the solution and for the auxiliary variable that approximates its first-order derivative. The order of convergence is proved to be p + 1, when piecewise polynomials of degree at most p ≥ 1 are employed. We further proved that the derivatives of the LDG solutions are superconvergent with order p + 1 toward the derivatives of Gauss-Radau projections of the exact solutions. Finally, we showed that the LDG solutions are superconvergent with order p + 3/2 toward Gauss-Radau projections of the exact solutions, while computational results show higher O(h p+2 ) convergence rate. However, a theoretical proof of this property remains open. The main purpose of our current work is to use a different approach to prove the (p + 2)th superconvergence rate and also the (2p + 1)th superconvergence rate at upwind and downwind points and for the domain average. To the best of our knowledge, these results are original.
In this paper, we design a superconvergent LDG method for nonlinear BVPs of the form (1.1). We prove several optimal L 2 error estimates for the LDG solutions. In particular, we prove that the LDG solutions to approximate u and u are (p + 1)th order convergent in the L 2 -norm, when the space of piecewise polynomials of degree p ≥ 1 is used. We further show that the derivatives of the LDG solutions are superconvergent with order p + 1 toward the derivatives of Gauss-Radau projections of the exact solutions. Moreover, we prove (p + 2)th order superconvergence of the LDG solutions toward Gauss-Radau projections of the exact solutions. Finally, we show that the errors between the LDG solutions and the exact solutions are (2p + 1)th order superconvergent at either the upwind point or downwind point in each element on regular meshes. Numerical experiments demonstrate that the theoretical orders of convergence and superconvergence are optimal. Our global error analysis is valid for any regular meshes and using piecewise polynomials of degree p ≥ 1 and for the classical set of boundary conditions. We would like to mention that the proposed LDG method has several advantages over the standard methods due to the following nice features: (i) it achieves arbitrary high order accuracy, (ii) it exhibits optimal convergence properties for the solution and for the auxiliary variables that approximate the derivatives, (iii) it can easily handle meshes using local spaces of different orders, and (iv) achieves superconvergence results that can be used to construct asymptot-ically exact a posteriori error estimates by solving a local problem on each element. This will be discussed in a separate paper. The rest of the paper is organized as follows: In section 2, we describe the LDG method for nonlinear second-order BVPs. We also present some preliminary results, which will be used in our error analysis. In section 3, we present a detailed proof of the optimal a priori error estimates of the LDG method. We state and prove our main superconvergence results in section 4. In section 5, we present several numerical examples to validate our theoretical results. Finally, we provide some concluding remarks in section 6.

The LDG scheme for nonlinear second-order BVPs
In order to define the LDG method, we introduce a new auxiliary variable q = u and convert (1.1a) into a first-order system of ODEs To obtain the LDG weak formulation, we partition the computational domain = [a, b] into a collection of non-overlapping elements We also define h = max 1≤i≤N h i and h min = min 1≤i≤N h i to be the lengths of the largest and smallest cells, respectively. We assume that the mesh is regular in the sense that there exists a constant K ≥ 1 independent of h such that h ≤ K h i , i = 1, . . . , N.
v(x i + s) to denote the left limit and the right limit of v at the discontinuity point Multiplying the two equations in (2.1) by arbitrary test functions v 1 and v 2 , integrating over the interval I i , and using integration by parts, we get (2.2b) Next, we introduce the following discontinuous finite element approximation space where P p (I i ) denotes the space of polynomials of degree at most p on I i . We would like to emphasize that polynomials in the finite element space V p h are allowed to be completely discontinuous at the mesh points. To obtain the LDG scheme, we replace the exact solutions u and q by piecewise polynomials of degree at most p and denote them by u h ∈ V p h and q h ∈ V p h . We also choose the test functions v 1 and v 2 to be piecewise polynomials of degree at most p. The LDG scheme can now be defined as: where u h and q h are the so-called numerical fluxes, which are, respectively, the discrete approximations to u and q at the nodes. These numerical fluxes must be designed based on different guiding principles for different differential equations to ensure stability and optimal error estimates. To complete the definition of the LDG scheme, we only need to define u h and q h on the boundaries of I i . It turns out that the following simple choices would guarantee the optimal convergence and superconvergence of our LDG scheme: For the mixed boundary conditions (1.1b), we take the following alternating numerical fluxes; see e.g., [4] where the stabilization parameter δ 1 for the LDG method is given by For the periodic boundary conditions (1.1d), we choose the following alternating fluxes (2.4c) Implementation: The LDG solution (u h , q h ) can be obtained using the following steps: to be a local basis of P p (I i ) and we express u h , q h as .
(3) We solve the nonlinear system for the unknown coefficients c 0,i , c 1,i . . . , c 2p+1,i , i = 1, . . . , N using e.g., Newton's method for nonlinear systems. Once we solve for the unknown coefficients, we get the LDG solutions u h and q h , which are piecewise discontinuous polynomials of degree ≤ p.

Norms:
We present some norms that will be used throughout the paper. Denote u 0, to be the standard L 2 -norm of the function u on I i . For any natural integer s, the Sobolev space H s (I i ) consists of functions that have generalized derivatives of order s in the space L 2 (I i ). It is defined by . We shall also use the following notation for the semi-norm |u| s,I i = D s u 0,I i . Finally, we define the norms on the whole computational domain as follows: For convenience, we use u , u s , and |u| s to denote u 0, , u s, , and |u| s, , respectively. We would like to mention that if u ∈ H s ( ), s = 1, 2, . . ., then the Sobolev norm u s, on the whole computational domain is the standard Sobolev Legendre polynomial: In our analysis we need the pth-degree Legendre polynomial defined by Rodrigues formula [1] which satisfies the following properties: L p (1) = 1, L p (−1) = (−1) p , and the orthogonality relation (2.5) Mapping the physical element Using the mapping (2.6) and the orthogonality relation (2.5), we obtain Gauss-Radau Projections: For p ≥ 1, we introduce two special Gauss-Radau projections P ± h . These projections are defined element-by-element as follows: For any integrable function u on , P ± h u ∈ V p h and the restrictions of P ± h u to I i are polynomials in P p (I i ) satisfying the conditions, see e.g [10] I i (2.9) By the scaling argument, we obtain the following projection results [12]: For any function u ∈ H p+1 ( ), there exists a positive constant C independent of the mesh size h, such that (2.10) Moreover, we recall the inverse properties of the finite element space V p h that will be used in our error analysis [26]: For 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 BVP (1.1). They also may have different values at different places.

A priori error estimates
In this section, we derive optimal L 2 error estimates for the LDG method. For convenience, we use e u and e q to denote the errors between the exact solutions of (2.1) and the LDG solutions defined in (2.3), i.e., As the usual treatment in finite element analysis, we divide the errors into the form e u = u +ē u , e q = q +ē q , (3.1) where the projection errors are defined by In our error analysis, we assume that the function f appearing in the right-hand side of (1.1a) is sufficiently differentiable function. More precisely, we assume that the f satisfies the following conditions: In the next theorem, we prove a priori error estimates for e u and e q in the L 2 -norm.
Proof. First, we derive some error equations which will be used repeatedly throughout this paper. Subtracting (2.3) from h , k = 1, 2 and using the numerical fluxes (2.4a), we obtain the error equations on Applying Taylor's Theorem with integral remainder and using (3.1), we write Under Assumption A2, we have (3.7) Using (3.6), we rewrite (3.5) as Adding the above equations, we get Summing over all elements gives . (3.11) On the other hand, performing integration by parts, we write (3.9) as We note that, with the numerical fluxes (2.4a), the jumps of e u and e q at an interior point x i are defined as Since e u (x − 0 ) = e q (x + N ) = 0, the jumps at the endpoints of the computational domain are given by Adding and subtracting P − h V 1 to V 1 and P + h V 2 to V 2 and using (3.8) with Combining (3.13) and (3.12) and using the properties of the projections P ± (3.14) By the property of the projection P ± h , we have (3.15) since w is a polynomial of degree at most p and thus w is a polynomial of degree at most p − 1. Substituting (3.1) into (3.14) and using (3.15) with w =ē u , ē q ∈ P p (I i ), we get Adding these two equations, we obtain Summing over all elements, we arrive at The main idea behind the proof of the theorem is to construct the following adjoint problem: find W 1 and W 2 such that The BVP (3.19) can be converted into the system of equations The solution to (3.20) can be expressed in terms of its fundamental matrix (3.21) where the 2 × 2 fundamental matrix M(x) satisfies the following initial-value problem which gives W k 1 ≤ C e u + e q , k = 1, 2.
Using (3.7) and applying the Cauchy-Schwarz inequality yields e q Applying the standard interpolation error estimate (2.10), we get e q 2 + e u 2 ≤ Ch p + K 1 e u + K 2 e q Ch |W 1 | 1 + Ch p + e q Ch |W 2 | 1 .
Applying the regularity estimate (3.23), we get e q 2 + e u 2 ≤ C 1 h p+1 e q + e u + C 2 h e q 2 + e u 2 .

Superconvergence error analysis
In this section, we investigate the superconvergence properties of the proposed LDG method. We prove that the derivatives of the LDG solutions are superconvergent with order p + 1 toward the derivatives of Gauss-Radau projections of the exact solutions. We further prove pointwise superconvergence results at the upwind and downwind points of each element.
, and applying (2.5) gives Applying the estimate (3.7), the Cauchy-Schwarz inequality, the inverse inequality, and the estimate (2.7) yields Consequently, we deduce ē q 0,I i ≤ C 1 e u 0,I i + e q 0,I i , ē u 0,I i ≤ C 2 e q 0,I i .

Pointwise superconvergence
First, we prove superconvergence results at the endpoints of the computational domain . We state them in the following theorem.
Proof. We construct the following adjoint problem: find U 1 and U 2 such that (4.7) The BVP (4.7) can be transformed into the system of equations   Taking V 1 = U 1 and V 2 = U 2 in (3.18) and applying (4.7), we get Using (3.7) and applying the Cauchy-Schwarz inequality yields Applying (3.3), (3.4), the standard interpolation error estimate (2.10) and using the regularity estimate (4.11), we get which completes the proof of (4.5).
Next, we will prove (4.6). We construct the following adjoint problem: find U k , k = 1 − 4 such that (4.12) As before, the BVP (4.12) satisfies the regularity estimate (4.13) Taking V k = U k , k = 1, 2 in (3.18) and applying (4.12), we get Following the steps used to prove (4.5), we establish (4.6). 2 In the next theorem, we state and prove superconvergence results at the downwind and upwind points. More precisely, we prove the following theorem.

Theorem 4.3. Suppose that the assumptions of Theorem 4.2 hold. Then there exists a constant C such that
Proof. For j = 1, 2, . . . , N, we consider the following terminal problem: find U 1 , U 2 ∈ H p+1 [a, x j ] such that Under the assumption of the theorem, one can easily verify that (4.16) satisfies the following regularity estimate (4.17) Taking V 1 = U 1 and V 2 = U 2 in (3.10), summing over the elements I i , i = 1, 2, . . . , j, using (4.16) and the fact that e u ( On the other hand, taking V k = U k , k = 1, 2 in (3.16) and summing over the elements (4.19) Combining the two formulas (4.18) and (4.19) yields Using (3.7) and applying the Cauchy-Schwarz inequality, we get Using the standard interpolation error estimates (2.10), the estimates (3.3), (3.4), (4.6), and the regularity estimate (4.17), we obtain which completes the proof of (4.14).
Next, we will show (4.15). We construct the following terminal problem: find (4.20) Under the assumptions of the theorem, one can easily verify that (4.20) satisfies the following regularity estimate (4.21) Taking V 1 = U 1 and V 2 = U 2 in (3.10), summing over the elements I i , i = j, . . . , N, using (4.20) and the fact that On the other hand, taking V k = U k , k = 1, 2 in (3.16) and summing over the elements I i , i = j, j + 1, . . . , N, we get Using (3.7) and applying the Cauchy-Schwarz inequality, we get Using the standard interpolation error estimates (2.10), the estimates (3.3), (3.4), (4.5), and the regularity estimate (4.21), we obtain which completes the proof of (4.15). 2 (4.24) Proof. Using (3.1) and the properties of the projections P ± h , we have e u (x − j ) =ē u (x − j ) and e q (x + j ) =ē q (x + j ). Invoking the estimates (4.14) and (4.15), we establish (4.24). 2

Superconvergence for average errors at downwind/upwind points
Next, we deduce the (2p + 1)th superconvergence rate for the average errors at downwind/upwind points. More specifically, we have the following superconvergence results.

Corollary 4.2. Suppose that the assumptions of Theorem 4.2 hold. Then there exists a constant C such that
(4.25) Proof. These results follow immediately from (4.14) and (4.15). 2

Superconvergence toward Gauss-Radau projections
Next, we prove that the LDG solutions are superconvergent with order p + 2 toward Gauss-Radau projections of the exact solutions.
Taking the square of both sides, applying the inequality (a + b) 2 ≤ 2a 2 + 2b 2 and applying the Cauchy-Schwartz inequality, we get Integrating these inequalities with respect to x, using the estimates in (4.24), and the fact that h i ≤ h, we get Summing over all elements and using the estimates in (4.1), we obtain (4.27b) where we used 4p + 2 ≥ 2p + 4 for p ≥ 1. This completes the proof of (4.26). 2

Numerical examples
In this section, we present numerical examples to verify our theoretical findings.
Example 5.1. We consider the following nonlinear second-order BVP (5.1a) subject to the periodic boundary conditions The exact solution of (5.1) is u = 2 + sin(x). We solve ( We display the all errors in Figures 5-8. These results suggest optimal convergence and superconvergence rates. Again these results are in full agreement with the theoretical results.
Finally, we solve the following BVP subject to the Dirichlet boundary conditions We present the errors in Figures 9-12. Again, these results suggest optimal convergence and superconvergence rates when the Dirichlet boundary conditions are used. The exact solution of (5.4) is given by u(x) = x 2 + 16/x. We solve (5.4) using the proposed LDG method with p = 1 − 4. We use uniform meshes obtained by subdividing [1,2] into N intervals with N = 4, 6, 8, 10, 12, 14. In Fig. 13 we display the         O(h p+1 ). Thus, the error estimates proved in this paper are optimal in the exponent of the parameter h. In Fig. 14, we report the L 2 -norm of the errors ē u and ē q and their orders of convergence. We observe that ē u = O(h p+2 ) and ē q = O(h p+2 ). Thus, the LDG solutions u h and q h are, respectively, superconvergent with order p + 2 to the particular projections P − h u and P + h q. Finally, we report the maximum errors at downwind and upwind points in Fig. 15. We observe that the convergence rate of max , u(1) = 0, u (2) = 4 + 12 ln(2), (5.5) which admits the unique solution u = x 3 ln(x) [6]. We solve this problem using the LDG method on uniform meshes having N = 5, 10, 15, 20, 25, 30, 35, 40 elements and using the spaces P p with p = 1, 2, 3, and 4. The L 2 errors as well as their order of convergence are shown in Fig. 16. These results indicate that the LDG method yields O(h p+1 ) convergent solutions. The rate of convergence is clearly optimal. The L 2 -norm of the errors ||ē u || and ||ē q || shown in Fig. 17            where the exact solution is u(x) = −2 ln(cos(x)). We use uniform meshes obtained by subdividing the computational domain [0, 1] into N intervals with N = 10, 20, 30, 40, 50. Fig. 19 shows the L 2 errors e u and e q and their orders of convergence. These results indicate that e u and e q are both O(h p+1 ). In Fig. 20, we report the L 2 -norm of the errors ē u and ē q and their orders of convergence. We observe that ē u = O(h p+2 ) and ē q = O(h p+2 ). Thus, the LDG solutions u h and q h are, respectively, superconvergent with order p + 2 to the particular projections P − h u and P + h q.

Concluding remarks
In this paper, we investigated the convergence and superconvergence of a local discontinuous Galerkin (LDG) finite element method for nonlinear second-order boundary-value problems (BVPs) for ordinary differential equations. We proved several L 2 error estimates, and superconvergence results toward special projections. To be more precise, we proved that the LDG solutions converge to the true solutions with order p + 1, when the space of piecewise polynomials of degree p ≥ 1 is used. Moreover, we showed that the derivatives of the LDG solutions converge with order p + 1 toward the derivatives of Gauss-Radau projections of the exact solutions. In addition, we also proved that the LDG solutions are superconvergent with order p + 2 to the Gauss-Radau projections of the exact solutions. Finally, we established superconvergence rate of order 2p + 1 for the maximum errors at the upwind or downwind points. Numerical experiments demonstrate that the error bounds are sharp. Our current and future works include two-dimensional elliptic, parabolic, and hyperbolic problems, which would be more challenging and interesting.

Declaration of Competing Interest
The author declares that he has no conflict of interest.