The errors and the global effectivity indices for Example 7.1 on uniform meshes having 5, 10, 20, 30, 40, 50 elements using , to 4.
We present an analysis of the discontinuous Galerkin (DG) finite element method for nonlinear ordinary differential equations (ODEs). We prove that the DG solution is $(p + 1) $th order convergent in the $L^2$-norm, when the space of piecewise polynomials of degree $p$ is used. A $ (2p+1) $th order superconvergence rate of the DG approximation at the downwind point of each element is obtained under quasi-uniform meshes. Moreover, we prove that the DG solution is superconvergent with order $p+2$ to a particular projection of the exact solution. The superconvergence results are used to show that the leading term of the DG error is proportional to the $ (p + 1) $-degree right Radau polynomial. These results allow us to develop a residual-based a posteriori error estimator which is computationally simple, efficient, and asymptotically exact. The proposed a posteriori error estimator is proved to converge to the actual error in the $L^2$-norm with order $p+2$. Computational results indicate that the theoretical orders of convergence are optimal. Finally, a local adaptive mesh refinement procedure that makes use of our local a posteriori error estimate is also presented. Several numerical examples are provided to illustrate the global superconvergence results and the convergence of the proposed estimator under mesh refinement.
- discontinuous Galerkin finite element method
- ordinary differential equations
- a priori error estimates
- a posteriori error estimates
- adaptive mesh refinement
In this chapter, we introduce and analyze the discontinuous Galerkin (DG) method applied to the following first‐order initial‐value problem (IVP)
where , and . We assume that the solution exists and is unique and we would like to approximate it using a discontinuous piecewise polynomial space. According to the ordinary differential equation (ODE) theory, the condition is sufficient to guarantee the existence and uniqueness of the solution to (1). We note that a general th‐order IVP of the form with initial conditions , can be converted into a system of equations in the form (1), where , , and .
The high‐order DG method considered here is a class of finite element methods (FEMs) using completely discontinuous piecewise polynomials for the numerical solution and the test functions. The DG method was first designed as an effective numerical method for solving hyperbolic conservation laws, which may have discontinuous solutions. Here, we will discuss the algorithm formulation, stability analysis, and error estimates for the DG method solving nonlinear ODEs. DG method combines the best proprieties of the classical continuous finite element and finite volume methods such as consistency, flexibility, stability, conservation of local physical quantities, robustness, and compactness. Recently, DG methods become highly attractive and popular, mainly because these methods are high‐order accurate, nonlinear stable, highly parallelizable, easy to handle complicated geometries and boundary conditions, and capable to capture discontinuities without spurious oscillations. The original DG finite element method (FEM) was introduced in 1973 by Reed and Hill  for solving steady‐state first‐order linear hyperbolic problems. It provides an effective means of solving hyperbolic problems on unstructured meshes in a parallel computing environment. The discontinuous basis can capture shock waves and other discontinuities with accuracy [2, 3]. The DG method can easily handle adaptivity strategies since the refinement (mesh refinement and coarsening) and the refinement (method order variation) can be done without taking into account the continuity restrictions typical of conforming FEMs. Moreover, the degree of the approximating polynomial can be easily changed from one element to the other . Adaptivity is of particular importance in nonlinear hyperbolic problems given the complexity of the structure of the discontinuities and geometries involved. Due to local structure of DG methods, physical quantities such as mass, momentum, and energy are conserved locally through DG schemes. This property is very important for flow and transport problems. Furthermore, the DG method is highly parallelizable [4, 5]. Because of these nice features, the DG method has been analyzed and extended to a wide range of applications. In particular, DG methods have been used to solve ODEs [6–9], hyperbolic [5, 6, 10–19] and diffusion and convection diffusion [20–23] partial differential equations (PDEs), to mention a few. For transient problems, Cockburn and Shu  introduced and developed the so‐called Runge‐Kutta discontinuous Galerkin (RKDG) methods. These numerical methods use DG discretizations in space and combine it with an explicit Runge‐Kutta time‐marching algorithm. The proceedings of Cockburn et al.  and Shu  contain a more complete and current survey of the DG method and its applications. Despite the attractive advantages mentioned above, DG methods have some drawbacks. Unlike the continuous FEMs, DG methods produce dense and ill‐conditioned matrices increasing with the order of polynomial degree.
Related theoretical results in the literature including superconvergence results and error estimates of the DG methods for ODEs are given in [7–9, 26–28]. In 1974, LaSaint and Raviart  presented the first error analysis of the DG method for the initial‐value problem (1). They showed that the DG method is equivalent to an implicit Runge‐Kutta method and proved a rate of convergence of
This chapter is organized as follows: In Section 2, we present the discrete DG method for the classical nonlinear initial‐value problem. In Section 3, we present a detailed proof of the optimal a priori error estimate of the DG scheme. We state and prove our main superconvergence results in Section 4. In Section 5, we present the a posteriori error estimation procedure and prove that these error estimates converge to the true errors under mesh refinement. In Section 6, we propose an adaptive algorithm based on the local a posteriori error estimates. In Section 7, we present several numerical examples to validate our theoretical results. We conclude and discuss our results in Section 8.
2. The DG scheme for nonlinear IVPs
The error analysis of nonlinear scalar and vector initial‐value problems (IVPs) having smooth solutions is similar. For this, we restrict our theoretical discussion to the following nonlinear initial‐value problem (IVP)
where is a sufficiently smooth function with respect to the variables and . More precisely, we assume that on the set , where is a positive constant. We note that the assumption is sufficient to ensure that satisfies a Lipschitz condition in on the convex set with Lipschitz constant
Next, we introduce the DG method for the model problem (2). Let be a partition of the interval . We denote the mesh by . We denote the length of by . We also denote and as the length of the largest and smallest subinterval, respectively. Here, we consider regular meshes, that is, , where is a constant (independent of ) during mesh refinement. If , then the mesh is uniformly distributed. In this case, the nodes and mesh size are defined by ,
Throughout this work, we define and to be the left limit and the right limit of the function at the discontinuity point , and . To simplify the notation, we denote by the jump of at the point .
If we multiply (2) by an arbitrary test function , integrate over the interval , and integrate by parts, we get the DG weak formulation
We denote by the finite element space of polynomials of degree at most in each interval , ,
where denotes the set of all polynomials of degree no more than on . We would like to emphasize that polynomials in are allowed to have discontinuities at the nodes .
Replacing the exact solution by a piecewise polynomial and choosing , we obtain the DG scheme: Find such that and ,
where is the so‐called numerical flux which is nothing but the discrete approximation at the node . We remark that is not necessarily continuous at the nodes.
To complete the definition of the DG scheme, we still need to define on the boundaries of . Since for IVPs, information travel from the past into the future, it is reasonable to take as the classical upwind flux
The DG solution can be efficiently obtained in the following order: first, we compute in the first element using (5a) and (5b) with since is known. Then, we can find in since in is already available. We can repeat the same process to compute in . More specifically, can be obtained locally for each using the following two steps: (i) express as a linear combination of orthogonal basis , where is the th degree Legendre polynomial on , i.e., , where is a local basis of , and (ii) choose the test functions . Thus, on each , we get a system of nonlinear algebraic equations, which can be solved for the unknown coefficients using, e.g., Newton's method for nonlinear systems. Once we obtain the DG solution on all elements , we get the DG solution which is a piecewise discontinuous polynomial of degree . We refer to [7–9] for more details about DG methods for ODEs as well as their properties and applications.
2.2. Linear stability for the DG method
Let us now establish a stability result for the DG method applied to the linear case, i.e., . Taking in the discrete weak formulation (5a), we get
which is equivalent to
Summing over all elements, we get the equality
Consequently, which gives the stability result provided that .
3. A priori error analysis
We begin by defining some norms that will be used throughout this work. We define the inner product of two integrable functions, and , on the interval as . Denote to be the standard norm of on . Moreover, the standard norm of on is defined by Let , where , denote the standard Sobolev space of square integrable functions on with all derivatives being square integrable on , i.e., and equipped with the norm The seminorm of a function on is given by . We also define the norms on the whole computational domain as follows:
The seminorm on the whole computational domain is defined as We note that if , , then the norm on the whole computational domain is the standard Sobolev norm For convenience, we use and to denote and , respectively.
For , we consider two special projection operators, , which are defined as follows: For a smooth function , the restrictions of and to are polynomials in satisfying
These two particular Gauss‐Radau projections are very important in the proofs of optimal error estimates and superconvergence results. We note that the special projections are mainly utilized to eliminate the jump terms at the cell boundaries in the error estimate in order to achieve the optimal order of accuracy .
For the projections mentioned above, it is easy to show that for any with , there exists a constant independent of the mesh size such that (see, e.g., )
Moreover, we recall the inverse properties of the finite element space that will be used in our error analysis: For any , there exists a positive constant independent of and , such that, ,
From now on, the notation etc. will be used to denote generic positive constants independent of , but may depend upon the exact solution of (1) and its derivatives. They also may have different values at different places.
Throughout this work, let us denote to be the error between the exact solution of (2) and the DG solution defined in (5a) and (5b), to be the projection error, and to be the error between the projection of the exact solution and the DG solution . We observe that the actual error can be written as
Now, we are ready to prove our optimal error estimates for in the and norms.
Theorem 3.1. Suppose that the exact solution of (2) is sufficiently smooth with bounded derivatives, i.e., is bounded. We also assume that on . Let and be the DG solution of (5a) and (5b), then, for sufficiently small , there exists a positive constant independent of such that,
Proof. We first need to derive some error equations which will be used repeatedly throughout this and the next sections. Subtracting (5a) from (4) with and using the numerical flux (5b), we obtain the following error equation: ,
By integration by parts, we get
Applying Taylor's series with integral remainder in the variable and using the relation , we write
To simplify the notation, we introduce the bilinear operator as
Thus, we can write (14) as
A direct calculation from integration by parts yields
On the other hand, if we add and subtract to then we can write (15) as
If is a polynomial of degree at most then is a polynomial of degree at most . Therefore, by the property of the projection , we immediately see
Now, we are ready to prove the theorem. We construct the following auxiliary problem: find such that
where . Clearly, the exact solution to (22) is given by the explicit formula
Next, we prove some regular estimates which will be needed in our error analysis. Using the assumption , we see that is bounded by
Using the definition of and the estimate (24a), we have
Similarly, we can easily estimate as follows
Squaring both sides and intergrading over yields
Squaring both sides, applying the inequality , integrating over the computational domain , and using (25a), we get
Applying the projection result and the estimate (3.20b) yields
Summing over the elements and using the fact that yields
On the other hand, if we choose in (21) then we get
Summing over all elements and applying the Cauchy‐Schwarz inequality, we get
Using the estimate (25c), we deduce that
Thus, , where is a positive constant independent of . Therefore, for sufficiently small , e.g., , we obtain , which yields for small. Thus, we completed the proof of (9).
We note that Applying (9) and the estimate yields
4. Superconvergence error analysis
In this section, we study the superconvergence properties of the DG method. We first show a th order superconvergence rate of the DG approximation at the downwind point of each element. Then, we apply this superconvergence result to show that the DG solution converges to the special projection of the exact solution at
First, we define some special polynomials. The th degree Legendre polynomial can be defined by Rodrigues formula 
It satisfies the following important properties: , , and the orthogonality relation
One can easily write the degree Legendre polynomial on as
The degree right Radau polynomial on is defined as . It has real distinct roots, .
Mapping into the reference element by the linear transformation , we obtain the shifted Legendre and Radau polynomials on :
Next, we define the monic Radau polynomial, , on as
Throughout this work the roots of are denoted by
In the next lemma, we recall the following results which will be needed in our error analysis .
Lemma 4.1. The polynomials and satisfy the following properties
where , , and .
Now, we are ready to prove the following superconvergence results.
Theorem 4.1. Suppose that the assumptions of Theorem 1 are satisfied. Also, we assume that is sufficiently smooth with respect to and (for example, is enough.). Then there exists a positive constant such that
Proof. To prove (33), we proceed by the duality argument. Consider the following auxiliary problem:
where and . The exact solution of this problem is , Using the assumption and the estimate (24a), we can easily show that there exists a constant such that
Summing over the elements , using , and the fact that , we obtain
Now, taking in (21) yields
Summing over all elements with and applying (39), we arrive at
Using (24a) and applying the Cauchy‐Schwarz inequality, we obtain
for all which completes the proof of (33).
Next, we will derive optimal error estimate for . By the property of , we have
By integration by parts on the first term, we obtain
Using (3) and applying the Cauchy‐Schwarz inequality gives
Consequently, . Taking the square of both sides, summing over all elements, and using (9), we conclude that
Finally, we will estimate . Using the fundamental theorem of calculus, we write
Taking the square of both sides, applying the inequality , and applying the Cauchy‐Schwartz inequality, we get
Integrating this inequality with respect to and using the estimate (34), we get
Summing over all elements and using the estimate (35) and the fact that , we obtain
where we used for . This completes the proof of the theorem.
The previous theorem indicates that the DG solution is closer to than to the exact solution . Next, we apply the results of Theorem 2 to prove that the actual error can be split into a significant part, which is proportional to the degree right Radau polynomial, and a less significant part that converges at
Next, we recall the following results from  which will be needed in our analysis.
Lemma 4.2. If , then interpolation error can be split as
and is the coefficient of in the degree polynomial . Furthermore,
Proof. The proof of this lemma can be found in , more precisely in its Lemma 2.1.
The main global superconvergence result is stated in the following theorem.
Theorem 4.2. Under the assumptions of Theorem 2, there exists a constant independent of such that
Moreover, the true error can be divided into a significant part and a less significant part as
Proof. Adding and subtracting to , we write Taking the norm and using the triangle inequality, we get
Next, we use the Cauchy‐Schwarz inequality and the inequality to write
Using the Cauchy‐Schwarz inequality and the inequality , we write
Using the inverse inequality (8), i.e., we obtain the estimate
5. A posteriori error estimation
In this section, we use the superconvergence results from the previous section to construct a residual‐based a posteriori error estimator which is computationally simple, efficient, and asymptotically exact. We will also prove its asymptotic exactness under mesh refinement. First, we present the weak finite element formulation to compute a posteriori error estimate for the nonlinear IVP (2).
In order to obtain a procedure for estimating the error , we multiply (2) by arbitrary smooth function and integrate over the element to obtain
Replacing by and choosing , we obtain
Using (32) and solving for , we obtain
Our error estimate procedure consists of approximating the true error on each element by the leading term as
where the coefficient of the leading term of the error, , is obtained from the coefficient defined in (56) by neglecting the terms and , i.e.,
We remark that our a posteriori error estimate is obtained by solving local problems with no boundary condition.
The global effectivity index, defined by , is an important criterion for evaluating the quality of an error estimator. The main results of this section are stated in the following theorem. In particular, we prove that the error estimate , in the norm, converges to the actual error . Moreover, we show that our a posterior error estimate is asymptotically exact by showing that the global effectivity index as .
Theorem 5.1. Suppose that the assumptions of Theorem 2 are satisfied. If where are defined in (57b), then
Thus, the post‐processed approximation yields superconvergent solution, i.e.,
Furthermore, then there exists a positive constant independent of such that
Finally, if there exists a constant independent of such that
then the global effectivity index in the norm converges to unity at
where we used the inequality . Summing over all elements yields
Using the Lipschitz condition (3) and applying the Cauchy‐Schwarz inequality yields
Applying the inequality , we obtain
Multiplying by and using (32), i.e., yields
where is a constant independent of the mesh size.
Summing over all elements and using , we arrive at
Next, we will prove (60). Using the reverse triangle inequality, we have
Remark 5.1. The previous theorem indicates that the computable quantity converges to at
Remark 5.2. The performance of an error estimator is typically measured by the global effectivity index which is defined as the ratio of the estimated error to the actual error . We say that the error estimator is asymptotically exact if as . The estimate (62) indicates that the global effectivity index in the norm converge to unity at
Remark 5.3. Our proofs are valid for any regular meshes and using piecewise polynomials of degree . If then (46) gives
Remark 5.4. The assumption (61), which is used to prove the convergence of to unity at
Remark 5. Our results readily extend to nonlinear systems of ODEs of the form
where , and . The DG method for this problem consists of finding such that: and ,
6. Application: adaptive mesh refinement (AMR)
A posteriori error estimates play an essential role in assessing the reliability of numerical solutions and in developing efficient adaptive algorithms. Adaptive methods based on a posteriori error estimates have become established procedures for computing efficient and accurate approximations to the solution of differential equations. The standard adaptive FEMs through local refinement can be written in the following loop
The local a posteriori errors estimator of Section 5 can be used to mark elements for refinement.
Next, we present a simple DG adaptive algorithm based on the local a posteriori error estimator proposed in the previous section. The adaptive algorithm that we propose has the following steps:
Select a tolerance and a maximum bound on the number of interval (say ). Put .
Construct an initial coarse mesh with nodes. For simplicity, we start with a uniform mesh having elements.
While and do
Solve the DG scheme to obtain the solution as described in Section 2.
For all elements
Choose a parameter . If the estimated global error then stop and accept the DG solution on the element .
Otherwise, reject the DG solution on and divide the element into two uniform elements by adding the coordinate of the midpoint of to the list of nodes.
Remark 6.1. There are many possibilities for selecting the elements to be refined given the local error indicator . In the above algorithm, we used the most popular fixed‐rate strategy which consists of refining the element if , where is a parameter provided by the user. Note that the choice gives uniform refinement, while gives no refinement. Also, there are other stopping criteria that may be used to stop the adaptive algorithm.
7. Computational results
In this section, we present several numerical examples to (i) validate our superconvergence results and the global convergence of the residual‐based a posteriori error estimates, and (ii) test the above local adaptive mesh refinement procedure that makes use of our local a posteriori error estimate.
Example 7.1. The test problem we consider is the following nonlinear IVP
Clearly, the exact solution is We use uniform meshes obtained by subdividing the computational domain into intervals with 5, 10, 20, 30, 40, 50. This example is tested by using polynomials with . Figure 1 shows the errors and with log‐log scale as well as their orders of convergence. These results indicate that
Next, we compute the maximum error at the shifted roots of the degree right Radau polynomial on each element and then take the maximum over all elements. For simplicity, we use to denote , where are the roots of . Similarly, we compute the true error at the downwind point of each element and then we denote to be the maximum over all elements , i.e., . In Figure 2, we present the errors , and their orders of convergence. We observe that
Next, we use (57a) and (57b) to compute the a posteriori error estimate for the DG solution. The global errors and their orders of convergence, using the spaces with , are shown in Figure 3. We observe that
In Table 1, we present the actual L2 errors and the global effectivity indices. These results demonstrate that the proposed a posteriori error estimates is asymptotically exact.
In Figure 4, we show the errors and . We see that
Example 7.2. In this example we test our error estimation procedure presented in Section 6 on adaptively refined meshes. We consider the following model problem
where the exact solution is simply . We apply our adaptive algorithm using (unstable), (stable), and (stiff). The DG solutions and the sequence of meshes obtained by applying our adaptive algorithm with for are shown in Figures 5–7 for , , and , respectively. As can be expected, the adaptive algorithm refines in the vicinity of the endpoint with coarser meshes for increasing polynomial degree . Furthermore, we observe that, when is closer to 0, we get more uniform refinement near the portion with high approximation error. When is near 1, we get less uniform refinement near the portion with high approximation error. We also observed that, both for and , the optimal convergence rates are achieved asymptotically and that the global effectivity indices converge to unity with increasing polynomial degree . Furthermore, we tested our adaptive algorithm on other problems and observed similar conclusions. These results are not included to save space.
8. Concluding remarks
In this chapter, we presented a detailed analysis of the original discontinuous Galerkin (DG) finite element method for the approximation of initial‐value problems (IVPs) for nonlinear ordinary differential equations (ODEs). We proved several optimal error estimates and superconvergence results. In particular, we showed that the DG solution converges to the true solution with order , when the space of piecewise polynomials of degree is used. We further proved the th superconvergence rate at the downwind points. Moreover, we proved that the DG solution is