## Abstract

We review the finite-element approximation of the classical obstacle problem in energy and max-norms and derive error estimates for both the solution and the free boundary. On the basis of recent regularity results, we present an optimal error analysis for the thin obstacle problem. Finally, we discuss the localization of the obstacle problem for the fractional Laplacian and prove quasi-optimal convergence rates.

## 1. Introduction

Obstacle and obstacle-type problems serve as a model for many phenomena and stand as the prototype for many theories: variational inequalities, constrained minimization and free boundary problems, to name a few. They are also the first truly nonlinear and non-smooth problem one encounters in the study of elliptic partial differential equations and their numerical approximation: linearization fails!

The purpose of this short note is to review the finite-element approximation of the solution to the classical, thin and fractional obstacle problems—the three basic prototypes. Some of the results in this note are classical but are scattered in the literature. Others, especially those about non-local operators, are completely new.

## 2. The classical obstacle problem

Let *Ω* be an open, bounded and connected domain of with boundary ∂*Ω*. Let *f*∈*H*^{−1}(*Ω*) and satisfying *ψ*≤0 on ∂*Ω*. We consider the following classical elliptic obstacle problem: find
2.1
where denotes the convex set of admissible displacements
2.2
(⋅,⋅) corresponds to the inner product in *L*^{2}(*Ω*) and 〈⋅,⋅〉 denotes the duality pairing between *H*^{−1}(*Ω*) and . It is known that under these conditions, (2.1) has a unique solution and there exists a non-negative Radon measure *μ* such that
and supp(*μ*)⊂*Λ*={*x*∈*Ω*:*u*(*x*)=*ψ*(*x*)}, where *Λ* is the so-called *coincidence* or *contact set*; see [1–3] for details. In addition, if *f*∈*L*^{2}(*Ω*) and *ψ*∈*H*^{2}(*Ω*), then *u*∈*H*^{2}(*Ω*) provided ∂*Ω* is *C*^{1,1} or *Ω* is a convex polyhedron. In this case, *μ* is absolutely continuous with respect to Lebesgue measure d*μ*=λ d*x* with 0≤λ∈*L*^{2}(*Ω*). With these assumptions, we can rewrite (2.1) as a *complementarity system* [1–3]:
2.3

We define the *non-coincidence set* *Ω*^{+}={*x*∈*Ω*:*u*(*x*)>*ψ*(*x*)} which, if *u* is continuous, is an open set, and the *free boundary* *Γ*=∂*Ω*^{+}∩*Ω*.

We remark that for the obstacle problem (2.1), the correspondence between the right-hand side *f* and *u* is nonlinear, non-differentiable and most notably not one-to-one as a change of *f* within the contact set *Λ* may yield no change in *u*.

### (a) Finite-element approximation

Let us now describe the discretization of problem (2.1). To avoid technical difficulties, let us assume that the boundary of *Ω* is polyhedral. Let be a mesh of *Ω* into cells *T* such that
where is a cell that is isoparametrically equivalent either to the unit cube [0,1]^{n} or the unit simplex in . The partition is assumed to be conforming or compatible, i.e. the intersection of any two cells *T* and *T*’ in is either empty or a common lower dimensional cell. We denote by the collection of all conforming meshes, which are refinements of a common mesh , and are *shape regular*, i.e. there exists a constant *σ*>1 such that
2.4
where *σ*_{T}=*h*_{T}/*ρ*_{T} is the shape coefficient of *T*. For simplicial elements, *h*_{T}=diam(*T*) and *ρ*_{T} is the diameter of the largest sphere inscribed in *T* [4]. For the definition of *h*_{T} and *ρ*_{T} in the case of cubes, we refer to [4]. The collection of meshes is *quasi-uniform* if for all and all , we have . In this case, the *mesh size* is .

Given a mesh , we define the following finite-element space:
2.5
where if *T* is a simplex , i.e. the set of polynomials of degree at most one, and if *T* is a cube, , that is, the set of polynomials of degree at most one in each variable.

Given a mesh and , we denote by N(*T*) the set of nodes of *T*. We then define and . Any discrete function is characterized by its nodal values on the set : , where is the canonical basis of , that is for .

The discrete admissible set , i.e. the discrete counterpart of , is
2.6
where is the Lagrange interpolation operator . The set is non-empty, convex, closed but in general is not a subset of . The classical literature would refer then to this setting as a non-conforming approximation. The *discrete problem* reads: find
2.7
It is well known that the solution exists and is unique [4,5].

### Remark 2.1 (positivity preserving operator)

In (2.6), the Lagrange interpolant can be replaced by the positivity preserving operator of Chen & Nochetto [6], which exhibits optimal approximation properties and preserves positivity, i.e. *w*≥0 implies . This allows us to weaken the regularity assumption on the obstacle to *ψ*∈*H*^{1}(*Ω*) (see remark 2.5).

### Remark 2.2 (discrete constraints)

As the discrete unilateral constraint in (2.6) is only enforced at the nodes of , its implementation is an easy task in practice. Alternative constraints such as *W*≥*ψ* or simplify the error analysis but complicate the implementation.

### (b) *A priori* error analysis in the energy norm

We now derive optimal *H*^{1}-error estimates under the assumption of full regularity. We follow the seminal paper by Brezzi *et al.* [5] (see also [8]).

### Theorem 2.4 (*H*^{1} error estimate)

*Let f*∈*L*^{2}(*Ω*) *and* *with ψ*≤0 *on* ∂*Ω, and let Ω be C*^{1,1} *or a convex polygon. Then*
2.8

### Proof.

Note, first of all, that the underlying assumptions imply that . Observe now that
2.9
and
Integrate by parts the second term in (2.9) and set in (2.7), to derive
We rewrite the right-hand side of the term above as follows:
2.10
In view of the complementarity relation in (2.3), we deduce that (λ,*u*−*ψ*)=0. The fact that , in conjunction with λ≥0, which also stems from (2.3), implies . We can now deal with the first term in (2.10) via a standard interpolation estimate
Combining the preceding estimates with , we obtain
This implies (2.8) and concludes the proof. ▪

### (c) Pointwise *a priori* error analysis

Let us now obtain pointwise error estimates which, as we will see in §2d, are essential to study the approximation of the free boundary *Γ*. This theory relies on two crucial ingredients, whose details are beyond the scope of this paper. The first one is a quasi-optimal *pointwise error estimate* for the so-called elliptic projection:
2.11

### Proposition 2.6 (pointwise error estimate)

*If* *then there exists a constant C independent of* *and u such that*
2.12

This classical result is due to Nitsche [9], Scott [10] and Schatz & Wahlbin [11]. For convenience, in what follows we denote the right-hand side of (2.12) by 2.13

The second ingredient is the *discrete maximum principle* (DMP): the mesh is *weakly acute* if
If *n*=2, this property holds over simplicial meshes provided that the sum of angles opposite to a side does not exceed *π* [12], §2.

### Proposition 2.7 (DMP)

*Let* *be weakly acute and* *be discrete subharmonic, i.e*.
*Then* .

With the help of the DMP, we prove pointwise error estimates between *u* and . This result is essentially due to Baiocchi [13] and Nitsche [9]. It extends to more general variational inequalities and lower than regularity. We begin with a simple but fundamental growth estimate.

### Lemma 2.8 (growth estimate)

*Assume that* *and let* . *If there is x*_{0}∈*T such that* (*u*−*ψ*)(*x*_{0})=0, *then there is a constant C proportional to* *such that*

### Proof.

The fact that *u*−*ψ*≥0 implies ∇(*u*−*ψ*)(*x*_{0})=0. As , we deduce . Applying the mean value theorem yields the asserted estimate. ▪

### Theorem 2.9 (pointwise error estimates)

*If* *is weakly acute and* *then there exists a constant C**>0 *depending on* *and* *such that, with* *as in* (*2.13*), *we have*
2.14

### Proof.

The proof relies on the construction of discrete super- and subsolutions to (2.7).

(1) *Supersolution.* Let with *C*_{1}>1 to be determined. We have
where denotes the Lagrange interpolation operator. The fact that implies the interpolation estimate which, together with *u*≥*ψ*, yields for *C*_{1} sufficiently large depending on .

For set in (2.1) to deduce
where we used that and that . To prove that , we argue by contradiction. Let
and note that because *C*_{1}>1. We thus have for all , whence for all , according to (2.7). Consequently,
Applying the DMP of proposition 2.7 to we infer that for all which contradicts the definition of . Hence,
2.15
where we used the pointwise estimate of proposition 2.6.

(2) *Subsolution*. Let for a suitable constant *C*_{2}>1. We again argue by contradiction to prove that . Let
and note that because *C*_{2}>1. We have that for each ,
This implies that *u*>*ψ* in supp *ϕ*_{z}, for otherwise we can apply the quadratic growth estimate of lemma 2.8 to obtain
which is a contradiction for *C*_{2} sufficiently large depending and .

Using (2.1), we get for all . Combining this with (2.7), we derive whence applying proposition 2.7 (DMP), we deduce for all . This contradicts the definition of and yields 2.16

(3) The desired estimate (2.14) is finally a consequence of (2.15) and (2.16). ▪

### (d) Error estimates for free boundaries

With the change of variables *v*=*u*−*ψ*, we can assume that the obstacle *ψ*=0. In this setting, let us address the approximation of the free boundary *Γ*=∂*Ω*^{+}∩*Ω*. We will follow the work of Nochetto [14], which is in turn inspired by the results of Brezzi & Caffarelli [15]. We present an elementary procedure to determine a discrete interface together with sharp convergence rates both in measure and in distance. The key ingredients are the so-called *non-degenerancy properties* (NDPs) of the obstacle problem (2.1) [1,16] and the definition of the discrete free boundary and non-coincidence set [14]: for a parameter to be properly selected, we set
2.17

Let us recall the NDP for the obstacle problem (2.1). Assume and that there is a negative constant *C*(*f*) for which
2.18
As , we have for and *s*<2+1/*p* [1,2], and pointwise regularity up to the fixed boundary ∂*Ω*, provided ∂*Ω*∈*C*^{2,α} [1,2].

We define *A*(*u*,*ϵ*)={*x*∈*Ω*:0<*u*(*x*)<*ϵ*^{2}}, which clearly depends on how the solution *u* behaves near the free boundary *Γ*, and . With this notation at hand, we present the following NDPs [1,16].

### Lemma 2.10 (local NDP)

*Let* . *If* (*2.18*) *holds, then the set A*(*u*,*ϵ*) *satisfies*
2.19
*where* . *In addition, we have*
2.20

A NDP thus prescribes a precise growth of *u* away from *Γ* and, as a result, *u* cannot be uniformly flat near *Γ*. A local NDP in measure was first established by Caffarelli in [16] under the qualitative property (2.18) and was exploited to derive regularity properties of the free boundary [1,16]. A global NDP in measure can be derived as soon as . On the other hand, a NDP in distance is more intricate in that it entails a pointwise behaviour of *u* near *Γ*. Let us now obtain an error estimate for the free boundaries [14].

### Theorem 2.11 (interface error estimates)

*Define* *with C** *given in* (*2.14*). *Then*
2.21
*If* ∂*Ω*∈*C*^{2,α}*, then one can set K*=*Ω. In addition, we have*
2.22

### Proof.

As , we analyse each set separately. We start by expressing as follows:
As , the NDP in measure (2.19) immediately implies for all . To estimate |*B*|, we first note that if *x*∈*B*, then . As the pointwise estimate (2.14) yields for almost every *x*, we deduce |*B*|=0. We thus arrive at
As , we argue as with *B* to get .

We now proceed to obtain (2.22). If , then , whence estimate (2.14) implies so that . The inclusion (2.22) now follows from the NDP in distance (2.20). ▪

## 3. The thin obstacle problem

To simplify the discussion, we assume that *Ω* is convex and replace the differential operators in (2.1) and (2.3) by −Δ+*I*. The admissible set is replaced by
If ∂*Ω*∈*C*^{1,1} and *g*∈*C*^{1,1/2}(∂*Ω*), then the corresponding variational inequality (2.1) has a unique solution ; the Hölder regularity is due to Athanasopoulos & Caffarelli [17]. With this regularity −Δ*u*+*u*=*f* a.e. in *Ω* and, if *z*=∂_{ν}*u* on ∂*Ω*, the following complementarity system is valid, which is also known as *Signorini complementarity conditions* (see [3], ch. 9):
3.1
As we do not assume *Ω* to be a polyhedral domain, we consider a family of triangulations of polyhedral domains that approximate *Ω* in such a way that , and . As *Ω* is convex, we have that and we extend continuously discrete functions to by a constant in the direction normal to . We define the discrete admissible set by
3.2
We present an optimal energy error estimate valid for any dimension and without assumptions on the free boundary. This improves upon the original result by Brezzi *et al.* [5]. Alternative results in this direction, for *n*=2, can be found in [18] and references therein.

### Theorem 3.1 (optimal energy error estimate)

*Let Ω be convex, C*_{Ω}>0 *denote the C*^{1,1}*-seminorm of* ∂*Ω and g*∈*C*^{1,1/2}(∂*Ω*). *Then there is C*>0, *that depends on* |*u*|_{H2(Ω)}, |*g*|_{C1,1/2(∂Ω)} *and C*_{Ω}*, for which*
3.3

### Proof.

The *C*^{1,1/2} regularity of *u* implies that 0≤*z*=∂_{ν}*u*∈*C*^{0,1/2}(∂*Ω*). We now proceed as in theorem 2.4 and write
To estimate I, we exploit that *u*∈*H*^{2}(*Ω*) and integrate by parts to obtain
where we have also used that solves the discrete problem. As −Δ*u*+*u*=*f* a.e. in *Ω*, the first term vanishes. For the second term, we write
The complementarity conditions (3.1) imply . As , we have that on ∂*Ω* which yields . For the remaining term, we only need to consider faces *S* on for which, on the corresponding subtended hypersurfaces , *u*−*g* is not identically zero nor strictly positive for otherwise or *z*=0, respectively. If *S* is one of these faces, then there exists such that the tangential gradient ∇_{Γ}(*u*−*g*)(*x*_{0})=0 because *u*−*g*≥0. As *u*−*g*∈*C*^{1,1/2}(∂*Ω*), employing an argument similar to lemma 2.8, we find *C*>0 that depends only on |*u*−*g*|_{C1,1/2(∂Ω)} for which
Conditions (3.1) imply there is with *z*(*x*_{1})=0. Now, as , we have that for all , with . Combine these estimates to obtain
We now split the term II into two contributions II and II) over the sets and , respectively. For the first one, we use interpolation theory to get We estimate the remaining term II as follows:
where we used that is stable in . As ∂*Ω* is *C*^{1,1}, we realize that
which, together with the estimate for I, yields the assertion (3.3). ▪

## 4. The fractional obstacle problem

We now consider the obstacle problem for fractional powers of the Laplace operator, which is a rather recent development. To be concrete, let *Ω* be an open, bounded and convex domain of (*n*≥1), with polyhedral boundary ∂*Ω*. The case of curved boundaries can be handled as in §3. For *s*∈(0,1), let be
4.1
which is a Hilbert space, and be its dual. Given and an obstacle satisfying *ψ*≤0 on ∂*Ω*, the fractional obstacle problem reads: find
4.2
where is the convex, closed and non-empty set of admissible displacements. Among the several definitions of (−Δ)^{s}, we adopt that in [19] which is based on the spectral theory of the Dirichlet Laplacian.

The problem (4.2) has a unique solution [1,2]. Moreover, if we assume that *Ω* and *f* are such that (−Δ)^{s}*u*∈*L*^{2}(*Ω*), we can rewrite (4.2) as a complementarity system:
4.3

### (a) Localization and truncation

The optimal Hölder regularity of the solution *u* to (4.3), namely *u*∈*C*^{1,s}(*Ω*), has been studied by Caffarelli *et al.* [20] using the extension proposed by Caffarelli & Silvestre in [21] to , which we now review. The problem
4.4
is equivalent to the Neumann-to-Dirichlet map associated with the following non-uniformly elliptic mixed boundary value problem:
4.5
where *α*=1−2*s*, , is the lateral boundary of , and *d*_{s} denotes a constant that depends only on *s* [19,21]. Define the weighted Sobolev space . The linear problem (4.5) can be formulated weakly as follows: find
For , we denote by tr_{Ω}*w* its trace onto *Ω*×{0}, and we recall that the trace operator *tr*_{Ω} satisfies [19], proposition 2.5
4.6

This extension problem is due to Caffarelli & Silvestre [21] for (see [19] and references therein for its modification to bounded domains). With its aid, we recast (4.2) as follows: find 4.7 where denotes the set of admissible displacements. If solves (4.7), then solves (4.2). The obstacle constraint is thus applied on part of .

We now exploit the exponential decay of the solution to (4.7) to truncate the cylinder . In what follows, we denote by λ_{1} the first eigenvalue of the Dirichlet Laplacian −Δ.

### Lemma 4.1 (exponential decay)

*Let* *denote the solution to* (*4.7*). *Then, we have*

### Proof.

Note that solves
The representation formula provided in [19], (2.24), together with the decay estimates of [19], proposition 3.1, applied to this problem yield
Stability of problem (4.7) in terms of *ψ* and *f* allows us to deduce the desired estimate. ▪

In view of such an exponential decay, we truncate the cylinder to a height *y*=Y, i.e. set (see [19], theorem 3.5 for the linear case). We then consider
4.8
where As in [19], we only incur an exponentially small error by considering the truncated version (4.8) of problem (4.7). In addition, satisfies a complementarity system such as (4.3):
4.9

Using the extended problem (4.7) to , Caffarelli *et al*. [20], theorem 6.7 showed that if *x*_{0}′∈*Ω* and , then
4.10
for a positive constant *C*. This implies that and because . As their techniques are *local*, they apply to our extension to the cylinder and truncation to as well, and give the optimal Hölder regularity for :
4.11

These results account only for the regularity of in *Ω*. However, we also need to know the regularity of over the cylinder , which is established by Allen *et al.* [22], theorem 6.4:
4.12

### (b) Discretization

In [19], we have used the local approach of Caffarelli & Silvestre [21] to approximate problem (4.4). After having truncated to , the next issue is to compensate for the rather singular behaviour of by *anisotropic* meshes of [0, Y] with mesh points:
4.13
This entails the development of a polynomial interpolation theory in weighted Sobolev spaces over anisotropic meshes with Muckenhoupt weights [19,23].

Let be a conforming and shape regular mesh of *Ω* as in §2a. The collection of these triangulations is denoted by . We construct the mesh of the cylinder as the tensor product triangulation of and , and denote the set of all such triangulations by . We assume that there is a constant *σ*_{Y} such that if *T*_{1}=*K*_{1}×*I*_{1} and have non-empty intersection, then , where *h*_{I}=|*I*|. This weak regularity condition on the mesh allows for anisotropy in the extended variable [19,23,24]. For , we define the finite-element space
4.14
where is the Dirichlet boundary and is defined as in §2a. The Galerkin approximation of the truncated problem is then given by
4.15

Note that and imply . Finally, if is shape regular and quasi-uniform, we have . All these considerations allow us to obtain the following result (see [19], theorem 5.4 and [23], corollary 7.11).

### Theorem 4.2 (a priori error estimate)

*Let* *be a tensor product grid, which is quasi-uniform in Ω and graded in the extended variable so that* (*4.13*) *holds. If* *is defined by* (*4.14*), *solves* (*4.15*) *and* *solves* (*4.5*), *then we have*
*where* *. Alternatively, if u solves* (*4.4*), *then*

The discretization of (4.2) is then carried out by a Galerkin approximation to (4.8), namely
4.16
for all . The discrete admissible set is defined by
4.17
where is the *α*-harmonic extension of *ψ* and is a positivity preserving variant of the quasi-interpolation operator on the mesh defined in [19,23], which is constructed by local averaging, is able to separate the variables *x*′∈*Ω* and *y*∈(0,Y), is stable in and and exhibits *quasi-local* optimal approximation properties for any .

### (c) Energy error analysis

We build on the preceding approach to obtain an almost optimal error estimate for the solution of (4.2). We start by quantifying the error due to truncation of .

### Proposition 4.3 (exponential error estimate)

*Let* *and* *solve problems* (*4.7*) *and* (*4.8*), *respectively. If Y*≥1, *then we have*

### Proof.

By definition, we have
We examine each term separately. For the first term, as , we set in (4.7) to get
4.18
For the second term, we cannot set in (4.8) since this is not a valid test function because . Instead, let be a smooth cut-off function such that *ρ*≡1 on [0,Y/2], *ρ*=(2/Y)(Y−*y*) on [Y/2,Y] and *ρ*=0 on [19], (3.7). Then, we write
4.19
Set in (4.8) to obtain
4.20
In light of (4.18)–(4.20), we see that the result will follow if we bound the last term in (4.19). For that we use lemma 4.1 and the arguments of Nochetto *et al*. [19], lemma 3.3 to obtain
which readily yields the asserted estimate. ▪

We now present an almost optimal *a priori* error estimate which relies on [5] and theorem 2.4.

### Theorem 4.4 (almost optimal error estimate)

*Let* *be the tensor product grid described in theorem 4.2. If* *and* *solve* (*4.7*) *and* (*4.16*), *respectively, we have*
4.21
*Alternatively, if u solves* (*4.2*), *then*
4.22
*where, in both inequalities, C depends on the Hölder moduli of smoothness of* *given by* (*4.11*) *and* (*4.12*), *and* .

### Proof.

We first compare the solution of the truncated problem and the discrete solution . To do so, we proceed as in theorem 2.4. We just need to examine
because in and is defined in (4.9). We can rewrite the preceding term as
4.23
Using (4.9), together with tr_{Ω}*Ψ*=*ψ* and , we infer the estimates and . We now express the last term on the right-hand side of (4.23) as follows:
We next examine separately the cells according to the value of in *S*_{K}, a discrete neighbourghood of *K*. The issue at stake is that the interpolation operator hinges on local averages over a discrete neighbourhood *S*_{T} of *T*=*K*×*I*. This leads to the following three cases.

Case 1: in *S*_{K}. In this situation, in *S*_{K} and thus vanishes.

Case 2: in *S*_{K}. The fact that is max-norm stable locally yields
4.24
The case follows immediately from (4.12). For , we use that solves an *α*-harmonic extension problem and then its conormal derivative is well defined. We realize that vanishes in *Ω* because , which together with (4.12) allows us to derive (4.24). This, combined with (4.11), implies the following bound:
As and *γ*>3/2*s*, we thus conclude .

Case 3: is not identically zero nor strictly positive in *S*_{K}. In view of (4.10), we have
4.25
On the other hand, using (4.11) we deduce and
4.26
Consequently, . The fact that there is a point *x*_{0}′∈*S*_{K} where , and an argument similar to the one that led us to (4.24) on the basis of (4.12), in conjunction with (4.25) and (4.26), yield
and . Hence, .

Collecting the estimates for the three cases, we thus conclude where Y accounts for the interpolation estimate of based on the mesh grading (4.13) [19]. The estimate (4.21) follows from proposition 4.3 and a suitable choice of the parameter Y in terms of [19], remark 5.5. Finally, (4.6) and (4.21) lead to (4.22). ▪

## 5. Conclusion and open problems

Several topics of interest were only mentioned in passing or not at all. Among them, mixed methods, monotone multigrid methods and *a posteriori* error estimation come to mind. Other discretization techniques such as non-conforming finite elements, virtual elements or mimetic finite differences have not been described. Let us also list some open problems of interest. The convergence properties of multigrid methods in higher dimensions are still not well understood. The lack of duality techniques makes obtaining error estimates in *L*^{2}(*Ω*) for the classical obstacle problem rather difficult. Pointwise error estimates for the thin and fractional obstacle problems are rather technical but appear to be an interesting open issue. Other obstacle-type problems, such as time-dependent problems and high-order equations, might require techniques other than those presented here and were not discussed.

## Competing interests

We declare we have no competing interests.

## Funding

R.H.N. and E.O. have been supported in part by NSF grant nos. DMS-1109325 and DMS-1411808. A.J.S. has been supported in part by NSF grant no. DMS-1418784.

## Footnotes

One contribution of 15 to a theme issue ‘Free boundary problems and related topics’.

- Accepted December 2, 2014.

- © 2015 The Author(s) Published by the Royal Society. All rights reserved.