Superconvergence Analysis of Discontinuous Galerkin Methods for Systems of Second-Order Boundary Value Problems

: In this paper, we present an innovative approach to solve a system of boundary value problems (BVPs), using the newly developed discontinuous Galerkin (DG) method, which eliminates the need for auxiliary variables. This work is the ﬁrst in a series of papers on DG methods applied to partial differential equations (PDEs). By consecutively applying the DG method to each space variable of the PDE using the method of lines, we transform the problem into a system of ordinary differential equations (ODEs). We investigate the convergence criteria of the DG method on systems of ODEs and generalize the error analysis to PDEs. Our analysis demonstrates that the DG error’s leading term is determined by a combination of speciﬁc Jacobi polynomials in each element. Thus, we prove that DG solutions are superconvergent at the roots of these polynomials, with an order of convergence of O ( h p + 2 ) .


Introduction
The discontinuous Galerkin (DG) technique was first introduced by Reed and Hill [1] in 1973, with the aim of addressing the hyperbolic equation of neutron transport.They developed a discontinuous Galerkin method, a locally conservative and parallelizable method that can effectively handle complex geometries and does not require continuity across element boundaries.These desirable features have contributed to the DG method's widespread use in solving various challenging real-world problems including, but not limited to, Troesch's problem [2], α-synuclein spreading in Parkinson's disease [3], fully coupled hydro-mechanical modeling of two-phase flow in deformable fractured porous media [4], the recovery of conductivity in electrical impedance tomography [5], wave propagation phenomena in thermo-poroelastic media [6], Maxwell's equations in optics and photonics [7], solidification problems in a semitransparent medium-filled cavity [8], and seismic wave propagation problems [9].
A notable refinement of the DG approach emerged with the introduction of the Runge-Kutta DG scheme by Cockburn et al.This advancement was detailed in a sequence of works [10][11][12][13], where they aimed to tackle nonlinear hyperbolic conservation laws.Over time, the DG methodology has been extended to handle scenarios involving higherorder derivatives.A comprehensive and current overview of the evolution and applications of the DG method was presented by Shu in his surveys [14,15].
In 2007, Adjerid and Temimi [16] developed a novel DG scheme for handling higherorder initial value problems without requiring any auxiliary variables like the local DG methods.Using p-degree piecewise polynomials, they proved that their method achieves a convergence rate of p + 1.In [17], Chen and Shu introduced a range of ultra-weak discontinuous Galerkin (UWDG) techniques.These methods were designed to address boundary value problems encompassing second-order to fifth-order ordinary differential equations (ODEs).Their strategy combined the DG method for spatial discretization with the TVD high-order Runge-Kutta method for temporal discretization.By carefully selecting numerical fluxes, they demonstrated the stability of their proposed schemes and numerically showed that their method achieves an optimal (p + 1)-th convergence rate despite the theoretical sub-optimal error estimates.
In [18], Adjerid and Temimi expanded upon their earlier research and employed it in the context of the wave equation through the utilization of the method of lines.This approach involves the integration of the standard finite element method and the discontinuous Galerkin method in the spatial and temporal dimensions, respectively.The researchers showcased the efficacy of this new technique by deriving a range of optimal error estimates in both space and time domains.In addition, they conducted a comparative analysis against existing methodologies.Furthermore, the study revealed the superconvergence of the DG solution within each space-time element, specifically at the intersecting points of the Lobatto polynomials in space and the Jacobi polynomials in time.Later, Temimi [19] solved the one-dimensional second-order BVPs using the developed DG scheme and demonstrated that a specific combination of Jacobi polynomials generates the leading term of the discretization error in each element.He proved that the p-degree DG solution achieves an O(h p+2 ) superconvergence rate at the roots of these particular polynomials.Baccouch and Temimi [20] further extended the DG error analysis to secondorder BVPs and showed that when using p-degree piecewise polynomials, the UWDG solution and its derivative exhibit a superconvergence rate of O(h 2p ) at the upwind and downwind endpoints.Moreover, Baccouch and Temimi [21] proposed a novel DG scheme for solving the wave equation based on the method of lines.They demonstrated that the DG solution achieves in the L 2 -norm an optimal rate of convergence.
Recently, several superconvergence studies have elaborated.In [22], Baccouch presented superconvergence results for the local discontinuous Galerkin method applied to the sine-Gordon nonlinear hyperbolic equation in one space dimension.Additionally, Baccouch [23,24] explored the convergence and superconvergence properties of a local discontinuous Galerkin method for nonlinear second-order two-point boundary value problems.Subsequently, Baccouch [25] delved into the convergence and superconvergence properties of an ultra-weak discontinuous Galerkin method for linear fourth-order boundary value problems.In a related vein, Ma [26] presented a comprehensive analysis of superconvergence properties for a wide class of mixed discontinuous Galerkin methods.In [27], Liu et al. investigated the superconvergence properties of local discontinuous Galerkin methods with generalized alternating fluxes for one-dimensional linear convection-diffusion equations.More recently, in 2023, Singh et al. [28] provided a detailed superconvergence error analysis of the discontinuous Galerkin method with interior penalties for 2D elliptic convection-diffusion-reaction problems.
To the best of our knowledge, none of the previous references have investigated a system of boundary value problems or provided a detailed form of the leading term of the DG errors, which serve as a foundation for interesting superconvergence criteria.In our current work, we developed a local error analysis of the DG method for systems of secondorder BVPs.Our findings revealed that, within each element, the error's leading term is a combination of specific Jacobi polynomials.Notably, we demonstrated that the DG solutions are O(h p+2 ) superconvergent at these combined Jacobi polynomial roots.These superconvergence results can be further applied to compute efficient a posteriori error estimates and to construct higher-order DG solutions.The results of our current analysis, combined with our previous work involving the method of lines, will be integrated in future phases to develop a fully DG scheme for solving higher-dimensional partial differential equations.We believe that the new scheme will exhibit superior accuracy and efficiency compared to existing methods.This manuscript is summarized in this way.We develop a DG scheme applied to systems of boundary value problems in Section 2. The local error analysis with DG superconvergence properties are provided in Section 3.Then, we carry out several computational simulations to exhibit the full agreement with theoretical findings in Section 4. Finally, some concluding remarks are stated in Section 5.

Model Problem
In this section, we develop the DG formulation of the system of BVPs given by ÿ subject to mixed boundary conditions The Dirichlet boundary conditions scenario is also explored: where where C and K are two matrices defined respectively by C = (C ij ) 1≤i,j≤n and K = (K ij ) 1≤i,j≤n .Also, a smooth exact vector solution is achieved by properly choosing the vector function f.To implement the discontinuous Galerkin method, we initially create a partition h = (b−a) N+1 .For k = 0, 1, 2, • • • , N, we let x k = a + k • h and I k = (x k , x k+1 ); we also construct a finite-dimensional space S N,p where P p designates the p-degree polynomial space.Next, we multiply (1a) by a vector test function We integrate the resulting system of equations over I k and we integrate it by parts twice to derive the weak discontinuous Galerkin (DG) formulation for (1).Then, for k = 0, 1, . . ., N we have In (4), replacing y by where and Ẏk (x k+1 ) denote the numerical fluxes defined by Cheng et al. [17] as follows: The numerical fluxes, for the case of Dirichlet boundary conditions, are given by and Therefore, the discrete formulation involves finding and When subjected to Dirichlet boundary conditions, (10c) is written as

Local Error Analysis
In order to perform a complete analysis of the method, we first define the Jacobi polynomials [29] using Rodgrigues' formula: Ref. [30] states some valuable properties of Jacobi polynomials that would be useful to formulate the error's leading term.
In order to investigate the local behavior of the error, we analyze the problem (1) in its reduced form defined on one reference element [0, h]: subject to Therefore, problem ( 12) can be written as, for i = 1, . . ., n, subject to This simplified version of the problem, even if it is not the original one, affords identifying the error's leading term.In this paper, we do not show the scenario of the Dirichlet boundary conditions given that it would be readily achievable by adopting almost identical analysis steps.
Therefore, the weak DG formulation of ( 13) is given by and the discrete formulation entails finding Next, we replace v by V ∈ P p in (14) and subtract (15) to achieve the discontinuous Galerkin condition of the local error e By implementing a linear mapping from [0, h] to the canonical element [−1, 1] and denoting by ê, the mapped local error on [−1, 1], (16) becomes Next, we present and demonstrate our major finding on the local discretization error.
Proof.Using Theorem 1, we have that, for each i = 1, . . ., n, the local discretization error satisfies using (18b), we obtain By letting xj be the roots of which can be written as Then, by mapping êi and shifting xj to each element I k , the proof is achieved.

Numerical Examples
In order to validate our theory, we consider two systems of second-order differential equations subject to Neumann and Dirichlet boundary conditions.The numerical rate of convergence is defined by ln( e (N 1 ) / e (N 2 ) ) ln(N 1 /N 2 ) , where e (N) denotes the error e = y i − Y i , for i = 1, 2, • • • , n using N elements.
Example 1.Let us consider the following system of three second-order differential equations subject to Neumann and Dirichlet boundary conditions We choose f 1 , f 2 and f 3 such that the exact solutions are given by We solve problem (37) using a uniform mesh with h = 0.05 and p = 2, 3, 4 and plot the errors u i − U i,DG for i = 1, 2, 3 versus x in Figures 1-3.These plots show that the errors are crossing the x-axis at the roots of the polynomial P 1,0 p (ξ) + p+1 p 2 P 1,0 p−1 (ξ) shifted to each element and at the downwind endpoints of each element.Therefore, these results confirm the established theory on the leading term of the errors for each of the three solutions.
Moreover, on uniform meshes with N = 10, 15, • • • , 30 steps, problem (37) is solved for p = 2, 3, 4, 5.We exhibit ||u i ( xj ) − U i,DG ( xj )|| ∞ for i = 1, 2, 3, where xj are the roots of the mapped polynomials Q i,p+1 , along with their convergence rates in Tables 1-3 As expected, there is full agreement between these results and the stated theory.Example 2. Next, let us consider the following system of ten second-order differential equations: subject to Dirichlet boundary conditions where We choose f i for i = 1, 2, • • • , 10 such that the exact solutions are given by We solve problem (39) using a uniform mesh with h = π 20 and p = 2, 3, 4. Without loss of generality, we plot the errors u 5 − U 5,DG and u 10 − U 10,DG versus x in Figures 4 and 5.These plots indicate that the errors intersect the x-axis at the roots of the polynomial P 1,0 p (ξ) + p+1 p 2 P 1,0 p−1 (ξ) shifted to each element and at the downwind endpoints of each element.Thus, these results perfectly support the established theory.

Conclusions
In this study, we conducted a local error analysis of the DG method for systems of second-order boundary value problems.We explored the superconvergence criteria and demonstrated that the leading error term on each element is a combination of specific Jacobi polynomials.We established that DG solutions exhibit superconvergence of O(h p+2 ) at the roots of these combined Jacobi polynomials.We further validated the established theory through some numerical simulations, and we showed the full agreement between the theory and the computational results.In our forthcoming paper, we will integrate the findings from this analysis with our earlier work involving the method of lines to create a comprehensive DG scheme for solving higher-dimensional partial differential equations.In our future research endeavors, we aim to expand the application of DG methods to address a variety of challenging real-world problems including, but not limited to, the two-dimensional hyperbolic partial differential equation of the Telegraph type [31], the integro-differential Beam problem [32], the nonlinear Bratu problem [33,34], delayed differential equations [35] and nonlinear Boussinesq and Klein-Gordon equations [36].