Error identities for variational problems with obstacles

TL;DR

This paper develops error identities for nonlinear free boundary problems solved via variational methods, linking a nonlinear measure of solution accuracy to the duality gap and providing computable bounds.

Contribution

It introduces specific error identities that relate a nonlinear measure of solution accuracy to the duality gap, enabling better error estimation in variational free boundary problems.

Findings

Error identities relate nonlinear measures to the duality gap.

Two-sided bounds provide efficient error estimates.

Different minimization sequences affect error component dominance.

Abstract

The paper is concerned with a class of nonlinear free boundary problems, which are usually solved by variational methods based on primal (or primal-dual) variational settings. We deduce and investigate special relations (error identities). They show that a certain nonlinear measure of the distance to the exact solution (specific for each problem) is equivalent to the respective duality gap, which minimization is a keystone of all variational numerical methods. Therefore, the identity defines the measure that contains maximal quantitative information on the quality of a numerical solution available through these methods. The measure has quadratic terms generated by the linear part of the differential operator and nonlinear terms associated with free boundaries. We obtain fully computable two sided bounds of this measure and show that they provide efficient estimates of the distance between the minimizer and any function from the corresponding energy space. Several examples show that for different minimization sequence the balance between different components of the overall error measure may be different and domination of nonlinear terms may indicate that coincidence sets are approximated incorrectly.

Tables

Table 1. Table 1. Terms in the primal error identity computed for various approximation v ϵ 1 subscript 𝑣 subscript italic-ϵ 1 v_{\epsilon_{1}} in case of the classical obstacle in 1D. The rightmost column shows the contribution of μ ϕ ​ ψ ​ ( v ϵ 1 ) subscript 𝜇 italic-ϕ 𝜓 subscript 𝑣 subscript italic-ϵ 1 \mu_{\phi\psi}(v_{\epsilon_{1}}) to the majorant value 𝕄 ​ ( u , v ϵ 1 ) 𝕄 𝑢 subscript 𝑣 subscript italic-ϵ 1 {\mathds{M}}(u,v_{\epsilon_{1}}) .

${ϵ}_1$	$\frac{1}{2}{\Vert \nabla \left(u-v_{{ϵ}_1}\right)\Vert }^2$	${\mu }_{\varphi \psi }\left(v_{{ϵ}_1}\right)$	$𝕄(u,v_{{ϵ}_1})$	$J\left(v_{{ϵ}_1}\right)-J\left(u\right)$	${\mu }_{\varphi \psi }\left(v_{{ϵ}_1}\right)$ [%]
0.1000	1.54e-01	4.09e-02	1.95e-01	1.95e-01	20.92
0.0500	4.82e-02	6.37e-03	5.45e-02	5.45e-02	11.68
0.0250	1.36e-02	8.98e-04	1.45e-02	1.45e-02	6.20
0.0125	3.62e-03	1.20e-04	3.73e-03	3.73e-03	3.20
0.0063	9.33e-04	1.54e-05	9.49e-04	9.49e-04	1.63

Table 2. Table 4. Terms in the primal error identity computed for various approximation v ϵ 1 subscript 𝑣 subscript italic-ϵ 1 v_{\epsilon_{1}} in case of the double obstacle in 1D. The rightmost column shows the contribution of μ ϕ ​ ψ ​ ( v ϵ 1 ) subscript 𝜇 italic-ϕ 𝜓 subscript 𝑣 subscript italic-ϵ 1 \mu_{\phi\psi}(v_{\epsilon_{1}}) to the majorant value 𝕄 ​ ( u , v ϵ 1 ) 𝕄 𝑢 subscript 𝑣 subscript italic-ϵ 1 {\mathds{M}}(u,v_{\epsilon_{1}}) .

${ϵ}_1$	$\frac{1}{2}{\Vert \nabla \left(u-v_{{ϵ}_1}\right)\Vert }^2$	${\mu }_{\omega }\left(v_{{ϵ}_1}\right)$	$𝕄(u,v_{{ϵ}_1})$	$J\left(v_{{ϵ}_1}\right)-J\left(u\right)$	${\mu }_{\omega }\left(v_{{ϵ}_1}\right)$ [%]
0.2000	2.18e-01	8.71e-02	3.05e-01	3.05e-01	28.57
0.1000	7.41e-02	1.48e-02	8.89e-02	8.89e-02	16.67
0.0500	2.20e-02	2.20e-03	2.42e-02	2.42e-02	9.09
0.0250	6.05e-03	3.02e-04	6.35e-03	6.35e-03	4.76
0.0125	1.59e-03	3.97e-05	1.63e-03	1.63e-03	2.44

Equations

$$\underset{\Omega }{\int }\nabla u\cdot \nabla wdx\ge \underset{\Omega }{\int }fwdx\forall w\in K:=\{w\in H^1(\Omega )\mid w=u_{D}on\partial \Omega ,w\ge \psi \}$$

$$\underset{v\in V}{\inf }J(v),J(v)=G(\Lambda v)+F(v),$$

$$I^{\ast }(y^{\ast }):=-G^{\ast }(y^{\ast })-F^{\ast }(-{\Lambda }^{\ast }{y}^{\ast })$$

$$D_F(v,v^{\ast }):=$$

$$\text{mathds}M(\{u,p^{\ast }\},\{v,y^{\ast }\}):=D_F(u,-{\Lambda }^{\ast }{y}^{\ast })+D_F(v,-{\Lambda }^{\ast }{p}^{\ast })+D_G(\Lambda u,y^{\ast })+D_G(\Lambda v,p^{\ast })\ge 0.$$

$$\Lambda v\in \partial G^{\ast }(p^{\ast }),y^{\ast }\in \partial G(\Lambda u),-{\Lambda }^{\ast }{y}^{\ast }\in \partial F(u),v\in \partial F^{\ast }(-{\Lambda }^{\ast }{p}^{\ast }).$$

$$\text{mathds}M(\{u,p^{\ast }\},\{v,y^{\ast }\})=J(v)-I^{\ast }(y^{\ast }).$$

$$\text{mathds}M(u,v)=$$

$$\text{mathds}M(p^{\ast },y^{\ast })=$$

$$K:=\{v\in V_0:=H_0^1(\Omega )\mid \varphi (x)\le v(x)\le \psi (x)a.e.\mathrm{in}\Omega \}.$$

$$\varphi (x)\le 0\mathrm{on}\partial \Omega ,\psi (x)\ge 0\mathrm{on}\partial \Omega ,\varphi (x)\le \psi (x),\forall x\in \Omega .$$

$$a(u,w-u)\ge (f,w-u)\forall w\in K$$

$$a(u,w):=\underset{\Omega }{\int }A\nabla u\cdot \nabla wdx.$$

$$A(x)\xi \cdot \xi \ge c_1|\xi {|}^2{c}_1>0,\forall \xi \in {\mathbb{R}}^d$$

$${\Omega }_{-}^u:=$$

$${\Omega }_{-}^v:=$$

$$p^{\ast }=A\nabla u$$

$$\mathrm{div}{p}^{\ast }+f\le 0$$

$$\Lambda v=\nabla v,$$

$$\displaystyle\chi_{K}(v):=\left\{\begin{array}[]{cc}0&{\rm if}\;\phi\leq v\leq\psi,\\ +\infty&{\rm else}.\end{array}\right.$$

$$G^{\ast }(y^{\ast })$$

$$D_G(\Lambda v,y^{\ast })$$

$$D_G(\Lambda v,p^{\ast })=\frac{1}{2}\underset{\Omega }{\int }A\nabla (u-v)\cdot \nabla (u-v)dx=:\frac{1}{2}||\nabla (u-v)|{|}_A^2,$$

$$D_G(\Lambda u,y^{\ast })=\frac{1}{2}\underset{\Omega }{\int }{A}^{-1}(p^{\ast }-y^{\ast })\cdot (p^{\ast }-y^{\ast })dx=:\frac{1}{2}||p^{\ast }-y^{\ast }|{|}_{A^{-1}}^2.$$

$$F^{\ast }(v^{\ast })=\underset{v\in K}{\sup }\underset{\Omega }{\int }v(v^{\ast }+f)dx=\underset{v\in K}{\sup }\underset{\Omega }{\int }(-v(v^{\ast }+f{)}_{-}+v(v^{\ast }+f{)}_{+})dx=\underset{\Omega }{\int }(-\varphi (v^{\ast }+f{)}_{-}+\psi (v^{\ast }+f{)}_{+})dx.$$

$$\begin{gathered} D_F(v,-{\Lambda }^{\ast }{y}^{\ast })=F(v)+F^{\ast }(-{\Lambda }^{\ast }{y}^{\ast })+⟨{\Lambda }^{\ast }{y}^{\ast },v⟩= \\ =\underset{\Omega }{\int }(-fv-\varphi (\mathrm{div}{y}^{\ast }+f{)}_{-}+\psi (\mathrm{div}{y}^{\ast }+f{)}_{+}-\mathrm{div}{y}^{\ast }v)dx= \\ =\underset{\Omega }{\int }((v-\varphi )(\mathrm{div}{y}^{\ast }+f{)}_{-}+(\psi -v)(\mathrm{div}{y}^{\ast }+f{)}_{+})dx. \end{gathered}$$

$$\begin{gathered} D_F(v,-{\Lambda }^{\ast }{p}^{\ast })=-\underset{{\Omega }_{-}^u}{\int }(v-\varphi )(\mathrm{div}{p}^{\ast }+f)dx+\underset{{\Omega }_{+}^u}{\int }(\psi -v)(\mathrm{div}{p}^{\ast }+f)dx= \\ =-\underset{{\Omega }_{-}^u}{\int }(v-\varphi )(\mathrm{div}A\nabla \varphi +f)dx+\underset{{\Omega }_{+}^u}{\int }(\psi -v)(\mathrm{div}A\nabla \psi +f)dx. \end{gathered}$$

$${\mu }_{\varphi \psi }(v):=\underset{{\Omega }_{-}^u}{\int }{\mathsf{W}}_{\varphi }(v-\varphi )dx+\underset{{\Omega }_{+}^u}{\int }{\mathsf{W}}_{\psi }(\psi -v)dx,$$

$${\mu }_{\varphi \psi }(v)=f\underset{{\Omega }_{+}^u}{\int }(\psi -v)dx=f\Vert \psi -v{\Vert }_{L^1({\Omega }_{+}^u)}=f\Vert \psi -v{\Vert }_{L^1({\Omega }_{+}^u\setminus {\Omega }_{+}^v)}.$$

$${\Omega }_{+}^u=({\Omega }_{+}^u\setminus {\Omega }_{+}^v)\cup ({\Omega }_{+}^u\cap {\Omega }_{+}^v)$$

$${\mu }_{\varphi \psi }(v)=-f\underset{{\Omega }_{-}^u}{\int }(v-\varphi )dx=-f\Vert v-\varphi {\Vert }_{L^1({\Omega }_{-}^u)}=-f\Vert v-\varphi {\Vert }_{L^1({\Omega }_{-}^u\cap {\Omega }_{-}^v)}.$$

Full text

1. Introduction

Variational inequalities form an important class of nonlinear models that describe free boundary phenomena arising in various applied problems (see, e.g., G. Duvaut and J. L. Lions [8] and other publications cited therein). Usually free boundaries separate regions where solutions possess quite different physical properties. Therefore, any reliable information on the shape and location of such a boundary is very important. Qualitative properties of free boundaries are studied by purely analytical (a priori) methods unlike quantitative information, which in the vast majority of cases can be obtained only by computational methods. In this context, it is necessary to know which quantitative information could be indeed extracted from a numerical solution.

In this paper, we are concerned with two classes of variational inequalities generated by obstacle type conditions. Differentiability properties of exact solutions to these problems are, in general, restricted even if all external data of a problem are smooth (e.g., see the works of H. Brezis [1], L.A. Caffarelli [7], D. Kinderlehrer and G. Stampacchia [16], A. Friedman [10], N. N. Uraltseva [28]). In [2] it was proved that there exists a unique solution $u\in W^{2,2}(\Omega)$ of an obstacle problem

[TABLE]

if $\psi\in W^{2,2}(\Omega)$, $f\in L^{2}$, the function $u_{D}$ (which defines the Dirichlét boundary condition) belongs to $W^{2,2}(\Omega)$ and satisfies the natural condition $u_{D}\geq\psi$ on $\partial\Omega$.

Many researches were focused on clarifying mathematical properties of the coincidence set. In particular, it was proved that if the domain $\Omega\subset R^{2}$ is strictly convex with a smooth boundary $\partial\Omega$ and if the obstacle $\psi\in C^{2}(\Omega)$ is strictly concave, then the coincidence set is connected and its boundary is smooth and homeomorphic to the unit circle (see, e.g. [16]). However, in general, the structure of a coincidence set can be very complicated and for any domain one can point out such an obstacle that this set has any number of disjoint subsets.

Numerical methods for problems with obstacles (and many other problems related to variational inequalities) were systematically studied in R. Glowinski, J.-L. Lions, and R. Tremolieres [12, 11]. Getting the respective a priori rate convergence estimates (in terms of the mesh size $h$) was the first question studied by many authors. In the context of finite element approximations such type estimates were derived by R. S. Falk [9] who proved the standard a priori convergence error estimates (with the rate $h$ for the $L^{2}$ norm of gradients and the rate $h^{2}$ for the $L^{2}$ norm of the functions) provided that $u\in W^{2,2}$. Convergence of mixed methods for problems with obstacles was established in F. Brezzi, W. W. Hager and P. A. Raviart [3] and numerical methods based on the augmented Lagrangian approach were studied in T. Kärkkäinen, K. Kunisch, and P. Tarvainen [15].

This paper is concerned with other important questions arising in quantitative analysis of nonlinear problems. One of them is which measure $\mathds{M}$ of the distance to the exact solution is adequate (natural) for a particular problem? (see a discussion in [23]). Furthermore, we must know which properties of a solution are controlled by $\mathds{M}$ and deduce explicitly computable bounds (minorants and majorants). In the paper, we study these questions in the context of obstacle type problems. Our analysis is based upon general type error identities derived in [17, 20, 21] for a wide class of convex variational problems. These identities establish equivalence of a certain nonlinear measure $\mathds{M}$ and the duality gap between the primal and dual energy functionals. Since variational methods are based on minimization of this gap, the measure $\mathds{M}$ shows limits of quantitative analysis for this class of methods.

For convenience of the reader we shortly recall the main items necessary for understanding of the material. Consider the class of variational problems

[TABLE]

where $\Lambda:Y^{*}\rightarrow{\mathbb{R}}$ is a bounded linear operator, $G:Y\rightarrow{\mathbb{R}}$ is a convex, coercive, and lower semicontinuous functional, $F:V\rightarrow{\mathbb{R}}$ is another convex lower semicontinuous functional, and $Y$ and $V$ are reflexive Banach spaces. The dual spaces are denoted by $Y^{*}$ and $V^{*}$, respectively, and the duality pairings are denoted by $(y^{*},y)$ and $\left<v^{*},v\right>$. The dual variational problem consists of finding $p^{*}\in Y^{*}$ maximizing the dual functional

[TABLE]

over the space $Y^{*}$. Here $G^{*}:Y^{*}\rightarrow{\mathbb{R}}$ and $F^{*}:V^{*}\rightarrow{\mathbb{R}}$ are the Young-Fenchel transforms (convex conjugates) of $G$ and $F$, respectively. Henceforth we use are the so called compound functionals

[TABLE]

generated by the convex functionals $F$ and $G$, respectively. These functionals are nonnegative and vanish if and only if $v$ and $v^{*}$ (resp. $y$ and $y^{*}$) are joined by special differential relations (see, e.g., [18]). Notice that in the simplest case where $V$ is a Hilbert space and $F(v)=\frac{1}{2}\|v\|^{2}$, the functional $D_{F}(v,v^{*})$ coincides with the norm $\frac{1}{2}\|v-v^{*}\|^{2}$. However, in general $D_{F}(v,v^{*})$ should be viewed as a nonlinear measure, which vanishes if and only if the pair $(v,v^{*})$ satisfies certain conditions.

Let $y^{*}\in Y^{*}$ and $v\in V$ be the functions compared with $p^{*}$ and $u$. Introduce the following (nonlinear) measure of the distance between $\{u,p^{*}\}$ and $\{v,y^{*}\}$:

[TABLE]

It vanishes if and only if

[TABLE]

The above conditions are satisfied if and only if $v=u$ and $y^{*}=p^{*}$ (i.e., if approximations coincide with the exact primal and dual solutions). In [20] and [17] (Section 7.2), it was proved that

[TABLE]

Hence ${\mathds{M}}\{(u,p^{*}),(v,y^{*})\}=0$ if and only if $J(v)=I^{*}(y^{*})$ (what means that $v$ is a minimizer of the problem ${\mathcal{P}}$ and $y^{*}$ is a maximizer of the problem ${\mathcal{P}}^{*}$).

Two particular forms of (4) arise if we set $v=u$ or $y^{*}=p^{*}$. They are ${\mathds{M}}(u,v):={\mathds{M}}(\{u,p^{*}\},\{v,p^{*}\})$ and ${\mathds{M}}(p^{*},y^{*}):={\mathds{M}}(\{u,p^{*}\},\{u,y^{*}\})$. In view of (4),

[TABLE]

Numerical methods are based either on minimization of the primal energy, or maximization of the dual energy, or on coupled minimization–maximization of both. The identities (4), (5), and (6) show that the functional ${\mathds{M}}(\{u,p^{*}\},\{v,y^{*}\})$ (and its particular forms ${\mathds{M}}(u,v)$ and ${\mathds{M}}(y^{*},p^{*})$) are in fact the error measures used by energy based numerical procedure designed to solve (1). Since the error measures are equal to the respective duality gaps, they present the strongest (and in a sense the most natural) measure for the class of problems considered.

Below we study these identities for two classes of nonlinear variational problems and show that they generated specific error measures containing two parts. The first part is presented by a norm equivalent to $H^{1}$ norm and the second one is a nonlinear measure, which controls (in a rather weak sense) how accurately an approximate solution recovers configuration of the free boundary. We deduce directly computable quantities which majorate the right hand sided of (4), (5), and (6). Furthermore, we prove that the majorants are sharp, i.e., they do not contain an irremovable gap between the left and right hand sides. The majorants possesses other important properties, namely, they need no a priori knowledge about the shape of a coincidence set, valid for any approximations of the admissible functional (energy) set, and do not contain unknown (e.g., interpolation) constants. In the last section of the paper, we collect computational results aimed to confirm theoretical analysis. They are mainly focused on two points. First we show that the measures correctly represent the quality of approximations for various minimizing sequences. Another observation is that for different sequences different parts of the measure may dominate, but their sum always correctly represent the error and can be efficiently estimated from above by the majorant.

2. Classical obstacle problem

2.1. Variational setting

We begin with the classical obstacle problem (see, e.g. [1, 10, 16]), where admissible functions belong to the set

[TABLE]

Here, $H^{1}_{0}(\Omega)$ denotes the Sobolev space of functions vanishing on $\partial\Omega$ (hence we consider the case $u_{D}=0$), $\Omega\subset\mathbb{R}^{d}$ ($d\in\{1,2,3\}$) is a bounded domain with a Lipschitz continuous boundary $\partial\Omega$ and $\phi,\psi\in H^{2}(\Omega)$ are two given functions (lower and upper obstacles) such that

[TABLE]

The problem is to find $u\in K$ satisfying the variational inequality

[TABLE]

for a given function $f\in L^{2}(\Omega)$ and a bilinear form

[TABLE]

It is assumed that $A$ is a symmetric matrix subject to the condition

[TABLE]

almost everywhere in $\Omega$. Under the assumptions made, the unique solution $u\in K$ exists. In general, the solution $u$ divides $\Omega$ into three sets:

[TABLE]

The sets $\Omega^{u}_{-}$ and $\Omega^{u}_{+}$ are the lower and upper coincidence sets and $\Omega^{u}_{0}$ is an open set, where $u$ satisfies the Poisson equation ${\rm div}(A\nabla u)+f=0$. Thus, the problem involves free boundaries, which are unknown a priori. Let $v$ be an approximation of $u$. It defines approximate sets

[TABLE]

Notice that unlike the sets in (2.1), the sets (2.1) are known.

Solution of the problem (7) can be represented in a mixed form, i.e., as a pair $(u,p^{*})$, where the flux

[TABLE]

satisfies the conditions

[TABLE]

The pair $(u,p^{*})\in K\times L^{2}(\Omega,\mathbb{R}^{d})$ is a saddle point of the respective minimax formulation. Under the above made assumptions it exists. Moreover, $p^{*}$ has square summable divergence and satisfies the relations (11) and (2.1) almost everywhere in $\Omega$.

2.2. Error measures

The variational inequality (7) is known to have the equivalent form (1) for

[TABLE]

where $\chi_{K}$ is the characteristic functional of the set $K$, i.e.,

[TABLE]

In this case, $V=V_{0}$, $V^{*}=H^{-1}(\Omega)$, $Y=L^{2}(\Omega,\mathbb{R}^{d})$

[TABLE]

and

[TABLE]

For $y^{*}=p^{*}$ and for $v=u$, we obtain

[TABLE]

Next, for $v^{*}\in L^{2}(\Omega)$,

[TABLE]

Here, $(z)_{-}$ and $(z)_{+}$ denote the negative and positive parts of the quantity $z$, i.e., $(z)_{-}:=-\min\{0,z\},(z)_{+}:=\max\{0,z\}$. They satisfy the relations $z=-(z)_{-}+(z)_{+}$ and $|z|=(z)_{-}+(z)_{+}$.

In view of (18), we deduce explicit form of the functional $D_{F}$ provided that $y^{*}\in Y^{*}_{\rm div}(\Omega):=\left\{y^{*}\in Y^{*}\,\mid\,{\rm div}y^{*}\in L^{2}(\Omega)\right\}$:

[TABLE]

Since $p^{*}$ belongs to $Y^{*}_{\rm div}(\Omega)$ and satisfies the relation (2.1), we find that

[TABLE]

This quantity can be viewed as a certain measure

[TABLE]

where ${\sf W}_{\phi}:=-({\rm div}A\nabla\phi+f)$, ${\sf W}_{\psi}:={\rm div}A\nabla\psi+f$ are two nonnegative weight functions generated by the source term $f$, the obstacles $\psi,\phi$ and the diffusion $A$. It is clear that $\mu_{\phi\psi}(v)=0$ if $\Omega^{v}_{-}\subset\Omega^{u}_{-}$ and $\Omega^{v}_{+}\subset\Omega^{u}_{+}$. In other words, if all points of approximate sets $\Omega^{v}_{-}$ and $\Omega^{v}_{+}$ indeed belong to the coincidence sets, then the measure is zero.

Remark 1**.**

Assume that $A={\mathbb{I}}$ (the identity matrix), obstacles $\phi,\psi$ are harmonic functions ($\triangle\phi=\triangle\psi=0$ in $\Omega$) satisfying $\phi<0<\psi$ almost everywhere in $\Omega$ and $f=const\not=0$. If $f>0$ then $\Omega^{u}_{-}=\emptyset$ (the lower obstacle $\phi$ is never active) and

[TABLE]

Here, we decomposed

[TABLE]

and applied the equality $\|\psi-v\|_{L^{1}(\Omega^{u}_{+})}=\|\psi-v\|_{L^{1}(\Omega^{u}_{+}\setminus\Omega^{v}_{+})}$ (which holds because $\|\psi-v\|_{L^{1}(\Omega^{u}_{+}\cap\Omega^{v}_{+})}={0}$). Analogously, if $f<0$ then $\Omega^{u}_{+}=\emptyset$ (the upper obstacle $\psi$ is neven active) and

[TABLE]

We see that $\mu_{\phi\psi}(v)$ represents a certain measure, which controls (in a weak integral sense) whether or not the function $v$ coincides with obstacles $\psi,\phi$ on true coincidence sets $\Omega^{u}_{-}$ and $\Omega^{u}_{+}$.

Analogously, the quantity

[TABLE]

forms another measure

[TABLE]

where the sets

[TABLE]

are approximations of $\Omega_{-}$, $\Omega_{0}$, and $\Omega^{u}_{-},\Omega^{u}_{+},\Omega^{u}_{0}$ generated on the basis of dual solution $y^{*}$. It is clear that this measure is zero if $\Omega^{y^{*}}_{-}\subset\Omega^{u}_{-}$ and $\Omega^{y^{*}}_{+}\subset\Omega^{u}_{+}$. Hence, the measure $\mu_{\phi\psi}^{*}(y^{*})$ is positive if the sets $\Omega^{y^{*}}_{-}$ and $\Omega^{y^{*}}_{+}$ contain parts which do not belong to true coincidence sets. We summarize properties of $\mu_{\phi\psi}(v)$ and $\mu_{\phi\psi}^{*}(y^{*})$ as follows:

[TABLE]

Now we use (4), (5), and (6) and deduce error identities for the obstacle problem.

Theorem 1** (energy identities for the classical obstacle problem).**

Let $v$ and $y^{*}$ be approximations of $u$ and $p^{*}$, respectively. Then,

[TABLE]

Theorem 1 establishes exact error identities for the classical obstacle problem in terms the primal and dual posings. In view of the relation between the primal and dual functionals, the identities (29) and (30) yield

[TABLE]

This error identity holds for the mixed nonlinear measure ${\mathds{M}}(\{u,p^{*}\},\{v,y^{*}\})$ (which decomposes additively to two primal nonlinear measures). It shows that the duality gap consists of four nonnegative quantities. Two of them are quadratic terms associated with energy errors. Two others are nonlinear measures $\mu_{\phi\psi}(v)$ and $\mu_{\phi\psi}^{*}(y^{*})$ defined by (21) and (25) Without taking them into account, only inequalities

[TABLE]

can be obtained.

2.3. Computable bounds of error measures

First we show that the measure ${\mathds{M}}(\{u,p^{*}\},\{v,y^{*}\})$ can be directly computed for any pair of approximate solutions $\{v,y^{*}\}$ provided that $y^{*}$ possesses an additional regularity.

Theorem 2**.**

Let $\{v,y^{*}\}\in K\times Y^{*}_{{\rm div}}(\Omega)$. Then,

[TABLE]

where

[TABLE]

Proof.

In view of (2)), (14), and (18), we have

[TABLE]

According to (31),

[TABLE]

Since

[TABLE]

and

[TABLE]

the substitution of last two equalities in (36) yields (32). ∎

Remark 2**.**

Assume that the right hand side of (32) is equal to zero. Then $y^{*}=A\nabla v$ and

[TABLE]

Hence, $\Omega^{y^{*}}_{-}\subset\Omega^{v}_{-}$ and $\Omega^{y^{*}}_{+}\subset\Omega^{v}_{+}$. The sets $\Omega^{v}_{+}$ and $\Omega^{v}_{-}$ do not intersect as well as the sets $\Omega^{y^{*}}_{+}$ and $\Omega^{y^{*}}_{-}$. Therefore, the set $\Omega^{v}_{0}=\Omega\setminus(\Omega^{v}_{+}\cup\Omega^{v}_{-})$ is contained in the set $\Omega^{y^{*}}_{0}=\Omega\setminus(\Omega^{y^{*}}_{+}\cup\Omega^{y^{*}}_{-})$. Thus, ${\rm div}y^{*}+f=0$ in $\Omega^{v}_{0}$. For any $w\in K$, we have

[TABLE]

The right hand side of the above relation is nonnegative. Indeed, the first two integrals are nonnegative and the last one is equal to zero. This means that $v$ satisfies the variational inequality and, consequently, the pair $\{v,y^{*}\}$ coincides with $\{u,p^{*}\}$.

Remark 3**.**

If approximations of the coincidence sets (constructed on the basis of $v$ and $y^{*}$) satisfy the relations $\Omega^{y^{*}}_{-}\subset\Omega^{v}_{-}$ and $\Omega^{y^{*}}_{+}\subset\Omega^{v}_{+}$, then (32) reads

[TABLE]

Moreover, if $\Omega^{y^{*}}_{-}\subset\Omega^{u}_{-}\subset\Omega^{v}_{-}$ and $\Omega^{y^{*}}_{+}\subset\Omega^{u}_{+}\subset\Omega^{v}_{+},$ then nonlinear terms of ${\mathds{M}}(\{u,p^{*}\},\{v,y^{*}\})$ vanish and we arrive at the equality

[TABLE]

However, the sets $\Omega^{u}_{-}$ and $\Omega^{u}_{+}$ are unknown, so that in practice it is impossible to verify the conditions that yield this simplest (hypercircle type) form of the error identity.

Theorem 2 provides a way to compute ${\mathds{M}}(\{u,p^{*}\},\{v,y^{*}\})$, which is the sum of error measures ${\mathds{M}}(u,v)$ and ${\mathds{M}}(p^{*},y^{*})$. These measures separately evaluate deviations of $v$ from $u$ and $y^{*}$ from $p^{*}$. It is desirable to have guaranteed bounds for them as well (notice that in view of (32) two sided bounds of ${\mathds{M}}(u,v)$ imply two sided bounds of ${\mathds{M}}(p^{*},y^{*})$ and vise versa). For this purpose, we require knowledge of the exact energy $J(u)$ (or $I^{*}(p^{*})$), which is generally unknown. However, their is a way to derive computable bounds of ${\mathds{M}}(u,v)$ without this knowledge (see [19, 22]). In this section, we briefly discuss some of them addressing the reader to a more systematic exposition and numerical tests to the above cited literature and [17].

The first bound of ${\mathds{M}}(u,v)$ has the form

[TABLE]

The majorant ${\mathds{M}}^{+}$ contains contains free variables: $\beta>0$, $y^{*}\in Y^{*}_{{\rm div}}(\Omega)$, and two nonnegative functions (Lagrange multipliers) $\lambda_{1},\lambda_{2}\in L^{2}(\Omega)$. The constant $C_{\Omega}>0$ is a minimal constant in a Friedrichs type inequality

[TABLE]

It is not difficult to show that for any $v$, there exist $\beta$, $\lambda_{1}$, $\lambda_{2}$, and $y^{*}$ such that (38) holds as the equality. Indeed, set $y^{*}=p^{*}$, and

[TABLE]

Then, the second term of ${\mathds{M}}^{+}$ vanishes (for any choice of $\beta$) and the third term is equal to $\mu_{\phi\psi}(v)$. By taking a limit $\beta\rightarrow+\infty$, the first term converges to

[TABLE]

The choice (40) of Lagrange multipliers is theoretically important since it depends on the exact solution $u$. It is replaced by different choices in practical computations. If we set alternatively

[TABLE]

the third term of (38) vanishes and we obtain another majorant (which is free of $\lambda_{1},\lambda_{2}$)

[TABLE]

where

[TABLE]

More accurate optimization of (38) with respect to $\lambda_{1},\lambda_{2}$ provides a sharper majorant [22] in the form

[TABLE]

where

[TABLE]

Practical computations of majorants for the classical obstacle problem are further explained in [6, 13, 14].

Remark 4**.**

Since $J(v)-J(w)\geq J(v)-J(u)$ holds for all $w\in K$, we always have a computable lower bound

[TABLE]

In practice, a suitable $w$ can be constructed by local (e.g., patch wise) improvement of $v$ and ideas of hierarchical basis methods.

3. Double obstacle problem

3.1. Variational setting

The following double-obstacle problem (also known as the two–phase obstacle problem), was studied in H. Shahgholian, N. N. Uraltseva, and G. S.Weiss [25], N.N. Uraltseva [29], G. S. Weiss [27] and some other papers cited therein. Here the variational (energy) functional $J(v)$ is defined by the relation

[TABLE]

The functional $J(v)$ is minimized on the set

[TABLE]

Here $u_{D}$ is a given bounded function that defines the boundary condition ($u_{D}$ may attain both positive and negative values on different parts of the boundary $\partial\Omega$). It is assumed that the coefficients $\alpha_{+},\alpha_{-}:\Omega\rightarrow\mathbb{R}$ are positive constants (without essential difficulties the consideration and main results can be extended to the case where they are positive Lipschitz continuous functions). Also, it is assumed that $f\in L^{\infty}(\Omega)$, $A\in L^{\infty}(\Omega,\mathbb{R}^{d\times d})$, and the condition (8) holds. Since the functional $J(v)$ is strictly convex and continuous on $V$, existence and uniqueness of a minimizer $u\in K$ is guaranteed by well known results of the calculus of variations (see, e.g., [18]). Analysis of the corresponding Euler-Lagrangian equation leads to the nonlinear problem ([25, 27, 29])

[TABLE]

where $\chi$ denotes the characteristic function of a set (attaining values 1 and 0 inside and outside the set, respectively). A physical interpretation of the problem (51) is presented by an elastic membrane touching the planar phase boundary between two liquid/gaseous phases (see, e.g., [25]).

We introduce two decompositions of $\Omega$ associated with the minimizer $u$ and an approximation $v$:

[TABLE]

and

[TABLE]

These decompositions generate exact and approximate free boundaries. Using the above notation we can rewrite (51) as follows

[TABLE]

3.2. Error measures

The problem is reduced to (1) if $V=V_{0}:=H^{1}_{0}(\Omega)$, $Y=L^{2}(\Omega,\mathbb{R}^{d})$, $\Lambda w=\nabla w$, $\Lambda^{*}y^{*}=-{\rm div}y^{*}$, and the functionals

[TABLE]

stand for $G$ and $F$, respectively. The problem is to find $u_{0}\in V_{0}$ such that the functional $\widehat{J}(v_{0})=\widehat{G}(\nabla v_{0})+\widehat{F}(v_{0})$ attains infimum on the space $V_{0}$.

[TABLE]

Hence,

[TABLE]

for any $v=v_{0}+u_{D}$. Computation of $\widehat{F}^{*}(v^{*})$ is more sophisticated.

Lemma 1**.**

Let $v^{*}\in L^{\infty}(\Omega)$. Then,

[TABLE]

Proof.

Assume that $v^{*}+f>\alpha_{+}$ on some open subset $\omega\subset\Omega$. Then this inequality holds on a ball $B\subset\omega$. Define two smooth cut off functions $\lambda^{\epsilon}_{1}$ and $\lambda^{\epsilon}_{2}$ such that

[TABLE]

Here $\epsilon$ is a positive quantity smaller than $\frac{1}{2}{\rm dist}(B,\partial\Omega)$. For any $\rho\in{\mathbb{R}}$, the function $v^{\epsilon}:=\lambda^{\epsilon}_{1}u_{D}+\rho\lambda^{\epsilon}_{2}$ belongs to $V_{0}+u_{D}$. It is not difficult to see that

[TABLE]

and

[TABLE]

Therefore,

[TABLE]

Let $\epsilon\rightarrow 0$ and $\rho\rightarrow+\infty$. Then the first integral in the right hand side vanishes, the second is positive and the third tends to $+\infty$. Hence, $F^{*}(v^{*})=+\infty$.

Quite analogously we prove that $F^{*}(v^{*})=+\infty$ if $v^{*}+f<\alpha_{-}$ on some open set $\omega\subset\Omega$. It remains to show that $F^{*}(v^{*})=-\int\limits_{\Omega}v^{*}u_{D}\,dx$ if $-\alpha_{-}\leq v^{*}+f\leq\alpha_{+}$. For this purpose, we define $v^{\epsilon}:=\lambda^{\epsilon}_{1}u_{D}$. In this case,

[TABLE]

We see that the first two integrals are nonpositive, so that $\widehat{F}^{*}(v^{*})\leq-\int\limits_{\Omega}v^{*}u_{D}\,dx$. On the other hand,

[TABLE]

as $\epsilon\rightarrow 0$ and we arrive at (59). ∎

Corollary 1**.**

If $v^{*}$ satisfies $-\alpha_{-}\leq v^{*}+f\leq\alpha_{+}$, then

[TABLE]

where $v=v_{0}+u_{D}$. Hence, if

[TABLE]

then

[TABLE]

To obtain error identities, we need to express (64) for two particular cases where $y^{*}=p^{*}$ and $v=u$. For the first case, we have

[TABLE]

Since $p^{*}=A\nabla u$ the relation (54) guarantees that ${\rm div}p^{*}+f\in[-\alpha_{-},-\alpha_{+}]$ almost everywhere in $\Omega$ and, therefore, $p^{*}\in Y^{*}_{{\rm div},[\alpha_{-},\alpha_{+}]}$. Introduce the sets

[TABLE]

which qualify the difference between exact coincidence sets and those formed by $v$ (see Fig. 1). The remaining part $\widehat{\Omega}:=\Omega\setminus\omega$ (where $\omega:=\omega_{+}\cup\omega_{-}\cup\omega_{\pm}$) contains the points of $\Omega$ which belong to $\Omega^{v}_{+}\cap\Omega^{u}_{+}$ or $\Omega^{v}_{0}\cap\Omega^{u}_{0}$. In view of (54), at these points integrands of (65) vanish and we obtain

[TABLE]

where

[TABLE]

The right hand side of (66) is a nonnegative functional (measure), which is equal to zero if $\Omega^{v}_{+}$ coincides with $\Omega^{u}_{+}$ and $\Omega^{v}_{-}$ coincides with $\Omega^{u}_{-}$.

For the second case, we have

[TABLE]

Again we may view the right hand side as a certain measure, which is zero if the sets

[TABLE]

coincide with the sets $\Omega^{u}_{-}$ and $\Omega^{u}_{+}$, respectively.

Now (4), (5), and (6) imply the following result.

Theorem 3**.**

Let $v\in V_{0}+u_{D}$ and $y^{*}\in Y^{*}_{{\rm div},[\alpha_{-},\alpha_{+}]}$ be approximations of $u$ and $y^{*}$, respectively. Then

[TABLE]

where

[TABLE]

is a nonnegative functional, which vanishes if $y^{*}=p^{*}$ and $v=u$.

Proof.

We apply (5) and (6). Notice that $\widehat{J}(v_{0})=G(\nabla v_{0})+F(v_{0})=J(v)$. Next,

[TABLE]

It is easy to see that for any $v=v_{0}+u_{D}\in V_{0}+u_{D}$, the functional

[TABLE]

coincides with $\widehat{J}(v_{0})$ and $\widehat{J}(u_{0})$ coincides with $J(u)$. Since

[TABLE]

we arrive at (73).

Since $u_{0}=u-u_{D}$ (where $u$ satisfies the relation $A\nabla u=p^{*}$), we use (56) and (72) and obtain

[TABLE]

Now (6) yields (74), where

[TABLE]

Finally, summation of (73) and (74) yields

[TABLE]

where

[TABLE]

∎

Corollary 2**.**

From (75) it follows that

[TABLE]

This inequality has a practical value because it provides a directly computable upper bound of the error.

Remark 5**.**

It is not difficult to show that $\Upsilon(v,y^{*})=0$ if and only if the set $\Omega^{y^{*}}_{<>}:=\Omega\setminus\left\{\Omega^{y^{*}}_{-}\cup\Omega^{y^{*}}_{+}\right\}$ (in this set $-\alpha_{-}<{\rm div}y^{*}+f<\alpha_{+}$) is a subset of $\Omega^{v}_{0}$, $v$ does not have positive values in $\Omega^{y^{*}}_{-}$ and negative values in $\Omega^{y^{*}}_{+}$. To prove this we represent $\Upsilon(v,y^{*})$ in the form

[TABLE]

where the terms are defined by the relations

[TABLE]

with the weights $W_{+}(y^{*})=(\alpha_{+}-f-{\rm div}y^{*})$ and $W_{-}(y^{*})=\alpha_{-}+f+{\rm div}y^{*}$. The term $\Upsilon_{1}(v,y^{*})$ vanishes if $v\leq 0$ in $\Omega^{y^{*}}_{-}$ and $v\geq 0$ in $\Omega^{y^{*}}_{+}$. In the set $\Omega^{y^{*}}_{<>}$ the weights $W_{+}(y^{*})$ and $W_{-}(y^{*})$ are positive. Therefore, $\Upsilon_{2}(v,y^{*})=\Upsilon_{3}(v,y^{*})=0$ implies $v=0$ almost everywhere in $\Omega^{y^{*}}_{<>}$, i.e., $\Omega^{y^{*}}_{<>}\subset\Omega^{v}_{0}$. If all the above conditions are satisfied, then $\Upsilon(v,y^{*})=0$ and we arrive at the identity

[TABLE]

It is clear that $\Upsilon(v,y^{*})=0$ if the set $\Omega^{v}_{-}$ coincides (up to a set of zero measure) with the set $\Omega^{y^{*}}_{-}$ and $\Omega^{v}_{+}$ coincides $\Omega^{y^{*}}_{+}$.

Remark 6**.**

Computable upper bound of the primal error measure ${\mathds{M}}(u,v)$ was first derived in [24]. It has the form

[TABLE]

The majorant ${\mathds{M}}^{+}$ contains contains free variables: $\beta>0$, $y^{*}\in Y^{*}_{{\rm div}}(\Omega)$, and two nonnegative functions (Lagrange multipliers) $\lambda_{+},\lambda_{-}\in L^{2}(\Omega)$ satisfying $\lambda_{+}(x),\lambda_{-}(x)\in[0,1]$ almost for all $x\in\Omega$. The constant $C_{\Omega}>0$ is given by (39). In practical computations [5] it is convenient to simplify ${\mathds{M}}^{+}(v;\beta,\lambda_{+},\lambda_{-},y^{*})$ to

[TABLE]

where only one multiplier $\lambda\in L^{2}(\Omega)$ satisfying $\lambda\in[-\alpha_{-},\alpha_{+}]$ almost for all $x\in\Omega$ is required.

4. Numerical verifications of the error identities

4.1. The classical obstacle problem

We consider an example from [13] with known exact solution. Here,

[TABLE]

and $u$ satisfies the homogeneous Dirichlet boundary conditions $u(0)=0,u(1)=0.$ The exact solution $u$ is in the form

[TABLE]

where $r=r_{f,\phi}:=\frac{1}{2}-\sqrt{\frac{2\phi}{f}}\in(0,\frac{1}{2})$. The parameter $r$ determines the radius of the exact lower coincidence set

[TABLE]

It is easy to show that the exact energy reads

[TABLE]

An approximation $v_{\epsilon_{1}}$ is considered in the form of $u$ corresponding to the same value of $\phi$ and a perturbed value $f$,

[TABLE]

for some small perturbation $\epsilon_{1}$. This choice ensures

[TABLE]

and in particular, $\Omega^{v_{\epsilon_{1}}}_{-}\subset\Omega^{u}_{-}$ for $\epsilon_{1}\in(0,r)$. An example of $u$ and $v_{\epsilon_{1}}$ is depicted in the top left picture of Figure 2. An approximation $y^{*}_{\epsilon_{2}}$ is taken as

[TABLE]

where $I$ denotes a piecewise linear nodal and continuous interpolation operator at nodes

[TABLE]

for some small positive perturbation $\epsilon_{2}$. The approximation $y^{*}_{\epsilon_{2}}$ differs from the exact flux $p^{*}$ only locally in $(\frac{1}{2}-r-\epsilon_{2},\frac{1}{2}-r+\epsilon_{2})\cup(\frac{1}{2}+r-\epsilon_{2},\frac{1}{2}+r+\epsilon_{2})$ and

[TABLE]

An example of $p^{*}$ and $y^{*}_{\epsilon_{2}}$ is shown in the top right picture of Figure 2 and corresponding equilibrium terms ${\rm div}p^{*}+f$ and ${\rm div}y^{*}_{\epsilon_{2}}+f$ in the bottom left picture.

For numerical verifications, we choose parameters

[TABLE]

resulting in $r\approx 0.1220$, $J(u)\approx-6.9446$ and approximations $v_{\epsilon_{1}},y^{*}_{\epsilon_{2}}$ corresponding to

[TABLE]

We first verify the primal error identity

[TABLE]

for all approximations $v_{\epsilon_{1}}$. Table 3 confirms that the primal error identity holds and both quadratic (gradient containing) and nonlinear parts of the primal error converge. For smaller values of $\epsilon_{1}$ the quadratic part dominates over the nonlinear part. This is due to the fact that the quadratic part of error is globally distributed over $\Omega$ and the nonlinear part $\mu_{\phi\psi}(v_{\epsilon_{1}})$ has a support in

[TABLE]

Table 3 verifies the dual error identity

[TABLE]

for all approximations $y^{*}_{\epsilon_{2}}$. Again, both quadratic and nonlinear parts converge. None of error parts dominates, since $y^{*}_{\epsilon_{2}}$ and $p^{*}_{\epsilon_{2}}$ differ only locally. The nonlinear part $\mu_{\phi\psi}^{*}(y^{*}_{\epsilon_{2}})$ has a support in

[TABLE]

An example of primal and dual nonlinear error functions is depicted in the bottom left picture of of Figure 2.

Table 3 verifies the majorant identity

[TABLE]

where the computable nonlinear majorant part $\Upsilon$ is given by (33). The majorant identity is valid for all considered approximations.

Remark 7**.**

Since the upper obstacle $\psi$ is not considered in this example,

[TABLE]

has to be satisfied. This condition is fulfilled for $y^{*}_{\epsilon_{2}}$ constructed above.

4.2. The double obstacle problem

We consider an example with known exact solution. Here,

[TABLE]

and $u$ satisfies Dirichlet boundary conditions $u(-1)=-1,u(1)=1.$ This example generalizes example of [4], in which $\alpha_{+}=\alpha_{-}=8$. It is possible to show the exact solution is given by a formula

[TABLE]

where $r_{-}:=\sqrt{\frac{2}{\alpha_{-}}}-1\in(-1,0)$ and $r_{+}:=1-\sqrt{\frac{2}{\alpha_{+}}}\in(0,1)$ determine exact coincidence sets

[TABLE]

The exact energy then reads

[TABLE]

An approximation $v_{\epsilon_{1}}$ is considered in the form of $u$ corresponding to perturbed values $\alpha_{+}$, $\alpha_{-}$,

[TABLE]

for some small perturbation $\epsilon_{1}$. This choice ensures

[TABLE]

and in particular, $\Omega^{v_{\epsilon_{1}}}_{0}\subset\Omega^{u}_{0}$ for $\epsilon_{1}\in(0,\min\{r_{+},-r_{-}\})$. An example of $u$ and $v_{\epsilon_{1}}$ is depicted in the top left picture of Figure 3. An approximation $y^{*}_{\epsilon_{2}}$ is taken as

[TABLE]

where $I$ denotes a piecewise linear nodal and continuous interpolation operator at nodes

[TABLE]

for some small positive perturbation $\epsilon_{2}$. The approximation $y^{*}_{\epsilon_{2}}$ differs from the exact flux $p^{*}$ only locally in $(r_{-}-\epsilon_{2},r_{-}+\epsilon_{2})\cup(r_{+}-\epsilon_{2},r_{+}+\epsilon_{2})$ and

[TABLE]

An example of $p^{*}$ and $y^{*}_{\epsilon_{2}}$ is shown in the top right picture of Figure 3 and corresponding equilibrium terms in the bottom left picture.

For numerical verifications, we choose parameters (identical to example of [4])

[TABLE]

resulting in

[TABLE]

and $J(u)=5\frac{1}{3}$ and approximations $v_{\epsilon_{1}},y^{*}_{\epsilon_{2}}$ corresponding to

[TABLE]

We again verify the primal error identity

[TABLE]

for all approximations $v_{\epsilon_{1}}$. Table 6 confirms that the primal error identity holds and both quadratic (gradient containing) and nonlinear parts of the primal error converge. For smaller values of $\epsilon_{1}$ the quadratic part dominates over the nonlinear part. This is due to the fact that the quadratic part of error is globally distributed over $\Omega$ and the nonlinear part $\mu_{\phi\psi}(v_{\epsilon_{1}})$ has a support in

[TABLE]

Table 6 verifies the dual error identity

[TABLE]

for all approximations $y^{*}_{\epsilon_{2}}$. Again, both quadratic and nonlinear parts converge. None of error parts dominates, since $y^{*}_{\epsilon_{2}}$ and $p^{*}_{\epsilon_{2}}$ differ only locally. The nonlinear part $\mu_{\phi\psi}^{*}(y^{*}_{\epsilon_{2}})$ has a support in

[TABLE]

An example of primal and dual nonlinear error functions is depicted in the bottom left picture of of Figure 3.

Table 6 verifies the majorant identity

[TABLE]

where the computable nonlinear majorant part $\Upsilon$ is given by (76). The majorant identity is valid for all considered approximations.

Acknowledgments

The first author acknowledges the support of the Johann Radon Institute for Computational and Applied Mathematics (RICAM) in Linz, Austria during Special Semester on Computational Methods in Science and Engineering in 2016. The second author has been supported by GA CR through the projects GF16-34894L and 17-04301S.