A posteriori Error Estimation

0 downloads 0 Views 5MB Size Report
Aug 24, 2007 - posteriori error estimates for discontinuous finite element solutions of scalar first- order hyperbolic partial ..... 1 + 455ϕ2. 2−. 54ϕ0. 1 + 108ϕ2.
The Discontinuous Galerkin Method for Two-dimensional Hyperbolic Problems Part II: A posteriori Error Estimation Slimane Adjerid and Mahboub Baccouch Department of Mathematics Virginia Polytechnic Institute and State University Blacksburg, VA 24061 August 24, 2007 Abstract In this manuscript we construct simple, efficient and asymptotically correct a posteriori error estimates for discontinuous finite element solutions of scalar firstorder hyperbolic partial differential problems on triangular meshes. We explicitly write the basis functions for the error spaces corresponding to several finite element spaces. The leading term of the discretization error on each triangle is estimated by solving a local problem. We also show global superconvergence for discontinuous solutions on triangular meshes. The a posteriori error estimates are tested on several linear and nonlinear problems to show their efficiency and accuracy under mesh refinement for smooth and discontinuous solutions. Keywords: Discontinuous Galerkin method; hyperbolic problems; superconvergence; a posteriori error estimation.

1

Introduction

The discontinuous Galerkin (DG) finite element method was first used to solve the neutron equation [18] and then studied for initial-value problems for ordinary differential equations [3, 8, 18, 20, 21]. Cockburn and Shu [12] extended the method to solve firstorder hyperbolic partial differential equations of conservation laws. They also developed the Local Discontinuous Galerkin (LDG) method for convection-diffusion problems [13]. Consult [11] and the references cited therein for more information on DG methods. DG methods allow discontinuous bases, which simplify both h-refinement (mesh refinement and coarsening) and p-refinement (method order variation). However, for DG methods to be used in an adaptive framework one needs a posteriori error estimates to guide adaptivity and stop the refinement process. The solution space consists of piecewise continuous polynomial functions relative to a structured or unstructured mesh. As 1

such, it can sharply capture solution discontinuities relative to the computational mesh. It maintains local conservation on an elemental basis. The success of the DG method is due to the following properties: (i ) does not require continuity across element boundaries, (ii ) is locally conservative, (iii) is well suited to solve problems on locally refined meshes with hanging nodes, (iv ) exhibits strong superconvergence that can be used to estimate the discretization error, (v ) has a simple communication pattern between elements with a common face that makes it useful for parallel computation and (vi ) it can handle problems with complex geometries to high order. A posteriori error estimates lie in the heart of every adaptive finite element algorithm for differential equations. The subject has attained a certain level of maturity for elliptic problems [7, 22]. However, a posteriori error estimation is much less developed for DG methods applied to hyperbolic problems. The first estimates were developed by Adjerid et al. [3] who constructed the first a posteriori DG error estimates for one-dimensional linear and nonlinear hyperbolic problems. Later, Adjerid and Massey [5, 6] showed how to construct accurate error estimates for multi-dimensional problems on rectangular meshes where they showed that the leading term of error is spanned by two (p+1)-degree Radau polynomials in the x and y directions, respectively. Krivodonova and Flaherty [17] showed that the leading term of the local discretization error on triangles having one outflow edge is spanned by a suboptimal set of orthogonal polynomials of degree p and p + 1 and computed DG error estimates by solving local problems involving numerical fluxes, thus requiring information from neighboring inflow elements. Adjerid and Klauser [4] constructed a posteriori error estimates for the local discontinuous Galerkin method applied to diffusion problems. Superconvergence properties for DG methods have been studied in [3, 14, 16, 18] for first-order ordinary differential equations, [3, 5, 6] for hyperbolic problems and [4, 9, 10] for diffusion and convection-diffusion problems. In this manuscript we write explicitly the form of the leading asymptotic term of the local error for each type of elements and polynomial space and derive optimal sets of orthogonal polynomials that span the leading term of the error. Moreover, we show that the local superconvergence results [2] extend to global DG solutions on general meshes when we use a corrected inflow boundary condition. We first classify triangular elements into three types: (i) type I with one inflow edge and two outflow edges, (ii) type II with two inflow edges and one outflow edge and (iii) type III with one inflow edge , one outflow edge , one edge parallel to the characteristics. We observe that the DG solution does not exhibit any superconvergence on meshes consisting of elements of type I. Thus, elements of type I prevent superconvergence from propagating through the whole mesh. In order to recover the superconvergence results [2] we use a modified DG method with corrected inflow boundary conditions and the augmented space Up that will be defined later. We present several numerical results confirming the O(hp+2 ) pointwise superconvergence rates on all elements for three polynomial spaces. As described above, we note that our error estimation procedure is closely tied to the global superconvergence properties of the DG solutions. We further use these new results to discover the optimal finite element space for the error on each element and solve a local Galerkin problem that does not require numerical fluxes for the error and thus are 2

much more efficient than the estimates presented in [17]. This paper is organized as follows, in §2 we present a model problem and recall the DG formulation and finite element spaces. In §3 we present optimal finite element spaces for the error. In §4 we discuss an error estimation procedure and present numerical examples in §5. We conclude and discuss our results in §6.

2

Discontinuous Galerkin formulation

We consider the model problem (x, y) ∈ Ω, a = [α, β]T ∈ R+2 ,

a.∇u + cu = f (x, y),

(2.1a)

subject to the boundary conditions u(x, 0) = g0 (x), u(0, y) = g1 (y),

(2.1b)

where the functions f (x, y), g0 (x) and g1 (y) are selected such that the exact solution u(x, y) ∈ C ∞ (Ω). Here we consider general meshes with triangles of type I, having one inflow edge and two outflow edges, triangles of type II, having two inflow edges and one outflow edge and triangles of type III, having one inflow edge, one outflow edge and one edge parallel to a. We approximate u(x, y) by a piecewise polynomial function U (x, y) whose restriction to ∆ is in Pp , Vp or Up defined as p m X X

i m−i cm }, i x y

(2.2a)

Vp = Pp ∪ {xi y p+1−i , 1 ≤ i ≤ p},

(2.2b)

Up = Pp ∪ {xp+1 , y p+1 },

(2.2c)

Pp = {q : q =

m=0 i=0

p

Pp = {q : q =

X

ci xi },

(2.2d)

i=0

where V0 = P0 and U0 = P1 . The discrete DG method consists of finding U ∈ Wp such that Z Z ZZ ZZ a.nUˆ V ds + a.nU V ds + (−a.∇V + cV )U dxdy = f V dxdy, Γ−

Γ+





∀ V ∈ Wp ,

(2.3)

where Wp is either Pp , Vp or Up . If n is the outer unit normal on the boundary of ∆, we let Γ− denote the inflow boundary of ∆ where a.n < 0 while Γ+ denotes the outflow boundary where a.n > 0. 3

In order to complete the definition of the DG method we need to select the flux Uˆ on Γ− . On meshes where all elements are of type III one may use the standard flux Uˆ = U − .

(2.4)

However, on meshes including elements of type I and II we will use the higher-order flux Uˆ = U − + E − ,

(2.5)

on all interior inflow edges and Uˆ = u on the physical boundary. Here E is an a posteriori error estimate that will be discussed later. Following the line of reasoning of [2], the DG orthogonality condition for the local error ² = u − U can be written on the canonical triangle ∆ as ZZ Z Z a.n²V ds + (−a.∇V + hcV )²dξdη = 0, ∀ V ∈ Wp . (2.6) a.nˆ²V ds + Γ−

Γ+



To enhance the conditioning of the finite element problem we use orthogonal polynomials [15] as shape functions given on the canonical triangle ∆ defined by the vertices A = (0, 0), B = (1, 0) and C = (0, 1) as ϕlk (ξ, η) = 2k Lk (

2ξ − 1)(1 − η)k Pl2k+1,0 (2η − 1), k, l = 0, 1, · · · 1−η

(2.7a)

where Pnα,β (x) and Ln are the n-th degree Jacobi and Legendre polynomials, respectively, on [−1, 1] [1]. Dubiner polynomials satisfy the L2 orthogonality Z

1 0

Z

1−ξ 0

ϕlk ϕqp dηdξ = clq kp δkp δlq .

(2.7b)

We also need Radau polynomials defined as Rp+1 (x) = Lp+1 (x) − Lp (x).

(2.8)

In the next Lemma we introduce a set of orthogonal polynomials in L2 with respect to the weight function w = 1 − ξ − η on the standard triangle. Lemma 2.1. The polynomials qij = Pj1,0 (

2ξ − 1)(1 − η)j Pi2j+2,0 (2η − 1), i, j = 0, 1, · · · 1−η

satisfy the orthogonality condition Z 1 Z 1−η (1 − ξ − η)qij qlm dξdη = cjm il δil δjm , 0

0

where δij is the Kronecker symbol equal to 1 if i = j and 0 otherwise. Thus, {qij , i, j = 0, 1, · · · , p, i + j ≤ p} is a basis for Pp . 4

(2.9a)

(2.9b)

Proof. Let us compute the inner product ZZ Z 1Z j m (1 − ξ − η)qi ql dξ dη = [ ∆

0

Z

1 0

1−η 0

(1 − ξ − η)qij qlm dξ] dη =

[Pi2j+2,0 (2η − 1)Pl2m+2,0 (2η − 1)(1 − η)j+m G(η)]dη,

where

Z

1−η

G(η) = 0

(1 − ξ − η)Pj1,0 (

2ξ 2ξ − 1)Pm1,0 ( − 1)dξ. 1−η 1−η

Applying the change of variables t = 2η − 1 and z = 1−ξ−η =

(1−η)(1−z) 2

2ξ 1−η

− 1, dt = 2dη , dz =

(2.10)

(2.11) 2dξ 1−η

and

to (2.10) we obtain ZZ (1 − ξ − η)qij qlm dξ dη = ∆

1

Z

1

Z [Pi2j+2,0 (t)Pl2m+2,0 (t)(1

1

(1 − z)Pj1,0 (z)Pm1,0 (z)dz] dt = 2j+m+5 −1 −1 Z 1 Z 1 1 ( (1 − t)j+m+2 Pi2j+2,0 (t)Pl2m+2,0 (t) dt)( (1 − z)Pj1,0 (z)Pm1,0 (z)dz). (2.12) j+m+5 2 −1 −1 j+m+2

− t)

Applying the orthogonality properties of Jacobi polynomials we prove that Z 1 Z 1 2m+2,0 j+m+2 2j+2,0 ( (1 − t) Pi (t)Pl (t) dt)( (1 − z)Pj1,0 (z)Pm1,0 (z)dz) = cjm il δil δjm . −1

−1

Thus, qij , 0 ≤ i + j ≤ p, are orthogonal and form a basis for Pp . For the sake of completeness, we include the following results from [2] which will be needed in our a posteriori error analysis. These theorems state the orthogonality conditions for the leading term of the local discontinuous Galerkin error. Theorem 2.1. Let u ∈ C ∞ (4) and U ∈ Pp (4) be the solutions of (2.1) and (2.3), respectively, with Uˆ |Γ− = u. If α ≥ 0 and β ≥ 0 such that 4 is either a triangle of type II or type III, then the local finite element error can be written as ²(ξ, η) =

∞ X

hk Qk (ξ, η),

(2.13)

k=p+1

where

ZZ Qp+1 V dξ dη = 0,

∀ V ∈ Pp−1 ,

(2.14)

a.nQp+1 V ds = 0,

∀ V ∈ Pp ,

(2.15)



Z

Γ+

Z

a.nQk ds = 0, Γ+

5

k ≥ p + 1,

(2.16)

Qp+1 (ξ, η) =

p X

cpi ϕip−i (ξ, η)

i=0

+

p+1 X

ϕip+1−i (ξ, η). cp+1 i

(2.17)

i=0

Furthermore, at the outflow boundary of the physical element ∆ the local error satisfies Z a.n²ds = O(h2p+2 ). (2.18) Γ+

Proof. Cf. Krivodonova and Flaherty [17]. Theorem 2.2. Under the same assumptions as in Theorem 2.1 there exist two constants C1 and C2 such that on the outflow edge we have Qp+1 (1 − η, η) = C1 Lp+1 (2η − 1).

(2.19a)

Qp+1 (ξ, 1 − ξ) = C2 Lp+1 (2ξ − 1).

(2.19b)

Furthermore, ZZ ∂Qp+1 i j ξ η dξ dη = 0, ∆ ∂ξ ZZ ∂Qp+1 i j ξ η dξ dη = 0, ∆ ∂η and

i = 1, ..., p, j = 0, ..., p − 1, i + j ≤ p,

(2.20)

i = 0, ..., p − 1, j = 1, ..., p, i + j ≤ p,

(2.21)

ZZ a.∇(1 − ξ − η)i Qp+1 dηdξ = 0,

i = 1, ..., p.

(2.22)

4

Proof. Cf. [2]. Theorem 2.3. Let p ≥ 1 and 4 = {(ξ, η), ξ, η ≥ 0, ξ + η ≤ 1}. Let u ∈ C ∞ (4) and U ∈ Up (4) be the solution of (2.1) and (2.3),respectively with Uˆ |Γ− = u. If α ≥ 0 and β ≥ 0 such that 4 is either of type II or type III, then the local finite element error can be written as in (2.13) where the leading term Qp+1 , satisfies the following conditions Z 1 Z 1 Z 1−ξ ∂Qp+1 (ξ, η) p+1 βQp+1 (ξ, 1 − ξ)ξ dξ + αξ p+1 dηdξ = 0, (2.23) ∂ξ 0 0 0 and

Z

Z

1

αQp+1 (1 − η, η)η 0

p+1

1

Z

1−η

dη + 0

βη p+1

0

∂Qp+1 (ξ, η) dξdη = 0. ∂η

(2.24)

If either α = 0, β > 0 or α > 0, β = 0, i.e., ∆ is of type III, then the leading term of the local error is zero on the outflow edge Qp+1 (1 − η, η) = 0,

0 ≤ η ≤ 1.

Furthermore, if α > 0, β = 0, then Z 1 Z 1−ξ ξ k Qp+1 (ξ, η)dηdξ = 0, 0

0

6

0 ≤ k ≤ p,

(2.25)

(2.26)

Z

1−ξ 0

∂ k Qp+1 (ξ, η) dη = 0, ∂ξ k

and

ZZ 4

k = 0, 1,

∂Qp+1 p+1 ξ dηdξ = 0, ∂ξ

Similarly, if α = 0, β > 0, then Z 1 Z 1−η η k Qp+1 (ξ, η)dξdη = 0, 0

Z

1−η 0

0 ≤ ξ ≤ 1.

(2.27)

(2.28)

0 ≤ k ≤ p,

(2.29)

0

∂ k Qp+1 (ξ, η) dξ = 0, ∂η k

k = 0, 1,

0 ≤ η ≤ 1.

(2.30)

Proof. Cf. [2]. Theorem 2.4. Under the assumption of Theorem 2.1 and using the polynomial space Vp the leading term in the local DG error on an element of type II or III satisfies ZZ a.∇Qp+1 ξ i η p+1−i dηdξ = 0, i = 1, ..., p. (2.31) 4

Moreover, on an element of type III with α > 0, β = 0 ZZ ∂Qp+1 i p+1−i ξη dηdξ = 0, i = 1, ..., p, 4 ∂ξ and Qp+1 (ξ, η) = C1 Lp+1 (2η − 1) + C2 (1 − η)p+1 Rp+1 ( Similarly, if α = 0, β > 0 ZZ

∂Qp+1 i p+1−i ξη dηdξ = 0, 4 ∂η

2ξ − 1). 1−η

i = 1, ..., p,

and Qp+1 (ξ, η) = C1 Lp+1 (2ξ − 1) + C2 (1 − ξ)p+1 Rp+1 (

2η − 1). 1−ξ

(2.32)

(2.33)

(2.34)

(2.35)

Proof. Cf. [2]. Theorem 2.5. Let p ≥ 1 and 4 denote the reference triangle defined by the vertices (1, 0), (1, 1) and (0, 1). Let u ∈ C ∞ (4) and U ∈ Up (4) be the solution of (2.1) and (2.3), respectively with Uˆ |Γ− = u. If α > 0 and β > 0, then the local error can be written as (2.13) where the leading term , Qp+1 , satisfies the following orthogonality conditions Z a · nQp+1 (βξ − αη)i ds = 0, i = 0, ..., p, (2.36) Γ+

ZZ

a.∇(ξ − 1)i (η − 1)j−i Qp+1 dηdξ = 0, 4

7

1 ≤ i < j ≤ p,

(2.37)

Z Z

1 0 1

ZZ βQp+1 (ξ, 1)(ξ − 1) dξ − i i

ZZ αQp+1 (1, η)(η − 1) dη − i

α(ξ − 1)i−1 Qp+1 dξ dη = 0,

i = 1, ..., p + 1, (2.38)

β(η − 1)i−1 Qp+1 dξ dη = 0,

i = 1, ..., p + 1. (2.39)

4

i

0

4

Furthermore, Z

1

βQp+1 (ξ, 1) + αQp+1 (ξ, 1 − ξ) + α 1−ξ

Z

1

αQp+1 (1, η) + βQp+1 (1 − η, η) + β 1−η

∂Qp+1 (ξ, η) dη = 0, ∂ξ

0 ≤ ξ ≤ 1,

(2.40)

∂Qp+1 (ξ, η) dξ = 0, ∂η

0 ≤ η ≤ 1.

(2.41)

Finally, we have the pointwise superconvergence Qp+1 (0, 1) = Qp+1 (1, 0) = 0.

(2.42)

Proof. Cf. [2].

3

Finite element spaces for the error

In order to construct efficient a posteriori error estimates for the DG method we will construct optimal finite spaces for the leading term, Qp+1 , of the error on each element. We apply several results from the previous section on the local error analysis and study each type of element separately for the three spaces Pp , Up and Vp . Finally, we present a technique to compute asymptotically correct a posteriori estimates of the DG error e = u − U on each element as well as on the whole domain.

3.1

Basis functions for element of type I

On triangles of type I there is no optimal set of functions that spans the leading term of the error for Pp and Vp . Therefore, we will only consider Up on the reference triangle of type I defined by the vertices (1, 0), (1, 1) and (0, 1) for α > 0 and β > 0. If we assume the conditions of Theorem 2.5 are satisfied, we show that the leading term Qp+1 of local DG finite element error can be written as Qp+1 (ξ, η) =

p X

cip+1−i χip+1−i (ξ, η),

(3.1)

i=1

where χip−i are computed using Mathematica in terms of ϕij−i , i, j = 0, ..., p, and are given in Table 1. Next, we outline a procedure to obtain the error basis functions χip−i . Since Dubiner polynomials form a basis for Pp , the finite element error the standard triangle of type I can be approximated by its leading term in Pp as e(ξ, η) ≈ E(ξ, η) = Qp+1 =

p+1 p+1 X X i=0 j=i

8

cji ϕij−i (ξ, η),

(3.2)

which is defined by (p+2)(p+3) parameters. 2 p2 +3p+6 Applying the error orthogonality conditions (2.36), (2.37), (2.39) and (2.42) from 2 Theorem 2.5 yields E(ξ, η) =

p X

cip+1−i χip+1−i (ξ, η),

(3.3)

i=1

χip−i

where the functions shown in Table 1 are given in terms of ϕij−i , i, j = 0, · · · , p, and s = α/β. We note that the leading term of the error in (3.3) is determined by p parameters. Table 1: Error basis functions for the spaces Up for p = 1 to 3 on elements of type I where s = α/β. p=1

p=2

χ11 = (12(ϕ01 + ϕ11 )s2 + 2(10ϕ00 − 2ϕ10 − 3ϕ01 − 2ϕ20 + 12ϕ11 + 5ϕ02 )s + 12ϕ10 + 6ϕ01 − 8ϕ20 + 6ϕ11 + 5ϕ02 )/10s + 5 χ12 = (75(2ϕ11 + ϕ02 + 2ϕ21 + ϕ12 )s4 + 5(28ϕ00 + 28ϕ10 − 42ϕ01 − 12ϕ20 + 54ϕ11 + 51ϕ02 − 12ϕ30 + 96ϕ21 + 51ϕ12 )s3 + (−140ϕ00 + 532ϕ10 + 126ϕ01 − 68ϕ20 + 36ϕ11 + 140ϕ02 − 180ϕ30 + 540ϕ21 + 315ϕ12 )s2 + (−140ϕ00 + 280ϕ10 + 168ϕ01 + 100ϕ20 − 12ϕ11 − 10ϕ02 − 180ϕ30 + 240ϕ21 + 165ϕ12 )s + 5(16ϕ20 + 6ϕ11 − ϕ02 − 12ϕ30 + 6ϕ21 + 6ϕ12 ))/15(s + 1)3 (5s + 2) χ21 = (−240(ϕ11 + ϕ21 )s5 + (−448ϕ00 − 112ϕ10 + 840ϕ01 + 128ϕ20 − 780ϕ11 − 350ϕ02 + 72ϕ30 − 1200ϕ21 + 105ϕ03 )s4 − 21(24ϕ00 + 56ϕ10 − 68ϕ01 − 24ϕ20 + 12ϕ11 + 30ϕ02 − 16ϕ30 + 108ϕ21 − 17ϕ03 )s3 + (616ϕ00 − 2072ϕ10 − 84ϕ01 + 408ϕ20 + 936ϕ11 + 576ϕ30 − 2004ϕ21 + 441ϕ03 )s2 + (392ϕ00 − 784ϕ10 − 336ϕ01 − 184ϕ20 + 624ϕ11 + 280ϕ02 + 432ϕ30 − 804ϕ21 + 231ϕ03 )s + 2(−80ϕ20 + 30ϕ11 + 35ϕ02 + 60ϕ30 − 54ϕ21 + 21ϕ03 ))/21(s + 1)3 (5s + 2) χ13 = (420(2ϕ31 + ϕ22 + 2ϕ21 + ϕ12 )s5 + 4(84ϕ00 + 228ϕ10 − 80ϕ40 + 870ϕ31 + 455ϕ22 − 54ϕ01 + 108ϕ20 − 324ϕ11 − 116ϕ30 + 600ϕ21 + 455ϕ12 )s4 + (−2352ϕ00 + 2976ϕ10 − 1280ϕ40 + 5520ϕ31 + 3080ϕ22 + 4608ϕ01 + 2736ϕ20 − 2592ϕ11 − 2520ϕ02 − 1352ϕ30 + 1596ϕ21 + 2450ϕ12 + 441ϕ03 )s3 + 6(8ϕ00 − 232ϕ10 − 320ϕ40 + 680ϕ31 + 420ϕ22 + 156ϕ01 + 888ϕ20 + 216ϕ11 − 240ϕ02 − 152ϕ30 − 112ϕ21 + 120ϕ12 + 63ϕ03 )s2 + (720ϕ00 − 1440ϕ10 − 1280ϕ40 + 1320ϕ31 + 980ϕ22 − 864ϕ01 + 2160ϕ20 + 1296ϕ11 + 520ϕ30 − 372ϕ21 − 190ϕ12 + 63ϕ03 )s − 20(16ϕ40 − 6ϕ31 − 7ϕ22 − 20ϕ30 − 6ϕ21 + 2ϕ12 ))/140(s + 1)4 (3s + 1)

p=3

χ22 = (−168(4ϕ31 − ϕ13 + 4ϕ21 − ϕ03 )s5 − 4(168ϕ00 + 240ϕ10 − 40ϕ40 + 708ϕ31 − 182ϕ13 − 216ϕ01 − 486ϕ11 + 135ϕ02 − 112ϕ30 + 438ϕ21 + 135ϕ12 − 182ϕ03 )s4 + 2(1512ϕ00 − 2088ϕ10 + 320ϕ40 − 2304ϕ31 + 616ϕ13 − 2916ϕ01 − 648ϕ20 + 2484ϕ11 + 1350ϕ02 + 752ϕ30 − 180ϕ21 − 540ϕ12 + 175ϕ03 )s3 + 3(48ϕ00 + 288ϕ10 + 320ϕ40 − 1184ϕ31 + 336ϕ13 − 576ϕ01 − 1392ϕ20 + 144ϕ11 + 600ϕ02 + 488ϕ30 + 724ϕ21 + 30ϕ12 − 21ϕ03 )s2 + (−528ϕ00 + 1056ϕ10 + 640ϕ40 − 1248ϕ31 + 392ϕ13 + 432ϕ01 − 1584ϕ20 − 648ϕ11 + 180ϕ02 − 8ϕ30 + 1236ϕ21 + 450ϕ12 − 49ϕ03 )s + 160ϕ40 − 144ϕ31 + 56ϕ13 − 200ϕ30 + 108ϕ21 + 90ϕ12 − 7ϕ03 )/56(s + 1)4 (3s + 1) χ31 = (3360(ϕ31 + ϕ21 )s6 + 14(384ϕ00 + 240ϕ10 − 64ϕ40 + 1440ϕ31 + 54ϕ04 − 648ϕ01 − 648ϕ11 + 540ϕ02 − 136ϕ30 + 1116ϕ21 + 270ϕ12 − 189ϕ03 )s5 + (−6720ϕ00 + 29280ϕ10 − 5216ϕ40 + 48240ϕ31 + 3276ϕ04 + 13968ϕ01 + 4320ϕ20 − 44712ϕ11 + 3780ϕ02 − 9680ϕ30 + 22320ϕ21 + 17640ϕ12 − 4410ϕ03 )s4 − 6(4984ϕ00 − 6248ϕ10 + 1984ϕ40 − 9760ϕ31 − 924ϕ04 − 10452ϕ01 − 4728ϕ20 + 9108ϕ11 + 3990ϕ02 + 2944ϕ30 − 196ϕ21 − 4200ϕ12 − 441ϕ03 )s3 + (144ϕ00 − 11232ϕ10 − 13376ϕ40 + 37440ϕ31 + 4536ϕ04 + 16416ϕ01 + 44208ϕ20 − 5184ϕ11 − 14400ϕ02 − 10856ϕ30 − 18900ϕ21 + 11610ϕ12 + 5985ϕ03 )s2 + (4848ϕ00 − 9696ϕ10 − 7424ϕ40 + 11520ϕ31 + 1764ϕ04 − 3600ϕ01 + 14544ϕ20 + 5400ϕ11 − 1260ϕ02 + 2008ϕ30 − 9180ϕ21 + 1530ϕ12 + 2835ϕ03 )s + 3(−544ϕ40 + 400ϕ31 + 84ϕ04 + 680ϕ30 − 188ϕ21 + 30ϕ12 + 147ϕ03 ))/252(s + 1)4 (3s + 1)

9

3.2

Basis functions for elements of type II

First, we state and prove a theorem for the space Pp . Theorem 3.1. Under the conditions of Theorem 2.1 the leading term Qp+1 of the local DG finite element error on an element of type II defined by the vertices (0, 0), (1, 0) and (0, 1) using the space Pp can be written as Qp+1 (ξ, η) = cLp+1 (2η−1)+

X

cji (1−ξ−η)Pi2j+2,0 (2η−1)(1−η)j Pj1,0 (

i,j≥0 i+j=p

2ξ −1). (3.4) 1−η

Proof. Applying Theorem 2.2 we write Qp+1 (ξ, η) = cLp+1 (2η − 1) + (1 − η − ξ)˜ qp (ξ, η),

(3.5)

where q˜p ∈ Pp . Applying Lemma 2.1 we can write q˜p =

X

cji Pi2j+2,0 (2η − 1)(1 − η)j Pj1,0 (

i,j≥0 i+j≤p

2ξ − 1). 1−η

The orthogonality condition (2.14) yields Z 1 Z 1−η (1 − ξ − η)˜ qp V dξdη = 0, ∀ V ∈ Pp−1 . 0

(3.6)

(3.7)

0

Combining (3.6) and (3.7) establishes (3.4) and completes the proof. Another basis for the leading term of the error on a triangle of type II with the space Pp may be obtained by applying Theorem 2.1 where the local finite error on an element ∆ of type II can be approximated by (2.17). On the outflow edge η = 1 − ξ, Qp+1 is orthogonal to ϕip−i , i = 0, · · · , p, i.,e., Z 1 Qp+1 (ξ, 1 − ξ)ϕip−i (ξ, 1 − ξ)dξ = 0, i = 0, ..., p. (3.8) 0

We note that (2.17) involves 2p + 3 parameters. Using (3.8), the p + 1 coefficients cpi , i = 0, · · · , p, can be expressed in terms of cp+1 , i = 0, · · · , p + 1, thus, Qp+1 can be i written as p+1 X Qp+1 = dp+1 χip+1−i (ξ, η), (3.9) i i=0

χip+1−i ,

computed using Mathematica, are given in Table 2 in terms where the functions i i of ϕp−i and ϕp+1−i . Thus, Qp+1 is determined by the p + 2 parameters, dp+1 , i = 0, · · · , p + 1. i In the following theorem we establish the asymptotic behavior of the DG error for the space Vp . 10

Table 2: Error basis functions for the spaces Pp for p = 0 to 3 on elements of type II. p=0

χ01 = −ϕ00 + ϕ01 χ10 = − 12 ϕ00 + ϕ10

p=1

χ02 = − 32 ϕ01 + 13 ϕ10 + ϕ02 χ11 = − 14 ϕ01 − 56 ϕ10 + ϕ11 χ20 = − 23 ϕ10 + ϕ20

p=2

χ03 χ12 χ21 χ30

= −10ϕ02 + 6ϕ11 − 83 ϕ20 + ϕ03 = −6ϕ02 + ϕ12 = − 52 ϕ11 + ϕ21 = − 43 ϕ20 + ϕ30

p=3

χ04 χ13 χ22 χ31 χ40

5 1 2 3 = − 74 ϕ03 + 14 ϕ2 − 17 ϕ21 + 35 ϕ0 + ϕ04 1 0 45 1 9 2 9 3 = − 8 ϕ3 − 28 ϕ2 + 14 ϕ1 − 35 ϕ0 + ϕ13 3 + ϕ2 = − 27 ϕ12 − 97 ϕ21 + 18 ϕ 2 35 0 = − 12 ϕ21 − 35 ϕ30 + ϕ31 2 4 = − 45 ϕ30 + 35 ϕ0

Theorem 3.2. Under the conditions of Theorem 2.1 the leading term Qp+1 of local DG finite element error on an the reference element of type II defined by the vertices (0, 0), (1, 0) and (0, 1) using the space Vp can be written as (3.4) such that C = [c1p−1 , ..., c0p ]T is solution of the linear system AC = B where ZZ aij = a.∇[(1 − ξ − η)qij ]ξ i η p+1−i dηdξ, i, j = 1, · · · , p, (3.10a) ∆

bi = (−cβωi − cp0 γi ), i = 1, · · · , p, with

(3.10b)

ZZ γi = ∆

a.∇[(1 − ξ − η)qi0 ]ξ i η p+1−i dηdξ, i = 1, · · · , p,

ZZ ωi = ∆

dLp+1 (2η − 1) i p+1−i ξη dηdξ, i = 1, · · · , p. dη

(3.10c) (3.10d)

Proof. Since Pp ⊂ Vp , Qp+1 satisfies (3.4), which is combined with (2.31) to obtain the following algebraic system cβωi +

p X

cp−j j aij = 0,

i = 1, ..., p.

(3.11)

j=0 0 Solving for cp−1 1 , ..., cp , completes the proof.

Another set of basis functions can be found by writing the leading term Qp+1 as p+1 p+1 0 Qp+1 (ξ, η) = cp+1 1 χp+1 (ξ, η) + c2 χ0 (ξ, η),

11

(3.12)

Table 3: Error basis functions for the spaces Vp for p = 1 to 4 on elements of type II, s = α/β.

p=1

χ02 = −sϕ10 − χ20 =

p=2

1 ((4s 6

7ϕ1 0 6



ϕ0 1 4

+

3sϕ2 0 2

+

ϕ2 0 2

+ ϕ11

− 6)ϕ10 − 15ϕ01 − 6sϕ20 + 14ϕ20 + 10ϕ02 )

1 χ03 = 30 (160ϕ30 s2 + 116ϕ30 s + 150ϕ21 s − 4(30s2 + 48s + 19)ϕ20 − 6(10s + 11)ϕ11 − 5ϕ02 + 64ϕ30 + 60ϕ21 + 30ϕ12 ) 1 χ30 = 21 (−32ϕ30 s2 + 92ϕ30 s − 30ϕ21 s + 4(6s2 − 12s − 25)ϕ20 + 6(2s − 5)ϕ11 − 0 35ϕ2 + 40ϕ30 + 96ϕ21 + 21ϕ03 )

p=3

1 χ04 = 56 (−1120ϕ30 s3 − 2320ϕ30 s2 − 840ϕ21 s2 + 490ϕ22 s − 1760ϕ30 s − 1320ϕ21 s− 140ϕ12 s + 40(35s3 + 41s2 + 37s + 11)ϕ40 + 12(140s2 + 115s + 43)ϕ31 + 210ϕ22 + 56ϕ13 − 568ϕ30 − 492ϕ21 − 150ϕ12 − 7ϕ03 ) 1 χ40 = 36 (160ϕ30 s3 − 400ϕ30 s2 + 120ϕ21 s2 − 70ϕ22 s − 992ϕ30 s − 360ϕ21 s + 20ϕ12 s− 8(25s3 − 85s2 − 65s − 63)ϕ40 − 60(4s2 − 15s − 7)ϕ31 + 290ϕ22 + 36ϕ04 − 504ϕ30 − 588ϕ21 − 70ϕ12 − 63ϕ03 )

p=4

1 χ05 = 270 (36288ϕ50 s4 + 59472ϕ05 s3 + 52920ϕ41 s3 − 7560ϕ22 s2 + 73248ϕ50 s2 + 67200ϕ41 s2 + 20160ϕ32 s2 − 11760ϕ22 s − 756ϕ13 s + 41672ϕ50 s + 45612ϕ41 s + 17500ϕ32 s + 3402ϕ23 s− 48(630s4 + 1610s3 + 1785s2 + 1055s + 251)ϕ40 − 24(1260s3 + 2590s2 + 1834s + 507)ϕ31 − 4260ϕ22 − 798ϕ13 − 27ϕ04 + 10128ϕ50 + 11592ϕ41 + 5640ϕ32 + 1512ϕ23 + 270ϕ14 ) 1 χ50 = 165 (−4032ϕ50 s4 + 15792ϕ50 s3 − 5880ϕ41 s3 + 840ϕ22 s2 + 19328ϕ50 s2 + 25200ϕ41 s2 − 2240ϕ32 s2 − 3360ϕ22 s + 84ϕ13 s + 27192ϕ50 s + 22932ϕ41 s + 10500ϕ32 s − 378ϕ23 s + 48(70s4 − 210s3 − 635s2 − 605s − 261)ϕ40 + 24(140s3 − 490s2 − 1074s − 477)ϕ31 − 4860ϕ22 − 378ϕ13 − 297ϕ04 + 9648ϕ50 + 12312ϕ41 + 5040ϕ32 + 1932ϕ23 + 165ϕ05 )

12

where χ0p+1 and χp+1 are given in Table 3 in terms of ϕip−i , ϕip+1−i and s = α/β. 0 These basis functions are obtained by combining (3.9) and (2.31). In the following theorem we establish the asymptotic behavior of the DG error for the space Up . Theorem 3.3. Under the conditions of Theorem 2.1 the leading term Qp+1 of local DG finite element error on an the reference element of type II defined by the vertices (0, 0), (1, 0) and (0, 1) using the space Up can be written as (3.4) where λp (cµ − (p + 1)s

p−1 P i=1

cp0 =

p−1 P i=1

cp−i i λi )

s(p + 1)(κ0 λp − κp λ0 )

−λ0 (cµ − (p + 1)s

p−1 P i=1

c0p = with

cp−i i κi ) − sκp (c(s + p + 2)ν − (p + 1)

cp−i i κi ) + sκ0 (c(s + p + 2)ν − (p + 1) s(p + 1)(κ0 λp − κp λ0 )

p−1 P i=1

,

(3.13a)

cp−i i λi ) ,

(3.13b)

Z 1 ZZ α p+1 (1 − ξ − η)qip−i ξ p dηdξ Lp+1 (1 − 2ξ)ξ dξ, κi = s= , µ= β ∆ 0 Z 1 ZZ ν= Lp+1 (2η − 1)η p+1 dη, λi = (1 − ξ − η)qip−i η p dηdξ. 0



Proof. First we note that the leading term Qp+1 of the error satisfies the orthogonality conditions [2] Z ZZ a · nQp+1 V ds − a · ∇V Qp+1 dξ dη = 0, ∀ V ∈ Up . (3.14) Γ+



Testing against V = ξ p+1 , η p+1 leads to Z 1 Z p+1 (α + β) Qp+1 (ξ, 1 − ξ)ξ dξ − α(p + 1) 0

and

Z

1 0

(α + β)

Qp+1 (1 − η, η)η

p+1

1−ξ

1

Z

1−ξ

dη − β(p + 1)

0

Qp+1 (ξ, η)ξ p dηdξ = 0,

(3.15)

Qp+1 (ξ, η)η p dηdξ = 0.

(3.16)

0

Z

1

Z

0

0

Combining these orthogonality conditions with (3.4) yields the algebraic system Z 1 ZZ p X p−i p+1 cβ Lp+1 (1 − 2ξ)ξ dξ − α(p + 1) ci (1 − ξ − η)qip−i ξ p dηdξ = 0, (3.17) 0

i=0



and Z

1

c(α + (p + 2)β)

Lp+1 (2η − 1)η

p+1

dη − β(p + 1)

0

Solving for

cp0

and

p X i=0

c0p

completes the proof. 13

ZZ cip−i



(1 − ξ − η)qip−i η p dηdξ = 0. (3.18)

We construct another set of shape functions for Qp+1 in terms of ϕip−i and ϕip+1−i . Noting that (3.9) still holds and using (2.23) and (2.24), the leading term of local finite error on a reference element of Type II can be written as Qp+1 =

p X

dp+1 χip+1−i (ξ, η), i

(3.19)

i=1

where the χip+1−i , computed using Mathematica, are given in Table 4 in terms of s = α/β, ϕip−i and ϕip+1−i . We note that Qp+1 in (3.19) is determined by p parameters. Table 4: Error basis functions for the spaces Up for p = 1 to 4 on elements of type II, s = α/β. p=1

p=2

χ11 = ((−24ϕ10 − 12ϕ01 + 16ϕ20 + 18ϕ11 + 5ϕ02 )s2 − 3(8ϕ10 + 12ϕ01 + 8ϕ20 − 18ϕ11 − 5ϕ02 )s − 24ϕ01 + 6ϕ11 + 15ϕ02 )5(s2 + 3s + 3) χ12 = −((32ϕ20 + 18ϕ11 + ϕ02 − 24ϕ30 − 24ϕ12 − 6ϕ03 )s2 + 4(8ϕ20 + 18ϕ11 + ϕ02 + 8ϕ30 − 24ϕ12 − 6ϕ03 )s + 6(10ϕ11 + ϕ02 − 4ϕ12 − 6ϕ03 ))6(s2 + 4s + 6) χ21 = ((320ϕ20 + 90ϕ11 − 35ϕ02 − 240ϕ30 − 204ϕ12 + 21ϕ21 )s2 + 4(80ϕ20 + 90ϕ11 − 35ϕ02 + 80ϕ30 − 204ϕ12 + 21ϕ21 )s + 6(10ϕ11 − 35ϕ02 − 4ϕ12 + 21ϕ21 ))21(s2 + 4s + 6) χ13 = ((32ϕ40 + 30ϕ31 + 7ϕ22 − 40ϕ30 − 24ϕ12 − 2ϕ03 )s2 − 5(8ϕ40 − 30ϕ31 − 7ϕ22 + 8ϕ30 + 24ϕ12 + 2ϕ03 )s + 10(6ϕ31 + 7ϕ22 − 2(6ϕ12 + ϕ03 )))7(s2 + 5s + 10)

p=3

χ22 = ((−800ϕ40 − 624ϕ31 + 56ϕ13 + 1000ϕ30 + 348ϕ12 − 90ϕ03 − 7ϕ21 )s2 + 5(200ϕ40 − 624ϕ31 + 56ϕ13 + 200ϕ30 + 348ϕ12 − 90ϕ03 − 7ϕ21 )s − 10(24ϕ31 − 56ϕ13 − 48ϕ12 + 90ϕ03 + 7ϕ21 ))56(s2 + 5s + 10) χ31 = ((2656ϕ40 + 2000ϕ31 + 84ϕ04 − 3320ϕ30 − 1012ϕ12 + 30ϕ03 − 147ϕ21 )s2 − 5(664ϕ40 − 2000ϕ31 − 84ϕ04 + 664ϕ30 + 1012ϕ12 − 30ϕ03 + 147ϕ21 )s+ 10(8ϕ31 + 84ϕ04 − 16ϕ12 + 30ϕ03 − 147ϕ21 ))84(s2 + 5s + 10) χ14 = ((−48ϕ40 − 30ϕ31 − 3ϕ22 + 40ϕ50 + 36ϕ41 + 8ϕ32 )s2 − 6(8ϕ40 + 30ϕ31 + 3ϕ22 + 8ϕ50 − 36ϕ41 − 8ϕ32 )s − 15(14ϕ31 + 3ϕ22 − 8(ϕ41 + ϕ32 )))8(s2 + 6s + 15)

p=4

χ23 = ((720ϕ40 + 282ϕ31 − 55ϕ22 − 8ϕ13 − 600ϕ50 − 444ϕ41 + 36ϕ23 )s2 + 6(120ϕ40 + 282ϕ31 − 55ϕ22 − 8ϕ13 + 120ϕ50 − 444ϕ41 + 36ϕ23 )s+ 15(42ϕ31 − 55ϕ22 − 8ϕ13 − 24ϕ41 + 36ϕ23 ))36(s2 + 6s + 15) χ32 = ((−4176ϕ40 − 1434ϕ31 + 55ϕ22 − 154ϕ13 − 9ϕ04 + 3480ϕ50 + 2460ϕ41 + 90ϕ14 )s2 − 6(696ϕ40 + 1434ϕ31 − 55ϕ22 + 154ϕ13 + 9ϕ04 + 696ϕ50 − 2460ϕ41 − 90ϕ14 )s− 15(42ϕ31 − 55ϕ22 + 154ϕ13 + 9ϕ04 − 24ϕ41 − 90ϕ14 ))90(s2 + 6s + 15) χ41 = ((47376ϕ40 + 15834ϕ31 − 55ϕ22 + 154ϕ13 − 891ϕ04 − 39480ϕ50 − 27660ϕ41 + 495ϕ05 )s2 + 6(7896ϕ40 + 15834ϕ31 − 55ϕ22 + 154ϕ13 − 891ϕ04 + 7896ϕ50 − 27660ϕ41 + 495ϕ05 )s + 15(42ϕ31 − 55ϕ22 + 154ϕ13 − 891ϕ04 − 24ϕ41 + 495ϕ05 ))495(s2 + 6s + 15)

14

3.3

Basis functions for elements of type III

In this section we state and prove results on the leading term of the local error by characterizing the optimal polynomial space and determining a basis for the leading term of the local error. Theorem 3.4. Under the conditions of Theorem 2.1 let Qp+1 be the leading term of local DG finite element error on an element of type III defined by the vertices (0, 0), (1, 0) and (0, 1) with α > 0 and β = 0. Then, For the space Pp we have X

Qp+1 = cLp+1 (2η − 1) +

cji (1 − ξ − η)Pi2j+2,0 (2η − 1)(1 − η)j Pj1,0 (

i,j≥0 i+j=p

2ξ − 1). (3.20) 1−η

For the space Up we have Qp+1 =

X

cji (1 − ξ − η)Pi2j+2,0 (2η − 1)(1 − η)j Pj1,0 (

i,j≥0 i+j=p

with

2ξ − 1), 1−η

(3.21a)

p

cp0

1 X = (−1)i+1 (p + 1 − i)cp−i i . p + 1 i=1

(3.21b)

For the space Vp we have Qp+1 = c1 Lp+1 (2η − 1) + c2 (1 − η)p+1 Rp+1 (

2ξ − 1). 1−η

(3.22)

Proof. The proof of (3.20) is the same as (3.4). Using Q(ξ, 1 − ξ) = 0 we obtain (3.21a). We apply (2.26) for k = p to prove (3.21b). Next, we note that for the space Vp (3.4) holds which in turn by the orthogonality condition Z 1 Z 1−η ∂V (1 − ξ − η)˜ qp (ξ, η) = 0, ∀ V ∈ Vp ∂ξ 0 0 leads to Qp+1 = c1 Lp+1 (2η − 1) + c˜2 (1 − ξ − η)(1 − η)p Pp1,0 (

2ξ − 1). 1−η

(3.23)

We complete the proof of (3.22) by noting that Rp+1 (x) = (1 − x)Pp1,0 (x). Finally, we note that the size of the error problem with an optimal basis is equal to the number of terms missing in the finite element space to have Pp+1 .

4

A posteriori error estimation

Here we present weak finite element formulations to compute a posteriori error estimates. 15

4.1

Error estimation for linear problems

We illustrate our procedure on meshes consisting of triangles of type II or III with Up where the finite element error on each element ∆ is approximated as e(x, y) ≈ E(x, y) =

p X

dp+1 χip+1−i (ξ(x, y), η(x, y)), i

(4.1)

i=1

where (ξ(x, y), η(x, y)) is the linear transformation from the physical triangle to the standard triangle. Replacing u in (2.1) by e + U and integrating on the physical element ∆ leads to ZZ ZZ (a.∇(U + e) + c(U + e))W dxdy = f W dxdy. (4.2) ∆



Replacing e by E in (4.2) leads to the p × p linear algebraic system for dp+1 , i = 1, ..., p i ZZ ZZ i (a.∇(U + E) + c(U + E))χp+1−i dxdy = f χip+1−i dxdy, i = 1, 2, · · · , p. (4.3) ∆



This local weak formulation is applied to estimate the error for other polynomial spaces and types of elements using the appropriate error basis functions given in §3. An accepted efficiency measure of a posteriori error estimates is the effectivity index. In this paper we use the local effectivity indices in the L2 norm θi =

||E||L2 (∆i ) , ||e||L2 (∆i )

(4.4)

θ =

||E||L2 (Ω) . ||e||L2 (Ω)

(4.5)

and the global effectivity index

Ideally, the effectivity indices should approach unity under mesh refinement.

4.2

Error estimation for nonlinear problems

We also consider nonlinear hyperbolic problems of the form uy + h(u)x = f (x, y),

(4.6)

subject to appropriate boundary conditions. The discrete DG method consists of finding U ∈ Up such that Z Z ZZ ZZ ˆ ˆ [h(U ), U ].nV ds + [h(U ), U ].nV ds − [h(U ), U ].∇V dxdy = f V dxdy, Γ−

Γ+



∀ V ∈ Up .



(4.7)

For smooth solutions the results of the previous section are still true. For instance, on an element of type II using Vp we have the following theorem. 16

Theorem 4.1. Under the conditions of Theorem 2.1 the leading term Qp+1 of local DG finite element error on an element of type III using the space Vp can be written as Qp+1 (ξ, η) = c1 Lp+1 (2η − 1) + c2 (1 − η)p+1 Rp+1 (

2ξ − 1). 1−η

(4.8)

Proof. First, we linearize (4.6) about the true solution u, find the DG orthogonality and expand the error into a power series. From this point we follow the same line of reasoning as for linear problems. The a posteriori error estimation problem on an element of type II and III using Pp consists determining p+1 X E= ckp+1−k χkp+1−k , k=0

such that ZZ ZZ 0 k [h (U ), 1] · ∇(U + E)χp+1−k dxdy = f χkp+1−k dxdy, ∆

k = 1, · · · , p.

(4.9)



We also use the solution of the linearized problem (4.9) as an initial guess for Newton iteration when solving the nonlinear finite element problem for E ZZ

ZZ 0

[h (U + E), 1] · ∇(U + ∆

E)χkp+1−k dxdy

= ∆

f χkp+1−k dxdy,

k = 1, · · · , p. (4.10)

For instance, on elements of type III using the space Vp , we determine p+1 + c0p+1 χ0p+1 , E = cp+1 0 χ0

such that

ZZ

ZZ 0

[h (U ), 1] · ∇(U + ZZ



= ZZ

0

[h (U ), 1] · ∇(U + ∆

E)χp+1 0 dxdy E)χ0p+1 dxdy



= ∆

f χp+1 0 dxdy,

(4.11a)

f χ0p+1 dxdy,

(4.11b)

where χp+1 0 (x, y) = Lp+1 (2η(x, y) − 1),

(4.11c)

χ0p+1 (x, y) = (1 − η(x, y))p+1 Rp+1 (

2ξ(x, y) − 1). 1 − η(x, y)

(4.11d)

We also use the solution of the linearized problem (4.9) as an initial guess for Newton iteration when solving the nonlinear finite element problem for E ZZ ZZ p+1 0 [h (U + E), 1] · ∇(U + E)χ0 dxdy = f χp+1 (4.12a) 0 dxdy, ∆

Z Z ∆



Z Z

[h0 (U + E), 1] · ∇(U + E)χ0p+1 dxdy = 17



f χ0p+1 dxdy.

(4.12b)

Table 5: L2 errors and global effectivity indices for Example 5.1 using uniform meshes having N = 32, 72, 128, 200, elements and the spaces Pp , Vp and Up . Pp N

5

p=0

p=1

32 72 128 200

||e||L2 6.8948e-3 3.1080e-3 1.7607e-3 1.1317e-3

θ 1.0256 1.0159 1.0116 1.0091

||e||L2 7.3089e-6 1.4556e-6 4.6249e-7 1.8991e-7

32 72 128 200

6.8948e-3 3.1080e-3 1.7607e-3 1.1317e-3

1.0256 1.0159 1.0116 1.0091

7.3089e-6 1.4556e-6 4.6249e-7 1.8991e-7

32 72 128 200

p=1 7.3089e-6 1.0162 1.4556e-6 1.0100 4.6249e-7 1.0072 1.8991e-7 1.0056

θ 1.0162 1.0100 1.0072 1.0056 Vp 1.0152 1.0089 1.0061 1.0046 Up

p=2 7.3089e-6 1.0161 1.4556e-6 1.0098 4.6249e-7 1.0071 1.8991e-7 1.0055

p=2 ||e||L2 3.2879e-9 2.9026e-10 5.1809e-11 1.3648e-11

θ 1.0127 1.0076 1.0054 1.0043

p=3 ||e||L2 3.7335e-13 1.4673e-14 1.4744e-15 2.4837e-16

θ 1.0138 1.0087 1.0064 1.0031

3.2879e-9 2.9026e-10 5.1809e-11 1.3648e-11

1.0117 1.0066 1.0044 1.0033

4.0145e-13 1.5778e-14 1.5854e-15 2.6706e-16

1.0128 1.0077 1.0054 1.0021

p=3 3.2879e-9 2.9026e-10 5.1809e-11 1.3648e-11

1.0126 1.0075 1.0053 1.0042

p=4 3.8490e-13 1.5127e-14 1.5200e-15 2.5605e-16

1.0137 1.0086 1.0063 1.0030

Computational examples

Example 5.1. We consider the linear problem ux + uy + u = f (x, y),

(x, y) ∈ [0, 1]2 ,

(5.1a)

subject to the boundary conditions u(x, 0) = g0 (x), u(0, y) = g1 (y).

(5.1b)

We select f (x, y), g0 and g1 such that the exact solution is u(x, y) = ex−y .

(5.1c)

We solve (5.1) on uniform meshes obtained by partitioning the domain into n × n squares with n = 4, 6, 8, 10, and dividing each square into two triangles using the diagonals parallel to y = x. Thus, the meshes have N = 32, 72, 128, 200 triangles of type III. We use the DG finite element method (2.3), (2.4) with the spaces Pp , Vp and Up , p = 0, 1, 2, 3. On each element we apply the error estimation procedure (4.3) to compute error estimates and present the local effectivity indices in Figure 1. The gobal effectivity indices shown in Table 5 indicate that, for smooth solutions, our a posteriori error estimates converge to the true error under both h and p refinements. As shown in [2] the DG solution is O(hp+2 ) superconvergent at Legengre points shifted to the outflow edges of the mesh. Example 5.2. We consider the linear hyperbolic problem 2ux + uy = f (x, y), 18

(x, y) ∈ [0, 1]2 ,

(5.2a)

1

1

1

0.99 0.99 0.99 0.99 0.99 0.97 0.98 0.98 0.98 0.98 0.98 0.99

0.9 0.8

0.9 0.8

0.99 0.99 0.99 0.99 0.97 0.97 0.98 0.98 0.98 0.98 0.98 1

0.7

0.99 0.99 0.99 0.97 0.97 0.98 0.98 0.98 0.98 1

0.6

0.99 0.99 0.99 0.99 0.98 0.99 0.99 0.99 0.99 0.99

0.99

0.99 0.99 0.99 0.98 0.98 0.99 0.99 0.99 0.99 1 0.99 0.99 0.98 0.98 0.99 0.99 0.99 1

0.98

0.99

0.99 0.98 0.98 0.99 0.99 1

0.98

0.98

0.99

0.98 0.98 0.99 1

0.98

0.99

1

0.7

0.97

0.99

1

0.6 1

1

0.5

0.5 0.99 0.99 0.97 0.97 0.98 0.98 0.98 1

0.4 0.3 0.2 0.1

0.99 0.97 0.97 0.98 0.98 1

0.97

0.97 0.97 0.98 1

0.97

0.97

0.97

0.97 1

1

0.4

1

0.3

0.97 1

0.97 1

1

0.97 1

1 0.97

1

1

0.2

0.98 1

1

1 0.98

1

1

0.98 0.98 0.98 0.98 0.98 0.98 0.99 1.01 1.01 1.01 1.01 1.01

0.1 1

0

0 0

0.2

0.4

0.6

0.8

0

1

0.2

0.4

0.6

0.8

1

1

1 1 0.9

1 0.99

0.8

1

1 0.99

0.7 1

0.6

1 0.99 1 0.99

1 0.99

1 0.99

0.99 0.99 0.99 0.99

0.99

0.98 0.98 0.99 1

1 0.99

1 0.99

0.99

0.99

0.98 0.98 0.99 1

1 0.9 0.8

0.99 0.99 0.99 0.99 0.98 0.99 0.99 0.99 0.99 1

0.99

0.99 0.99 0.99 0.97 0.97 0.99 0.99 0.99 0.99 1.01

0.99

0.99 0.99 0.97 0.97 0.97 0.99 0.99 0.99 1.01 1.01

0.99

0.99 0.97 0.97 0.97 0.97 0.99 0.99 1.01 1.01 1.01

0.99

0.97 0.97 0.97 0.97 0.97 0.99 1.02 1.02 1.02 1.02

1

0.7

0.98

0.99

1

0.6 1

0.5

0.5 1 0.4

1 0.99

0.3

1 0.99

0.2

0.99

0.98 0.98 0.99 1

0.98

1

0.98 1

0.4

1

0.3

0.98 0.98 0.98 0.98 0.98 0.99 1.01 1.01 1.01 1.01

0.98 0.98 0.98 0.98 0.98 0.98 0.99 1.01 1.01 1.01 1.01 1.01

0.1

1

0.2

0.96 0.96 0.96 0.96 0.96 0.96 0.99 1.02 1.02 1.02 1.02 1.02

0.1 0

0 0

0.2

0.4

0.6

0.8

0

1

0.2

0.4

0.6

0.8

1

Figure 1: Local effectivity indices on a uniform mesh having N = 72 elements with Pp , p = 0 to 3 (upper left to lower right) for Example 5.1.

19

subject to the boundary conditions u(x, 0) = g0 (x), u(0, y) = g1 (y).

(5.2b)

We select f (x, y), g0 and g1 such that the exact solution is u(x, y) = ex+y .

(5.2c)

We create uniform triangular meshes obtained by partitioning the domain into n × n squares with n = 4, 6, 8, 10, 15, 20, 25 squares and dividing each square into two triangles by connecting the upper-left and lower-right vertices. These meshes consist of N = 32, 72, 128, 200, 450, 800, 1250 triangles of type I and II. First, we solve (5.2) using the standard DG method (2.3), (2.4) with the spaces Pp , Up and Vp and present the maximum errors at the superconvergence points over all elements in Table 6 which show only O(hp+1 ) convergence rates. These results indicate that the local superconvergence results of [2] do not hold globally on meshes having elements of type I and II since the convergence rates at the superconvergence points are the same as the global L2 convergence rates. Table 6: Maximum errors and orders of convergence at the superconvergence points for Example 5.2 on uniform meshes having N = 32, 72, 128, 200 elements. Pp N

p=1

32 72 128 200

||e||∗∞ 3.3166e-2 1.5485e-2 8.9288e-3 5.8001e-3

32 72 128 200

4.9642e-2 2.3084e-2 1.3269e-2 8.5970e-3

32 72 128 200

7.3522e-3 3.3148e-3 1.8793e-3 1.2089e-3

order 1.8784 1.9139 1.9333 Vp 1.8884 1.9247 1.9449 Up 1.9646 1.9726 1.9771

p=2 ||e||∗∞ order 8.3482e-4 2.5927e-4 2.8840 1.1198e-4 2.9183 5.8186e-5 2.9339 1.4481e-3 4.5006e-4 1.9888e-4 1.0424e-4

2.8822 2.8389 2.8948

2.6398e-5 7.8662e-6 3.3211e-6 1.7008e-6

2.9860 2.9974 2.9989

In order to maintain the O(hp+2 ) superconvergence globally on meshes consisting of elements of type I and II, we used the modified DG method (2.3), (2.5) where the inflow boundary conditions are corrected using our a posteriori error estimate. Now, we solve (5.2) using the same parameters as above with the space Up , p = 1, 2, 3, 4, N = 32, 72, 128, 200 and present the maximum errors and orders of convergence in Table 7. In Figure 2 we plot the zero-level curves of the true error and mark the superconvergence points by x for N = 32 and p = 1, 2, 3, 4, to show that the level curves pass close to the superconvergence points on most elements. These results show that the modified DG method (2.3), (2.5) on Up yields O(hp+2 ) superconvergent solutions at the Legendre points on the outflow edges of elements of type II and at the endpoints of the inflow edges of elements of type I. 20

The local effectivity indices of Figure 3 for N = 32 and Up , p = 1, 2, 3, 4, are close to unity on most elements. Moreover, the true L2 errors and the corresponding global effectivity indices shown in Table 8 indicate that the a posteriori error estimates are asymptotically exact under mesh refinement. 1

1

0.9

0.9

0.8

0.8

0.7

0.7

0.6

0.6

0.5

0.5

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

0

1

1

1

0.9

0.9

0.8

0.8

0.7

0.7

0.6

0.6

0.5

0.5

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

0

1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Figure 2: Zero-level curves of the true error for Example 5.2 on a mesh having N = 32 elements using the spaces Up , p = 1, 2, 3, 4 (upper left to lower right).

Example 5.3. We solve problem (5.2) using the modified DG method (2.3), (2.5) on a uniform mesh having 72 elements and the spaces Up with the nonuniform degree distributions shown in Figure 4, where on the left the degree on each element is higher or equal to the degree on its downwind neighbors and on the right the degree on each element is lower or equal to the degree on its downwind neighbors. In both situations, we use exact boundary conditions. We show the zero-level curves of the true error in Figure 5 which pass close to the superconvergence points on most elements on the left while this not true for the mesh on the right. Consequently, the local effectivity indices shown in left of Figure 6 are close to unity on most elements while on the right they are smaller than unity on elements near the downwind end of the domain. The purpose of this example is to show 21

1

1 0.98

0.9

0.98

0.8 0.7

0.98 0.98

0.98

0.6

0.98

0.98 0.98

0.98 0.98

0.98 0.98

0.98 0.98

0.99

0.9

0.99

0.8 0.7

0.98

0.99

0.6

0.98

0.99 0.99

0.99

0.99 0.99

0.99 0.99

0.99 0.99

0.99 0.99

0.99 0.99

0.5

0.5 0.98

0.4

0.98

0.3 0.2

0.98 0.98

0.97

0.1

0.96

0.98 0.98

0.97 0.97

0.98

0.98 0.97

0.99

0.4

0.98

0.99

0.3 0.2

0.98

0.98

0.1

0.97

0.99 0.99

0.99

0.99 0.99

0.99 0.98

0.99 0.99

0.99 0.98

0.99 0.98

0

0 0

0.2

0.4

0.6

0.8

0

1

0.2

0.4

0.6

0.8

1

1

1 0.98

0.9

0.98

0.8 0.7

0.98 0.97

0.98

0.6

0.98

0.98 0.97

0.98 0.97

0.98

0.98 0.98

0.98

0.9

0.97

0.98

0.8 0.7

0.98

0.98

0.6

0.98

0.98 0.98

0.98

0.98 0.98

0.98 0.98

0.98 0.98

0.98 0.98

0.98 0.98

0.5

0.5 0.98

0.4

0.99

0.3 0.2

0.99 0.98

0.99

0.1

1

0.98 0.98

0.99 1

0.98

0.99 1

0.98

0.4

0.98

0.98

0.3 0.2

0.99

0.97

0.1

1

0.98 0.98

0.98

0.98 0.97

0.98 0.98

0.98 0.98

0.98 0.97

0.98 0.97

0

0 0

0.2

0.4

0.6

0.8

0

1

0.2

0.4

0.6

0.8

1

Figure 3: Local effectivity indices for Example 5.2 on a mesh having N = 32 elements using the spaces Up , p = 1, 2, 3, 4 (upper left to lower right).

Table 7: Maximum errors and orders of convergence at the superconvergence points (marked by × in Figure 2) on meshes having N = 32, 72, 128, 200 elements with the spaces Up , p = 1, 2, 3, 4, for Example 5.2. N 32 72 128 200

p=1 ||e||∗∞ order 2.3849e-3 7.3721e-4 2.8955 3.1584e-4 2.9464 1.6315e-4 2.9603

p=2 ||e||∗∞ order 9.9233e-5 2.0906e-5 3.8411 6.9342e-6 3.8361 2.9342e-6 3.8542

22

p=3 ||e||∗∞ order 1.1025e-5 1.4639e-6 4.9796 3.5114e-7 4.9627 1.1634e-7 4.9505

p=4 ||e||∗∞ order 1.9016e-6 1.7836e-7 5.8369 3.3221e-8 5.8419 9.0300e-9 5.8376

Table 8: L2 errors and global effectivity indices for Example 5.2 on uniform meshes having 32, 72, 128, 200, 450, 800 and 1250 elements using Up , p = 1, 2, 3, 4. N

p=1 ||e||L2 θ 1.7511e-2 0.9813 4.4195e-3 0.9843 1.9707e-3 0.9843 1.1103e-3 0.9850 7.1128e-4 0.9850 4.9426e-4 0.9850 3.8579e-4 0.9850

32 72 128 200 450 800 1250

p=2 ||e||L2 θ 3.2864e-4 0.9883 4.1188e-5 0.9834 1.2216e-5 0.9845 5.1564e-6 0.9848 2.6410e-6 0.9850 1.5287e-6 0.9850 1.1932e-6 0.9850

p=3 ||e||L2 θ 4.5666e-6 0.9711 2.8433e-7 0.9932 5.6123e-8 0.9940 1.7754e-8 0.9946 7.2711e-9 0.9950 3.5064e-9 0.9950 2.7369e-9 0.9952

p=4 ||e||L2 5.0213e-08 1.5488e-09 2.0345e-10 4.8243e-11 1.5801e-11 6.3491e-12 4.9557e-12

θ 0.9914 0.9883 0.9932 0.9938 0.9941 0.9942 0.9943

that nonuniform polynomial degree distributions shown on the left of Figure 4 lead to accurate error estimates while degree distributions shown on the right lead to poor error estimates. 1

1 1

0.9

1

0.8

1 1

2 2

0.7

1 1

2

0.6 2

1

1

1

1

1

1

1

2 2

1 1

1

1

1 1

1

1

1

1

1 0.9

1

0.8

1 1

1 1

2

1 1

1

0.6

1

3 3

2

1

3 3

2

2

1

1

0.5

2 2

1

1

0.7 1

2

1

2

1

3 3

2

1

2

3 3

3 3

2 2

3 3

0.5 3

0.4

3

0.3

2 2

3 3

0.2

3 3

3

0.1 3

2

1

2

1

2

2

2

3 3

1 1

2

3

1 1

2

3

2

1

1 0.4

1

0.3

1 1

2 2

1

1 1

1

0.1

1

1 1

1

1

2 2

1

1

1

1

0

1 1

1

1

0.2 1

1

1

1

1

1 1

1

1

1

2 2

2 2

1 1

1 1

0 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Figure 4: Two uniform meshes for Example 5.3 with nonuniform polynomial degree distributions.

Example 5.4. In this example we show that the modified DG solution is also O(hp+2 ) superconvergent on triangular meshes consisting of all three types of elements. Let us consider the problem (x, y) ∈ [0, 1]2 ,

(5.3a)

u(x, 0) = g0 (x), u(0, y) = g1 (y).

(5.3b)

uy = f (x, y), subject to the boundary conditions

We select f (x, y), g0 and g1 such that the exact solution is u(x, y) = e3x−y . 23

(5.3c)

1

1

0.9

0.9

0.8

0.8

0.7

0.7

0.6

0.6

0.5

0.5

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

0

1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Figure 5: Zero-level curves of the error for Example 5.3 with the p-distributions of Figure 4.

1

1 1.18 0.9

1.29

0.8

0.89 0.79

0.99 1.05

0.7

1.25

1.15

1.02

1 0.86

0.99

0.6

1.04

0.83

0.85

0.67

0.95

1.05

0.78

0.71

1.07

1.06

1

0.78

0.9

0.98

1.13 0.83

0.86

0.87

0.98

1.07

0.8

1.07 0.96

1.05 1.06

0.99

0.98 0.99

0.7

1.13

0.21 0.5

0.99 0.99

0.98

0.6

0.99

0.83

0.03

0.03

0.59

0.11

0.22

0.03

0.51

0.99 0.99

0.02 0.02

0.59

0.99

0.03 0.11

0.22

0.99

0.51

0.01 0.02

0.02 0.02

0.03 0.59

0.03 0.11

0.5

0.5 0.93 0.4

0.97

0.3

0.66 0.95

0.94 0.97

0.2

0.99

1

0.71

0.96 0.97

0.96

0.1

0.68 0.87

0.66

0.97

0.78

0.69

0.95

0.98

0.86

0.67

0.97

0.88

0.4

1.06

1 0.87

0.62

0.98

0.98

1.04

0.3

0.87 0.83

0.98 0.97

0.99

0.99 0.99

0.2

1.16

0.99 0.99

0.99 0.99

0.98

0.1

0.98

0.73

0.99

0.99

0.99

0.99

0.99

0.99

0.99

0.98 0.98

0.22 0.51

0.99

0.99

0.99 0.99

0.99

0.98

0.98

0.03 0.59

0.22 0.51

0.99 0.98

0.99 0.98

0

0 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

0

1

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Figure 6: Local effectivity indices for Example 5.3 with the p-distributions of Figure 4.

24

Following [19] we construct meshes consisting of elements of type I, II and III, obtained by partitioning the domain Ω = [0, 1]2 into N = 2n × (2n + 1) triangles where h = 1/n as shown in Figure 7. 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0

0

0.2

0.4

0.6

0.8

1

Figure 7: Triangulation of Ω for h = 1/6. We solve (5.3) using the modified DG method (2.3), (2.5) on meshes having N = 42, 72, 110, 156 elements with the spaces Up , p = 1, 2, 3, 4. We present the maximum errors and the orders of convergence at the shifted roots of (p + 1)-degree Legendre polynomial on the outflow edge of type II and III in Table 9. We also present the maximum errors and orders of convergence at the endpoints of the inflow edges of elements of type I in Table 10. The results of Tables 9 and 10 indicate that the p−degree modified DG solution is O(hp+2 ) superconvergent at Legendre points on the outflow edges for each element of type II and at the endpoints of inflow edges of elements of type I. In Figure 8 and Table 11 we present the global L2 effectivity indices versus N and p which show that the error estimates converge to the true error under h-refinement. Table 9: Maximum errors and orders of convergence at the roots of (p+1)-degree Legendre polynomial on the outflow edges over all elements of type II for Example 5.4 on meshes having N = 42, 72, 110, 156 elements using Up with p = 1, 2, 3, 4. N 42 72 110 156

p=1 ||e||∗∞ order 7.6186e-2 3.2586e-2 2.9521 1.6913e-2 2.9390 9.8807e-3 2.9482

p=2 ||e||∗∞ order 8.5507e-4 2.7507e-4 3.9424 1.1629e-4 3.8583 5.7284e-5 3.8835

Example 5.5. 25

p=3 ||e||∗∞ order 1.5323e-4 3.7432e-5 4.8991 1.2694e-5 4.8463 5.2571e-5 4.8350

p=4 ||e||∗∞ order 1.8248e-8 3.3048e-9 5.9395 8.6990e-10 5.9815 2.9259e-10 5.9763

Table 10: Maximum errors and rates of convergence at the endpoints of the inflow edge on all elements of type I for Example 5.4 on uniform meshes having N = 42, 72, 110, 156 elements using Up with p = 1, 2, 3, 4. N

p=1 ||e||∗∞

42 72 110 156

3.9097e-2 1.6601e-2 8.7853e-3 5.1994e-3

p=2 order

||e||∗∞

2.9775 2.8519 2.8770

8.1613e-4 2.6449e-4 1.1182e-4 5.5081e-5

order 3.9167 3.8583 3.8835

3.1591e-4 7.5874e-5 2.5232e-5 1.0205e-5

p=1 p=2 p=3 p=4

1.15

1.1

1.1

1.05

1.05

1

0.95

0.9

0.9

0.85

0.85

60

80

100 Number of elements

120

140

order

order

4.9581 4.9338 4.9654

1.4332e-5 2.6136e-6 6.9288e-7 2.2986e-7

5.9155 5.9497 5.9572

N=42 N=72 N=110 N=156

1

0.95

0.8 40

p=4 ||e||∗∞

1.15

Effectivity Index

Effectivity Index

p=3 ||e||∗∞

160

0.8

1

1.5

2

2.5 p

3

3.5

4

Figure 8: Global effectivity indices for Example 5.4 on uniform meshes having N = 42, 72, 110, 156 elements using the spaces Up , p = 1, 2, 3, 4 versus N (left) and p (right).

Table 11: ||e||L2 and global effectivity indices for Example 5.4 versus h and p using Up . p N 42 72 110 156

1 ||e||L2 6.801e-2 6.808e-2 4.425e-2 3.109e-2

2 θ 0.989 0.980 0.986 0.989

||e||L2 2.721e-4 3.626e-4 1.570e-4 7.872e-5

3 θ 1.023 1.019 1.016 1.014

26

||e||L2 3.317e-6 1.251e-6 3.500e-7 1.224e-7

4 θ 1.052 1.049 1.040 1.034

||e||L2 7.294e-9 2.252e-9 4.035e-10 9.856e-11

θ 1.066 1.061 1.048 1.038

We consider the following problem with a contact discontinuity ux + uy = 0,

(x, y) ∈ [0, 1]2 ,

(5.4a)

0 ≤ x ≤ 1,

(5.4b)

subject to the boundary conditions u(x, 0) = e−x , u(0, y) = ey + .25, The exact solution is

( u(x, y) =

0 < y ≤ 1.

(5.4c)

e−x+y + .25 if x < y . e−x+y if x ≥ y

(5.4d)

The true solution has a contact discontinuity along y = x. Therefore, the smoothness assumption of Theorem 2.1 is violated and as a result we expect the a posteriori error estimate to perform poorly near the discontinuity. We solve (5.4) on meshes having 1250, 5000 and 20000 elements with Up , p = 1, 2 and present the local effectivity indices in Figure 9. These computational results indicate that the local effectivity indices on elements away from the discontinuity converge to unity under mesh refinement while they perform poorly on elements near the discontinuity. Since we are not using limiting to suppress spurious oscillations near the discontinuity, the region around the discontinuity where the error is underestimated gets wider as p increases. We will consider a nonlinear example to validate our superconvergence results for nonlinear problems with smooth solutions. Example 5.6. Let us consider the inviscid Burger’s equation (x, y) ∈ [0, 1]2 ,

(5.5a)

u(x, 0) = g0 (x), and u(y, 0) = g1 (x).

(5.5b)

uy + uux = f (x, y), subject to the boundary conditions

We select f , g0 and g1 such that the exact solution is p u(x, y) = 1 + x2 + 5y 2 .

(5.5c)

We solve (5.5) using the modified DG method (4.7), (2.5) on uniform meshes having N = 32, 72, 128, 200, 450, 800, 1250 triangular elements. Since the true solution is positive over the whole domain, our meshes will consist of elements of type I and II. In Table 12 we present the maximum errors and their order of convergence at the shifted roots of (p + 1)-degree Legendre polynomial on the outflow edges on elements of type II and at the endpoints of the inflow edges for elements of type I. These results indicate that the p-degree DG solution is O(hp+2 ) superconvergent at Legendre points on the outflow edges of elements of type II and at the endpoints of the inflow edges of elements of type 27

Figure 9: Local effectivity indices for 5.5 with (N, p) = (1250, 1), (1250,2), (5000, 1), (5000,2), (20000,1), and (20000,2) (upper left to lower right).

28

II. Thus, our local superconvergence results [2] on elements of type I and II extend to the whole mesh. The a posteriori error estimates are computed by solving the linear problem (4.9) and show the local effectivity indices in Figure 10. We also show the global effectivity indices in Table 13. Here, we also applied Newton’s iteration to solve the nonlinear problem (4.10) where we used linear error estimate (4.9) as an initial guess and present the effectivity indices in Table 14. While the effectivity indices for both estimators are within 2% from unity for most meshes and polynomial degrees, the error estimates obtained by solving the linear problem (4.9) are more efficient. Thus, for the remaining nonlinear computational examples we will present numerical results for the linear error estimator (4.9) only. Table 12: Maximum error and its order of convergence at the superconvergence points for problem (5.5) on meshes having N = 32, 72, 128 elements with the spaces Up , p = 1, 2, 3, 4. N 32 72 128 200

p=1 ||e||∗∞ order 8.4554e-4 2.6084e-4 2.9005 1.1330e-4 2.8986 5.9373e-5 2.8959

p=2 ||e||∗∞ order 1.5292e-4 3.0934e-5 3.9413 9.9362e-6 3.9477 4.1362e-6 3.9275

N 32 72 128 200 450 800 1225

p=1 0.9911 0.9941 0.9946 0.9948 0.9949 0.9949 0.9949

p=2 0.9982 0.9932 0.9943 0.9946 0.9948 0.9948 0.9949

p=3 ||e||∗∞ order 6.7141e-8 9.1431e-9 4.9173 2.2028e-9 4.9474 7.3728e-10 4.9050

p=3 0.9808 1.0031 1.0039 1.0045 1.0049 1.0050 1.0052

p=4 ||e||∗∞ order 4.8521e-8 4.4483e-9 5.8932 8.0178e-10 5.9560 2.1447e-10 5.9095

p=4 1.0013 0.9982 1.0031 1.0037 1.0040 1.0041 1.0042

Table 13: Effectivity indices for problem (5.5) using the linearized error estimator (4.9) on meshes having 32, 72, 128, 200, 450, 800 and 1250 elements and Up , p = 1, 2, 3, 4.

Example 5.7. Here we consider the homogeneous inviscid Burger’s equation uy + uux = 0, (x, y) ∈ [0, 1] × [0, b], b > 0,

(5.6a)

subject to the initial conditions 1 g0 (x, 0) = 1 + sin(2πx). 2 29

(5.6b)

1

1 1.01

0.9 1.02

1.02

1.02

1.02

1.02

1.02

1

0.9

1

1.02

1.01

1.01

1.02

1.01

1.01 1.01

0.8

0.8

0.7

0.7 1.01 0.6

1.02

1.02

1.02

1.02

1.01

1

1.02

0.6

1.01

1

1.02

1

1.02

1.01

1.01

1.01

0.5

0.5 1.02

0.4

1.03

1.01

1.03

1

1.03

1

1

0.4

1

1.02

1

1

1.02

1

1.02

1.01

0.3

0.3

0.2

0.2 1.05 0.1

0.93

1.05

1.05

0.92

0.92

1

1.05

0.1

0.91

1

1.02

1

1.02

1

1.02

1.02

0

0 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

0

1

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

1

1 1

0.9

1

1.01

1.01

1.01

1.01

1.01

0.99

0.9 1.01

1

1

1

1.01

1

1.01

1.01

0.8

0.8

0.7

0.7 1 0.6

1.01

1.01

1.01

1.01

0.99

1.01

1

0.6

1

1.01

0.99

1

1.01

1

1.01

1.01

0.5

0.5 1.01

0.4

1.01

1

1.01

1

1

1.02

0.99

0.4 1.02

0.99

0.99

1

1.01

1

1.01

1.01

0.3

0.3

0.2

0.2 1.03 0.1

0.96

1.04

1.04

0.96

0.95

0.99

1.04

0.1

0.95

1.02

0.99

0.99

1.01

1.01

1 1.01

0

0 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

0

1

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Figure 10: Local effectivity indices for problem (5.5) on a uniform mesh having 32 elements with Up , p = 1, 2, 3, 4, (upper left to lower right).

30

Table 14: Global effectivity indices for problem (5.5) using the nonlinear error estimator (4.10) on meshes having 32, 72, 128, 200, 450, 800 and 1250 elements and Up , p = 1, 2, 3, 4. N 32 72 128 200 450 800 1250

p=1 0.9925 0.9945 0.9948 0.9949 0.9949 0.9950 0.9950

p=2 0.9997 0.9938 0.9945 0.9947 0.9948 0.9949 0.9949

p=3 0.9832 1.0037 1.0041 1.0047 1.0050 1.0051 1.0052

p=4 1.0141 0.9986 1.0032 1.0038 1.0040 1.0042 1.0043

and select g1 (0, y) such that the true solution is periodic and forms a shock discontinuity at the point ( π1 , π1 ) which propagates along y = x. First, we apply the modified DG method (4.7), (2.5) to solve this problem on [0, 1]×[0, 0.3] with a smooth solution on uniform meshes having N = 32, 72, 128, 200, 450, 800, 1250 elements of type I and II with Up , p = 1, 2, 3, 4. We compute error estimates by solving (4.9) and present global effectivity indices in Table 15. The computational results show that the effectivity indices converge to unity under h refinement. Table 15: Effectivity indices for problem (5.6) on [0, 1] × [0, 0.3] with the error estimate (4.9) on meshes with N = 32, 72, 128, 200, 450, 800, 1250 elements and Up , p = 1, 2, 3, 4. N 32 72 128 200 450 800 1250

p=1 0.9281 1.0423 1.0157 1.0214 1.0224 1.0230 1.0197

p=2 1.0916 0.9371 0.9636 1.0052 1.0033 1.0136 1.0067

p=3 0.9505 0.8564 0.9561 0.9315 0.9833 0.9878 0.9867

p=4 0.9048 0.8394 0.8499 0.9681 0.9458 0.9958 0.9969

Next, we solve (5.6) on [0, 1] × [0, 0.999] using uniform meshes with N = 32, 72, 128, 200, 450, 800, 1250 and plot the zero-level curves of the discontinuous Galerkin error for N = 32 and p = 1, 2 in Figure 11 with Legendre points for every element of type II and endpoints of inflow edges of elements of type I marked by ×. We observe that the level curves pass close to the superconvergence points on most triangles away from the shock discontinuity. These computational results indicate that the maximum errors at the shifted roots of (p + 1)-degree Legendre polynomial on the outflow edge(for type II) and at the endpoints (for type I) on elements away from the discontinuity converge as O(hp+2 ) under 31

mesh refinement as shown in Table 16. The results shown in Table 17 indicate that superconveregnce is lost near the shock discontinuity. Table 16: Maximum errors and orders of convergence at the superconvergence points in the strip [0, 1] × [0, 0.3] for Example 5.7 on meshes having N = 32, 72, 128, 200 elements with the spaces Up , p = 1, 2, 3, 4. N 32 72 128 200

p=1 ||e||∗∞ 9.9732e-2 3.0301e-2 1.3142e-2 6.8637e-3

p=2 order

||e||∗∞

2.9381 2.9039 2.9108

5.8855e-4 1.1989e-4 3.9004e-5 1.6256e-5

p=3 order

||e||∗∞

3.9240 3.9034 3.9221

9.8014e-6 1.3545e-6 3.3397e-7 1.1273e-7

p=4 order

||e||∗∞

order

4.8810 4.8671 4.8670

5.9773e-7 5.3416e-8 9.7535e-9 2.5949e-9

5.9562 5.9110 5.9338

Table 17: Maximum errors and orders of convergence at the superconvergence points in [0, 1] × [0, 0.999] for Example 5.7 on meshes having N = 32, 72, 128, 200 elements with the spaces Up , p = 1, 2, 3, 4. N 32 72 128 200

p=1 ||e||∗∞ 1.0176e-2 6.5140e-3 3.3835e-3 2.2358e-3

p=2 order

||e||∗∞

1.1002 2.2770 1.8567

1.5465e-3 5.1831e-4 2.4891e-4 1.2656e-4

p=3 order

||e||∗∞

2.6960 2.5496 3.0311

1.9711e-3 5.8724e-4 2.1983e-4 1.0760e-4

p=4 order

||e||∗∞

order

2.9865 3.4155 3.2016

1.4709e-3 2.9563e-4 9.4182e-5 3.5710e-5

3.9572 3.9762 4.3461

Thus, our superconvergence results hold only on elements with a smooth solution and away from the discontinuity. We conclude by solving the previous problem on [0, 1] × [0, 0.999] on two uniform meshes having N = 1800, 20000 elements and Up , p = 1, 2. We plot the local effectivity indices in Figure 12. These computational results indicate that our theoretical results are valid for nonlinear problems in regions where the solution is smooth and away from discontinuities. Thus, the local effectivity indices converge to unity under mesh refinement in regions where the solution is smooth.

6

Conclusion

In this manuscript we have constructed a new discontinuous method for first-order hyperbolic problems in two dimensions with a new finite element space. For smooth solutions the new method exhibits O(hp+2 ) superconvergence at the roots of p + 1-degree Legendre polynomials for elements of type II and III and at the endpoints of inflow edges of elements of type I. We applied these results to construct efficient and accurate a posteriori 32

1

1

0.9

0.9

0.8

0.8

0.7

0.7

0.6

0.6

0.5

0.5

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0

0

0.2

0.4

0.6

0.8

0

1

0

0.2

0.4

0.6

0.8

1

Figure 11: Zero level curves for of the error using p = 1, 2 approximations with 32 elements (left to right) for Example 5.7 on [0, 1] × [0, 0.999].

Figure 12: Local effectivity indices for the homogeneous Burger’s equation 5.7 with initial condition (5.6b) on [0, 1] × [0, 0.999] using meshes having (N, p) = (1800,1), (1800, 2),(20000,1), (20000, 2) (upper left to lower right).

33

estimates of the multi-dimensional discontinuous Galerkin finite element error for hyperbolic problems on triangular meshes. The error estimation procedure presented in this manuscript yields error estimates in several norms and does not require communication across neighboring elements which makes useful for parallel computations. We plan to extend these results to locally refined meshes and unstructured meshes in two and three dimensions for hyperbolic systems. Finally, we note that our error estimates are not accurate near discontinuities.

Acknowledgement This research was partially supported by the National Science Foundation (Grant Number DMS 0511806).

References [1] M. Abramowitz and I. A. Stegun, editors. Handbook of Mathematical Functions. Dover, New York, 1965. [2] S. Adjerid and M. Baccouch. The discontinuous Galerkin method for twodimensional hyperbolic problems Part I: Superconvergence error analysis. Journal of Scientific Computing, 2007. to appear. [3] S. Adjerid, K. D. Devine, J. E. Flaherty, and L. Krivodonova. A posteriori error estimation for discontinuous Galerkin solutions of hyperbolic problems. Computer Methods in Applied Mechanics and Engineering, 191:1097–1112, 2002. [4] S. Adjerid and A. Klauser. Superconvergence of discontinuous finite element solutions for transient convection-diffusion problems. Journal of Scientific computing, 22:5–24, 2005. [5] S. Adjerid and T. C. Massey. A posteriori discontinuous finite element error estimation for two-dimensional hyperbolic problems. Computer Methods in Applied Mechanics and Engineering, 191:5877–5897, 2002. [6] S. Adjerid and T. C. Massey. Superconvergence of discontinuous finite element solutions for nonlinear hyperbolic problems. Computer Methods in Applied Mechanics and Engineering, 195:3331–3346, 2006. [7] M. Ainsworth and J. T. Oden. A Posteriori Error Estimation in Finite Element Analysis. John Wiley, New York, 2000. [8] K. Bottcher and R. Rannacher. Adaptive error control in solving ordinary differential equations by the discontinuous Galerkin method. Tech. report, University of Heidelberg, 1996.

34

[9] P. Castillo. A superconvergence result for discontinuous Galerkin methods applied to elliptic problems. Computer Methods in Applied Mechanics and Engineering, 192:4675–4685, 2003. [10] F. Celiker and B. Cockburn. Superconvergence of the numerical traces for discontinuous Galerkin and hybridized methods for convection-diffusion problems in one space dimension. Math. Comp., 2006. to appear. [11] B. Cockburn, G. E. Karniadakis, and C. W. Shu, editors. Discontinuous Galerkin Methods Theory, Computation and Applications, Lecture Notes in Computational Science and Engineering, volume 11. Springer, Berlin, 2000. [12] B. Cockburn and C. W. Shu. TVB Runge-Kutta local projection discontinuous Galerkin methods for scalar conservation laws II: General framework. Mathematics of Computation, 52:411–435, 1989. [13] B. Cockburn and C. W. Shu. The local discontinuous Galerkin finite element method for convection-diffusion systems. SIAM Journal on Numerical Analysis, 35:2240– 2463, 1998. [14] M. Delfour, W. Hager, and F. Trochu. Discontinuous Galerkin methods for ordinary differential equation. Math. Comp., 154:455–473, 1981. [15] M. Dubiner. Spectral methods on triangles and other domains. Journal of Scientific Computing, 6:345–390, 1991. [16] O. Karakashian and C. Makridakis. A space-time finite element method for the nonlinear Shroedinger equation: The discontinuous Galerkin method. Preprint #969, Department of Mathematics, University of Crete, 71409 Heraklion-Crete, Greece, 1996. [17] L. Krivodonova and J. E. Flaherty. Error estimation for discontinuous Galerkin solutions of two-dimensional hyperbolic problems. Advances in Computational Mathematics, 19:57–71, 2003. [18] P. Lesaint and P. Raviart. On a finite element method for solving the neutron transport equations. In C. de Boor, editor, Mathematical Aspects of Finite Elements in Partial Differential Equations, pages 89–145, New York, 1974. Academic Press. [19] T. Peterson. A note on the convergence of the discontinuous Galerkin method for scalar hyperbolic equation. SIAM Journal on Numerical Analysis, 28:133–140, 1991. [20] D. Sch¨otzau and C. Schwab. An hp a-priori error analysis of the DG time-stepping method for initial value problems. Calcolo, 37:207–232, 2000. [21] D. Sch¨otzau and C. Schwab. Time discretization of parabolic problems by the hpversion of the discontinuous Galerkin finite element method. SIAM Journal on Numerical Analysis, 38:837–875, 2000. 35

[22] R. Verfurth. A Review of a Posteriori Error Estimation and Adaptive Mesh Refinement Techniques. Teubner-Wiley, Teubner, 1996.

36