The Discontinuous Galerkin Finite Element Method for Ordinary Differential Equations

Abstract

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.

1. Introduction

In this chapter, we introduce and analyze the discontinuous Galerkin (DG) method applied to the following first‐order initial‐value problem (IVP)

$\frac{d\overrightarrow{u}}{dt}=\overrightarrow{f}(t,\overrightarrow{u}),t\in [0,T],\overrightarrow{u}(0)={\overrightarrow{u}}_0,$E1

where $\overrightarrow{u}:[0,T]\to {ℝ}^n$, ${\overrightarrow{u}}_0\in {ℝ}^n,$ and $\overrightarrow{f}:[0,T]\times {ℝ}^n\to {ℝ}^n$. 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 $\overrightarrow{f}\in C^1([0,T]\times {ℝ}^n)$ is sufficient to guarantee the existence and uniqueness of the solution to (1). We note that a general $n$th‐order IVP of the form $y^{(n)}=g(t,y,y^{\prime },\dots ,y^{(n-1)})$ with initial conditions $y(0)=a_0$, $y^{\prime }(0)=a_1,\dots ,y^{(n-1)}(0)=a_{n-1}$ can be converted into a system of equations in the form (1), where $\overrightarrow{u}=[y,y^{\prime },\dots ,y^{(n-1)}{]}^t$, $\overrightarrow{f}(t,\overrightarrow{u})=[u_2,u_3,\dots ,u_n,g(t,u_1,\dots ,u_n{)]}^t$, and ${\overrightarrow{u}}_0=[a_0,a_1,\dots ,a_{n-1}{]}^t$.

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 [1] 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 $h$ refinement (mesh refinement and coarsening) and the $p$ 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 [3]. 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 [17] 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. [24] and Shu [25] 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[23].

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 [9] 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

(h^p) for general triangulations and of

(h^p+1) for Cartesian grids. Delfour et al. [7] investigated a class of Galerkin methods which lead to a family of one‐step schemes generating approximations up to order $2p+2$ for the solution of an ODE, when polynomials of degree $p$ are used. In their proposed method, the numerical solution $u_h$ at the discontinuity point $t_n$ is defined as an average across the jump, i.e., ${\alpha }_n{u}_h(t_n^{-})+(1-{\alpha }_n)u_h(t_n^{+})$. By choosing special values of ${\alpha }_n$, one can obtain the original DG scheme of LeSaint and Raviart [9] and Euler's explicit, improved, and implicit schemes. Delfour and Dubeau [27] introduced a family of discontinuous piecewise polynomial approximation schemes. They presented a more general framework of one‐step methods such as implicit Runge‐Kutta and Crank‐Nicholson schemes, multistep methods such as Adams‐Bashforth and Adams‐Moulton schemes, and hybrid methods. Later, Johnson [8] proved new optimal a priori error estimates for a class of implicit one‐step methods for stiff ODEs obtained by using the discontinuous Galerkin method with piecewise polynomials of degree zero and one. Johnson and Pitkaränta [29] proved a rate of convergence of

(h^p+1/2) for general triangulations and Peterson [19] confirmed this rate to be optimal. Richter [30] obtained the optimal rate of convergence

(h^p+1) for some structured two‐dimensional non‐Cartesian grids. We also would like to mention the work of Estep [28], where the author outlined a rigorous theory of global error control for the approximation of the IVP (1). In [6], Adjerid et al. showed that the DG solution of one‐dimensional hyperbolic problems exhibit an

(h^p+2) superconvergence rate at the roots of the right Radau polynomial of degree $p+1$. Furthermore, they obtained a $(2p+1)$th order superconvergence rate of the DG approximation at the downwind point of each element. They performed a local error analysis and showed that the local error on each element is proportional to a Radau polynomial. They further constructed implicit residual‐based a posteriori error estimates but they did not prove their asymptotic exactness. In 2010, Deng and Xiong [31] investigated a DG method with interpolated coefficients for the IVP (1). They proved pointwise superconvergence results at Radau points. More recently, the author [12, 15, 26, 32–39] investigated the global convergence of the several residual‐based a posteriori DG and local DG (LDG) error estimates for a variety of linear and nonlinear problems.

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)

$u^{\prime }=f(t,u),t\in [0,T],u(0)=u_0,$E2

where $f(t,u):[\mathrm{0,}T]\times ℝ\to ℝ$ is a sufficiently smooth function with respect to the variables $t$ and $u$. More precisely, we assume that $|f_u(t,u)|\le M_1$ on the set $D=[\mathrm{0,}T]\times ℝ\subset {ℝ}^2$, where $M_1$ is a positive constant. We note that the assumption $|f_u(t,u)|\le M_1$ is sufficient to ensure that $f(t,u)$ satisfies a Lipschitz condition in $u$ on the convex set $D$ with Lipschitz constant $M_1$

$|f(t,u)-f(t,v)|\le M_1|u-v|,\text{for any (t,u) and (t,v)}\in \text{D}.$E3

Next, we introduce the DG method for the model problem (2). Let $0=t_0<t_1<\cdots <t_N=T$ be a partition of the interval $\Omega =[\mathrm{0,}T]$. We denote the mesh by $I_j=[t_{j-1},t_j],j=1,\dots ,N$. We denote the length of $I_j$ by $h_j=t_j-t_{j-1}$. We also denote $h=\underset{1\le j\le N}{\max }{h}_j$ and $h_{\text{min}}=\underset{1\le j\le N}{\min }{h}_j$ as the length of the largest and smallest subinterval, respectively. Here, we consider regular meshes, that is, $\frac{h}{h_{\text{min}}}\le \lambda$, where $\lambda \ge 1$ is a constant (independent of $h$) during mesh refinement. If $\lambda =1$, then the mesh is uniformly distributed. In this case, the nodes and mesh size are defined by $t_j=jh,j=0,1,\dots ,N$, $h=T/N.$

Throughout this work, we define $v(t_j^{-})$ and $v(t_j^{+})$ to be the left limit and the right limit of the function $v$ at the discontinuity point $t_j$, $\text{i}\text{.e}.,$ $v(t_j^{-})=\underset{s\to 0^{-}}{{\displaystyle \lim }}v(t_j+s)$ and $v(t_j^{+})=\underset{s\to 0^{+}}{{\displaystyle \lim }}v(t_j+s)$. To simplify the notation, we denote by $[v](t_j)=v(t_j^{+})-v(t_j^{-})$ the jump of $v$ at the point $t_j$.

If we multiply (2) by an arbitrary test function $v$, integrate over the interval $I_j$, and integrate by parts, we get the DG weak formulation

${\displaystyle {\int }_{I_j}}{v}^{\prime }udt+{\displaystyle {\int }_{I_j}}f(t,u)vdt-u(t_j)v(t_j)+u(t_{j-1})v(t_{j-1})=0.$E4

We denote by $V_h^p$ the finite element space of polynomials of degree at most $p$ in each interval $I_j$, $\text{i}\text{.e}\text{.}$,

$V_h^p=\{v:v{|}_{I_j}\in P^p(I_j),j=1,\dots ,N\},$E101

where $P^p(I_j)$ denotes the set of all polynomials of degree no more than $p$ on $I_j$. We would like to emphasize that polynomials in $V_h^p$ are allowed to have discontinuities at the nodes $t_j$.

Replacing the exact solution $u(t)$ by a piecewise polynomial $u_h(t)\in V_h^p$ and choosing $v\in V_h^p$, we obtain the DG scheme: Find $u_h\in V_h^p$ such that $\forall v\in V_h^p$ and $j=1,\dots ,N$,

${\displaystyle {\int }_{I_j}}{v}^{\prime }{u}_{h}dt+{\displaystyle {\int }_{I_j}}f(t,u_h)vdt-{\widehat{u}}_h(t_j)v(t_j^{-})+{\widehat{u}}_h(t_{j-1})v(t_{j-1}^{+})=\mathrm{0,}$E5a

where ${\widehat{u}}_h(t_j)$ is the so‐called numerical flux which is nothing but the discrete approximation $u$ at the node $t=t_j$. We remark that $u_h$ is not necessarily continuous at the nodes.

To complete the definition of the DG scheme, we still need to define ${\widehat{u}}_h$ on the boundaries of $I_j$. Since for IVPs, information travel from the past into the future, it is reasonable to take ${\widehat{u}}_h$ as the classical upwind flux

${\widehat{u}}_h(t_0)=u_0,\text{and}{\widehat{u}}_h(t_j)=u_h(t_j^{-}),j=1,\dots ,N.$E5b

2.1. Implementation

The DG solution $u_h(t)$ can be efficiently obtained in the following order: first, we compute $u_h(t)$ in the first element $I_1$ using (5a) and (5b) with $j=1$ since $u_h(t_0^{-})=u_0$ is known. Then, we can find $u_h(t)$ in $I_2$ since $u_h(t)$ in $I_1$ is already available. We can repeat the same process to compute $u_h(t)$ in $I_3,\dots ,I_N$. More specifically, $u_h(t)$ can be obtained locally for each $I_j$ using the following two steps: (i) express $u_h(t)$ as a linear combination of orthogonal basis $L_{i,j}(t),i=0,\dots ,p$, where $L_{i,j}$ is the $i$th degree Legendre polynomial on $I_j$, i.e., $u_h(t)={\displaystyle {\sum }_{i=0}^p}{c}_{i,j}{L}_{i,j}(t),t\in I_j$, where ${\{L_{i,j}(t)\}}_{i=0}^{i=p}$ is a local basis of $P^p(I_j)$, and (ii) choose the test functions $v=L_{k,j}(t),k=0,\dots ,p$. Thus, on each $I_j$, we get a $(p+1)\times (p+1)$ system of nonlinear algebraic equations, which can be solved for the unknown coefficients $c_{\mathrm{0,}j},\dots ,c_{p,j}$ using, e.g., Newton's method for nonlinear systems. Once we obtain the DG solution on all elements $I_j,j=1,\dots ,N$, we get the DG solution which is a piecewise discontinuous polynomial of degree $\le p$. 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., $f(t,u)=\lambda u$. Taking $v=u_h$ in the discrete weak formulation (5a), we get

$\frac{u_h^2(t_j^{-})}{2}+\frac{u_h^2(t_{j-1}^{+})}{2}-u_h(t_{j-1}^{-})u_h(t_{j-1}^{+})=\lambda {\displaystyle {\int }_{I_j}}{u}_h^{2}dt,$E201

which is equivalent to

$\frac{u_h^2(t_j^{-})}{2}-\frac{u_h^2(t_{j-1}^{-})}{2}+\frac{1}{2}{\left(u_h(t_{j-1}^{-})-u_h(t_{j-1}^{+})\right)}^2=\lambda {\displaystyle {\int }_{I_j}}{u}_h^{2}dt.$E202

Summing over all elements, we get the equality

$\frac{u_h^2(T^{-})}{2}-\frac{u_0^2}{2}+\frac{1}{2}{\displaystyle \sum _{j=1}^N}{\left(u_h(t_{j-1}^{-})-u_h(t_{j-1}^{+})\right)}^2=\lambda {\displaystyle {\int }_{\Omega }}{u}_h^{2}dt.$E203

Consequently, $\frac{u_h^2(T^{-})}{2}-\frac{u_0^2}{2}\le \lambda {\displaystyle {\int }_{\Omega }}{u}_h^{2}dt,$ which gives the stability result $u_h^2(T^{-})\le u_0^2$ provided that $\lambda \le 0$.

3. A priori error analysis

We begin by defining some norms that will be used throughout this work. We define the $L^2$ inner product of two integrable functions, $u$ and $v$, on the interval $I_j=[t_{j-1},t_j]$ as ${(u,v)}_{I_j}={\displaystyle {\int }_{I_j}}u(t)v(t)dt$. Denote ${\Vert u\Vert }_{\mathrm{0,}{I}_j}=(u,u{)}_{I_j}^{1/2}$ to be the standard $L^2$ norm of $u$ on $I_j$. Moreover, the standard $L^{\infty }$ norm of $u$ on $I_j$ is defined by ${\Vert u\Vert }_{\infty ,I_j}=\underset{t\in I_j}{\sup }|u(t)|.$ Let $H^s(I_j)$, where $s=0,1,\dots$, denote the standard Sobolev space of square integrable functions on $I_j$ with all derivatives $u^{(k)},k=0,1,\dots ,s$ being square integrable on $I_j$, i.e., $H^s(I_j)=\left\{u:{\displaystyle {\int }_{I_j}}|u^{(k)}(t{)|}^{2}dt<\infty ,0\le k\le s\right\},$ and equipped with the norm ${\Vert u\Vert }_{s,I_j}={\left({\displaystyle {\sum }_{k=0}^s}{\Vert u^{(k)}\Vert }_{\mathrm{0,}{I}_j}^2\right)}^{1/2}.$ The $H^s(I_j)$ seminorm of a function $u$ on $I_j$ is given by ${\left|u\right|}_{s,I_j}={\Vert u^{(s)}\Vert }_{\mathrm{0,}{I}_j}$. We also define the norms on the whole computational domain $\Omega$ as follows:

${\Vert u\Vert }_{\mathrm{0,}\Omega }={\left({\displaystyle \sum _{j=1}^N}{\Vert u\Vert }_{\mathrm{0,}{I}_j}^2\right)}^{1/2},{\Vert u\Vert }_{\infty ,\Omega }=\underset{1\le j\le N}{{\displaystyle \max }}{\Vert u\Vert }_{\infty ,I_j},{\Vert u\Vert }_{s,\Omega }={\left({\displaystyle \sum _{j=1}^N}{\Vert u\Vert }_{s,I_j}^2\right)}^{1/2}.$E204

The seminorm on the whole computational domain $\Omega$ is defined as ${\left|u\right|}_{s,\Omega }={\left({\displaystyle {\sum }_{j=1}^N}{\left|u\right|}_{s,I_j}^2\right)}^{1/2}.$ We note that if $u\in H^s(\Omega )$, $s=1,2,\dots$, then the norm ${\Vert u\Vert }_{s,\Omega }$ on the whole computational domain is the standard Sobolev norm ${\left({\displaystyle {\sum }_{k=0}^s}{\Vert u^{(k)}\Vert }_{\mathrm{0,}\Omega }^2\right)}^{1/2}.$ For convenience, we use ${\Vert u\Vert }_{I_j}$ and $\Vert u\Vert$ to denote ${\Vert u\Vert }_{\mathrm{0,}{I}_j}$ and ${\Vert u\Vert }_{\mathrm{0,}\Omega }$, respectively.

For $p\ge 1$, we consider two special projection operators, $P_h^{\pm }$, which are defined as follows: For a smooth function $u$, the restrictions of $P_h^{+}u$ and $P_h^{-}u$ to $I_j$ are polynomials in $P^p(I_j)$ satisfying

${\displaystyle {\int }_{I_j}}(P_h^{-}u-u)vdt=0,\forall v\in P^{p-1}(I_j),\text{and}(P_h^{-}u-u)(t_j^{-})=\mathrm{0,}$E6a

${\displaystyle {\int }_{I_j}}(P_h^{+}u-u)vdt=0,\forall v\in P^{p-1}(I_j),\text{and}(P_h^{+}u-u)(t_{j-1}^{+})=0.$E6b

These two particular Gauss‐Radau projections are very important in the proofs of optimal $L^2$ error estimates and superconvergence results. We note that the special projections $P_h^{\pm }u$ 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 [22].

For the projections mentioned above, it is easy to show that for any $u\in H^{p+1}(I_j)$ with $j=1,\dots ,N$, there exists a constant $C$ independent of the mesh size $h$ such that (see, e.g., [40])

${\Vert u-P_h^{\pm }u\Vert }_{I_j}\le C{h}_j^{p+1}{\left|u\right|}_{p+\mathrm{1,}{I}_j},{\Vert (u-P_h^{\pm }u{)}^{\prime }\Vert }_{I_j}\le C{h}_j^p{\left|u\right|}_{p,I_j}.$E7

Moreover, we recall the inverse properties of the finite element space $V_h^p$ that will be used in our error analysis: For any $v_h\in V_h^p$, there exists a positive constant $C$ independent of $v_h$ and $h$, such that, $\forall j=1,\dots ,N$,

${\Vert v_h^{(k)}\Vert }_{I_j}\le C{h}_j^{-k}\left|\right|v_h{\left|\right|}_{I_j},k\ge \mathrm{1,}\left|v_h(t_{j-1}^{+})\right|+\left|v_h(t_j^{-})\right|\le C{h}_j^{-1/2}\left|\right|v_h{\left|\right|}_{I_j}.$E8

From now on, the notation $C,C_1,C_2,$ etc. will be used to denote generic positive constants independent of $h$, 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 $e=u-u_h$ to be the error between the exact solution of (2) and the DG solution defined in (5a) and (5b), $\epsilon =u-P_h^{-}u$ to be the projection error, and $\overline{e}=P_h^{-}u-u_h$ to be the error between the projection of the exact solution $P_h^{-}u$ and the DG solution $u_h$. We observe that the actual error can be written as $e=(u-P_h^{-}u)+(P_h^{-}u-u_h)=\epsilon +\overline{e}.$

Now, we are ready to prove our optimal error estimates for $e$ in the $L^2$ and $H^1$ norms.

Theorem 3.1. Suppose that the exact solution of (2) is sufficiently smooth with bounded derivatives, i.e., ${\Vert u\Vert }_{p+\mathrm{1,}\Omega }$ is bounded. We also assume that $|f_u(t,u)|\le M_1$ on $D=[\mathrm{0,}T]\times ℝ$. Let $p\ge 0$ and $u_h$ be the DG solution of (5a) and (5b), then, for sufficiently small $h$, there exists a positive constant $C$ independent of $h$ such that,

$\Vert e\Vert \le C{h}^{p+1},$E9

${\displaystyle \sum _{j=1}^N}{\Vert e^{\prime }\Vert }_{I_j}^2\le C{h}^{2p},{\Vert e\Vert }_{\mathrm{1,}\Omega }\le C{h}^p.$E10

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 $v\in V_h^p$ and using the numerical flux (5b), we obtain the following error equation: $\forall v\in V_h^p$,

${\displaystyle {\int }_{I_j}}{v}^{\prime }edt+{\displaystyle {\int }_{I_j}}(f(t,u)-f(t,u_h))vdt+e(t_{j-1}^{-})v(t_{j-1}^{+})-e(t_j^{-})v(t_j^{-})=0.$E11

By integration by parts, we get

${\displaystyle {\int }_{I_j}}{e}^{\prime }vdt-{\displaystyle {\int }_{I_j}}\left(f(t,u)-f(t,u_h)\right)vdt+[e](t_{j-1})v(t_{j-1}^{+})=0.$E12

Applying Taylor's series with integral remainder in the variable $u$ and using the relation $u-u_h=e$, we write

$f(t,u)-f(t,u_h)=\theta (u-u_h)=\theta e,\text{where}\theta ={\displaystyle {\int }_0^1}{f}_u(t,u+s(u_h-u))ds={\displaystyle {\int }_0^1}{f}_u(t,u-se)ds.$E13

Substituting (13) into (12), we arrive at

${\displaystyle {\int }_{I_j}}\left(e^{\prime }-\theta e\right)vdt+[e](t_{j-1})v(t_{j-1}^{+})=\mathrm{0,}\forall v\in V_h^p.$E14

To simplify the notation, we introduce the bilinear operator ${\mathcal{A}}_j(e;V)$ as

${\mathcal{A}}_j(e;V)={\displaystyle {\int }_{I_j}}(e^{\prime }-\theta e)Vdt+[e](t_{j-1})V(t_{j-1}^{+}).$E15

Thus, we can write (14) as

${\mathcal{A}}_j(e;v)=\mathrm{0,}\forall v\in V_h^p.$E16

A direct calculation from integration by parts yields

${\mathcal{A}}_j(e;V)={\displaystyle {\int }_{I_j}}(-V^{\prime }-\theta V)edt+e(t_j^{-})V(t_j^{-})-e(t_{j-1}^{-})V(t_{j-1}^{+}).$E17

On the other hand, if we add and subtract $P_h^{+}V$ to $V$ then we can write (15) as

${\mathcal{A}}_j(e;V)={\mathcal{A}}_j(e;V-P_h^{+}V)+{\mathcal{A}}_j(e;P_h^{+}V).$E18

Combining (18) and (16) with $v=P_h^{+}V\in P^p(I_j)$ and applying the property of the projection $P_h^{+}$, i.e., $(V-P_h^{+}V)(t_{j-1}^{+})=0$, we obtain

${\mathcal{A}}_j(e;V)={\displaystyle {\int }_{I_j}}(e^{\prime }-\theta e)(V-P_h^{+}V)dt+[e](t_{j-1})(V-P_h^{+}V)(t_{j-1}^{+})={\displaystyle {\int }_{I_j}}(e^{\prime }-\theta e)(V-P_h^{+}V)dt.$E19

If $v$ is a polynomial of degree at most $p$ then $v^{\prime }$ is a polynomial of degree at most $p-1$. Therefore, by the property of the projection $P_h^{+}$, we immediately see

${\displaystyle {\int }_{I_j}}{v}^{\prime }(V-P_h^{+}V)dt=0,\forall v\in P^p(I_j).$E20

Substituting the relation $e=\epsilon +\overline{e}$ into (19) and invoking (20) with $v=\overline{e}$, we get

${\mathcal{A}}_j(e;V)={\displaystyle {\int }_{I_j}}({\epsilon }^{\prime }-\theta e)(V-P_h^{+}V)dt+{\displaystyle {\int }_{I_j}}\overline{e}\text{'}(V-P_h^{+}V)dt={\displaystyle {\int }_{I_j}}({\epsilon }^{\prime }-\theta e)(V-P_h^{+}V)dt.$E21

Now, we are ready to prove the theorem. We construct the following auxiliary problem: find $\phi$ such that

$-{\phi }^{\prime }-\theta \phi =e,t\in [0,T]\text{subject to}\phi (T)=0.$E22

where $\theta =\theta (t)={\displaystyle {\int }_0^1}{f}_u(t,u(t)-se(t))ds$. Clearly, the exact solution to (22) is given by the explicit formula

$\phi (t)=\frac{1}{\Theta (t)}{\displaystyle {\int }_t^T}\Theta (y)e(y)dy,\text{where}\Theta (t)=\exp \left(-{\displaystyle {\int }_t^T}\theta (s)ds\right).$E23

Next, we prove some regular estimates which will be needed in our error analysis. Using the assumption $|f_u(t,u)|\le M_1$, we see that $\theta (t),t\in [0,T]$ is bounded by $M_1$

$|\theta (t)|\le {\displaystyle {\int }_0^1}|f_u(t,u(t)-se(t))|ds\le {\displaystyle {\int }_0^1}{M}_{1}ds=M_1,\forall t\in [0,T].$E24a

Using the definition of $\Theta$ and the estimate (24a), we have

$0\le \Theta (t)\le \exp \left({\displaystyle {\int }_0^T}|\theta (s)|ds\right)\le \exp \left({\displaystyle {\int }_0^T}{M}_{1}ds\right)=\exp \left(M_{1}T\right)=C_1.$E24b

Similarly, we can easily estimate $\frac{1}{\Theta (t)}$ as follows

$0\le \frac{1}{\Theta (t)}=\exp \left({\displaystyle {\int }_t^T}\theta (s)ds\right)\le \exp \left({\displaystyle {\int }_0^T}|\theta (s)|ds\right)\le \exp \left({\displaystyle {\int }_0^T}{M}_{1}ds\right)=\exp \left(M_{1}T\right)=C_1.$E24c

Applying the estimates (24b), (24c), and the Cauchy‐Schwarz inequality, we get

$\left|\phi (t)\right|\le \frac{1}{\Theta (t)}{\displaystyle {\int }_t^T}\Theta (y)|e(y)|dy\le C_1{\displaystyle {\int }_t^T}{C}_1|e(y)|dy\le C_1^2{\displaystyle {\int }_0^T}|e(y)|dy\le C_1^2{T}^{1/2}\Vert e\Vert ,t\in [0,T].$E205

Squaring both sides and intergrading over $\Omega$ yields

${\Vert \phi \Vert }^2\le C_1^4{T}^2{\Vert e\Vert }^2=C_2{\Vert e\Vert }^2.$E25a

We also need to obtain an estimate of ${\left|\phi \right|}_{\mathrm{1,}\Omega }$. Using (22) and (24a) gives

$\left|{\phi }^{\prime }\right|=|\theta \phi +e|\le M_1\left|\phi \right|+\left|e\right|,t\in [0,T].$E206

Squaring both sides, applying the inequality ${(a+b)}^2\le 2a^2+2b^2$, integrating over the computational domain $\Omega$, and using (25a), we get

${\left|\phi \right|}_{\mathrm{1,}\Omega }^2={\displaystyle {\int }_0^T}|{\phi }^{\prime }{|}^{2}dt\le 2(M_1^2{\Vert \phi \Vert }^2+{\Vert e\Vert }^2)\le 2(M_1^2{C}_2+1){\Vert e\Vert }^2\le C_3{\Vert e\Vert }^2.$E25b

Applying the projection result and the estimate (3.20b) yields

$\Vert \phi -P_h^{+}\phi \Vert \le C_{4}h{\left|\phi \right|}_{\mathrm{1,}\Omega }\le C_{5}h\Vert e\Vert .$E25c

Now, we are ready to show (9). Using (17) with $V=\phi$ and (22), we obtain

${\mathcal{A}}_j(e;\phi )={\displaystyle {\int }_{I_j}}(-{\phi }^{\prime }-\theta \phi )edt-e(t_{j-1}^{-})\phi (t_{j-1})+e(t_j^{-})\phi (t_j)={\displaystyle {\int }_{I_j}}{e}^{2}dt-e(t_{j-1}^{-})\phi (t_{j-1})+e(t_j^{-})\phi (t_j).$E207

Summing over the elements and using the fact that $\phi (T)=e(t_0^{-})=0$ yields

${\displaystyle \sum _{j=1}^N}{\mathcal{A}}_j(e;\phi )={\Vert e\Vert }^2-e(t_0^{-})\phi (t_0^{-})+e(T^{-})\phi (T^{-})={\Vert e\Vert }^2.$E26

On the other hand, if we choose $V=\phi$ in (21) then we get

${\mathcal{A}}_j(e;\phi )={\displaystyle {\int }_{I_j}}({\epsilon }^{\prime }-\theta e)(\phi -P_h^{+}\phi )dt.$E27

Summing over all elements and applying the Cauchy‐Schwarz inequality, we get

${\displaystyle \sum _{j=1}^N}{\mathcal{A}}_j(e;\phi )\le (\Vert {\epsilon }^{\prime }\Vert +M_1\Vert e\Vert )\Vert \phi -P_h^{+}\phi \Vert \mathrm{.}$E208

Using the estimate (25c), we deduce that

${\displaystyle \sum _{j=1}^N}{\mathcal{A}}_j(e;\phi )\le (C_0{h}^p{\left|u\right|}_{p+\mathrm{1,}\Omega }+M_1\Vert e\Vert )C_{1}h\Vert e\Vert \le C(h^{p+1}+h\Vert e\Vert )\Vert e\Vert .$E28

Combining (26) and (28), we conclude that

$\Vert e\Vert \le C{h}^{p+1}+Ch\Vert e\Vert .$E29

Thus, $(1-Ch)\Vert e\Vert \le C{h}^{p+1}$, where $C$ is a positive constant independent of $h$. Therefore, for sufficiently small $h$, e.g., $h\le \frac{1}{2C}$, we obtain $\frac{1}{2}\Vert e\Vert \le (1-Ch)\Vert e\Vert \le C{h}^{p+1}$, which yields $\Vert e\Vert \le 2C{h}^{p+1}$ for $h$ small. Thus, we completed the proof of (9).

To show (10), we use $e=\overline{e}+\epsilon$, the classical inverse inequality (8), the estimate (9), and the projection result (7) to obtain

$\Vert e^{\prime }\Vert =\Vert \overline{e}\text{'}+{\epsilon }^{\prime }\Vert \le \Vert \overline{e}\text{'}\Vert +\Vert {\epsilon }^{\prime }\Vert \le C_1{h}^{-1}\Vert \overline{e}\Vert +C_2{h}^p\le C_1{h}^{-1}\left(\Vert e\Vert +\Vert \epsilon \Vert \right)+C_2{h}^p\le C_3{h}^p+C_2{h}^p\le C{h}^p.$E209

We note that ${\Vert e\Vert }_{\mathrm{1,}\Omega }^2={\Vert e\Vert }^2+{\Vert e^{\prime }\Vert }^2.$ Applying (9) and the estimate $\Vert e^{\prime }\Vert \le C{h}^p$ yields ${\Vert e\Vert }_{\mathrm{1,}\Omega }^2\le C_1{h}^{2p+2}+C_2{h}^{2p}=$

(h^2p), which completes the proof of the theorem.

4. Superconvergence error analysis

In this section, we study the superconvergence properties of the DG method. We first show a $(2p+1)$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 $P_h^{-}u$ at

(h^p+2). This result allows us to prove that the leading term of the DG error is proportional to the $(p+1)$ degree right Radau polynomial.

First, we define some special polynomials. The $p$th degree Legendre polynomial can be defined by Rodrigues formula [41]

${\tilde{L}}_p(\xi )=\frac{1}{2^{p}p!}\frac{d^p}{d{\xi }^p}\left({({\xi }^2-1)}^p\right),-1\le \xi \le 1.$E210

It satisfies the following important properties: ${\tilde{L}}_p(1)=1$, ${\tilde{L}}_p(-1)=(-{1)}^p$, and the orthogonality relation

${\displaystyle {\int }_{-1}^1}{\tilde{L}}_p(\xi ){\tilde{L}}_q(\xi )d\xi =\frac{2}{2p+1}{\delta }_{pq},\text{where }{\delta }_{pq}\text{ is the Kronecker symbol}.$E30

One can easily write the $(p+1)$ degree Legendre polynomial on $[-\mathrm{1,1}]$ as

${\tilde{L}}_{p+1}(\xi )=\frac{(2p+2)!}{2^{p+1}{((p+1)!)}^2}{\xi }^{p+1}+{\tilde{q}}_p(\xi ),\text{where}{\tilde{q}}_p\in P^p([-\mathrm{1,1}]).$E211

The $(p+1)$ degree right Radau polynomial on $[-\mathrm{1,1}]$ is defined as ${\tilde{R}}_{p+1}(\xi )={\tilde{L}}_{p+1}(\xi )-{\tilde{L}}_p(\xi )$. It has $p+1$ real distinct roots, $-1<{\xi }_0<\cdots <{\xi }_p=1$.

Mapping $I_j$ into the reference element $[-\mathrm{1,1}]$ by the linear transformation $t=\frac{t_j+t_{j-1}}{2}+\frac{h_j}{2}\xi$, we obtain the shifted Legendre and Radau polynomials on $I_j$:

$L_{p+\mathrm{1,}j}(t)={\tilde{L}}_{p+1}\left(\frac{2t-t_j-t_{j-1}}{h_j}\right),R_{p+\mathrm{1,}j}(t)={\tilde{R}}_{p+1}\left(\frac{2t-t_j-t_{j-1}}{h_j}\right).$E301

Next, we define the monic Radau polynomial, ${\psi }_{p+\mathrm{1,}j}(t)$, on $I_j$ as

${\psi }_{p+\mathrm{1,}j}(t)=\frac{h_j^{p+1}{[(p+1)!]}^2}{(2p+2)!}{R}_{p+\mathrm{1,}j}(t)=c_p{h}_j^{p+1}{R}_{p+\mathrm{1,}j}(t),\text{where}{c}_p=\frac{{((p+1)!)}^2}{(2p+2)!}.$E31

Throughout this work the roots of $R_{p+\mathrm{1,}j}(t)$ are denoted by $t_{j,i}=\frac{t_j+t_{j-1}}{2}+\frac{h_j}{2}{\xi }_i,i=0,1,\dots ,p.$

In the next lemma, we recall the following results which will be needed in our error analysis [32].

Lemma 4.1. The polynomials $L_{p,j}$ and ${\psi }_{p+\mathrm{1,}j}$ satisfy the following properties

${\Vert L_{p,j}\Vert }_{I_j}^2=\frac{h_j}{2p+1},{\displaystyle {\int }_{I_j}}{{\psi }^{\prime }}_{p+\mathrm{1,}j}{\psi }_{p+\mathrm{1,}j}dt=-k_1{h}_j^{2p+2},{\Vert {\psi }_{p+\mathrm{1,}j}\Vert }_{I_j}^2=(2p+2)k_2{h}_j^{2p+3},$E32

where $k_1=2c_p^2$, $k_2=\frac{k_1}{(2p+1)(2p+3)}$, and $c_p=\frac{{((p+1)!)}^2}{(2p+2)!}$.

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 $f_u$ is sufficiently smooth with respect to $t$ and $u$ (for example, $h(t)=f_u(t,u(t))\in C^p([0,T])$ is enough.). Then there exists a positive constant $C$ such that

$\left|e(t_k^{-})\right|\le C{h}^{2p+1},k=1,\dots ,N,$E33

$\left|\overline{e}(t_k^{-})\right|\le C{h}^{2p+1},k=1,\dots ,N,$E34

$\Vert \overline{e^{\prime }}\Vert \le C{h}^{p+1},$E35

$\Vert \overline{e}\Vert \le C{h}^{p+2}.$E36

Proof. To prove (33), we proceed by the duality argument. Consider the following auxiliary problem:

$W^{\prime }+\theta W=0,t\in [0,t_k]\text{subject to}W(t_k)=\mathrm{1,}$E37

where $1\le k\le N$ and $\theta =\theta (t)={\displaystyle {\int }_0^1}{f}_u(t,u(t)-se(t))ds$. The exact solution of this problem is $W(t)=\exp \left({\displaystyle {\int }_t^{t_k}}\theta (s)ds\right)$, $t\in {\Omega }_k=[\mathrm{0,}{t}_k].$ Using the assumption $h(t)=f_u(t,u(t))\in C^p([0,T])$ and the estimate (24a), we can easily show that there exists a constant $C$ such that

${\Vert W\Vert }_{p+\mathrm{1,}{\Omega }_k}\le C.$E38

Using (17) and (37), we get

${\mathcal{A}}_j(e;W)={\displaystyle {\int }_{I_j}}(-W^{\prime }-\theta W)edt+-e(t_{j-1}^{-})W(t_{j-1})+e(t_j^{-})W(t_j)=-e(t_{j-1}^{-})W(t_{j-1})+e(t_j^{-})W(t_j).$E302

Summing over the elements $I_j,j=1,\dots ,k$, using $W(t_k)=1$, and the fact that $e(t_0^{-})=0$, we obtain

${\displaystyle \sum _{j=1}^k}{\mathcal{A}}_j(e;W)=-e(t_0^{-})W(t_0)+e(t_k^{-})W(t_k)=e(t_k^{-}).$E39

Now, taking $V=W$ in (21) yields

${\mathcal{A}}_j(e;W)={\displaystyle {\int }_{I_j}}({\epsilon }^{\prime }-\theta e)(W-P_h^{+}W)dt.$E303

Summing over all elements $I_j,j=1,\dots ,k$ with $k=1,\dots ,N$ and applying (39), we arrive at

$e(t_k^{-})={\displaystyle \sum _{j=1}^k}{\displaystyle {\int }_{I_j}}({\epsilon }^{\prime }-\theta e)(W-P_h^{+}W)dt.$E304

Using (24a) and applying the Cauchy‐Schwarz inequality, we obtain

$\left|e(t_k^{-})\right|\le ({\Vert {\epsilon }^{\prime }\Vert }_{\mathrm{0,}{\Omega }_k}+M_1{\Vert e\Vert }_{\mathrm{0,}{\Omega }_k}){\Vert W-P_h^{+}W\Vert }_{\mathrm{0,}{\Omega }_k}\le (\Vert {\epsilon }^{\prime }\Vert +M_1\Vert e\Vert ){\Vert W-P_h^{+}W\Vert }_{\mathrm{0,}{\Omega }_k}.$E305

Invoking the estimates (7), (9), and (38), we conclude that

$\left|e(t_k^{-})\right|\le (C_0{h}^p{\left|u\right|}_{p+\mathrm{1,}\Omega }+M_1{C}_1{h}^{p+1})C_2{h}^{p+1}{\left|W\right|}_{p+\mathrm{1,}{\Omega }_k}\le C(h^p+h^{p+1})h^{p+1}=\mathcal{O}(h^{2p+1}),$E40

for all $k=1,\dots ,N,$ which completes the proof of (33).

In order to prove (34), we use the relation $e=\overline{e}+\epsilon$, the property of the projection $P_h^{-}$, i.e., $\epsilon (t_k^{-})=0$, and the estimate (33) to get

$\left|\overline{e}(t_k^{-})\right|=\left|e(t_k^{-})-\epsilon (t_k^{-})\right|=\left|e(t_k^{-})\right|=\mathcal{O}(h^{2p+1}).$E306

Next, we will derive optimal error estimate for $\Vert \overline{e}\text{'}\Vert$. By the property of $P_h^{-}$, we have

${\displaystyle {\int }_{I_j}}\epsilon v^{\prime }dt=0,\forall v\in P^p(I_j),\text{and}\epsilon (t_j^{-})=\mathrm{0,}j=1,\dots ,N.$E41

Using the relation $e=\epsilon +\overline{e}$, applying (41) and (11) yields

${\displaystyle {\int }_{I_j}}{v}^{\prime }\overline{e}dt+{\displaystyle {\int }_{I_j}}(f(t,u)-f(t,u_h))vdt+\overline{e}(t_{j-1}^{-})v(t_{j-1}^{+})-\overline{e}(t_j^{-})v(t_j^{-})=0.$E307

By integration by parts on the first term, we obtain

${\displaystyle {\int }_{I_j}}\left(\overline{e}\text{'}-f(t,u)+f(t,u_h)\right)vdt+[\overline{e}](t_{j-1})v(t_{j-1}^{+})=0.$E42

Choosing $v(t)=\overline{e}\text{'}(t)-{(-1)}^p\overline{e}\text{'}(t_{j-1}^{+})L_{p,j}(t)\in P^p(I_j)$ in (42), we have, by the property ${\tilde{L}}_p(-1)=(-{1)}^p$ and the orthogonality relation (30), $v(t_{j-1}^{+})=0$ and

$\begin{array}{c}{\displaystyle {\int }_{I_j}}{(}^{\overline{e}}dt=(-1^{)}\overline{e}\text{'}(t_{j-1}^{+}){\displaystyle {\int }_{I_j}}{L}_{p,j}\overline{e}\text{'}dt+{\displaystyle {\int }_{I_j}}(f(t,u)-f(t,u_h))(\overline{e}\text{'}-{(}^{-}\overline{e}\text{'}(t_{j-1}^{+})L_{p,j})dt\\ ={\displaystyle {\int }_{I_j}}(f(t,u)-f(t,u_h))(\overline{e}\text{'}-{(}^{-}\overline{e}\text{'}(t_{j-1}^{+})L_{p,j})dt.\end{array}$E43

Using (3) and applying the Cauchy‐Schwarz inequality gives

$\begin{array}{c}{\Vert {\overline{e}}^{\prime }\Vert }_{I_j}^2\le {\displaystyle {\int }_{I_j}}\left|f(t\mathrm{,}u)-f(t,u_h)\right|\left(\left|{\overline{e}}^{\prime }\right|+\left|{\overline{e}}^{\prime }(t_{j-1}^{+})\right|\left|L_{p\mathrm{,}j}\right|\right)dt\le M_1{\displaystyle {\int }_{I_j}}\left|e\right|\left(\left|{\overline{e}}^{\prime }\right|+\left|{\overline{e}}^{\prime }(t_{j-1}^{+})\right|\left|L_{p\mathrm{,}j}\right|\right)dt\\ \le M_1{\Vert e\Vert }_{I_j}\left({\Vert {\overline{e}}^{\prime }\Vert }_{I_j}+\left|{\overline{e}}^{\prime }(t_{j-1}^{+})\right|{\Vert L_{p\mathrm{,}j}\Vert }_{I_j}\right)\mathrm{.}\end{array}$E44

Combining (44) with (8) and (32), we obtain

${\Vert {\overline{e}}^{\prime }\Vert }_{I_j}^2\le M_1{\Vert e\Vert }_{I_j}\left({\Vert {\overline{e}}^{\prime }\Vert }_{I_j}+\left(C_1{h}_j^{-1/2}{\Vert {\overline{e}}^{\prime }\Vert }_{I_j}\right)\left(\frac{h_j^{1/2}}{{(2p+1)}^{1/2}}\right)\right)\le C{\Vert e\Vert }_{I_j}{\Vert {\overline{e}}^{\prime }\Vert }_{I_j}.$E367

Consequently, ${\Vert {\overline{e}}^{\prime }\Vert }_{I_j}\le C{\Vert e\Vert }_{I_j}$. Taking the square of both sides, summing over all elements, and using (9), we conclude that

${\Vert \overline{e}\text{'}\Vert }^2\le C{\Vert e\Vert }^2\le C{h}^{2p+2}.$E45

Finally, we will estimate $\Vert \overline{e}\Vert$. Using the fundamental theorem of calculus, we write

$|\overline{e}(t)|=|\overline{e}(t_j^{-})+{\displaystyle {\int }_{t_j}^t}{\overline{e}}^{\prime }(s)ds|\le |\overline{e}(t_j^{-})|+{\displaystyle {\int }_{I_j}}|{\overline{e}}^{\prime }(s)|ds,\forall t\in I_j.$E308

Taking the square of both sides, applying the inequality ${(a+b)}^2\le 2a^2+2b^2$, and applying the Cauchy‐Schwartz inequality, we get

$|\overline{e}(t{)|}^2\le 2|\overline{e}(t_j^{-}{)|}^2+2{\left({\displaystyle {\int }_{I_j}}|{\overline{e}}^{\prime }(s)|ds\right)}^2\le 2|\overline{e}(t_j^{-}{)|}^2+2h_j{\displaystyle {\int }_{I_j}}|{\overline{e}}^{\prime }(s{)|}^{2}ds=2|\overline{e}(t_j^{-}{)|}^2+2h_j{\Vert {\overline{e}}^{\prime }\Vert }_{I_j}^2.$E309

Integrating this inequality with respect to $t$ and using the estimate (34), we get

${\Vert \overline{e}\Vert }_{I_j}^2\le 2h_j|\overline{e}(t_j^{-}{)|}^2+2h_j^2{\Vert {\overline{e}}^{\prime }\Vert }_{I_j}^2\le 2C{h}_j^{4p+3}+2h_j^2{\Vert {\overline{e}}^{\prime }\Vert }_{I_j}^2.$E310

Summing over all elements and using the estimate (35) and the fact that $h=\max h_j$, we obtain

${\Vert \overline{e}\Vert }^2\le C_1{h}^{4p+2}+2h^2{\Vert \overline{e}\text{'}\Vert }^2\le C_1{h}^{4p+2}+2C_2{h}^{2p+4}=\mathcal{O}(h^{2p+4}),$E46

where we used $4p+2\ge 2p+4$ for $p\ge 1$. This completes the proof of the theorem.

The previous theorem indicates that the DG solution $u_h$ is closer to $P_h^{-}u$ than to the exact solution $u$. Next, we apply the results of Theorem 2 to prove that the actual error $e$ can be split into a significant part, which is proportional to the $(p+1)$ degree right Radau polynomial, and a less significant part that converges at $$

(h^p+2) rate in the $L^2$ norm. Before we prove this result, we introduce two interpolation operators $\pi$ and $\widehat{\pi }$. The interpolation operator $\pi$ is defined as follows: For smooth $u=u(t)$, $\pi u{|}_{I_j}\in P^p(I_j)$ and interpolates $u$ at the roots of the $(p+1)$ degree right Radau polynomial shifted to $I_j$, i.e., at $t_{j,i},i=0,1,\dots ,p,$. The interpolation operator $\widehat{\pi }$ satisfies $\widehat{\pi }u{|}_{I_j}\in P^{p+1}(I_j)$ and $\widehat{\pi }u{|}_{I_j}$ interpolates $u$ at $t_{j,i},i=0,1,\dots ,p,$ and at an additional point ${\overline{t}}_j$ in $I_j$ with ${\overline{t}}_j\ne t_{j,i},i=0,1,\dots ,p$. The choice of the additional point is not important and can be chosen as ${\overline{t}}_j=t_{j-1}$.

Next, we recall the following results from [12] which will be needed in our analysis.

Lemma 4.2. If $u\in H^{p+2}(I_j)$, then interpolation error can be split as

$u-\pi u={\varphi }_j+{\gamma }_j,on{I}_j,$E47a

where

${\varphi }_j(t)={\alpha }_j{\psi }_{p+\mathrm{1,}j}(t),{\psi }_{p+\mathrm{1,}j}(t)={\displaystyle \prod _{i=0}^p}(t-t_{j,i}),{\gamma }_j=u-\widehat{\pi }u,$E47b

and ${\alpha }_j$ is the coefficient of $t^{p+1}$ in the $(p+1)$ degree polynomial $\widehat{\pi }u$. Furthermore,

${\Vert {\varphi }_j\Vert }_{s,I_j}\le C{h}_j^{p+1-s}{\Vert u\Vert }_{p+\mathrm{1,}{I}_j},s=0,\dots ,p,$E47c

${\Vert {\gamma }_j\Vert }_{s,I_j}\le C{h}_j^{p+2-s}{\Vert u\Vert }_{p+\mathrm{2,}{I}_j},s=0,\dots ,p+1.$E47d

Finally,

${\Vert \pi u-P_h^{-}u\Vert }_{I_j}\le C{h}_j^{p+2}{\Vert u\Vert }_{p+\mathrm{2,}{I}_j}.$E48

Proof. The proof of this lemma can be found in [12], 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 $C$ independent of $h$ such that

$\Vert u_h-\pi u\Vert \le C{h}^{p+2}.$E49

Moreover, the true error can be divided into a significant part and a less significant part as

$e(t)={\alpha }_j{\psi }_{p+\mathrm{1,}j}(t)+{\omega }_j(t),on{I}_j,$E50a

where

${\omega }_j={\gamma }_j+\pi u-u_h,$E50b

and

${\displaystyle \sum _{j=1}^N}{\Vert {\omega }_j\Vert }_{I_j}^2\le C{h}^{2(p+2)},{\displaystyle \sum _{j=1}^N}{\Vert {\omega }_{j^{\prime }}\Vert }_{I_j}^2\le C{h}^{2(p+1)}.$E50c

Proof. Adding and subtracting $P_h^{-}u$ to $u_h-\pi u$, we write $u_h-\pi u=u_h-P_h^{-}u+P_h^{-}u-\pi u=-\overline{e}+P_h^{-}u-\pi u.$ Taking the $L^2$ norm and using the triangle inequality, we get

$\Vert u_h-\pi u\Vert \le \Vert \overline{e}\Vert +\Vert P_h^{-}u-\pi u\Vert .$E311

Applying the estimates (36) and (48), we deduce (49). Next, adding and subtracting $\pi u$ to $e$, we write $e=u-\pi u+\pi u-u_h$. Moreover, one can split the interpolation error $u-\pi u$ on $I_j$ as in (47a) to get

$e={\varphi }_j+{\gamma }_j+\pi u-u_h={\varphi }_j+{\omega }_j,\text{where}{\omega }_j={\gamma }_j+\pi u-u_h.$E51

Next, we use the Cauchy‐Schwarz inequality and the inequality $\left|ab\right|\le \frac{1}{2}(a^2+b^2)$ to write

$\begin{array}{c}{\Vert {\omega }_j\Vert }_{I_j}^2={\left({\gamma }_j+\pi u-u_h,{\gamma }_j+\pi u-u_h\right)}_{I_j}={\Vert {\gamma }_j\Vert }_{I_j}^2+2{\left({\gamma }_j,\pi u-u_h\right)}_{I_j}+{\Vert \pi u-u_h\Vert }_{I_j}^2\\ \le 2\left({\Vert {\gamma }_j\Vert }_{I_j}^2+{\Vert \pi u-u_h\Vert }_{I_j}^2\right).\end{array}$E312

Summing over all elements and applying (47d) with $s=0$ and (49) yields the first estimate in (50c).

Using the Cauchy‐Schwarz inequality and the inequality $\left|ab\right|\le \frac{1}{2}(a^2+b^2)$, we write

$\begin{array}{l}{\Vert {\omega }_{j^{\prime }}\Vert }_{I_j}^2={\left({\gamma }_{j^{\prime }}+(\pi u-u_h{)}^{\prime },{\gamma }_{j^{\prime }}+(\pi u-u_h{)}^{\prime }\right)}_{I_j}\le 2\left({\Vert {\gamma }_{j^{\prime }}\Vert }_{I_j}^2+{\Vert (\pi u-u_h{)}^{\prime }\Vert }_{I_j}^2\right).\hfill \end{array}$E52

Using the inverse inequality (8), i.e., ${\Vert (\pi u-u_h{)}^{\prime }\Vert }_{I_j}\le C{h}^{-1}{\Vert (\pi u-u_h)\Vert }_{I_j},$ we obtain the estimate

${\Vert {\omega }_{j^{\prime }}\Vert }_{I_j}^2\le C\left({\Vert {\gamma }_{j^{\prime }}\Vert }_{I_j}^2+h^{-2}{\Vert \pi u-u_h\Vert }_{I_j}^2\right).$E313

Summing over all elements and applying (49) and the estimate (47d) with $s=1$, we establish the second estimate in (50c).

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 $e$, we multiply (2) by arbitrary smooth function $v$ and integrate over the element $I_j$ to obtain

${\displaystyle {\int }_{I_j}}{u}^{\prime }vdt={\displaystyle {\int }_{I_j}}f(t,u)vdt.$E53

Replacing $u$ by $u_h+e$ and choosing $v={\psi }_{p+\mathrm{1,}j}(t)$, we obtain

${\displaystyle {\int }_{I_j}}{e}^{\prime }{\psi }_{p+\mathrm{1,}j}dt={\displaystyle {\int }_{I_j}}\left(f(t,u_h+e)-{u^{\prime }}_h\right){\psi }_{p+\mathrm{1,}j}dt.$E54

Substituting (50a), i.e., $e(t)={\alpha }_j{\psi }_{p+\mathrm{1,}j}(t)+{\omega }_j(t)$, into the left‐hand side of (54) yields

${\alpha }_j{\displaystyle {\int }_{I_j}}{{\psi }^{\prime }}_{p+\mathrm{1,}j}{\psi }_{p+\mathrm{1,}j}dt={\displaystyle {\int }_{I_j}}\left(f(t,u_h+e)-{u^{\prime }}_h-{{\omega }^{\prime }}_j\right){\psi }_{p+\mathrm{1,}j}dt.$E55

Using (32) and solving for ${\alpha }_j$, we obtain

${\alpha }_j=-\frac{1}{k_1{h}_j^{2p+2}}{\displaystyle {\int }_{I_j}}\left(f(t,u_h+e)-{u^{\prime }}_h-{{\omega }^{\prime }}_j\right){\psi }_{p+\mathrm{1,}j}dt.$E56

Our error estimate procedure consists of approximating the true error on each element $I_j$ by the leading term as

$e(t)\approx E(t)=a_j{\psi }_{p+\mathrm{1,}j}(t),t\in I_j,$E57a

where the coefficient of the leading term of the error, $a_j$, is obtained from the coefficient ${\alpha }_j$ defined in (56) by neglecting the terms ${\omega }_j$ and $e$, i.e.,

$a_j=-\frac{1}{k_1{h}_j^{2p+2}}{\displaystyle {\int }_{I_j}}\left(f(t,u_h)-u_{h^{\prime }}\right){\psi }_{p+\mathrm{1,}j}dt.$E57b

We remark that our a posteriori error estimate is obtained by solving local problems with no boundary condition.

The global effectivity index, defined by $\sigma =\frac{\Vert E\Vert }{\Vert e\Vert }$, 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 $E$, in the $L^2$ norm, converges to the actual error $e$. Moreover, we show that our a posterior error estimate is asymptotically exact by showing that the global effectivity index $\sigma \to 1$ as $h\to 0$.

Theorem 5.1. Suppose that the assumptions of Theorem 2 are satisfied. If $E(t)=a_j{\psi }_{p+\mathrm{1,}j}(t),t\in I_j,$ where $a_j,j=1,\dots ,N,$ are defined in (57b), then

${\Vert e-E\Vert }^2\le C{h}^{2p+4}.$E58

Thus, the post‐processed approximation $u_h+E$ yields $\mathcal{O}(h^{p+2})$ superconvergent solution, i.e.,

${\Vert u-(u_h+E)\Vert }^2={\displaystyle \sum _{j=1}^N}{\Vert u-(u_h+a_j{\psi }_{p+\mathrm{1,}j})\Vert }_{I_j}^2\le C{h}^{2p+4}.$E59

Furthermore, then there exists a positive constant $C$ independent of $h$ such that

$\left|{\Vert e\Vert }^2-{\Vert E\Vert }^2\right|\le C{h}^{2p+4}.$E60

Finally, if there exists a constant $C=C(u)>0$ independent of $h$ such that

$\Vert e\Vert \ge C{h}^{p+1},$E61

then the global effectivity index in the $L^2$ norm converges to unity at

(h) rate, i.e.,

$\frac{\Vert E\Vert }{\Vert e\Vert }=1+\mathcal{O}(h).$E62

Proof. First, we will prove (58) and (59). Since $e={\alpha }_j{\psi }_{p+\mathrm{1,}j}+{\omega }_j$ and $E=a_j{\psi }_{p+\mathrm{1,}j}$ on $I_j$, we have

${\Vert e-E\Vert }_{I_j}^2={\Vert ({\alpha }_j-a_j){\psi }_{p+\mathrm{1,}j}+{\omega }_j\Vert }_{I_j}^2\le 2({\alpha }_j-a_j{)}^2{\Vert {\psi }_{p+\mathrm{1,}j}\Vert }_{I_j}^2+2{\Vert {\omega }_j\Vert }_{I_j}^2,$E401

where we used the inequality ${(a+b)}^2\le 2a^2+2b^2$. Summing over all elements yields

${\Vert e-E\Vert }^2={\displaystyle \sum _{j=1}^N}{\Vert e-E\Vert }_{I_j}^2\le 2{\displaystyle \sum _{j=1}^N}{({\alpha }_j-a_j)}^2{\Vert {\psi }_{p+\mathrm{1,}j}\Vert }_{I_j}^2+2{\displaystyle \sum _{j=1}^N}{\Vert {\omega }_j\Vert }_{I_j}^2.$E63

Next, we will derive upper bounds for ${\displaystyle \sum _{j=1}^N}{({\alpha }_j-a_j)}^2{\Vert {\psi }_{p+\mathrm{1,}j}\Vert }_{I_j}^2$. Subtracting (56) from (57b), we obtain

$a_j-{\alpha }_j=\frac{1}{k_1{h}_j^{2p+2}}{\displaystyle {\int }_{I_j}}\left(f(t,u_h+e)-f(t,u_h)-{\omega }_{j^{\prime }}\right){\psi }_{p+\mathrm{1,}j}dt.$E64

Thus,

$\left|a_j-{\alpha }_j\right|\le \frac{1}{k_1{h}_j^{2p+2}}{\displaystyle {\int }_{I_j}}\left(\left|f(t,u_h+e)-f(t,u_h)\right|+\left|{{\omega }^{\prime }}_j\right|\right)\left|{\psi }_{p+\mathrm{1,}j}\right|dt.$E65

Using the Lipschitz condition (3) and applying the Cauchy‐Schwarz inequality yields

$\left|a_j-{\alpha }_j\right|\le \frac{1}{k_1{h}_j^{2p+2}}{\displaystyle {\int }_{I_j}}\left(M_1\left|e\right|+\left|{{\omega }^{\prime }}_j\right|\right)\left|{\psi }_{p+\mathrm{1,}j}\right|dt\le \frac{{\Vert {\psi }_{p+\mathrm{1,}j}\Vert }_{I_j}}{k_1{h}_j^{2p+2}}\left(M_1{\Vert e\Vert }_{I_j}+{\Vert {{\omega }^{\prime }}_j\Vert }_{I_j}\right).$E66

Applying the inequality ${(a+b)}^2\le 2(a^2+b^2)$, we obtain

${(a_j-{\alpha }_j)}^2\le \frac{2{\Vert {\psi }_{p+\mathrm{1,}j}\Vert }_{I_j}^2}{k_1^2{h}_j^{4p+4}}\left(M_1^2{\Vert e\Vert }_{I_j}^2+{\Vert {{\omega }^{\prime }}_j\Vert }_{I_j}^2\right).$E67

Multiplying by ${\Vert {\psi }_{p+\mathrm{1,}j}\Vert }_{I_j}^2$ and using (32), i.e., ${\Vert {\psi }_{p+\mathrm{1,}j}\Vert }_{I_j}^2=(2p+2)k_2{h}_j^{2p+3}$ yields

${(a_j-{\alpha }_j)}^2{\Vert {\psi }_{p+\mathrm{1,}j}\Vert }_{I_j}^2\le \frac{2{\Vert {\psi }_{p+\mathrm{1,}j}\Vert }_{I_j}^4}{k_1^2{h}_j^{4p+4}}\left(M_1^2{\Vert e\Vert }_{I_j}^2+{\Vert {{\omega }^{\prime }}_j\Vert }_{I_j}^2\right)\le k_3{h}_j^2\left({\Vert e\Vert }_{I_j}^2+{\Vert {{\omega }^{\prime }}_j\Vert }_{I_j}^2\right),$E68

where $k_3=\frac{2(2p+{2)}^2{k}_2^2}{k_1^2}\max (M_1^2,1)$ is a constant independent of the mesh size.

Summing over all elements and using $h=\underset{1\le j\le N}{{\displaystyle \max }}{h}_j$, we arrive at

${\displaystyle \sum _{j=1}^N}{(a_j-{\alpha }_j)}^2{\Vert {\psi }_{p+\mathrm{1,}j}\Vert }_{I_j}^2\le k_3{h}^2\left[{\Vert e\Vert }^2+{\displaystyle \sum _{j=1}^N}{\Vert {{\omega }^{\prime }}_j\Vert }_{I_j}^2\right].$E402

Combining this estimate with (9) and (50c), we establish

${\displaystyle \sum _{j=1}^N}{\left(a_j-{\alpha }_j\right)}^2{\Vert {\psi }_{p+\mathrm{1,}j}\Vert }_{I_j}^2\le C{h}^{2p+4}.$E69

Now, combining (63) and the estimates (50c) and (69) yields

${\Vert e-E\Vert }^2\le 2C_1{h}^{2p+4}+2C_2{h}^{2p+4}=C{h}^{2p+4},$E403

which completes the proof of (58). Using the relation $e=u-u_h$ and the estimate (58), we obtain

${\displaystyle \sum _{j=1}^N}{\Vert u-(u_h+a_j{\psi }_{p+\mathrm{1,}j})\Vert }_{I_j}^2={\Vert u-(u_h+E)\Vert }^2={\Vert e-E\Vert }^2\le C{h}^{2p+4}.$E404

Next, we will prove (60). Using the reverse triangle inequality, we have

$|\Vert E\Vert -\Vert e\Vert |\le \Vert E-e\Vert ,$E70

which, after applying the estimate (58), completes the proof of (60).

In order to show (62), we divide (70) by $\Vert e\Vert$ to obtain $|\sigma -1|\le \frac{\Vert E-e\Vert }{\Vert e\Vert }$. Applying the estimate (58) and the inverse estimate (61), we arrive at

$|\sigma -1|\le Ch.$E405

Therefore, $\sigma =\frac{\Vert E\Vert }{\Vert e\Vert }=1+$

(h), which establishes (62).

Remark 5.1. The previous theorem indicates that the computable quantity $\Vert E\Vert$ converges to $\Vert e\Vert$ at

(h^p+2) rate. This accuracy enhancement is achieved by adding the error estimate $E$ to the DG solution $u_h$.

Remark 5.2. The performance of an error estimator $\sigma$ is typically measured by the global effectivity index which is defined as the ratio of the estimated error $\Vert E\Vert$ to the actual error $\Vert e\Vert$. We say that the error estimator is asymptotically exact if $\sigma \to 1$ as $h\to 0$. The estimate (62) indicates that the global effectivity index in the $L^2$ norm converge to unity at

(h) rate. Therefore, the proposed estimator $\Vert E\Vert$ is asymptotically exact. We would like to emphasize that $E$ is a computable quantity since it only depends on the DG solution $u_h$ and $f$. It provides an asymptotically exact a posteriori estimator on the actual error $\Vert e\Vert$. Finally, we would like to point out that our procedure for estimating the error $e$ is computationally simple. Furthermore, our DG error indicator is obtained by solving a local problem with no boundary condition on each element. This makes it useful in adaptive computations. We demonstrate this in Section 6.

Remark 5.3. Our proofs are valid for any regular meshes and using piecewise polynomials of degree $p\ge 1$. If $p=0$ then (46) gives $\Vert \overline{e}\Vert =$

(h) which is the same as $\Vert e\Vert =$

(h). Thus, our superconvergence results are not valid when using $p=0$. Also, our error estimate procedure does not apply.

Remark 5.4. The assumption (61), which is used to prove the convergence of $\sigma$ to unity at

(h), requires that terms of order

(h^p+1) are present in the error. If not, $E$ might not be a good approximation of $e$. We note that the exponent of $h$ in the estimate (9) is optimal in the sense that it cannot be improved. In fact, for the $h$version finite element method one may show that provided that the $(p+1)$th order derivatives of the exact solution $u$ do not vanish identically over the whole domain, then an inverse estimate of the form $\Vert e\Vert \ge C(u)h^{p+1}$ is valid for some positive constant $C(u)$ depending only on $u$ [42–44].

Remark 5. Our results readily extend to nonlinear systems of ODEs of the form

$\frac{d\overrightarrow{u}}{dt}=\overrightarrow{f}(t,\overrightarrow{u}),t\in [0,T],\overrightarrow{u}(0)={\overrightarrow{u}}_0,$E406

where $\overrightarrow{u}=[u_1,\dots ,u_n{]}^t:[0,T]\to {ℝ}^n$, ${\overrightarrow{u}}_0\in {ℝ}^n,$ and $\overrightarrow{f}=[f_1,\dots ,f_n{]}^t:[0,T]\times {ℝ}^n\to {ℝ}^n$. The DG method for this problem consists of finding ${\overrightarrow{u}}_h\in {\overrightarrow{V}}_h^p=\{\overrightarrow{v}:\overrightarrow{v}{|}_{I_j}\in {(P^p(I_j))}^n,j=1,\dots ,N\}$ such that: $\forall \overrightarrow{v}\in {\overrightarrow{V}}_h^p$ and $j=1,\dots ,N$,

${\displaystyle {\int }_{I_j}}{(\overrightarrow{v}\text{'})}^t{\overrightarrow{u}}_{h}dt+{\displaystyle {\int }_{I_j}}{(\overrightarrow{v})}^t\overrightarrow{f}(t,{\overrightarrow{u}}_h)dt-{(\overrightarrow{v})}^t(t_j^{-})\overrightarrow{u}(t_j^{-})+{(\overrightarrow{v})}^t(t_{j-1}^{+})\overrightarrow{u}(t_{j-1}^{-})=0.$E407

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

$\text{SOLVE}\to \text{ESTIMATE}\to \text{MARK}\to \text{REFINE}.$E410

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:

1. Select a tolerance $Tol$ and a maximum bound on the number of interval (say $N_{\text{max}}=1000$). Put $\Vert E\Vert =1$.

2. Construct an initial coarse mesh with $N+1$ nodes. For simplicity, we start with a uniform mesh having $N=2$ elements.

3. While $N+1\le N_{\text{max}}$ and $\Vert E\Vert \ge Tol$ do

   1. (a)

      Solve the DG scheme to obtain the solution $u_h$ as described in Section 2.

   2. (b)

      For each element, use (57a) and (57b) to compute the local error estimators ${\Vert E\Vert }_{I_j},j=1,\dots ,N$ as described in Section 5 and the global error estimator $\Vert E\Vert ={\left({\displaystyle {\sum }_{j=1}^N}{\Vert E\Vert }_{I_j}^2\right)}^{1/2}$.

   3. (c)

      For all elements $I_j$

      1. Choose a parameter $0\le \lambda \le 1$. If the estimated global error ${\Vert E\Vert }_{I_j}<\lambda \underset{j=1,\dots ,N}{{\displaystyle \max }}{\Vert E\Vert }_{I_j}$ then stop and accept the DG solution on the element $I_j$.

      2. Otherwise, reject the DG solution on $I_j$ and divide the element $I_j$ into two uniform elements by adding the coordinate of the midpoint of $I_j$ to the list of nodes.

4. Endwhile.

Remark 6.1. There are many possibilities for selecting the elements to be refined given the local error indicator ${\Vert E\Vert }_{I_j}$. In the above algorithm, we used the most popular fixed‐rate strategy which consists of refining the element $I_j$ if ${\Vert E\Vert }_{I_j}>\lambda \underset{j=1,\dots ,N}{{\displaystyle \max }}{\Vert E\Vert }_{I_j}$, where $0\le \lambda \le 1$ is a parameter provided by the user. Note that the choice $\lambda =0$ gives uniform refinement, while $\lambda =1$ 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

$u^{\prime }=-u-u^2,t\in [0,1],u(0)=1.$E408

Clearly, the exact solution is $u(t)=\frac{1}{2e^t-1}.$ We use uniform meshes obtained by subdividing the computational domain $[\mathrm{0,1}]$ into $N$ intervals with $N=$ 5, 10, 20, 30, 40, 50. This example is tested by using $P^p$ polynomials with $p=0-4$. Figure 1 shows the $L^2$ errors $\Vert e\Vert$ and $\Vert \overline{e}\Vert$ with log‐log scale as well as their orders of convergence. These results indicate that $\Vert e\Vert =$

(h^p+1) and $\Vert \overline{e}\Vert =$

(h^p+2). This example demonstrates that our theoretical convergence rates are optimal.

Figure 1.

Log‐log plots of $\Vert e\Vert$ (left) and $\Vert \overline{e}\Vert$ (right) versus mesh sizes $h$ for Example 7.1 on uniform meshes having $N=$ 5, 10, 20, 30, 40, 50 elements using $P^p$, $p=0$ to 4.

Figure 2.

Log‐log plots of $\Vert e{\Vert }^{*}$ (left) versus $h$ for Example 7.1 using $N=$ 5, 10, 20, 30, 40, 50 and $P^p$, $p=0$ to 4. Log‐log plots of ${\left|e\right|}^{*}$ (right) versus $h$ using $N=$ 5, 10, 20, 30 elements using $P^p$, $p=0$ to 3.

Next, we compute the maximum error at the shifted roots of the $(p+1)$ degree right Radau polynomial on each element $I_j$ and then take the maximum over all elements. For simplicity, we use ${\Vert e\Vert }^{*}$ to denote ${{\displaystyle \max }}_{1\le j\le N}\left({{\displaystyle \max }}_{0\le i\le p}|e(t_{j,i}^{-})|\right)$, where $t_{j,i}$ are the roots of $R_{p+\mathrm{1,}j}(t)$. Similarly, we compute the true error at the downwind point of each element and then we denote ${\left|e\right|}^{*}$ to be the maximum over all elements $I_j,j=1,\dots ,N$, i.e., ${\left|e\right|}^{*}={{\displaystyle \max }}_{1\le j\le N}|e(t_j^{-})|$. In Figure 2, we present the errors ${\Vert e\Vert }^{*}$, ${\left|e\right|}^{*}$ and their orders of convergence. We observe that ${\Vert e\Vert }^{*}=$

(h^p+2) and ${\left|e\right|}^{*}=$

(h^2p+1) as expected. Thus, the error at right Radau points converges at

(h^p+2). Similarly, the error at the downwind point of each element converge with an order $2p+1$. This is in full agreement with the theory.

Next, we use (57a) and (57b) to compute the a posteriori error estimate for the DG solution. The global errors $\Vert e-E\Vert$ and their orders of convergence, using the spaces $P^p$ with $p=1-4$, are shown in Figure 3. We observe that $\Vert e-E\Vert =$

(h^p+2). This is in full agreement with the theory. This example demonstrates that the convergence rate proved in this work is sharp. Since $\Vert e-E\Vert =\Vert u-(u_h+E)\Vert =$

(h^p+2), we conclude that the computable quantities $u_h+E$ converges to the exact solution $u$ at

(h^p+2) rate in the $L^2$ norm. We would like to emphasize that this accuracy enhancement is achieved by adding the error estimate $E$ to the DG solution $u_h$ only once at the end of the computation. This leads to a very efficient computation of the postprocessed approximation $u_h+E$.

Figure 3.

The errors $\Vert e-E\Vert$ and their orders of convergence for Example 1 on uniform meshes having $N=$ 5, 10, 20, 30, 40, 50 elements using $P^p$, $p=1$ to 4.

In Table 1, we present the actual L² errors and the global effectivity indices. These results demonstrate that the proposed a posteriori error estimates is asymptotically exact.

$N$	$p=1$		$p=2$		$p=3$		$p=4$
	$\Vert e\Vert$	$\sigma$	$\Vert e\Vert$	$\sigma$	$\Vert e\Vert$	$\sigma$	$\Vert e\Vert$	$\sigma$
5	4.7637e-03	1.0362	2.7867e-04	1.0531	1.6847e-05	1.0637	1.0386e-06	1.0705
10	1.2750e-03	1.0179	3.7805e-05	1.0271	1.1742e-06	1.0326	3.7481e-08	1.0363
20	3.2849e-04	1.0089	4.8747e-06	1.0136	7.6227e-08	1.0164	1.2290e-09	1.0182
30	1.4736e-04	1.0059	1.4568e-06	1.0090	1.5201e-08	1.0109	1.6369e-10	1.0122
40	8.3262e-05	1.0044	6.1698e-07	1.0068	4.8296e-09	1.0082	3.9026e-11	1.0091
50	5.3429e-05	1.0035	3.1660e-07	1.0054	1.9827e-09	1.0066	1.2820e-11	1.0073

Table 1.

The errors $\Vert e\Vert$ and the global effectivity indices for Example 7.1 on uniform meshes having $N=$ 5, 10, 20, 30, 40, 50 elements using $P^p$, $p=1$ to 4.

In Figure 4, we show the errors $\delta e=\text{|‖}e\Vert -\Vert E\text{‖|}$ and $\delta \sigma =|\sigma -1|$. We see that $\delta e=$

(h^p+2) and $\delta \sigma =$

(h) as the theory predicts.

Figure 4.

Convergence rates for $\delta e$ (left) and $\delta \sigma$ (right) for Example 1 on uniform meshes having $N=$5, 10, 20, 30, 40, 50 elements using $P^p$, $p=1$ to 4.

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

$u^{\prime }=\beta u,t\in [\mathrm{0,5}],u(0)=\mathrm{1,}$E409

where the exact solution is simply $u(t)=e^{\beta t}$. We apply our adaptive algorithm using $\beta =1$ (unstable), $\beta =-1$ (stable), and $\beta =-20$ (stiff). The DG solutions and the sequence of meshes obtained by applying our adaptive algorithm with $Tol{=10}^{-2}$ for $p=1-4$ are shown in Figures 5–7 for $\beta =1$, $\beta =-1$, and $\beta =-20$, respectively. As can be expected, the adaptive algorithm refines in the vicinity of the endpoint $t=5$ with coarser meshes for increasing polynomial degree $p$. Furthermore, we observe that, when $\lambda$ is closer to 0, we get more uniform refinement near the portion with high approximation error. When $\lambda$ is near 1, we get less uniform refinement near the portion with high approximation error. We also observed that, both for $\lambda =0.2$ and $\lambda =0.9$, the optimal convergence rates are achieved asymptotically and that the global effectivity indices converge to unity with increasing polynomial degree $p$. Furthermore, we tested our adaptive algorithm on other problems and observed similar conclusions. These results are not included to save space.

Figure 5.

$u$, $u_h$, and final meshes for Example 7.2 with $\beta =-20$ using $P^p$, $p=1$ to 4, and $Tol{=10}^{-2}$.

Figure 6.

u, u_h, and final meshes for Example 7.2 with β = −1 using P^p, p = 1 to 4, and Tol = 10⁻².

Figure 7.

u, u_h, and final meshes for Example 7.2 with β = −20 using P^p, p = 1 to 4, and Tol = 10⁻².

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 $p+1$, when the space of piecewise polynomials of degree $p$ is used. We further proved the $(2p+1)$th superconvergence rate at the downwind points. Moreover, we proved that the DG solution is

(h^p+2) superconvergent toward a particular projection of the exact solution. We used these results and showed that the leading term of the DG error is proportional to the $(p+1)$ degree right Radau polynomial. This result allowed us to construct computationally simple, efficient, and asymptotically exact a posteriori error estimator. It is obtained by solving a local residual problem with no boundary condition on each element. Furthermore, we proved that the proposed a posteriori error estimator converges to the actual error in the $L^2$ norm. The order of convergence is proved to be $p+2$. All proofs are valid for regular meshes and for $P^p$ polynomials with $p\ge 1$. Finally, we presented a local adaptive procedure that makes use of our local a posteriori error estimate. Future work includes the study of superconvergence of DG method for nonlinear boundary‐value problems.

Abbreviations

Symbols