POINTWISE A POSTERIORI ERROR ANALYSIS FOR AN ADAPTIVE PENALTY FINITE ELEMENT METHOD FOR THE

OBSTACLE PROBLEM

DONALD A. FRENCH^{1}, STIG LARSSON^{2}, AND RICARDO H. NOCHETTO^{3}

Abstract. Finite element approximations based on a penalty formulation of the elliptic obstacle problem are analyzed in the maximum norm. A posteriori error estimates, which involve a residual of the approximation and a spatially variable penalty parameter, are derived in the cases of both smooth and rough obstacles. An adaptive algorithm is suggested and implemented in one dimension.

1. Introduction

We consider finite element approximations of the obstacle problem

(1.1)  \begin{aligned}−Δu(x) + β(u(x) − ψ(x)) & ∋ f(x), & x ∈ Ω, \\ u(x) & = 0, & x ∈ ∂Ω, \end{aligned}

where ψ and f are given functions with ψ ≤ 0 on ∂Ω, and where β is the maximal monotone graph defined by

(1.2)  β(s) = \begin{cases}\{0\}, & s > 0, \\ (−∞, 0], & s = 0, \\ ∅, & s < 0. \end{cases}

Our analysis and our finite element method are based on the following penalized (or regularized) form of (1.1): find u_{ϵ} such that

(1.3)  \begin{aligned}−Δu_{ϵ}(x) + ϵ(x)^{−1}(u_{ϵ}(x) − ψ(x))^{−} & = f(x), & x ∈ Ω, \\ u_{ϵ}(x) & = 0, & x ∈ ∂Ω, \end{aligned}

where ϵ = ϵ(x) is a positive function on Ω and s^{−} = min{0, s}. This can also be written

(1.4)  \begin{aligned}−Δu_{ϵ} + β_{ϵ}(u_{ϵ} − ψ) & = f, & in Ω, \\ u_{ϵ} & = 0, & in ∂Ω, \end{aligned}

Date: September 17, 2000. 1991 Mathematics Subject Classification. 65N30, 35J85. Key words and phrases. elliptic obstacle problem, a posteriori error estimate, residual, maxi-

mum norm, adaptive, finite element, penalty method, duality argument.

Published in Comput. Methods Appl. Math. 1 (2001), 18–38,

http://www.emis.de/journals/CMAM .

1Partially supported by the Taft Foundation at the University of Cincinnati. 2Partially supported by the Swedish Research Council for Engineering Sciences (TFR). 3Partially supported by the NSF grants DMS–9623394 and DMS–9971450.

where β_{ϵ} is an approximation of β defined by

(1.5)  β_{ϵ}(s) = \begin{cases}0, & s ≥ 0, \\ s/ϵ, & s ≤ 0. \end{cases}

The weak formulation of (1.4) reads

(1.6)  u_{ϵ} ∈ H_{0}^{1}(Ω) : (∇u_{ϵ}, ∇v) + (β_{ϵ}(u_{ϵ} − ψ), v) = (f, v), ∀v ∈ H_{0}^{1}(Ω).

Our finite element method is based on discretizing (1.6). Let T = {K} be a member of a regular family of partitions of Ω into simplices K and let V_{h} ⊂ H_{0}^{1}(Ω) be the space of continuous piecewise polynomials of degree < r (r ≥ 2) with respect to the mesh T . Our finite element problem is:

(1.7)  U_{ϵ} ∈ V_{h} : (∇U_{ϵ}, ∇χ) + (β_{ϵ}(U_{ϵ} − ψ), χ) = (f, χ), ∀χ ∈ V_{h}.

Detailed assumptions about the domain Ω, the mesh T , the penalty function ϵ, and the data f and ψ are stated in Section 2 below.

The obstacle problem is often considered as a prototype for a class of problems that involve free boundaries, modelling many phenomena such as phase transitions, jet flow, and gas expansion in a porous medium; see Friedman [8]. Thus this is a natural beginning problem to consider when extending the rapidly growing body of research on a posteriori error estimates for adaptive finite element algorithms. There are already several papers on this subject; see Ainsworth et al. [1], Hoppe and Kornhuber [12], Kornhuber [15], [16], [17], Johnson [13], Chen and Nochetto [3], Nochetto, Siebert, and Veeser [22], and Veeser [25]. We note also the related work by Nochetto, Paolini, and Verdi [21] on the Stefan problem.

The a posteriori estimate that we provide is based on the splitting

U_{ϵ} − u = (U_{ϵ} − u_{ϵ}) + (u_{ϵ} − u),

where the first term is the discretization error and the second is the penalty error.

We have chosen to conduct our analysis using the penalty formulation (1.4). This approach is motivated both by the partial differential equations analysis (see Friedman [8]) and by the desire to regularize the original non-smooth problem before computation. Finite element error analysis in the penalty formulation was first done by Scholz [23], who provided an a priori estimate in the energy norm using constant ϵ. Johnson [13] proved an a posteriori estimate in the energy norm and introduced the possibility of letting ϵ vary with x. Due to the monotonicity of the nonlinearity, the results obtained in the energy norm are essentially the same as for the corresponding linear problem. However, the standard Aubin-Nitsche duality argument does not go through here, because the linearized adjoint problem lacks the required regularity in L_{2}. Thus one does not obtain the usual rates of convergence that would be expected in the L_{2}-norm.

The aim of the present work is twofold. We first exploit the fact that the linearized adjoint problem does essentially have the necessary smoothing property in L_{1}, leading by duality to an error bound in the maximum norm. This bound is of the form

∥U_{ϵ} − u_{ϵ}∥_{L_{∞}(Ω)} ≤ C| log h_{max}| ∥h^{2}R_{∞}∥_{L_{∞}(Ω)},

where h is the piecewise constant mesh function defined by h|_{K} = diam(K) and R_{∞} is a computable function derived from the residual R = −ΔU_{ϵ} + β_{ϵ}(U_{ϵ} − ψ) − f ∈ H^{−1}(Ω) of the computed solution; see Theorem 3.4.

The second important feature of the present work is that we allow the penalty parameter ϵ to vary with x. Our result for the penalty error in the case of a smooth obstacle function, ψ ∈ W_{∞}^{2}(Ω), is

∥u − u_{ϵ}∥_{L_{∞}(Ω)} ≤ ∥ϵ(f + Δψ)∥_{L_{∞}(Ω̂)},

where Ω̂ = {x ∈ Ω : u(x) − ψ(x) = 0, u_{ϵ}(x) − ψ(x) ≤ 0} is the “contact set”; see Lemma 4.1. We also consider the case of a rough obstacle where ψ is merely Hölder continuous, but the resulting estimate is more involved; see Lemma 5.3.

The time-dependent obstacle problem is analyzed in a similar way by Boman [2]. For a posteriori error estimates in the maximum norm for linear elliptic problems we refer to Eriksson [6], Nochetto [20], and Dari, Durán, and Padra [4].

This paper is organized as follows. In Section 2 we state our assumptions on the continuous problem (1.1) and the finite element method and we introduce some notation. In Section 3 we provide the a posteriori estimate of the discretization error. In Section 4 we analyze the penalty error in the case of a smooth obstacle. Section 5 is devoted to the case of a rough obstacle. Finally in Section 6 we present the results of our numerical experiments with an adaptive algorithm in one dimension.

2. Notation and assumptions

We assume that Ω ⊂ R^{d}, d = 1, 2, 3, is a bounded polyhedral domain. Let ω ⊂ Ω with |ω| = meas(ω) and write

(u, v)_{ω} = ∫_{ω} uv dx, (u, v) = (u, v)_{Ω}.

We use the standard Lebesgue spaces L_{p}(ω) with the convention that L_{p} = L_{p}(Ω), and the corresponding Sobolev spaces W_{p}^{k}(ω), W_{p}^{k} = W_{p}^{k}(Ω), H^{k} = W_{2}^{k}. For k = 1, 2 we let W̊_{p}^{k} = {v ∈ W_{p}^{k} : v|_{∂Ω} = 0} with dual space W_{p}^{−k} and H_{0}^{1} = W̊_{2}^{1} with dual space H^{−1}. Moreover, for f ∈ (L_{p})^{d} and v ∈ W_{p}^{k} we write

(2.2)  ∥f∥_{L_{p}} = ( ∑_{i=1}^{d} ∥f_{i}∥_{L_{p}}^{p})^{1/p}, ∥D^{k}v∥_{L_{p}} = ( ∑_{|α|=k} ∥D^{α}v∥_{L_{p}}^{p})^{1/p}.

Then ∥D^{k}v∥_{L_{p}} is equivalent to the standard norm in W̊_{p}^{k} and we define the dual norm

(2.3)  ∥f∥_{W_{p}^{−k}} = sup_{v∈W̊_{p′}^{k}} \frac{|⟨f, v⟩|}{∥D^{k}v∥_{L_{p′}}} ,

where ⟨f, v⟩ is the duality pairing, and p, p′ are dual exponents, \frac{1}{p} + \frac{1}{p′} = 1.

We assume further that Ω has the property that there are constants C and p^{∗} ∈ (1, 2] such that

(2.4)  ∥D^{2}v∥_{L_{p}} ≤ \frac{C}{p − 1}∥Δv∥_{L_{p}}, ∀v ∈ W̊_{p}^{2}, 1 < p ≤ p^{∗}.

This holds with p^{∗} = 2 if Ω is bounded and convex. To see this, let v = T f be the solution of the Dirichlet problem −Δv = f in Ω, v = 0 on ∂Ω, and let D_{ij}^{2} denote a partial derivative of second order. It is well known [10] that the operator D_{ij}^{2}T is bounded on L_{2}, i.e., it is strong type (2,2); this is the case p = 2 of (2.4). Moreover, D_{ij}^{2}T is weak type (1,1); this is an unpublished result of Dahlberg, Verchota, and Wolff, but a proof can be found in [9] and a generalization in [7]. An application of the Marcinkiewicz interpolation theorem now yields (2.4). The inequality (2.4) also holds if ∂Ω is smooth by the Calderón-Zygmund theory of singular integrals.

It is plausible that (2.4) is true even for nonconvex polyhedral domains for some p^{∗} near 1, because the regularity implied by it then holds for p < 4/3; see [5], [10], [11], [18]; however, we do not know if the constant behaves like (p−1)^{−1} as p → 1 in this case. Related estimates are derived in [4], [20] for general polyhedral domains under the additional assumption that Δv has a small support. We formulate (2.4) as an assumption rather than a statement about convex domains, because we do not want to rule out this possibility; in any case we are not primarily interested in dealing with corner singularities in this work.

We further assume that the data satisfy f ∈ L_{∞}(Ω) and ψ ∈ H^{1}(Ω) ∩ C^{α}(Ω̄) is Hölder continuous with some exponent α ∈ (0, 1) and also ψ ≤ 0 on ∂Ω. Then (1.1) has a unique bounded weak solution u ∈ H_{0}^{1}(Ω)∩L_{∞}(Ω), which is also Hölder continuous, perhaps with a different exponent α. In our theorems below we make additional regularity assumptions about ψ.

Let F = {T } denote a regular family of partitions T = {K} of Ω into simplices K with diameters h_{K}, i.e., there is C_{0} > 0 such that, for all T ∈ F,

(2.5)  max_{K∈T} \frac{h_{K}}{ρ_{K}} ≤ C_{0},

where ρ_{K} denotes the radius of the largest ball contained in K̄. Let ϵ ∈ C(Ω̄) be a positive function and set

h_{min} = min_{K∈T} h_{K}, h_{max} = max_{K∈T} h_{K}, ϵ_{min} = min_{x∈Ω̄} ϵ(x), ϵ_{max} = max_{x∈Ω̄} ϵ(x).

We also define for each T ∈ F a piecewise constant mesh function h = h(x) by h|_{K} = h_{K}. We assume that there are constants γ ≥ 1 and C_{1} > 0 such that, for all T ∈ F,

(2.7)  h_{min} ≥ C_{1}h_{max}^{γ}, ϵ_{min} ≥ C_{1}h_{max}^{γ}.

Let r ≥ 2 be an integer and for each T ∈ F let V_{h} ⊂ H_{0}^{1}(Ω) be the space of continuous piecewise polynomials of degree < r with respect to the mesh T .

Under these assumptions it follows that problems (1.6) and (1.7) have unique solutions u_{ϵ} and U_{ϵ}, respectively.

3. The Discretization Error

In this section we compare the solutions u_{ϵ} and U_{ϵ}. We first define the errors

e = U_{ϵ} − u_{ϵ}, e_{β} = β_{ϵ}(U_{ϵ} − ψ) − β_{ϵ}(u_{ϵ} − ψ),

and the residual R = −ΔU_{ϵ} + β_{ϵ}(U_{ϵ} − ψ) − f ∈ H^{−1}, i.e.,

(3.2)  ⟨R, v⟩ = (∇U_{ϵ}, ∇v) + (β_{ϵ}(U_{ϵ} − ψ), v) − (f, v), ∀v ∈ H_{0}^{1},

where ⟨·, ·⟩ denotes the duality between H_{0}^{1} and H^{−1}. In view of (1.6) we thus have

(3.3)  e ∈ H_{0}^{1} : (∇e, ∇v) + (e_{β}, v) = ⟨R, v⟩, ∀v ∈ H_{0}^{1},

and (1.7) means that R is orthogonal to V_{h}, i.e.,

(3.4)  ⟨R, χ⟩ = 0, ∀χ ∈ V_{h}.

As a consequence of this orthogonality we have the following a posteriori bound of the residual with respect to the dual norms (2.3).

A POSTERIORI ERROR ANALYSIS FOR THE OBSTACLE PROBLEM 5

Lemma 3.1. For each p ∈ [1, ∞] and l = 1, 2 there is a constant C such that,

|⟨R, v⟩| ≤ C∥h^{l}R_{p}∥_{L_{p}}∥D^{l}v∥_{L_{p′}} , ∀v ∈ W̊_{p′}^{l},

where the functions h = h(x) and R_{p} = R_{p}(x) are defined piecewise by

(3.6) \begin{aligned}h|_{K} & = h_{K}, \\ R_{p}|_{K} & = | − ΔU_{ϵ} + β_{ϵ}(U_{ϵ} − ψ) − f| + h_{K}^{−1/p′}|K|^{−1/p}∥[∂_{n}U_{ϵ}]∥_{L_{p}(∂K\∂Ω)}. \end{aligned}

Here [∂_{n}U_{ϵ}] is the jump across ∂K in the exterior normal derivative ∂_{n}U_{ϵ} = n·∇U_{ϵ}. Proof. Write β_{ϵ} = β_{ϵ}(U_{ϵ} − ψ). Elementwise integration by parts in (3.2) gives

\begin{aligned}⟨R, v⟩ & = ∑_{K} ((∇U_{ϵ}, ∇v)_{K} + (β_{ϵ} − f, v)_{K}) \\ & = ∑_{K} ((−ΔU_{ϵ} + β_{ϵ} − f, v)_{K} + (∂_{n}U_{ϵ}, v)_{∂K}) \\ & = ∑_{K} ((−ΔU_{ϵ} + β_{ϵ} − f, v)_{K} − \frac{1}{2}([∂_{n}U_{ϵ}], v)_{∂K\∂Ω}), ∀v ∈ H_{0}^{1}. \end{aligned} ; (3.7)

The orthogonality (3.4) implies that ⟨R, v⟩ = ⟨R, v − χ⟩ for any χ ∈ V_{h}, so that

\begin{aligned}|⟨R, v⟩| & ≤ ∑_{K} (∥ − ΔU_{ϵ} + β_{ϵ} − f∥_{L_{p}(K)}∥v − χ∥_{L_{p′} (K)} \\ & + \frac{1}{2}∥[∂_{n}U_{ϵ}]∥_{L_{p}(∂K\∂Ω)}∥v − χ∥_{L_{p′} (∂K)}). \end{aligned} ; (3.8)

We choose χ = Πv, where Π : H_{0}^{1} → V_{h} is an interpolation operator such that

(3.9)  ∥D^{m}(v − Πv)∥_{L_{p′} (K)} ≤ Ch_{K}^{l−m}∥D^{l}v∥_{L_{p′} (S_{K})}, m = 0, 1, l = 1, 2,

where S_{K} is the union of all simplices adjacent to K; see [24]. Using also the (scaled) trace inequality

(3.10)  ∥w∥_{L_{p′} (∂K)} ≤ C(h_{K}^{−1/p′}∥w∥_{L_{p′} (K)} + h_{K}^{1−1/p′}∥Dw∥_{L_{p′} (K)}),

we obtain

∥v − Πv∥_{L_{p′} (∂K)} ≤ Ch_{K}^{l−1/p′}∥D^{l}v∥_{L_{p′} (S_{K})}, l = 1, 2.

Hence

\begin{aligned}|⟨R, v⟩| & ≤ C ∑_{K} h_{K}^{l}(∥ − ΔU_{ϵ} + β_{ϵ} − f∥_{L_{p}(K)} \\ & + h_{K}^{−1/p′}∥[∂_{n}U_{ϵ}]∥_{L_{p}(∂K\∂Ω)})∥D^{l}v∥_{L_{p′} (S_{K})} \\ & ≤ C∥h^{l}R_{p}∥_{L_{p}}∥D^{l}v∥_{L_{p′}} , \end{aligned} ; (3.12)

where the constant depends only on the constants in (3.9) and (3.10), which in turn depend only on the constant in (2.5). □

As a result of the inequality

(3.13)  (s^{−} − t^{−})(s − t) ≥ (s^{−} − t^{−})^{2}, ∀s, t ∈ R,

we have

(3.14)  (β_{ϵ}(v − ψ) − β_{ϵ}(w − ψ), v − w) ≥ ∥ϵ^{1/2}(β_{ϵ}(v − ψ) − β_{ϵ}(w − ψ))∥_{L_{2}}^{2}.

Taking v = e in (3.3) we immediately conclude from (3.14) and Lemma 3.1 with p = 2, l = 1, that

∥∇e∥_{L_{2}}^{2} + 2∥ϵ^{1/2}e_{β}∥_{L_{2}}^{2} ≤ C∥hR_{2}∥_{L_{2}}^{2},

which is Johnson’s result [13].

In order to estimate ∥e∥_{L_{p}} we argue by duality. We first note that (3.3) may be written as

(3.16)  e ∈ H_{0}^{1} : (∇e, ∇v) + (be, v) = ⟨R, v⟩, ∀v ∈ H_{0}^{1},

where

b(x) = \begin{cases}e_{β}(x)/e(x), & if e(x) ≠ 0, \\ 0, & if e(x) = 0. \end{cases}

Clearly, cf. (3.13), we have

(3.18)  0 ≤ b(x) ≤ ϵ(x)^{−1} ≤ ϵ_{min}^{−1}, ∀x ∈ Ω̄.

This suggests the introduction of the adjoint problem

(3.19)  G ∈ H_{0}^{1} : (∇w, ∇G) + (w, bG) = (w, g), ∀w ∈ H_{0}^{1},

with data g ∈ L_{p′}. Combining (3.16) and (3.19), and using Lemma 3.1, we get

(3.20)  |(e, g)| = |(∇e, ∇G) + (be, G)| = |⟨R, G⟩| ≤ C∥h^{2}R_{p}∥_{L_{p}}∥D^{2}G∥_{L_{p′}} .

Here we need to bound ∥D^{2}G∥_{L_{p′}} /∥g∥_{L_{p′}}. Trying first p = p′ = 2, we are only able to show, using (2.4), assuming temporarily that Ω is convex,

∥D^{2}G∥_{L_{2}} ≤ C∥ΔG∥_{L_{2}} ≤ C(∥g∥_{L_{2}} + ∥bG∥_{L_{2}}) ≤ C(1 + ϵ_{min}^{−1/2})∥g∥_{L_{2}};

see (3.32) below. This bound is too large to be useful, reflecting the fact that the nonlinearity is “too strong” to be controlled in the L_{2}-norm. On the other hand we have ∥bG∥_{L_{1}} ≤ ∥g∥_{L_{1}}, see (3.32), and we therefore use p = ∞, p′ = 1, in (3.20). However, ∥D^{2}G∥_{L_{1}}/∥g∥_{L_{1}} cannot be estimated by (2.4) for general data g and we therefore proceed as in [20] and choose g of a special form.

Let x_{0} ∈ Ω and let g = g_{x_{0}} be a regularized δ-function such that

(3.22)  ∫_{R^{d}} g dx = 1; supp g ⊂ B(x_{0}; ρ); 0 ≤ g(x) ≤ Cρ^{−d}.

Here B(x_{0}; ρ) denotes the closed ball with center x_{0} and small radius ρ to be chosen. Lemma 3.2. Let e satisfy (3.16) and x_{0} ∈ Ω be such that ∥e∥_{L_{∞}} = |e(x_{0})|. Let g = g_{x_{0}} be as in (3.22). There are constants h^{∗} > 0 and σ^{∗} > 0, depending only on Ω and the constants in Section 2, such that, if ρ ≤ h_{max}^{σ} with σ > σ^{∗} and h_{max} ≤ h^{∗}, then

(3.23)  ∥e∥_{L_{∞}} ≤ |(e, g)| + ∥h^{2}R_{∞}∥_{L_{∞}}.

Proof. In this proof we extend all functions by 0 outside Ω. By the mean value theorem there is x_{1} ∈ B(x_{0}; ρ) ∩ Ω̄ such that (e, g) = e(x_{1}). In order to estimate e(x_{0}) − e(x_{1}) we use a classical Hölder estimate of DeGiorgi and Nash: there are constants K > 0 and α ∈ (0, 1) such that if v ∈ H_{0}^{1} is a weak solution of −Δv = ∇ · F + G, where F ∈ (L_{p})^{d}, G ∈ L_{p/2} for some p > d, then for all x_{0} ∈ Ω, ρ > 0

sup_{x,y∈B(x_{0};ρ)∩Ω̄} |v(x) − v(y)| ≤ Kρ^{α}(∥F ∥_{L_{p}} + ∥G∥_{L_{p/2}}).

The constants K and α depend only on p and Ω; see [14, Theorem C.2].

We want to apply this result to equation (3.16) written in the form −Δe = R−be. By a well known characterization of the dual space W_{p}^{−1} we have

∥R∥_{W_{p}^{−1}} = inf_{F} ∥F ∥_{L_{p}},

where the norms are defined in (2.2) and (2.3), and the infimum is taken over all F ∈ (L_{p})^{d} such that R = ∇ · F, i.e., ⟨R, v⟩ = − ∫_{Ω} F · ∇v dx for all v ∈ W̊_{p′}^{1}. Hence, with p = 4 > d,

|e(x_{0}) − e(x_{1})| ≤ Kρ^{α}(∥R∥_{W_{4}^{−1}} + ∥be∥_{L_{2}}).

Inequality (3.18) and a standard energy argument, based on taking v = e in (3.16), imply

∥be∥_{L_{2}} ≤ ϵ_{min}^{−1/2}∥b^{1/2}e∥_{L_{2}} ≤ ϵ_{min}^{−1/2}∥R∥_{W_{2}^{−1}}.

Using Hölder’s inequality ∥Dv∥_{L_{1}} ≤ |Ω|^{1/p}∥Dv∥_{L_{p′}} in (2.3) and then Lemma 3.1 we get

∥R∥_{W_{p}^{−1}} ≤ |Ω|^{1/p}∥R∥_{W_{∞}^{−1}} ≤ C∥hR_{∞}∥_{L_{∞}} ≤ Ch_{min}^{−1}∥h^{2}R_{∞}∥_{L_{∞}},

and we conclude

|e(x_{0}) − e(x_{1})| ≤ Cρ^{α}ϵ_{min}^{−1/2}h_{min}^{−1}∥h^{2}R_{∞}∥_{L_{∞}}.

Recalling assumption (2.7) and taking ρ ≤ h_{max}^{σ} we finally get

|e(x_{0}) − e(x_{1})| ≤ Ch_{max}^{σα−3γ/2}∥h^{2}R_{∞}∥_{L_{∞}}

and the assertion of the lemma follows with σ^{∗} = \frac{3γ}{2α}. □

It remains to bound ∥D^{2}G∥_{L_{1}}/∥g∥_{L_{1}} = ∥D^{2}G∥_{L_{1}}. Lemma 3.3. Let G = G_{x_{0}} be the solution of (3.19) with data g = g_{x_{0}} as in (3.22). For any κ > 0 there is C such that, if ρ ≤ ϵ_{min}^{κ}, then, for all x_{0} ∈ Ω,

∥D^{2}G∥_{L_{1}} ≤ C| log ρ|.

Proof. Let first g be arbitrary. We shall prove that, for 1 ≤ p ≤ 2,

(3.32)  ∥bG∥_{L_{p}} ≤ (Cϵ_{min})^{−1/p′}∥g∥_{L_{p}}, ∀g ∈ L_{p}.

From (2.4) and (3.19) we then obtain, for 1 < p ≤ 2, (note that p−1 = p/p′ > 1/p′)

(3.33) ∥D^{2}G∥_{L_{p}} ≤ \frac{C}{p − 1}∥ΔG∥_{L_{p}} ≤ Cp′(∥g∥_{L_{p}} + ∥bG∥_{L_{p}}) ≤ Cp′ϵ_{min}^{−1/p′}∥g∥_{L_{p}}.

Now take g as in (3.22). A direct calculation gives

∥g∥_{L_{p}} ≤ Cρ^{−d(1−1/p)} = Cρ^{−d/p′},

and since ρ ≤ ϵ_{min}^{κ} we conclude from (3.33) with p′ = | log ρ| that

∥D^{2}G∥_{L_{1}} ≤ C∥D^{2}G∥_{L_{p}} ≤ Cp′ρ^{−c/p′} = C| log ρ|.

We now prove (3.32). Taking w = G/\sqrt{G^{2} + δ} with δ > 0 in (3.19), and passing to the limit as δ ↓ 0, yields

(3.36)  ∥bG∥_{L_{1}} ≤ ∥g∥_{L_{1}}, ∀g ∈ L_{1}.

Taking w = G in (3.19) shows ∥b^{1/2}G∥_{L_{2}} ≤ C∥g∥_{L_{2}}, so that in view of (3.18),

(3.37)  ∥bG∥_{L_{2}} ≤ ϵ_{min}^{−1/2}∥b^{1/2}G∥_{L_{2}} ≤ (Cϵ_{min})^{−1/2}∥g∥_{L_{2}}, ∀g ∈ L_{2}.

The Riesz-Thorin theorem applied to the linear operator g ↦ bG now yields

sup_{g} \frac{∥bG∥_{L_{p}}}{∥g∥_{L_{p}}} ≤ ( sup_{g} \frac{∥bG∥_{L_{1}}}{∥g∥_{L_{1}}} )^{1−2/p′}( sup_{g} \frac{∥bG∥_{L_{2}}}{∥g∥_{L_{2}}} )^{2/p′} ≤ (Cϵ_{min})^{−1/p′},

which is (3.32). □

We now state the main result of this section. Recall that R_{∞} is defined in (3.6). Theorem 3.4. There are C > 0 and h^{∗} > 0 such that, if h_{max} ≤ h^{∗}, then

∥U_{ϵ} − u_{ϵ}∥_{L_{∞}} ≤ C| log h_{max}| ∥h^{2}R_{∞}∥_{L_{∞}}.

Proof. Take ρ = h_{max}^{σ} with σ as in Lemma 3.2. In view of (2.7) we may find κ such that h_{max}^{σ} = ρ ≤ ϵ_{min}^{κ} and Lemma 3.3 applies. The error bound now follows from (3.20) in conjunction with these lemmas. □

Remark 3.1. The constant C of Lemma 3.3 enters as a “stability factor” in the error bound. The essence of Lemma 3.3 is that this stability factor is of moderate size; in particular, it is almost independent of ϵ.

4. The Penalization Error: Smooth Obstacle

We now compare the solutions u and u_{ϵ} of problems (1.1) and (1.4), respectively. In the following lemma we assume that the obstacle ψ is smooth. Lemma 4.1. Let u and u_{ϵ} be the solutions of (1.1) and (1.4), respectively. If ψ ∈ W_{∞}^{2}(Ω), then

(4.1)  ∥u − u_{ϵ}∥_{L_{∞}} ≤ ∥ϵ(f + Δψ)∥_{L_{∞}(Ω̂)},

where Ω̂ = Ω̂_{ψ} = {x ∈ Ω : u(x) − ψ(x) = 0, u_{ϵ}(x) − ψ(x) ≤ 0} is the “contact set”.

Proof. We define

Ω^{−} = {x ∈ Ω : u(x) − ψ(x) = 0},  Ω^{+} = {x ∈ Ω : u(x) − ψ(x) > 0},
Ω_{ϵ}^{−} = {x ∈ Ω : u_{ϵ}(x) − ψ(x) ≤ 0},  Ω_{ϵ}^{+} = {x ∈ Ω : u_{ϵ}(x) − ψ(x) > 0},

so that Ω̂ = Ω^{−} ∩ Ω_{ϵ}^{−}. Let v = u − u_{ϵ}. We shall show that there is a constant C such that

(4.3)  ∥v∥_{L_{q}} ≤ (Cq′)^{1/q}∥ϵ^{1/q′}(f + Δψ)∥_{L_{q}(Ω̂)},

for all even integers q ≥ 2 (q and q′ are conjugate exponents). Letting q → ∞ we then obtain (4.1).

In order to prove (4.3) we define

\begin{aligned}B(x) & := Δu(x) + f(x) ∈ β(u(x) − ψ(x)), \\ B_{ϵ}(x) & := Δu_{ϵ}(x) + f(x) = β_{ϵ(x)}(u_{ϵ}(x) − ψ(x)), \end{aligned} ; (4.4)

so that, for any even integer q ≥ 2,

(4.5)  (∇v, ∇v^{q−1}) = (B_{ϵ} − B, v^{q−1}),

where the left side is equal to (note that q − 1 = q/q′)

(4.6)  (∇v, ∇v^{q−1}) = (q − 1)∥v^{−1+q/2}∇v∥_{L_{2}}^{2} = \frac{4}{qq′} ∥∇(v^{q/2})∥_{L_{2}}^{2}.

We now turn to the right side of (4.5). We first show that

(4.7)  (B_{ϵ} − B, v^{q−1}) ≤ (B_{ϵ} − B, v^{q−1})Ω̂.

This follows from (4.4) and the monotonicity of the graphs β and β_{ϵ(x)}. More precisely, if x ∈ Ω^{+}, then B(x) = 0 = β_{ϵ(x)}(u(x) − ψ(x)), so that

(B_{ϵ} − B)v^{q−1} = −(β_{ϵ}(u_{ϵ} − ψ) − β_{ϵ}(u − ψ))(u_{ϵ} − u)v^{q−2} ≤ 0 in Ω^{+}.

Here we used the monotonicity (β_{ϵ}(s) − β_{ϵ}(t))(s − t) ≥ 0, and the assumption that q is even. Similarly, if x ∈ Ω_{ϵ}^{+}, then B_{ϵ}(x) = 0 ∈ {0} = β(u_{ϵ}(x) − ψ(x)), so that

(B_{ϵ} − B)v^{q−1} ∈ −(β(u_{ϵ} − ψ) − β(u − ψ))(u_{ϵ} − u)v^{q−2} ⊂ R^{−} in Ω_{ϵ}^{+},

since (β(s) − β(t))(s − t) ⊂ R^{+}. Therefore, (B_{ϵ} − B, v^{q−1})_{Ω^{+}∪Ω_{ϵ}^{+}} ≤ 0 and (4.7) follows.

It now remains to bound the right side of (4.7). In order to do so, we note that in Ω^{−} we have u = ψ, so that v = ψ − u_{ϵ}, and Δu = Δψ a. e., so that B = f + Δψ a. e. in Ω^{−}. (Note that this is where the smoothness of ψ is required; it implies that u is smooth, whence u = ψ yields Δu = Δψ.) Moreover, in Ω_{ϵ}^{−} we have B_{ϵ} = −ϵ^{−1}(ψ − u_{ϵ}) in view of (4.4).

Summing up: in Ω̂ = Ω^{−} ∩ Ω_{ϵ}^{−} we have B = f + Δψ, B_{ϵ} = −ϵ^{−1}v, so that

\begin{aligned}(B_{ϵ} − B, v^{q−1})Ω̂ & = (−ϵ^{−1}v − (f + Δψ), v^{q−1})Ω̂ \\ & = −∥ϵ^{−1/q}v∥_{L_{q}(Ω̂)}^{q} − (f + Δψ, v^{q−1})Ω̂. \end{aligned} ; (4.10)

Using Hölder’s and Young’s inequalities

|(f, g)| ≤ ∥f∥_{L_{q}}∥g∥_{L_{q′}} ≤ \frac{1}{q}∥f∥_{L_{q}}^{q} + \frac{1}{q′} ∥g∥_{L_{q′}}^{q′} ,

we get

\begin{aligned}|(f + Δψ, v^{q−1})Ω̂| & = |(ϵ^{1/q′}(f + Δψ), (ϵ^{−1/q}v)^{q/q′})Ω̂| \\ & ≤ \frac{1}{q}∥ϵ^{1/q′}(f + Δψ)∥_{L_{q}(Ω̂)}^{q} + (1 − \frac{1}{q})∥ϵ^{−1/q}v∥_{L_{q}(Ω̂)}^{q}. \end{aligned} ; (4.12)

Combining (4.5), (4.6), (4.7), (4.10), and (4.12), we conclude

(4.13)  \frac{4}{q′} ∥∇(v^{q/2})∥_{L_{2}}^{2} + ∥ϵ^{−1/q}v∥_{L_{q}(Ω̂)}^{q} ≤ ∥ϵ^{1/q′}(f + Δψ)∥_{L_{q}(Ω̂)}^{q}.

Using also Poincaré’s inequality

∥∇(v^{q/2})∥_{L_{2}}^{2} ≥ c∥v^{q/2}∥_{L_{2}}^{2} = c∥v∥_{L_{q}}^{q},

we finally obtain (4.3). □

Remark 4.1. Using q = 2 in (4.13) we get

(4.15)  2∥∇(u − u_{ϵ})∥_{L_{2}}^{2} + ∥ϵ^{−1/2}(ψ − u_{ϵ})∥_{L_{2}(Ω̂)}^{2} ≤ ∥ϵ^{1/2}(f + Δψ)∥_{L_{2}(Ω̂)}^{2},

which is (a slight improvement of) the corresponding result derived in [23] and used in [13]. More precisely, in [23] it is assumed that −Δψ ≤ 0, so that in (4.10) (−Δψ, v^{q−1})Ω̂ = (−Δψ, (ψ − u_{ϵ})^{q−1})Ω̂ ≤ 0, and this term may be dropped. The right side of (4.13) then becomes ∥ϵ^{1/q′}f∥_{L_{q}(Ω̂)}^{q}. A similar comment applies to [13], where it is assumed that ψ = 0 after the change of dependent variable u ← u − ψ. Note also that we localize the error bound to the contact set; this was not done in either [23] or [13]. Remark 4.2. Comparing with a similar analysis in [19], we note that we allow a variable penalty parameter ϵ = ϵ(x), and we do not truncate the nonlinearities β and β_{ϵ}.

Remark 4.3. Lemma 4.1 immediately generalizes to the “double obstacle problem;” cf. [23], where it is required that ψ(x) ≤ u(x) ≤ ϕ(x), and where the nonlinear terms are replaced by

β(u(x) − ψ(x)) − β(ϕ(x) − u(x)), β_{ϵ}(u_{ϵ}(x) − ψ(x)) − β_{ϵ}(ϕ(x) − u_{ϵ}(x)).

The result is

∥u − u_{ϵ}∥_{L_{∞}} ≤ ∥ϵ(f + Δψ)∥_{L_{∞}(Ω̂_{ψ})} + ∥ϵ(f + Δϕ)∥_{L_{∞}(Ω̂_{ϕ})},

with the obvious definitions of the contact sets Ω̂_{ψ}, Ω̂_{ϕ}.

To conclude this section we further point out that, in view of Theorem 3.4 and Lemma 4.1, we readily obtain an almost a posteriori error bound for the entire approximation procedure. Theorem 4.2. Let u be the solution of (1.1) and U_{ϵ} that of (1.7). There are C > 0 and h^{∗} > 0 such that if ψ ∈ W_{∞}^{2}(Ω) and h_{max} ≤ h^{∗}, then

(4.18)  ∥u − U_{ϵ}∥_{L_{∞}} ≤ C| log h_{max}| ∥h^{2}R_{∞}∥_{L_{∞}} + ∥ϵ(f + Δψ)∥_{L_{∞}(Ω̂)},

where Ω̂ = {x ∈ Ω : u(x) − ψ(x) = 0, u_{ϵ}(x) − ψ(x) ≤ 0} is the “contact set”. Remark 4.4. The error bound in Theorem 4.2 is not computable because the set Ω̂ is defined in terms of functions that are not known. To remedy this, we define a computable set

Ω_{τ} = {x ∈ Ω : U_{ϵ}(x) − ψ(x) ≤ τ},

where τ > 0 is a given tolerance. We claim that if

(4.19)  ∥u_{ϵ} − U_{ϵ}∥_{L_{∞}} ≤ τ,

then Ω̂ ⊂ Ω_{τ}. In fact, if x ∈ Ω̂, then

(U_{ϵ} − ψ)(x) = (U_{ϵ} − u_{ϵ})(x) + (u_{ϵ} − ψ)(x) ≤ τ + 0 = τ.

Remark 4.5. The basis for the corresponding algorithm of Section 6 is as follows. Produce an appropriate mesh so that (4.19) holds. Then evaluate the second term on the right side of (4.18), replacing Ω̂ by Ω_{τ}. Remark 4.6. Since β_{ϵ}(U_{ϵ} − ψ) = 0 on Ω \ Ω̂, the size of ϵ in the domain Ω \ Ω̂ does not affect the choice of meshsize h in Ω \ Ω̂. Therefore, singularities of ψ in Ω \ Ω̂ will not yield unnecessary refinement. Remark 4.7. In light of Theorem 4.2 we may expect ϵ(x) and h(x) to satisfy the local relation

(4.20)  ϵ(x) ≈ h^{2}(x), ∀x ∈ Ω.

A similar, but global, relation with ϵ constant was derived in [19] as a result of an a priori analysis.

5. The penalization error: rough obstacle

The proof of Lemma 4.1 does not generalize to the case when the obstacle function ψ has less than two derivatives in L_{∞}. This situation is analyzed in the following lemmas by means of auxiliary problems with a regularized obstacle. We stress that the auxiliary problems are only used in the analysis and thus do not alter the definitions of u_{ϵ} and U_{ϵ}. We assume the existence of μ ∈ (0, 2] such that

(5.1)  ψ ∈ C^{μ}(Ω̄),

which, in case μ is not an integer, means that the derivatives of order [μ] are Hölder continuous with exponent μ − [μ].

We start by examining the additional regularization of ψ. Let δ be a smooth function satisfying, for some constant C^{∗},

(5.2)  |D^{k}δ(x)| ≤ C^{∗}δ(x)^{1−k}, k = 0, 1, 2, ∀x ∈ Ω.

This function will be used as a space dependent regularization parameter, and can thus be viewed as a distance function. Its existence is guaranteed by the following lemma. Lemma 5.1. For every mesh T ∈ F there exists a positive function δ ∈ C^{∞}(Ω) satisfying

(5.3)  ch(x)^{1−k} ≤ |D^{k}δ(x)| ≤ C_{k}h(x)^{1−k}, ∀k ≥ 0, x ∈ Ω.

Proof. The case d = 2 is proved in [21, Lemma 5.1]. The proof, which is based on a C^{∞} partition of unity and the fact that any regular mesh is locally quasi-uniform, is readily adapted to any dimension d. □

We next introduce a measure of the regularity of ψ. We set, for ν ∈ (0, 2],

(5.4)  (D_{δ}^{ν}ψ)(x) = max_{y∈B(x;δ(x))} \frac{|ψ(y) − ψ(x) − V (x) · (y − x)|}{δ(x)^{ν}},

(5.5)  V (x) = \begin{cases}0, & if ν ∈ (0, 1], \\ ∇ψ(x), & if ν ∈ (1, 2]. \end{cases}

We are going to define, via convolution, a regularization ψ^{δ} of ψ with parameter δ = δ(x) as in (5.2). Let ζ ∈ C^{∞}(R^{d}) be positive, radially symmetric, with support in B(0; 1) and ∫ ζ dz = 1. We define

(5.6)  φ^{δ}(y, x) = \frac{1}{δ(x)^{d}} ζ(\frac{y − x}{δ(x)}),
(5.7)  ψ^{δ}(x) = ∫_{R^{d}} ψ(y)φ^{δ}(y, x) dy,

where ψ is extended outside Ω as a C^{μ}-function. Lemma 5.2. There exists C > 0, depending only on C^{∗}, such that, for all x ∈ Ω,

|ψ^{δ}(x) − ψ(x)| ≤ δ(x)^{ν}(D_{δ}^{ν}ψ)(x), (5.8)

|Δψ^{δ}(x)| ≤ Cδ(x)^{ν−2}(D_{δ}^{ν}ψ)(x). (5.9)

Proof. The following properties of φ^{δ} are obvious from its definition:

(5.10)  ∫_{R^{d}} φ^{δ}(y, x) dy = 1,
(5.11)  ∫_{R^{d}} \frac{∂}{∂x_{i}} φ^{δ}(y, x) dy = 0,
(5.12)  ∫_{R^{d}}(y − x)φ^{δ}(y, x) dy = 0.

Taking the derivative of the i^{th} component of (5.12) with respect to x_{j} and using (5.10) yields

(5.13)  ∫_{R^{d}}(y − x)_{i} \frac{∂}{∂x_{j}} φ^{δ}(y, x) dy = δ_{ij} (Kronecker’s delta).

Taking the derivative of this with respect to x_{k} and using (5.11) yields

(5.14)  ∫_{R^{d}}(y − x)_{i} \frac{∂^{2}}{∂x_{k}∂x_{j}} φ^{δ}(y, x) dy = 0.

From the definition (5.7) of ψ^{δ}, (5.10) and (5.12), we immediately obtain

\begin{aligned}|ψ^{δ}(x) − ψ(x)| & = | ∫_{R^{d}} (ψ(y) − ψ(x) − V (x) · (y − x))φ^{δ}(y, x) dy| \\ & ≤ δ(x)^{ν} ∫_{R^{d}} \frac{|ψ(y) − ψ(x) − V (x) · (y − x)|}{δ(x)^{ν}}φ^{δ}(y, x) dy \\ & ≤ δ(x)^{ν}D_{δ}^{ν}ψ(x), \end{aligned} ; (5.15)

and (5.8) follows. In a similar way, using also (5.14), we get

\begin{aligned}|Δψ^{δ}(x)| & = | ∫_{R^{d}} (ψ(y) − ψ(x) − V (x) · (y − x))Δ_{x}φ^{δ}(y, x) dy| \\ & ≤ δ(x)^{ν} ∫_{R^{d}} \frac{|ψ(y) − ψ(x) − V (x) · (y − x)|}{δ(x)^{ν}}|Δ_{x}φ^{δ}(y, x)| dy \\ & ≤ δ(x)^{ν}D_{δ}^{ν}ψ(x) ∫_{R^{d}} |Δ_{x}φ^{δ}(y, x)| dy, \end{aligned} ; (5.16)

and (5.9) follows once we have shown

(5.17)  |Δ_{x}φ^{δ}(y, x)| ≤ Cδ(x)^{−2}N(y, x), with ∫_{R^{d}} N(y, x) dy ≤ C.

A tedious but straightforward calculation starting from (5.6) gives

(5.18) |Δ_{x}φ^{δ}(y, x)| ≤ δ(x)^{−2}(1 + |∇δ(x)| + |∇δ(x)|^{2} + δ(x)|Δδ(x)|)N(y, x),

where N(y, x) is a linear combination of

(5.19)  \frac{1}{δ(x)^{d}} (\frac{|y − x|}{δ(x)})^{j}|D^{k}ζ(\frac{y − x}{δ(x)})|, j, k = 0, 1, 2,

and |D^{k}ζ(z)| denotes the norm of the k^{th} order gradient of ζ at z. Together with (5.2) this implies (5.17), and concludes the proof. □

Let u^{δ} and u_{ϵ}^{δ} denote the solutions of (1.1) and (1.4), respectively, with ψ replaced by ψ^{δ} and with Dirichlet boundary data max(0, ψ^{δ}); note that the convolution may cause ψ^{δ} > 0 on ∂Ω. The following lemma provides a localized error estimate of the effect of penalization. To state it, we need the enlarged contact set

(5.20)  Ω̃_{0} = ⋃_{i=1}^{4} Ξ_{i}

defined by

Ξ_{1} := {x ∈ Ω : u^{δ}(x) − ψ^{δ}(x) = 0}, Ξ_{2} := {x ∈ Ω : u(x) − ψ(x) = 0},
Ξ_{3} := {x ∈ Ω : u_{ϵ}^{δ}(x) − ψ^{δ}(x) ≤ 0}, Ξ_{4} := {x ∈ Ω : u_{ϵ}(x) − ψ(x) ≤ 0}.

Lemma 5.3. Let ν ∈ (0, 2] and let δ satisfy (5.2). Let u and u_{ϵ} be the solutions of (1.1) and (1.4), respectively. There exist C and δ^{∗} such that if (5.1) is valid and ∥δ∥_{L_{∞}} ≤ δ^{∗}, then

(5.21) ∥u − u_{ϵ}∥_{L_{∞}} ≤ ∥ϵf∥_{L_{∞}(Ω̃_{0})} + C∥ϵδ^{ν−2}D_{δ}^{ν}ψ∥_{L_{∞}(Ω̃_{0})} + 2∥δ^{ν}D_{δ}^{ν}ψ∥_{L_{∞}(Ω̃_{0})}.

Proof. We have

(5.22)  ∥u − u_{ϵ}∥_{L_{∞}} ≤ ∥u − u^{δ}∥_{L_{∞}} + ∥u^{δ} − u_{ϵ}^{δ}∥_{L_{∞}} + ∥u_{ϵ}^{δ} − u_{ϵ}∥_{L_{∞}}.

Lemma 4.1 gives

(5.23) ∥u^{δ} − u_{ϵ}^{δ}∥_{L_{∞}} ≤ ∥ϵ(f + Δψ^{δ})∥_{L_{∞}(Ω̂_{ψ^{δ}})} ≤ ∥ϵf∥_{L_{∞}(Ω̃_{0})} + ∥ϵΔψ^{δ}∥_{L_{∞}(Ω̃_{0})},

because Ω̂_{ψ^{δ}} = {x ∈ Ω : u^{δ}(x) − ψ^{δ}(x) = 0, u_{ϵ}^{δ}(x) − ψ^{δ}(x) ≤ 0} is contained in Ω̃_{0}; note that u^{δ} − u_{ϵ}^{δ} = 0 on ∂Ω, whence the proof of Lemma 4.1 is still valid. In view of (5.9) we readily get

(5.24)  ∥ϵΔψ^{δ}∥_{L_{∞}(Ω̃_{0})} ≤ C∥ϵδ^{ν−2}D_{δ}^{ν}ψ∥_{L_{∞}(Ω̃_{0})}.

We now assert that

(5.25)  ∥u − u^{δ}∥_{L_{∞}} ≤ ∥ψ − ψ^{δ}∥_{L_{∞}(Ξ_{1}∪Ξ_{2})},
(5.26)  ∥u_{ϵ} − u_{ϵ}^{δ}∥_{L_{∞}} ≤ ∥ψ − ψ^{δ}∥_{L_{∞}(Ξ_{3}∪Ξ_{4})},

which combined with (5.8) yield (5.21). We resort to the maximum principle to prove (5.25) and (5.26). In order to prove (5.25) we define

σ = ∥ψ − ψ^{δ}∥_{L_{∞}(Ξ_{1}∪Ξ_{2})},

set w = σ + u, and show that w ≥ u^{δ}. Suppose by contradiction that the open set G = {x ∈ Ω : w(x) < u^{δ}(x)} is nonempty. Since σ ≥ ψ − ψ^{δ} in Ξ_{1}, we deduce that G ∩ Ξ_{1} is empty, because

u^{δ} − ψ^{δ} > w − ψ^{δ} = σ + u − ψ^{δ} ≥ σ + ψ − ψ^{δ} ≥ 0 in G ∩ Ξ_{1}

which is inconsistent with the definition of Ξ_{1}. In addition, u^{δ} − ψ^{δ} > 0 outside Ξ_{1}, so that β(u^{δ} − ψ^{δ}) = {0} in G.

Also the set G ∩ ∂Ω is empty, because in G ∩ ∂Ω we have

0 > u + σ − u^{δ} = σ − max{0, ψ^{δ}},

which implies u^{δ} = ψ^{δ} and hence G ∩ ∂Ω ⊂ Ξ_{1}. But then, since σ ≥ ψ^{δ} − ψ, we have 0 > σ − ψ^{δ} ≥ −ψ ≥ 0 in G ∩ ∂Ω, which is a contradiction.

By the variational inequalities defining u^{δ} and u, we thus have

−Δu^{δ} = f ≤ −Δu = −Δw in G.

Hence −Δ(w − u^{δ}) ≥ 0 in G, w − u^{δ} = 0 on ∂G, so that w − u^{δ} ≥ 0 in G, which is a contradiction. This proves w ≥ u^{δ}, i.e., u^{δ} − u ≤ σ. The inequality u − u^{δ} ≤ σ is proved by interchanging the roles of ψ and ψ^{δ}, setting w = σ + u^{δ}, and making use of Ξ_{2}.

The proof of (5.26) is accomplished in a similar way: we define σ = ∥ψ − ψ^{δ}∥_{L_{∞}(Ξ_{3}∪Ξ_{4})}, w = σ + u_{ϵ} and show that w ≥ u_{ϵ}^{δ}. Suppose by contradiction that the open set G = {x ∈ Ω : w(x) < u_{ϵ}^{δ}(x)} is nonempty. Since σ ≥ ψ − ψ^{δ} in Ξ_{3} we have in G ∩ Ξ_{3}

u_{ϵ}^{δ} − ψ^{δ} ≥ w − ψ^{δ} = σ + u_{ϵ} − ψ^{δ} ≥ u_{ϵ} − ψ,

whence, exploiting the monotonicity of β_{ϵ},

β_{ϵ}(u_{ϵ}^{δ} − ψ^{δ}) ≥ β_{ϵ}(u_{ϵ} − ψ).

On the other hand, we have u_{ϵ}^{δ} − ψ^{δ} > 0 outside Ξ_{3} and consequently

0 = β_{ϵ}(u_{ϵ}^{δ} − ψ^{δ}) ≥ β_{ϵ}(u_{ϵ} − ψ).

Thus β_{ϵ}(u_{ϵ}^{δ} − ψ^{δ}) ≥ β_{ϵ}(u_{ϵ} − ψ) in G and therefore, from the equations defining u_{ϵ}^{δ} and u_{ϵ} we infer that

−Δu_{ϵ}^{δ} = f − β_{ϵ}(u_{ϵ}^{δ} − ψ^{δ}) ≤ f − β_{ϵ}(u_{ϵ} − ψ) = −Δu_{ϵ} = −Δw in G.

We check that G ∩ ∂Ω is empty in the same way as above. This implies that −Δ(w − u_{ϵ}^{δ}) ≥ 0 in G, and since w − u_{ϵ}^{δ} = 0 on ∂G, that w − u_{ϵ}^{δ} ≥ 0 in G, which is a contradiction. This proves w ≥ u_{ϵ}^{δ}, i.e., u_{ϵ}^{δ} − u_{ϵ} ≤ σ. The inequality ue − u_{ϵ}^{δ} ≤ σ is proved by interchanging the roles of ψ and ψ^{δ}, setting w = σ + u_{ϵ}^{δ}, and making use of Ξ_{4}. □

We see from (5.21) that a natural relation between δ and ϵ is

(5.35)  ϵ(x) ≈ δ(x)^{2}, ∀x ∈ Ω̃_{0}.

Coupling Theorem 3.4 with Lemma 5.3, we readily obtain an almost a posteriori error estimate for the whole approximation process. Theorem 5.4. Let ν ∈ (0, 2] and let δ satisfy (5.2). Let u and u_{ϵ} be the solutions of (1.1) and (1.4), respectively. There exist C, h^{∗} > 0, and δ^{∗} > 0 such that if (5.1) is valid, h_{max} ≤ h^{∗}, and ∥δ∥_{L_{∞}} ≤ δ^{∗}, then (with Ω̃_{0} defined in (5.20))

∥u − U_{ϵ}∥_{L_{∞}} ≤ C| log h_{max}| ∥h^{2}R_{∞}∥_{L_{∞}}

+ C(∥ϵf∥_{L_{∞}(Ω̃_{0})} + ∥ϵδ^{ν−2}D_{δ}^{ν}ψ∥_{L_{∞}(Ω̃_{0})} + ∥δ^{ν}D_{δ}^{ν}ψ∥_{L_{∞}(Ω̃_{0})}). (5.36)

Remark 5.1. As in Theorem 4.2 for smooth obstacles, the error bound (5.36) is not computable because the enlarged contact set Ω̃_{0} is defined in terms of functions that are not known. We now show that Ω̃_{0} may be replaced by a slightly larger, but computable, set Ω_{τ}. For a given tolerance τ > 0 let

Ω_{τ} = {x ∈ Ω : U_{ϵ}(x) − ψ(x) ≤ 3τ}.

We recall from (5.20) that Ω̃_{0} = ⋃_{i=1}^{4} Ξ_{i} and suppose x ∈ Ξ_{1} so that (u^{δ} −ψ^{δ})(x) = 0. From now on we assume that the full right side of (5.36) is less than or equal to τ. From Lemma 5.3 and its proof it follows that the quantities ∥u^{δ} − u∥_{L_{∞}} and ∥ψ − ψ^{δ}∥_{L_{∞}} are each less than or equal to τ. Then

U_{ϵ} − ψ = (U_{ϵ} − u) + (u − u^{δ}) + (u^{δ} − ψ^{δ}) + (ψ^{δ} − ψ) ≤ 3τ;

thus Ξ_{1} ⊂ Ω_{τ}. The remaining inclusions Ξ_{i} ⊂ Ω_{τ} are obtained in the same way. Thus Ω̃_{0} ⊂ Ω_{τ} if the right side of (5.36) is less than or equal to τ. Remark 5.2. An important consequence of the localized estimate (5.36) is that if the obstacle ψ were rough outside the enlarged contact set Ω̃_{0} it would not lead to unnecessary mesh refinement in Ω \ Ω̃_{0}. On the other hand, the size of ϵ in Ω \ Ω̃_{0} does not affect the choice of meshsize h in Ω \ Ω̃_{0} because β_{ϵ}(U_{ϵ} − ψ) = 0. This is consistent with both Remark 4.6 for smooth obstacles and Example 6.5 below.

6. Numerical experiments

In this section we discuss the results of preliminary computations in 1D using adaptive algorithms with h and ϵ refinement. We examine three distinct examples: one has a smooth obstacle and the other two have rough obstacles with jump discontinuities in the first derivative. In the first rough obstacle example the discontinuity occurs inside the contact set and we observe the algorithm setting ϵ to a very small number and refining the mesh heavily near this point. In the second rough obstacle example the discontinuity is outside the contact set. The algorithm sets ϵ to a small value near this point but does not refine the mesh heavily there, since the solution does not touch the rough obstacle and hence is smooth. The small value of ϵ outside the contact set does not degrade the performance of the iteration scheme, since β_{ϵ} is zero in that region. Our preliminary conclusions are that the a posteriori estimates provide useful upper bounds on the true errors and the possibility to vary ϵ provides more accuracy in the rough obstacle case.

6.1. Implementation. Let τ > 0 be a given tolerance, and let Ω_{τ} be the computable set of either Remark 4.4 or Remark 5.1. We use the a posteriori estimate

(6.1)  ∥u − U_{ϵ}∥_{L_{∞}(Ω)} ≤ C_{1} max_{K∈T} h_{K}^{2}∥R_{∞}∥_{L_{∞}(K)} + C_{2}∥ϵ(f + Δψ)∥_{L_{∞}(Ω_{τ})}

for both smooth obstacles (Theorem 4.2) and rough obstacles (Theorem 5.4), in the latter case with the second term on the right side replaced by

C_{2} (∥ϵf∥_{L_{∞}(Ω_{τ})} + ∥ϵD_{h}^{2}ψ∥_{L_{∞}(Ω_{τ})}) .

This keeps both algorithms almost identical and reflects the local nature of the a posteriori estimates (4.18) and (5.36). The term D_{h}^{2}ψ is defined by a difference at nodes on the mesh.

We use continuous piecewise linear and quadratic finite element functions in our approximations. In order to linearize we use Newton’s method and thus solve the following problem on each iteration: Find U_{ϵ}^{j+1} ∈ V_{h} such that

(∇U_{ϵ}^{j+1}, ∇χ) + (β_{ϵ}(U_{ϵ}^{j} − ψ) + (β′_{ϵ}(U_{ϵ}^{j} − ψ)(U_{ϵ}^{j+1} − U_{ϵ}^{j}), χ)_{h} = (f, χ)_{h},

for all χ ∈ V_{h}. The inner products with subscript h are Gauss quadrature approximations of the corresponding integrals.

6.2. Adaptive Algorithm. We now describe the full adaptive algorithm for the smooth obstacle example. \underline{Choose} the functions ϵ_{0}, U_{ϵ}^{−1}, and U_{ϵ}^{0}; define a mesh T_{0} (we started with a uniform mesh with N element domains); set index j = 0; and define tolerances τ (for adaptivity) and Tol (for the nonlinear iteration). \underline{Do Until}

C_{1} max_{K∈T_{j}}h_{K}^{2}∥R_{∞}∥_{L_{∞}(K)} + C_{2}∥ϵ(f + Δψ)∥_{L_{∞}(Ω_{τ})} ≤ τ and ∥U_{ϵ}^{j} − U_{ϵ}^{j−1}∥_{L_{∞}(Ω)} ≤ Tol:

• Compute U_{ϵ}^{j+1} on mesh T_{j} with function ϵ_{j} and using function U_{ϵ}^{j} to linearize.

• Define new mesh T_{j+1} by refining any elements, K, in the mesh T_{j} where h_{K}^{2}∥R_{∞}(U_{ϵ}^{j+1}, ϵ_{j})∥_{L_{∞}(K)} > τ/2.

• Define Ω_{j+1} = {x ∈ Ω : U_{ϵ}^{j+1} − ψ ≤ τ}.

• Choose ϵ_{j+1} = τ/(2 max{1, |f + Δψ|}) on the set Ω_{j+1}.

• Increment j. \underline{End Do Until Loop}.

Note that the mesh and ϵ refinement all happen at the same time. Moreover, there is only one Newton step for each pass through the loop.

Note that, in view of the discussions after Theorems 4.2 and 5.4, once the algorithm has converged, except for the possibility that the constants were chosen incorrectly, we expect that the estimate (6.1) with U_{ϵ} defined as the current iterate is an upper bound for the error.

N Estimated Error ∥u − U_{ϵ}∥_{L_{∞}(Ω)} Rate
12  3.0(-2)  2.6(-2)  —
24  3.2(-3)  2.6(-3)  3.3
71  8.7(-5)  7.6(-5)  3.3
154  5.3(-6)  5.0(-6)  3.5

Table 1. Example 6.3: Smooth obstacle.

6.3. Example: Smooth obstacle. We take Ω = (−1, 1), ψ(x) = 1 − 4x^{2}, f(x) = −8, and

u(x) = \{ \begin{matrix}4x^{2} − 16bx − (1 + 16b), & if x > b, \\ ψ(x), & if −b ≤ x ≤ b, \\ 4x^{2} + 16bx − (1 + 16b), & if x < −b, \end{matrix}

where b = 1 − \sqrt{6}/4 ≅ 0.3876275 . . .. Note that b is an irrational number. This ensures that the typical mesh does not have any nodes at the free boundary.

Since quadratics usually give L_{∞} a priori estimates with an N^{−3} convergence rate assuming the function being approximated has three derivatives this example is interesting since the true solution has only W_{∞}^{2} regularity with jumps in the second derivative at the free boundary.

Table 1 shows the results of a series of test runs on this problem, where we explore the rate of convergence obtained as a result of the mesh refinement and variable epsilon. We find that C_{1} = 1/50 is an upper bound for the first constant and that C_{2} = 1 is sufficient for the second. (Other estimates of the constant in front of the ∥u − u_{ϵ}∥_{L_{∞}} term have shown that the 6 may be replaced by a 1 with the addition of a very small term on the right side of the inequality).

Examining the table we note that the estimated error was above the true error as was generally the case and that the rate was very close to 3 as one would hope for quadratics.

6.4. Example: Rough obstacle. We first consider the specific case where the roughness is inside the contact set and thus the true solution, u, is also rough. Let Ω = (−1, 1), f(x) = −8,

ψ̃(x) = \{ \begin{matrix}3 − 4x, & if x > 0, \\ 3 + 4x, & if x ≤ 0, \end{matrix} ũ(x) = \{ \begin{matrix}4(x − 1)^{2}, & if x > 1/2, \\ ψ(x), & if −1/2 ≤ x ≤ 1/2, \\ 4(x + 1)^{2}, & if x < −1/2, \end{matrix}

and then define ψ(x) = ψ̃(x − δ) as well as u(x) = ũ(x − δ), where δ = \sqrt{2}/20 ≅ .0707106 …, so that the point of discontinuity is not likely to occur at a mesh point.

Below we list some results from the experiments with this example. We use continuous piecewise linear approximation functions in this case with C_{1} = C_{2} = 2 in the a posteriori error estimate.

Again we note that the mesh and ϵ refinement helps us recover a convergence rate of N^{−2} and that the estimator is an upper bound. The variable ϵ is particularly useful in this example. Near the jump discontinuity it becomes very small to compensate for the large D_{h}^{2}ψ term in the a posteriori estimate.

N Estimated Error ∥u − U_{ϵ}∥_{L_{∞}(Ω)} Rate
74  9.5(-2)  1.2(-2)  —
96  3.3(-2)  5.0(-3)  3.3
168  8.2(-3)  1.3(-3)  2.4
243  4.5(-3)  6.3(-4)  2.0

Table 2. Example 6.4: Rough obstacle.

Figure 1. Example 6.4: Comparison of the true solution and the refined approximations with variable and constant ϵ.

To demonstrate this we show the results of two computations, one with and one without allowing ϵ to vary in x. In the “Refine” case shown in Figure 6.4 the mesh has 78 elements, the L_{∞} error is .243, and ϵ = .003. In the “Refine/Epsilon” graph in Figure 6.4 the mesh has 74 elements, the L_{∞} error is .012, the mesh is shown in Figure 6.4, and the function ϵ is shown in Figure 6.4. The (apparent) payoff for the variable ϵ is seen in Figure 6.4, which shows how much closer and sharper the variable ϵ approximation is to the true solution. Note that the results displayed in Figure 6.4 are consistent with the solution function involved. On [−1, −0.5] and [0.5, 1] the size of the mesh elements is small, since u is quadratic, while it is piecewise linear on (−0.5, 0.5). The refinement is heavy at δ, where u′ has its only jump.

6.5. Example: Rough obstacle. We finally turn to the case, where the rough part of the obstacle is outside the contact set. Again, we take Ω = (−1, 1), and let f(x) = −9/2,

ψ̃(x) = \{ \begin{matrix}−3 + 3x, & if x > 0, \\ −3 − 3x, & if x ≤ 0, \end{matrix} ũ(x) = \{ \begin{matrix}ψ(x), & if x > 2/3, \\ 9x^{2}/4 − 2, & if −2/3 ≤ x ≤ 2/3, \\ ψ(x) & if x < −2/3, \end{matrix}

Figure 2. Example 6.4: The mesh function h versus x.

Figure 3. Example 6.4: The penalty function ϵ versus x.

and define ψ(x) = ψ̃(x − δ) as well as u(x) = ũ(x − δ), where δ = \sqrt{2}/20 ≅ .0707106 …, so that the point of discontinuity is not likely to occur at a mesh point. We use continuous piecewise linear approximation functions, and take the constants in the estimators to each be 2.

Since the jump discontinuity of ψ′ is not in the contact set, and so u is smooth, one would expect that the algorithm does not refine such a singularity of ψ and this is what we observe in the experiments. In Table 3 we display some results for

N Estimated Error ∥u − U_{ϵ}∥_{L_{∞}(Ω)} Rate
38  9.5(-2)  1.2(-2)  —
68  2.4(-2)  3.1(-3)  2.3
136  5.8(-3)  7.5(-4)  2.0
275  1.4(-3)  1.9(-4)  2.0

Table 3. Example 6.5: Rough obstacle with jump discontinuity in ψ′ outside the contact set.

this problem: we see that the estimator bounds the true error and the convergence rate is proportional to N^{−2}.

Acknowledgement. We thank A. Veeser for a suggestion leading to the present form of Lemma 3.2.

References

[1] M. Ainsworth, J. T. Oden, and C.-Y. Lee, Local a posteriori error estimators for variational inequalities, Numer. Methods Partial Differential Equations 9 (1993), 23–33.

[2] M. Boman, A posteriori error analysis in the maximum norm for a penalty finite element method for the time-dependent obstacle problem, (2000), Preprint, Department of Mathematics, Chalmers University of Technology and Göteborg University.

[3] Z. Chen and R. H. Nochetto, Residual type a posteriori error estimates for elliptic obstacle problems, Numer. Math. 84 (2000), no. 4, 527–548.

[4] E. Dari, R. G. Durán, and C. Padra, Maximum norm error estimators for three-dimensional elliptic problems, SIAM J. Numer. Anal. 37 (2000), no. 2, 683–700 (electronic).

[5] M. Dauge, Neumann and mixed problems on curvilinear polyhedra, Integral Equations Oper. Theory. 15 (1992), 227–261.

[6] K. Eriksson, An adaptive finite element method with efficient maximum norm error control for elliptic problems, Math. Models Methods Appl. Sci. 4 (1994), 313–329.

[7] M. Fassihi, L^{p} integrability of the second order derivatives of green potentials in convex domains, (1998), Preprint 1998–22, Department of Mathematics, Chalmers University of Technology and Göteborg University.

[8] A. Friedman, Variational Principles and Free-Boundary Problems, Wiley, New York, 1982.

[9] S. J. Fromm, Potential space estimates for Green potentials in convex domains, Proc. Amer. Math. Soc. 119 (1993), no. 1, 225–233.

[10] P. Grisvard, Elliptic problems in nonsmooth domains, Pitman, Boston, 1985.

[11] , Singular behavior of elliptic problems in non-Hilbertian Sobolev spaces, J. Math. Pures Appl. (9) 74 (1995), no. 1, 3–33.

[12] R. H. W. Hoppe and R. Kornhuber, Adaptive multilevel methods for obstacle problems, SIAM J. Numer. Anal. 31 (1994), 301–323.

[13] C. Johnson, Adaptive finite element methods for the obstacle problem, Math. Models Methods Appl. Sci. 2 (1992), 483–487.

[14] D. Kinderlehrer and G. Stampacchia, An Introduction to Variational Inequalities and Their Applications, Academic Press Inc., New York, 1980.

[15] R. Kornhuber, Monotone multigrid methods for elliptic variational inequalities. I, Numer. Math. 69 (1994), 167–184.

[16] , Monotone multigrid methods for elliptic variational inequalities. II, Numer. Math. 72 (1996), 481–499.

[17] R. Kornhuber, A posteriori error estimates for elliptic variational inequalities, Comput. Math. Appl. 31 (1996), 49–60.

[18] V. G. Maz’ya and B. A. Plamenevskii, L^{p} estimates of solutions of elliptic boundary value problems in a domain with edges, Trans. Moscow Math. Soc. 1 (1980), 49–97.

[19] R. H. Nochetto, Sharp L^{∞}-error estimates for semilinear elliptic problems with free boundaries, Numer. Math. 54 (1988), 243–255.

[20] , Pointwise a posteriori error estimates for elliptic problems on highly graded meshes, Math. Comp. 64 (1995), 1–22.

[21] R. H. Nochetto, M. Paolini, and C. Verdi, An adaptive finite element method for two-phase Stefan problems in two space dimensions. I. Stability and error estimates, Math. Comp. 57 (1991), 73–108, S1–S11.

[22] R. H. Nochetto, K. G. Siebert, and A. Veeser, Pointwise aposteriori error control for elliptic obstacle problems, in preparation.

[23] R. Scholz, Numerical solution of the obstacle problem by the penalty method, Computing 32 (1984), 297–306.

[24] L. R. Scott and S. Zhang, Finite element interpolation of nonsmooth functions satisfying boundary conditions, Math. Comp. 54 (1990), 483–493.

[25] A. Veeser, Efficient and reliable a posteriori error estimators for elliptic obstacle problems, (2000), preprint 02/2000, Mathematische Fakultät, Universität Freiburg.

Department of Mathematical Sciences, University of Cincinnati, Cincinnati, Ohio 45221

   E-mail address: donald.french@math.uc.edu

Department of Mathematics, Chalmers University of Technology and Göteborg University, SE–412 96 Göteborg, Sweden

   E-mail address: stig@math.chalmers.se

Department of Mathematics and Institute for Physical Science and Technology, University of Maryland, College Park, Maryland 20742

   E-mail address: rhn@math.umd.edu