A posteriori error estimators for locally conservative methods of nonlinear elliptic problems

Applied Numerical Mathematics 57 (2007) 1065–1080 www.elsevier.com/locate/apnum

A posteriori error estimators for locally conservative methods of nonlinear elliptic problems Kwang Y. Kim Department of Applied Mathematics, Sejong University, Seoul 143-747, South Korea Available online 13 November 2006

Abstract In this article we present a unified analysis of some implicit and explicit error estimators for a general class of locally conservative methods applied to nonlinear elliptic problems. The implicit error estimator is computed by solving local Neumann problems, like the Bank–Weiser type, where the Neumann data are provided in a natural way by the locally conservative methods. This can be considered as the equilibrated residual method extended to locally conservative methods which are typically based on discontinuous approximation spaces. We also derive some explicit error estimators of residual type by bounding the implicit error estimator further from above. Our analysis is based on a proper decomposition of the error itself, and seems more direct than other approaches based on the Helmholtz decomposition of the gradient of the error. © 2006 IMACS. Published by Elsevier B.V. All rights reserved. MSC: 65N30; 65N15 Keywords: A posteriori error estimator; Locally conservative method; Nonlinear elliptic problem

1. Introduction In this work we are concerned with developing some implicit and explicit error estimators for a class of discretization schemes known as locally conservative methods which include the P 1 nonconforming finite element method and the discontinuous Galerkin methods. Nowadays a posteriori error estimators and adaptive algorithms based on them are considered to be an indispensable tool for efficient numerical approximation of partial differential equations. This is particularly important for numerical methods employing discontinuous finite element spaces, like the two methods mentioned above, since they involve much more unknowns than standard conforming methods. There are several different approaches to the development of a posteriori error estimators; see the books [4,19] for an overview on this topic. The simplest one is the explicit or residual-type error estimator which evaluates the norms of residuals arising in elements and interfaces. Some related results can be found, for example, in [12,13,18,20] for the P 1 nonconforming finite element method, and we refer to [7,11,14,17] for the discontinuous Galerkin methods. See also [16] for the unified analysis of explicit error estimators of locally conservative mixed methods. While explicit error estimators are most easily computable among others, the most suitable type of an error estimator for locally conservative methods is the one proposed by Bank and Weiser [5] and further studied by Ainsworth and Oden [3]. This E-mail address: [email protected] (K.Y. Kim). 0168-9274/$30.00 © 2006 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.apnum.2006.09.010

1066

K.Y. Kim / Applied Numerical Mathematics 57 (2007) 1065–1080

technique, also called the equilibrated residual method, is based on solving local Neumann problems with Neumann data carefully constructed on the interfaces for conforming finite element spaces. Locally conservative methods, which are originally designed to satisfy the local conservation of mass and thus work well for heterogeneous problems, provide these Neumann data in a natural way from their formulations. An interesting feature of the equilibrated residual method is that it yields a guaranteed computable upper bound on the error measured in the energy norm without involving any generic constants. In this direction some significant results were obtained by Ainsworth [1,2] for the P 1 nonconforming method and the discontinuous Galerkin method of linear elliptic problems. The main objective of this work is to present a unified analysis of some implicit and explicit error estimators for a general class of locally conservative methods applied to nonlinear elliptic problems. Similarly as in [1,2], the implicit error estimator is obtained by extending the equilibrated residual method to locally conservative methods which are typically based on discontinuous approximation spaces. Our analysis, however, takes a more direct approach used in [16] and makes use of a proper decomposition of the error itself, instead of relying on the Helmholtz decomposition of the gradient of the error as is usually done before. We also derive new explicit error estimators by bounding the implicit error estimator further from above which turn out to be simpler than the previous ones. These results are applied to three popular locally conservative methods, namely, the P 1 nonconforming finite element method, the interior penalty discontinuous Galerkin method and the local discontinuous Galerkin method. We remark that a residual error estimator for the third method has been recently obtained in [11] in a different manner. The rest of the paper is organized as follows. We introduce the model problem in the next section, and define the locally conservative methods to approximate its solution in Section 3 along with some notation to be used throughout the paper. In Section 4 we present some implicit and explicit error estimators and establish their reliability within a unified framework. In Section 5 the results of Section 4 are applied to the P 1 nonconforming finite element method, the interior penalty discontinuous Galerkin method and the local discontinuous Galerkin method. 2. Model problem 1

For a given data f ∈ L2 (Ω), uD ∈ H 2 (ΓD ) and gN ∈ L2 (ΓN ), we consider the following nonlinear second order elliptic problem  − div a(x, ∇u) = f in Ω, u = uD a(x, ∇u) · n = gN

on ΓD , on ΓN ,

(1)

where Ω is a bounded domain in R2 with boundary ∂Ω = Γ D ∪ Γ N , ΓD ∩ ΓN = ∅, and n is the unit normal outward to Ω. The standard weak formulation of problem (1) consists of finding u ∈ H 1 (Ω) such that u|ΓD = uD and    a(x, ∇u) · ∇v dx = f v dx + gN v ds ∀v ∈ HD1 (Ω), (2) Ω

Ω

ΓN

where HD1 (Ω) consists of functions in H 1 (Ω) with vanishing traces on ΓD . From now on we use the standard notation H m (G) for the Sobolev space on a domain G consisting of all functions whose derivatives of order up to m are square-integrable, with its norm and seminorm denoted by  · m,G and | · |m,G . The vector-valued function a : Ω × R2 → R2 is assumed to be strongly monotone and Lipschitz continuous, i.e., there exist constants C0 , C1 > 0 such that, for all x ∈ Ω and ξ , η ∈ R2 , 

 a(x, ξ ) − a(x, η) · (ξ − η)  C0 |ξ − η|2 ,   a(x, ξ ) − a(x, η)  C1 |ξ − η|,

(3) (4)

where the dot stands for the Euclidean inner product in R2 and | · | the corresponding vector norm. If a(·, ξ ) has the Jacobian matrix Da(·, ξ ) with respect to the variable ξ , then these inequalities are equivalent to   Da(x, ξ )η  C1 |η|. Da(x, ξ )η · η  C0 |η|2 , (5)

K.Y. Kim / Applied Numerical Mathematics 57 (2007) 1065–1080

1067

We are also interested in the special case a(·, ξ ) = κ(·, |ξ |)ξ for some scalar function κ : Ω × (R+ ∪ {0}) → R+ . In this case we suppose that the following conditions are satisfied: there exist constants C0 , C1 > 0 such that, for all (x, s) ∈ Ω × (R+ ∪ {0}),  ∂  κ(x, s)s  C1 . (6) ∂s It is not difficult to see that the second condition implies the inequalities (3) and (4) with the same upper and lower bounds. C0  κ(x, s)  C1 ,

C0 

3. Locally conservative methods In order to define a locally conservative method  for problem (1), we consider a family {Th }h>0 of triangulations composed of triangular elements such that Ω = T ∈Th T for all h > 0, and the nonempty intersection of any two elements is either a vertex or a complete edge. As usual, nT stands for the unit normal outward to T , and the index h denotes the mesh size of Th defined by h = max hT , T ∈Th

hT = diam(T ).

The family {Th }h>0 is supposed to be shape-regular in the usual sense, i.e., there exists a constant R > 0, hereafter called the shape-regularity constant of {Th }h>0 , satisfying for all h > 0 max

T ∈Th

hT  R, ρT

where ρT is the diameter of the largest ball inscribed in T . Let Eh,I , Eh,D and Eh,N be the collections of all edges of Th belonging to Ω, ΓD and ΓN , respectively, and let Eh = Eh,I ∪ Eh,D ∪ Eh,N . The notation ET is used to denote the set of edges of a given element T . With each edge E ∈ Eh we associate a fixed unit normal nE which is taken outward to Ω for E ∈ Eh,D ∪ Eh,N . We then define the average and the jump of a function v across an interior edge E shared by two elements T + and T − to be {v}E =

v+ + v− , 2

[[v]]E = v + − v − ,

where nE is directed from T + to T − , and v ± denotes the trace of v from the interior of T ± . We also define the jump from T to T across E = ∂T ∩ ∂T

[[v]]T = v|T − v|T . Note that [[v]]T = [[v]]E if nE = nT , and [[v]]T = −[[v]]E otherwise. For the boundary edge E ∈ Eh,D ∩ ET , we set {v}E = [[v]]E = [[v]]T = v|T . For later use it is convenient to define the broken Sobolev space  H m (Th ) := v ∈ L2 (Ω): v|T ∈ H m (T ) ∀T ∈ Th with similar definitions for subsets of Th and Eh as well. We also define the local Sobolev space  HD1 (T ) := v ∈ H 1 (T ): v|∂T ∩ΓD = 0 . Notice that HD1 (T ) = H 1 (T ) if ∂T ∩ ΓD = ∅. Now suppose that we are given a discretization method on a triangulation Th which produces a possibly discontinuous approximation uh ∈ H 1 (Th ) to the solution u of problem (1) or (2). The method is said to be locally conservative if there exist flux functions {

gT ∈ L2 (∂T )}T ∈Th (designed to approximate the normal flux a(x, ∇u) · nT on ∂T ) which are easily computable from the data f, uD , gN and uh , and fulfill the following three conditions for all T , T ∈ Th : (A1)

gT +

gT = 0 on E = ∂T ∩ ∂T . (A2)

gT = gN on E = ∂T ∩ ΓN .

1068

K.Y. Kim / Applied Numerical Mathematics 57 (2007) 1065–1080



 f dx +

(A3) T

gT ds = 0 if ∂T ∩ ΓD = ∅.

∂T

Let us point out that

gT may be undefined on the Dirichlet part ∂T ∩ ΓD . Before closing this section, it is worthwhile to emphasize that the main results presented in the next section do not require any information about the nature of the approximation space for uh nor the formulation of the discretization method, but depend only on the flux functions {

gT ∈ L2 (∂T )}T ∈Th satisfying assumptions (A1)–(A3). 4. A posteriori error estimators The aim of this section is to control the error u − uh in the broken H 1 seminorm ∇ · 0,Th defined by

1/2  ∇v0,Th := ∇v20,T , T ∈Th

and derive some computable error estimators. The main idea is to decompose the error into two components which can be dealt with in a separate way. 4.1. Decomposition of error Following the idea in [16], we decompose the error u − uh as u − uh = (u − φ) + (φ − uh ), where φ ∈ H 1 (Ω) is the function such that φ|ΓD = uD and   a(x, ∇φ) · ∇v dx = a(x, ∇uh ) · ∇v dx ∀v ∈ HD1 (Ω).

(7)

T ∈Th T

Ω

For reasons that become clear soon, we will refer to the first component u − φ as conforming and to the second one φ − uh as nonconforming. The following lemma indicates that we can derive an error estimator for u − uh by considering the two components u − φ and φ − uh separately. Lemma 1. Let φ ∈ H 1 (Ω) be defined by (7). Then there exists a constant C > 0, depending only on C0 , C1 in (3) and (4), such that         ∇(u − uh )  ∇(u − φ) + ∇(φ − uh )  C ∇(u − uh ) . 0,Th

0,Th

0,Ω

0,Th

Proof. The left inequality is a direct consequence of the triangle inequality. For the proof of the right inequality, we note that       a(x, ∇u) − a(x, ∇φ) · ∇(u − φ) dx = a(x, ∇u) − a(x, ∇uh ) · ∇(u − φ) dx, T ∈Th T

Ω

by taking v = u − φ in (7). Hence it follows from the properties (3)–(4) that     ∇(u − φ)  C0−1 C1 ∇(u − uh )0,T . 0,Ω h

Similarly, from the equality       a(x, ∇φ) − a(x, ∇uh ) · ∇(φ − uh ) dx = a(x, ∇φ) − a(x, ∇uh ) · ∇(u − uh ) dx, T ∈T h T

T ∈Th T

we obtain   ∇(φ − uh )

0,Th

   C0−1 C1 ∇(u − uh )0,T .

Thus the proof is completed by taking C

h

= 2C0−1 C1 .

2

K.Y. Kim / Applied Numerical Mathematics 57 (2007) 1065–1080

1069

Now we present two lemmas which provide an important step in deriving error estimators for the conforming and the nonconforming components. Lemma 2. Suppose that {

gT ∈ L2 (∂T )}T ∈Th satisfy assumptions (A1)–(A3), and RT : HD1 (T ) → R is the local residual functional defined by    RT (v) = f v dx +

gT v ds − a(x, ∇uh ) · ∇v dx. (8) T

∂T

T

Let C0 > 0 be the constant in (3). Then we have    T ∈Th RT (v) −1 ∇(u − φ)  C0 sup . 0,Ω ∇v0,Ω v∈H 1 (Ω) D

Proof. By using (A1), (A2) and then (2), (7) one can easily check that, for all v ∈ HD1 (Ω),    RT (v) = f v dx + gN v ds − a(x, ∇uh ) · ∇v dx T ∈Th

Ω

T ∈Th T

ΓN



  a(x, ∇u) − a(x, ∇φ) · ∇v dx.

= Ω

Thus we get the desired result by taking v = u − φ and applying (3).

2

Lemma 3. Let C0 , C1 > 0 be the constants in (3) and (4). Then we have     ∇(φ − uh )  C0−1 C1 inf ∇(χ − uh )0,T . 0,T h

h

χ∈H 1 (Ω)

χ|ΓD =uD

Proof. Let χ ∈ H 1 (Ω) be any function such that χ|ΓD = uD . Then we obtain by taking v = φ − χ in (7)       a(x, ∇φ) − a(x, ∇uh ) · ∇(φ − uh ) dx = a(x, ∇φ) − a(x, ∇uh ) · ∇(χ − uh ) dx, T ∈Th T

T ∈Th T

which proves the assertion by applying (3) and (4).

2

From Lemma 3 we see that the component φ − uh measures how “nonconforming” the approximate solution uh can be; more precisely, the distance from the solution space {v ∈ H 1 (Ω): v|ΓD = uD }. This is exactly why we called it the nonconforming component. 4.2. Error estimator for conforming component Now we are going to derive computable upper bounds for the conforming component u − φ. Let us define ψT ∈ H 1 (T ) to be the solution of the local problem    a(x, ∇ψT ) · ∇v dx = f v dx +

gT v ds ∀v ∈ HD1 (T ) T

T

(9)

∂T

subject to the Dirichlet boundary condition ψT |∂T ∩ΓD = uh if ∂T ∩ ΓD = ∅. Otherwise, (9) is a Neumann problem and admits a unique solution up to an additive constant (which will not affect subsequent results) thanks to assumption (A3). In practical computations, we solve the local problem

1070

K.Y. Kim / Applied Numerical Mathematics 57 (2007) 1065–1080

(9) in an approximate way, e.g., by replacing the space H 1 (T ) by Pk (T ) for sufficiently high k. As usual, Pk (T ) denotes the space of polynomials restricted to T whose degrees are at most k. The following theorem establishes an implicit error estimator of Bank–Weiser type based on the solution of the local problem (9). Theorem 4. Let ψT ∈ H 1 (T ) be defined by (9), and let C0 , C1 > 0 be the constants in (3) and (4). Then, under the assumptions of Lemma 2, we have 

1/2  2   −1   ∇(u − φ) ∇(ψT − uh ) 0,T  C0 C 1 . 0,Ω T ∈Th

Proof. By the definitions (8) and (9), we immediately have for v ∈ HD1 (T )    a(x, ∇ψT ) − a(x, ∇uh ) · ∇v dx. RT (v) =

(10)

T

Applying the Lipschitz continuity (4) and then the Cauchy–Schwarz inequality, we obtain   ∇(ψT − uh ) ∇v0,T RT (v)  C1 0,T T ∈T h

T ∈Th



1/2   ∇(ψT − uh )2  C1 ∇v0,Ω , 0,T T ∈Th

2

which yields the desired result by Lemma 2.

It would be computationally more convenient to establish an error estimator based on a linearized version of the nonlinear local problem (9). For example, one may propose the following linear local problem which corresponds to one iteration of Newton’s method:      Da(x, ∇uh )∇(ψT − uh ) · ∇v dx = f v dx +

gT v ds − a(x, ∇uh ) · ∇v dx, (11) T

T

∂T

T

where the Dirichlet boundary condition T |∂T ∩ΓD = uh , ψ is again imposed if ∂T ∩ ΓD = ∅. The unique solvability of this problem is a consequence of the inequalities (5) and the Lax–Milgram lemma. It turns out that the implicit error estimator based on the linear problem (11) is locally equivalent to that of the nonlinear problem (9), as proved in the following theorem. T ∈ H 1 (T ) be defined by (9) and (11), respectively, and let C0 , C1 > 0 be the constants in (3) Theorem 5. Let ψT , ψ and (4). Then we have       T − uh )  C −1 C1 ∇(ψT − uh ) . C0 C1−1 ∇(ψT − uh )0,T  ∇(ψ 0 0,T 0,T Proof. By the definitions (9) and (11), it is easy to see that, for all v ∈ HD1 (T ),     T − uh ) · ∇v = a(x, ∇ψT ) − a(x, ∇uh ) · ∇v. Da(x, ∇uh )∇(ψ T

T

Now we get the left inequality by taking v = ψT − uh and then applying the strong monotonicity of a in (3) and T − uh and applying the the upper bound for Da in (5). Similarly, the right inequality is obtained by taking v = ψ Lipschtiz continuity of a in (4) and the lower bound for Da in (5). 2

K.Y. Kim / Applied Numerical Mathematics 57 (2007) 1065–1080

1071

Now we derive an upper bound of residual type for the local error estimator ∇(ψT − uh )0,T established in Theorem 4. For this sake, let Π0 v denote the piecewise constant approximation of v ∈ H 1 (Th ) such that   1 Π0 v|T = |T | T v dx if ∂T ∩ ΓD = ∅, 0 if ∂T ∩ ΓD = ∅. It is well known that Π0 v satisfies the approximation estimate  1/2 v − Π0 v20,T + hT v − Π0 v20,∂T  C2 hT ∇v0,T ,

(12)

with C2 depending only on the shape-regularity constant of {Th }h>0 . Theorem 6. Let C0 , C2 > 0 be the constants in (3) and (12). Then we have    2    1/2 ∇(ψT − uh )  C −1 C2 h2 f + div a(·, ∇uh )2 + hT 

gT − a(·, ∇uh ) · nT 0,∂T \Γ . T 0 0,T 0,T D

Proof. The equality (10) obviously gives for v ∈ HD1 (T )    a(x, ∇ψT ) − a(x, ∇uh ) · ∇v dx = RT (v − Π0 v). T

On the other hand, using the integration by parts in (8), one finds that       RT (v) = f + div a(x, ∇uh ) v dx +

gT − a(x, ∇uh ) · nT v ds, T

∂T

from which it follows by Hölder’s inequality and (12) that     gT − a(·, ∇uh ) · nT 0,∂T \Γ v − Π0 v0,∂T RT (v − Π0 v)  f + div a(·, ∇uh )0,T v − Π0 v0,T + 

D 2  2  2 1/2     gT − a(·, ∇uh ) · nT 0,∂T \Γ  C2 hT f + div a(·, ∇uh ) 0,T + hT

∇v0,T . D

Thus the assertion is proved by taking v = ψT − uh and then using the strong monotonicity (3).

2

Remark 7. If we can define the flux function

gT on Dirichlet edges Eh,D so that assumption (A3) is verified on all elements T ∈ Th , then we may consider the Neumann problem    a(x, ∇ψT ) · ∇v dx = f v dx +

gT v ds ∀v ∈ H 1 (T ) T

T

∂T

even on those elements T with ∂T ∩ ΓD = ∅. A careful study of the proofs shows that Theorems 4 and 5 are still gT − a(·, ∇uh ) · nT 20,∂T ∩ΓD should be included in Theorem 6. valid, but the additional contribution hT 

A different form of a residual error estimator can be derived if we can find a vector field σh ∈ H (div; Ω) such that σh · nT ∈ L2 (∂T ), and  (f + div σh ) dx = 0 ∀T ∈ Th , (13) T

 (gN − σh · nE ) ds = 0 ∀E ∈ Eh,N ,

(14)

E

where

  2 H (div; Ω) := τ ∈ L2 (Ω) : div τ ∈ L2 (Ω) .

The constraint σh ∈ H (div; Ω) requires that σh should have continuous normal components across interior edges. Note that assumptions (A1)–(A3) are fulfilled if we define the flux functions by  σh · nT |E if E ∈ Eh,I ∪ Eh,D ,

gT |E := (15) if E ∈ Eh,N . gN

1072

K.Y. Kim / Applied Numerical Mathematics 57 (2007) 1065–1080

In the following theorem we derive a residual estimator which is very similar to the one proposed by Braess and Verfürth [8]. gT defined by (15), and let Theorem 8. Let ψT ∈ H 1 (T ) be the solution of problem (9) with the flux function

C0 , C2 > 0 be the constants in (3) and (12). Then we have      1/2 ∇(ψT − uh )  C −1 σh − a(·, ∇uh ) + C −1 C2 h2 f + div σh 2 + hT gN − σh · nT 2 . T 0,T 0,∂T ∩ΓN 0 0 0,T 0,T Proof. By the definitions of ψT and

gT , we obtain for v ∈ HD1 (T )    a(x, ∇ψT ) − a(x, ∇uh ) · ∇v dx T



 f v dx +

= T

 σh · nT v ds −

∂T

 a(x, ∇uh ) · ∇v dx +

(gN − σh · nT )v ds.

∂T ∩ΓN

T

It then follows by (13)–(14) that    a(x, ∇ψT ) − a(x, ∇uh ) · ∇v dx T



 (f + div σh )(v − Π0 v) dx +

= T

  σh − a(x, ∇uh ) · ∇v dx +

T

Thus the proof is done by taking v = ψT − uh and then using (3), (12).

 (gN − σh · nT )(v − Π0 v) ds.

∂T ∩ΓN

2

At this point let us summarize all the results derived so far. We established an upper bound for the conforming component

1/2    2 ∇(u − φ)  Cη := C η , c c,T 0,Ω T ∈Th

where the local contribution ηc,T is computed by • implicit error estimator based on the nonlinear local problem (9)   ηc,T := ∇(ψT − uh ) ; 0,T

• implicit error estimator based on the linear local problem (11)   T − uh ) ; ηc,T := ∇(ψ 0,T

• explicit error estimator of the classical type 2  2 1/2   gT − a(·, ∇uh ) · nT 0,∂T \Γ ; ηc,T := h2T f + div a(·, ∇uh )0,T + hT 

D

• explicit error estimator of the Braess–Verfürth type 2  1/2 . ηc,T := σh − a(·, ∇uh )0,T + h2T f + div σh 20,T + hT gN − σh · nT 20,∂T ∩ΓN It is worthwhile to mention that the upper bound constant C depends only on the constants C0 , C1 in (3) and (4) for the implicit error estimators, whereas it depends also on the shape-regularity constant of {Th }h>0 through the constant C2 in (12) for the explicit error estimators. One notable thing for the explicit error estimator of the Braess–Verfürth type is that the dependence of C on the shape-regularity constant of {Th }h>0 can be removed if the second and the third terms are of higher order than the first one (cf. Theorem 8), in which case we get a simpler estimator   ηc,T := σh − a(·, ∇uh ) . 0,T

We will elaborate on this fact for the P 1 nonconforming finite element method and the interior penalty discontinuous Galerkin method in Section 5.

K.Y. Kim / Applied Numerical Mathematics 57 (2007) 1065–1080

1073

4.3. Error estimator for nonconforming component Lemma 3 suggests that we can get a computable upper bound for the nonconforming component ∇(φ − uh )0,Th by choosing χ to be a convenient function in H 1 (Ω) satisfying χ|ΓD = uD . A common choice is to construct a continuous piecewise polynomial function by averaging the values of the discontinuous approximation uh at appropriate nodes. For clarity of exposition, let Pk (Th ) and Pk (ED ) denote the spaces of piecewise polynomials on Th and ED , respectively, whose local degrees are at most k, and suppose that uD ∈ C 0 (Γ D ) and uh ∈ Pk (Th ) for some k  1. Let N (S) be the set of all Lagrangian nodes associated with the space Pk (Th ) ∩ C 0 (Ω) which are contained in a domain S. We then define  uh to be the unique function in Pk (Th ) ∩ C 0 (Ω) interpolating  T ∈ωz αz,T uh |T (z) for z ∈ N (Ω) \ N (Γ D ), (16)  uh (z) := for z ∈ N (Γ D ), uD (z) where ωz is the set of all elements in Th sharing z, and the weights {αz,T } form a partition of unity αz,T = 1, 0  αz,T  1. T ∈ωz

For example, one can choose the simple averaging  uh (z) :=

1 uh |T (z) card(ωz ) T ∈ωz

or the weighted averaging  T ∈ωz |T |uh |T (z)   uh (z) := , T ∈ωz |T | where card(S) is the number of elements in a finite set S, and |T | is the measure of T . Although the function  uh defined above does not exactly satisfy the Dirichlet boundary condition χ|ΓD = uD , the following theorem shows that this will cause only a higher order perturbation. Theorem 9. Let  uh ∈ Pk (Th ) ∩ C 0 (Ω) be constructed from uh ∈ Pk (Th ) by interpolating (16), and let  uD ∈ Pk (ED ) ∩ 0 C (Γ D ) be the standard Lagrange interpolant of uD ∈ C 0 (Γ D ). Then there exists a constant C > 0, depending only on the shape-regularity constant of {Th }h>0 , such that     ∇(φ − uh ) uh − uh )0,T + C1 ,  C0−1 C1 ∇( 0,T h

h

where the extra term 1 is defined by

1/2  2 1 := hE |uD −  uD |1,E . E∈ED

Proof. The proof is essentially provided in Lemma 4.7 of [16].

2

Remark 10. The standard approximation estimate yields for uD ∈ H k+1 (ED )

1/2  2 h2k+1 u  , 1  C D k+1,E E E∈ED

which implies that 1 is of higher order than ∇(u − uh )0,Th ∼ hk , and thus negligible. As a corollary of Theorem 9, we can obtain a computable upper bound which does not involve computation of  uh .

1074

K.Y. Kim / Applied Numerical Mathematics 57 (2007) 1065–1080

Corollary 11. Under the assumptions of Theorem 9, there exists a constant C > 0, depending only on the shaperegularity constant of {Th }h>0 and the degree k, such that 

1/2   2  −1  −1 2 2 ∇(φ − uh )  [[u  C h ]] + h u − u  +  . h 0,E h D 0,E 1 E E 0,T h

E∈EI

E∈ED

Proof. The result follows directly from the estimates uD −  uD 0,E  ChE |uD −  uD |1,E and   ∇( uh − uh )

0,Th

C

 E∈EI

∀E ∈ ED

2    h−1 E [[uh ]] 0,E

+



h−1 E uh

1/2 − uD 20,E

.

E∈ED

The first estimate can be readily proved by the one-dimensional integration (notice that uD −  uD vanishes at the two endpoints of each E ∈ ED ). The second estimate is a generalization of Theorem 2.2 in [14] to the cases of nonhomogeneous Dirichlet data lying in Pk (Th ) ∩ C 0 (Ω) and EN = ∅, which can be done in a trivial way. 2 5. Applications In this section we apply the theoretical results obtained in the previous section to three well-known discretization methods. As the general strategy for estimating the nonconforming component described in Section 4.3 applies to all the methods, we mainly focus on estimating the conforming component. To begin with, we briefly review some results on the Raviart–Thomas vector space (cf. [9]). The Raviart–Thomas space of order k is defined as  RTk := τ ∈ H (div; Ω): τ |T ∈ RTk (T ) ∀T ∈ Th , where the local space is given by  2 RTk (T ) := Pk (T ) + xPk (T )

  x = (x1 , x2 ) .

It is now well known that we have for ξh ∈ RTk (T ) div ξh ∈ Pk (T )

and ξh · nT ∈ Pk (ET ),

and that the degrees of freedom for ξh are provided by the moments of order up to k of ξh · nT on ∂T   ξh · nT μ ds: μ ∈ Pk (ET ) ∂T

and the moments of order up to k − 1 of ξh on T (k  1)    2 ξh · τ dx: τ ∈ Pk−1 (T ) . T

For the lowest order k = 0, we can readily construct the basis functions dual to this set of degrees of freedom, yielding the identity ξh =

E∈ET

ξh · nT |E

|E| (x − xE ) 2|T |

∀ξh ∈ RT0 (T ),

(17)

where xE is the vertex of T opposite to the edge E. Throughout the section we will denote by Phk and Qkh the L2 -projections onto Pk (Th ) and Pk (Eh ), respectively.

K.Y. Kim / Applied Numerical Mathematics 57 (2007) 1065–1080

1075

5.1. P 1 nonconforming finite element method With the definition of the Crouzeix–Raviart space    nc P1 (Th ) := vh ∈ P1 (Th ): [[vh ]] ds = 0 ∀E ∈ EI , E

the P 1 nonconforming finite element method for problem (1) reads as follows: find uh ∈ Pnc 1 (Th ) such that uD ) ds = 0 on E ∈ ED and    a(x, ∇uh ) · ∇vh dx = f vh dx + gN vh ds T ∈Th T

Ω



vh ∈ Pnc 1 (Th )



E (uh

−

(18)

ΓN

for all with E vh ds = 0 on E ∈ ED . To see that the method (18) is locally conservative, we define the flux function

gT on E ∈ ET \ EN to be the constant 

 1 a(x, ∇uh ) · ∇φE dx − f φE dx , (19)

gT |E := |E| T

T

φE ∈ Pnc 1 (Th )

where 

is the basis function associated with E, i.e.,  φE ds = |E| and φE ds = 0 if E = E. E

E

gT |E := gN . For a Neumann edge E ∈ ET ∩ EN , we simply set

Now we check assumptions (A1)–(A3). It is easy to see that (A1) follows directly by setting vh = φE in (18) for each E ∈ EI , and that (A2) is obvious by the definition of

gT on ΓN . To check out assumption (A3), we note that, for every E ∈ ET ,   

gT ds = a(x, ∇uh ) · ∇φE dx − f φE dx. (20) E

T

T

For E ∈ ET \ EN , this is obvious by the definition (19). For E ∈ ET ∩ EN , the result is obtained by setting vh = φE in (18) and using the fact that φE ≡ 1 on E. Consequently, it follows that

   

gT ds = a(x, ∇uh ) · ∇φE dx − f φE dx = − f dx ∂T



E∈ET

T

T

T

since E∈ET φE ≡ 1 on T . We can also construct a vector field σh ∈ RT0 whose normal component on E ∈ ET is given by  1 σh · nT |E :=

gT ds, |E|

(21)

E

which, due to (20), can be written as 

 1 σh · nT |E = a(x, ∇uh ) · ∇φE dx − f φE dx . |E| T

(22)

T

It is a simple matter to verify the conditions (13) and (14), which means that div σh = −Ph0 f,

σh · n|ΓN = Q0h gN .

Combining this with Theorems 4 and 8, we can derive the following theorem. Note that if the data f and gN are piecewise smooth, then the extra term 2 below becomes a higher order perturbation and thus negligible.

1076

K.Y. Kim / Applied Numerical Mathematics 57 (2007) 1065–1080

Theorem 12. Let uh ∈ Pnc 1 (Th ) be the solution of (18) and let σh ∈ RT0 be constructed by (21) or (22). Then there exists a constant C > 0 depending only on the shape-regularity constant of {Th }h>0 such that 

1/2 2    −2   ∇(u − φ)  C C − a(·, ∇u ) + C2 , σ 1 h h 0,T 0 0,Ω T ∈Th

where the extra term 2 is due to the data oscillation defined by

1/2     2 2 0 2 0   2 := hT f − Ph f 0,T + hE gN − Qh gN 0,E . T ∈Th

E∈EN

Now we derive an explicit formula for σh which helps to compute the error estimator of Theorem 12 directly from ∇uh and f . Proposition 13. Let σh |T ∈ RT0 (T ) be constructed by (22). Then we have (x − xE )  σh |T = a(·, ∇uh ) − f φE dx, 2|T | E∈ET

(23)

T

where the overline denotes the integral mean value over T . Proof. Since we have |E| nT |E , ∇φE |T = |T | the formula (22) yields    1 σh − a(·, ∇uh ) · nT |E = − f φE dx. |E| T

Hence it follows from (17) that   |E| (x − xE ) σh − a(·, ∇uh ) · nT |E σh − a(·, ∇uh ) = 2|T | E∈ET (x − xE )  =− f φE dx. 2|T | E∈ET

This completes the proof.

T

2

Using the formula (23), we can even find an error estimator for the conforming component which depends only on the source term f , provided a(x, ∇uh ) is piecewise constant on Th . The same result was arrived at by Ainsworth [1] for the linear problem. Theorem 14. Let uh ∈ Pnc 1 (Th ) be the solution of (18) and let σh ∈ RT0 be constructed by (23). Suppose that a(x, ∇uh ) is piecewise constant on Th . Then there exists a constant C > 0 depending only on the shape-regularity constant of {Th }h>0 such that 2 1/2   0   Ph f   −2   ∇(u − φ)  C C ) + C2 , (x − x 1 T 0   0,Ω 2 0,T T ∈Th

where xT is the barycenter of T , and 2 was defined in Theorem 12. Proof. The assertion is proved if we show that  0        σh − a(·, ∇uh )   Ph f (x − xT ) + ChT f − P 0 f  . h  2  0,T 0,T 0,T

K.Y. Kim / Applied Numerical Mathematics 57 (2007) 1065–1080

1077

To do this, we put the following equality into (23)     |T | f − Ph0 f φE dx + Ph0 f f φE dx = , 3 T

T

which gives

(x − xE )    P 0f f − Ph0 f φE dx − h (x − xT ). σh − a(x, ∇uh ) = − 2|T | 2 E∈ET

T

It is elementary to show that the L2 -norm of the first term on the right-hand side is bounded by ChT f − Ph0 f 0,T . 2 We stress again that the explicit error estimators in Theorems 12 and 14 yield an upper bound that depends only on the constants C0 , C1 in (3) and (4). For the nonconforming component one can proceed in a standard way by choosing  uh to be the continuous piecewise linear interpolant whose values at the vertices of Th are computed by averaging the nodal values of uh around those vertices. An alternative choice is to consider the piecewise quadratic interpolant by further specifying the common values of uh at the midpoints of Th . Numerical experiments indicate that the latter choice leads to much better results in some cases. Remark 15. For the efficiency of the error estimators derived above, one may follow the argument of [1,16] to show that for all T ∈ Th ,

1/2  2    2 0 2   ∇(u − uh ) 0,T + hT f − Ph f 0,T ηT  C , (24) T ∈ωT

with C > 0 depending only on the shape-regularity constant of {Th }h>0 . Here ωT denotes the set of all elements sharing an edge with T , and the local error estimator ηˆ T is defined by  0 2  Ph f  2  1 2 2    (x − xT ) + h−1 h−1 (25) ηT :=  E [[uh ]] 0,E + E uh − uD 0,E .  2 2 0,T E∈ET ∩EI

E∈ET ∩ED

5.2. Interior penalty discontinuous Galerkin method For a given k  1 the interior penalty discontinuous Galerkin method consists of finding uh ∈ Pk (Th ) such that, for all vh ∈ Pk (Th ),     a(x, ∇uh ) · nE E [[vh ]]E ds + a(x, ∇uh ) · ∇vh dx − γ h−1 [[uh ]]E [[vh ]]E ds E T ∈Th T

f vh dx +

= Ω

E∈EI ∪ED E





gN vh ds +



γ h−1 E

E∈ED

ΓN

E∈EI ∪ED

 uD vh ds,

E

(26)

E

where γ > 0 is a stabilization parameter to be chosen sufficiently large but independently of the mesh size. The existence and uniqueness of a solution uh for problem (26) as well as some a priori error estimates in the energy norm can be established analogously to [15]. Now we define the flux function as follows: ⎧ if E ∈ ET ∩ EI , ⎨ {a(x, ∇uh ) · nT }E − γ h−1 E [[uh ]]T (27)

gT |E = a(x, ∇uh ) · nT − γ h−1 (uh − uD ) if E ∈ ET ∩ ED , E ⎩ if E ∈ ET ∩ EN . gN Assumptions (A1) and (A2) are trivial, and (A3) is an immediate consequence of the equality    a(x, ∇uh ) · ∇vh dx = f vh dx +

gT vh ds ∀vh ∈ Pk (T ), (28) T

T

∂T

1078

K.Y. Kim / Applied Numerical Mathematics 57 (2007) 1065–1080

which is nothing but the formulation (26) combined with (27) when the support of vh is restricted to an element T . We can also construct a vector field σh ∈ RTk−1 locally on each element T by (see also [6])     2 σh − a(x, ∇uh ) · τ dx = 0 ∀τ ∈ Pk−2 (T ) (k  2), (29) T

 (σh · nT −

gT )μ ds = 0 ∀μ ∈ Pk−1 (ET ).

(30)

∂T

The conditions (13) and (14) are checked out in an obvious way. Moreover, by letting vh ∈ Pk−1 (T ) in (28) we get    a(x, ∇uh ) · ∇vh dx = f vh dx + σh · nT vh ds, (31) T

T

∂T

and thus (by using (29) in the case k  2)  (f + div σh )vh dx = 0,

(32)

T

which implies that div σh = −Phk−1 (f ). We also get σh · nT |∂T = Qk−1 gT from (30). This establishes the following h

theorem in the same way as Theorem 12. Theorem 16. Let uh ∈ Pk (Th ) be the solution of (26) and let σh ∈ RTk−1 be constructed by (29)–(30) with the flux function

gT defined by (27). Then there exists a constant C > 0 depending only on the shape-regularity constant of {Th }h>0 such that 

1/2     σh − a(·, ∇uh )2 ∇(u − φ)  C −2 C1 + C3 , 0,Ω

0

0,T

T ∈Th

where the extra term 3 is due to the data oscillation defined by

1/2     2 k−1 2 k−1 2   3 := hT f − Ph f 0,T + hE gN − Qh gN 0,E . T ∈Th

E∈EN

It is even possible to make the term 3 one order higher by constructing a vector field σh ∈ RTk by     2 σh − a(x, ∇uh ) · τ dx = 0 ∀τ ∈ Pk−1 (T ) , T

 (σh · nT −

gT )μ ds = 0 ∀μ ∈ Pk (ET ).

∂T

It is not difficult to see that Eqs. (31)–(32) now hold for all v ∈ Pk (T ), and thus we can take

1/2  2 2   h2T f − Phk (f )0,T + hE gN − Qkh (gN )0,E . 3 := T ∈Th

E∈EN

In the lowest-order case k = 1, we can derive the same explicit formula (23) for the vector field σh and thus the same result of Theorem 14, which is stated in the following theorem. Theorem 17. Let uh ∈ P1 (Th ) be the solution of (26) in the lowest-order case k = 1 and let σh ∈ RT0 be constructed by (21) with the flux function

gT defined by (27). Suppose that a(x, ∇uh ) is piecewise constant

K.Y. Kim / Applied Numerical Mathematics 57 (2007) 1065–1080

1079

on Th . Then there exists a constant C > 0 depending only on the shape-regularity constant of {Th }h>0 such that 2 1/2   0   Ph f   −2  ∇(u − φ)  (x − x  C C ) + C2 , 1 T  0  2 0,Ω 0,T T ∈Th

where 2 was defined in Theorem 12. Proof. Recall that Theorem 14 was based solely on the formula (23), which was in turn derived from Eq. (22) or (20). Thus we only have to prove Eq. (20) for the flux function (27). This is indeed done by taking vh = φE in (28) and using the fact that φE ≡ 1 on E. 2 Remark 18. In the lowest-order case k = 1, one can also establish a similar result to (24) with the same local error estimator (25) for the mesh-dependent energy norm typically used in a priori error analysis of discontinuous Galerkin methods. We refer to [16] for more details. 5.3. Local discontinuous Galerkin method In the local discontinuous Galerkin (LDG, for short) method for problem (1), one introduces the gradient θ = ∇u and the flux σ = a(x, θ ) as additional unknowns, and seeks the approximate solution (θh , σh , uh ) in the discontinuous finite element space (Pk−1 (Th ))2 × (Pr (Th ))2 × Pk (Th ) with r = k − 1 or k. We refer to [10] for a detailed discussion on this method. For our purpose here we only mention that this method is locally conservative with the flux function  {σ · n } − (β · n )[[σ · n ]] − α[[u ]] if E ∈ E ∩ E ,

gT |E =

h

T E

T

σh · nT − α(uh − uD ) gN

h

T T

h T

T

I

if E ∈ ET ∩ ED , if E ∈ ET ∩ EN ,

where the auxiliary functions α (scalar) and β (vector) are single-valued on each E ∈ Eh and are to be chosen appropriately. A residual error estimator for the LDG method was derived in Theorem 4.1 of [11], which can be stated as

1/2     2 2 ∇(u − uh ) ∇(  u  C η + − u ) , h h 0,T c,T 0,T h

T ∈T h

h

where the local error estimator for the conforming component is given by 2   2 2 gT − a(·, θh ) · nT 0,∂T \Γ + ∇uh − θh 20,T . ηc,T = h2T f + div a(·, θh )0,T + hT 

D

On the other hand, our result gives a simpler one 2   2 2 gT − a(·, ∇uh ) · nT 0,∂T \Γ . ηc,T = h2T f + div a(·, ∇uh )0,T + hT 

D

In particular, we go back to the same error estimator of Theorem 3.2 in [11] for the linear problem, which is not true for the error estimator above. It is possible to construct a vector field σh by (29)–(30) which satisfies the conditions (13) and (14) but not the equalities div σh = −Phk−1 (f ),

σh · nT |∂T = Qk−1 gT . h

Thus the second and third terms in the error estimator of the Braess–Verfürth type cannot be neglected, except in the lowest-order case k = 1. References [1] M. Ainsworth, Robust a posteriori error estimation for nonconforming finite element approximation, SIAM J. Numer. Anal. 42 (6) (2005) 2320–2341.

1080

K.Y. Kim / Applied Numerical Mathematics 57 (2007) 1065–1080

[2] M. Ainsworth, A synthesis of a posteriori error estimation techniques for conforming, non-conforming and discontinuous Galerkin finite element methods, in: Z. Shi, Z. Chen, T. Tang, D. Yu (Eds.), Recent Advances in Adaptive Computation, Contemporary Mathematics, vol. 383, American Mathematical Society, Providence, RI, 2005, pp. 1–14. [3] M. Ainsworth, J.T. Oden, A unified approach to a posteriori error estimation using element residual methods, Numer. Math. 65 (1) (1993) 23–50. [4] M. Ainsworth, J.T. Oden, A Posteriori Error Estimation in Finite Element Analysis, John Wiley & Sons, New York, 2000. [5] R.E. Bank, A. Weiser, Some a posteriori error estimators for elliptic partial differential equations, Math. Comp. 44 (170) (1985) 283–301. [6] P. Bastian, B. Riviére, Superconvergence and H (div) projection for discontinuous Galerkin methods, Internat. J. Numer. Methods Fluids 170 (10) (2003) 1043–1057. [7] R. Becker, P. Hansbo, M.G. Larson, Energy norm a posteriori error estimation for discontinuous Galerkin methods, Comput. Methods Appl. Mech. Engrg. 192 (5–6) (2003) 723–733. [8] D. Braess, R. Verfürth, A posteriori error estimators for the Raviart–Thomas element, SIAM J. Numer. Anal. 33 (6) (1996) 2431–2444. [9] F. Brezzi, M. Fortin, Mixed and Hybrid Finite Element Methods, Springer, Berlin, 1991. [10] R. Bustinza, G.N. Gatica, A local discontinuous Galerkin method for nonlinear diffusion problems with mixed boundary conditions, SIAM J. Sci. Comput. 26 (1) (2004) 152–177. [11] R. Bustinza, G.N. Gatica, B. Cockburn, An a posteriori error estimate for the local discontinuous Galerkin method applied to linear and nonlinear diffusion problems, J. Sci. Comput. 22–23 (2005) 147–185. [12] C. Carstensen, S. Bartels, S. Jansche, A posteriori error estimates for nonconforming finite element methods, Numer. Math. 92 (2) (2002) 233–256. [13] E. Dari, R. Duran, C. Padra, V. Vampa, A posteriori error estimators for nonconforming finite element methods, RAIRO Modél. Math. Anal. Numér. 30 (4) (1996) 385–400. [14] O. Karakashian, F. Pascal, A posteriori error estimates for a discontinuous Galerkin approximation of second-order elliptic problems, SIAM J. Numer. Anal. 41 (6) (2003) 2374–2399. [15] Kwang Y. Kim, Mixed finite volume method for nonlinear elliptic problems, Numer. Methods Partial Differential Equations 21 (4) (2005) 791–809. [16] Kwang Y. Kim, A posteriori error analysis for locally conservative mixed methods, Math. Comp., posted on October 4, 2006, PII S00255718(06)01903-X, in press. [17] R. Riviére, M.F. Wheeler, A posteriori error estimates for a discontinuous Galerkin method applied to elliptic problems, Comput. Math. Appl. 46 (1) (2003) 141–163. [18] F. Schieweck, A posteriori error estimates with post-processing for nonconforming finite elements, M2AN Math. Model. Numer. Anal. 36 (3) (2002) 489–503. [19] R. Verfürth, A Review of a Posteriori Error Estimation and Adaptive Mesh-Refinement Techniques, Wiley–Teubner, Stuttgart, 1996. [20] B. Wohlmuth, A residual based error estimator for mortar finite element discretizations, Numer. Math. 84 (1) (1999) 143–171.