Boundary element methods for variational inequalities

In this paper we present a priori error estimates for the Galerkin solution of variational inequalities which are formulated in fractional Sobolev trace spaces, i.e. in H˜1/2(Γ)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\widetilde{H}^{1/2}(\Gamma )$$\end{document}. In addition to error estimates in the energy norm we also provide, by applying the Aubin–Nitsche trick for variational inequalities, error estimates in lower order Sobolev spaces including L2(Γ)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L_2(\Gamma )$$\end{document}. The resulting discrete variational inequality is solved by using a semi-smooth Newton method, which is equivalent to an active set strategy. A numerical example is given which confirms the theoretical results.


Introduction
In this paper we are interested in the numerical analysis of the Galerkin boundary element approximation of first kind variational inequalities to find In memoriam Christof Eck (1968Eck ( -2011 We assume that is either a (n − 1)-dimensional Lipschitz manifold in R n , n = 2, 3, or = ∂ is the Lipschitz boundary of a bounded domain ⊂ R n , n = 2, 3. By H 1/2 ( ) we denote the Sobolev space of functions which can be extended by zero when ⊂ is embedded in a closed Lipschitz surface , i.e.
Unique solvability of the first kind variational inequality (1.2) follows by applying standard arguments, see, e.g., [2,11,20,21]. Boundary element error estimates in the energy norm are discussed, e.g., in [13,22,27], related finite element error estimates for Galerkin approximations of variational inequalities formulated in H 1 ( ) are given, for example, in [3,10]. Although the energy error estimate for the boundary element approximation of the variational inequality (1.2) follows similar as for the finite element approximation of a variational inequality in H 1 ( ), we provide a proof for completeness. However, the boundary element error estimate given here differs from the related finite element error estimate due to the different approximation properties of functions defined in a bounded domain , or on its boundary = ∂ . Note that the latter also requires an increased regularity of the function to be approximated on the boundary, as compared to an approximation defined in . The proof as given here is also different as presented, e.g., in [13,27]. In particular we present a generalisation of Cea's lemma in the case of variational inequalities, and we prove a related approximation property.
The main interest of this paper is to provide an error estimate in L 2 ( ). In the case of variational equations these results are due to the well known Aubin-Nitsche trick, see, e.g., [1,16,28]. It seems that related results in the case of variational inequalities are not so well known, and to the best of our knowledge, not available for the problem class as considered in this paper. Note that finite element error estimates in L 2 ( ) for the solution of variational inequalities in H 1 ( ) are given in [24], see also [30].
This paper is organized as follows: in Sect. 2 we describe related complementarity conditions and discuss a rather general regularity result. Moreover, we introduce the boundary element discretization of the variational inequality (1.2). In Sect. 3 we provide an error estimate for the approximate solution in the energy norm · H 1/2 ( ) . The Nitsche trick for variational inequalities to derive an error estimate in L 2 ( ) is considered in Sect. 4. For the solution of the discrete variational inequality we describe a semi-smooth Newton approach in Sect. 5, which is equivalent to an active set strategy, and in Sect. 6 we discuss some applications and provide a numerical example.

Complementarity conditions and discretization of variational inequalities
The aim of this section is to describe the Galerkin discretization of the variational inequality (1.2) by using boundary element methods and to present an equivalent characterization of the unique solution of the discrete variational inequality by means of some discrete complementarity conditions. But first we consider related complementarity conditions in the continuous case.
For u ∈ K being the unique solution of the variational inequality (1.2) we introduce the active and inactive boundary parts as and we define Now we are in a position to state a result which is well known in convex analysis and nonlinear analysis of complementarity problems, see, e.g., [18]. However, we provide a proof for the particular application as considered here. Lemma 2.1 Let u ∈ K be the unique solution of the variational inequality (1.2), and let λ ∈ H −1/2 ( ) be defined as in (2.1). Then there hold the complementarity conditions (

2.2)
Proof We first consider the variational inequality (1.2), i.e. for v = g ∈ K we have i.e. λ ≤ 0 in the sense of H −1/2 ( ). In particular for w := g − u ≥ 0 on we have λ, g − u ≤ 0, and therefore, Note that the complementarity conditions (2.2) are nothing than the Karush-Kuhn-Tucker conditions which describe the saddle point (u, λ) of the Lagrange functional Now we are in a position to state some regularity result for the solution of the variational inequality (1.2).

Theorem 2.2 Let
Proof We consider the complementarity conditions (2.2) on the inactive boundary part in , i.e.
The unique solution u ∈ K of the variational inequality (1.2) can be written as u = g + w for w ∈ H 1/2 ( in ) (and satisfying w < 0 on in ). Hence we find that w ∈ H 1/2 ( in ) is the unique solution of the operator equation Note that A : H 1/2 ( in ) → H −1/2 ( in ) is bounded and H 1/2 ( in )-elliptic, and hence invertible. From the mapping properties of A and the assumptions made on g and f it follows that w ∈ H 1/2+s ( in ). To conclude w ∈ H σ ( ) we have to take into account that w is zero on act , and therefore the first derivative of w is in general discontinuous when changing from the inactive to the active part. This gives the restriction σ < 3 2 .
Remark 2.1 Depending on the particular application in mind, i.e. on the particular choice of A, one may conclude some higher values of σ than given in Theorem 2.2.
Next we consider the Galerkin discretization of the variational inequality (1.2) by using boundary element methods. Let S 1 h ( ) = span{ϕ k } M k=1 be the space of piecewise linear and continuous nodal basis functions ϕ k which are defined with respect to an admissible and quasi-uniform boundary element mesh h = ∪ N =1 τ of mesh size h, and with nodal points x k , k = 1, . . . , M. In the three-dimensional case n = 3 we assume that the triangular boundary elements τ are shape regular. Let be the piecewise linear interpolation of the barrier function g which now is assumed to be continuous. Then we define and we consider the variational inequality to find u h ∈ K h such that As in the continuous case we conclude unique solvability of the discrete variational inequality (2.3), see, e.g., [11]. Note that (2.3) is equivalent to the discrete variational inequality to find u ∈ R M ↔ u h ∈ K h such that Then there hold the discrete complementarity conditions (2.5) Note that Lemma 2.3 is the discrete counterpart of Lemma 2.1, and that (2.5) are the Karush-Kuhn-Tucker conditions which are related to the discrete variational inequality (2.4). Since the proof of Lemma 2.3 uses the same arguments as the proof of Lemma 2.1, we omit the proof of Lemma 2.3.

Remark 2.2
The discrete Lagrange multiplier λ = A h u − f ∈ R M is in general not an approximation of the continuous Lagrange multiplier λ = Au − f ∈ H −1/2 ( ). i.e.
with the mass matrix defined by Hence, only when using bi-orthogonal basis functions [17] satisfying ϕ k , ψ = δ k we conclude λ = λ. However, in our approach as presented in this paper we do not consider any approximation of the continuous Lagrange parameter λ, we just introduce the discrete counterpart λ.

Error estimates in the energy norm
In this section we present an a priori error estimate in the energy norm u − u h H 1/2 ( ) for the approximate solution u h ∈ K h of the discrete variational inequality (2.3). While the general idea is similar to the related proof in the case of a finite element approximation [3], the handling of the inequality constraints is rather different. An alternative proof for a particular application, following [10], is discussed in [27]. For energy error estimates for hp boundary element methods in the case of Signorini problems, see also [22].
The first result is the extension of Cea's lemma to variational inequalities.
where the constants c A 1 and c A 2 are the ellipticity and boundedness constants of A as given in (1.3).
Proof From the variational inequality (2.3) we first have where we have used λ ≤ 0 and u h ≤ g h to ensure Taking into account v h = g h on act and λ = 0 on in = \ act we have

Now the assertion follows from the boundedness of A.
It remains to construct v h = u * h ∈ K h in a suitable way to be able to derive an approximation property in H 1/2 ( ). In particular, we define u * h ∈ K h by, see also Note that u * h coincides with g h in all boundary elements τ which include some nonzero measure part of the active boundary act . In all other nodes, u * h is the piecewise linear interpolation of u. We start to give an error estimate for Then there holds the error estimate (i) For τ ⊂ act we have u = g and u * h = g h = I h g. Hence we have, by using a local interpolation error estimate, (ii) Next we consider all boundary elements τ where u * h = I h u is the piecewise linear interpolation of u. As in the first case we obtain It remains to consider two additional cases as depicted in Fig. 1. (iii) We now consider the case of elements τ 1 with τ 1 ⊂ act but meas(τ 1 where the second part again corresponds to the standard local interpolation error. Due to u(x) = g(x) for x ∈ U ε (x * ) we conclude that also the surface gradients of u and g coincide in x * . Hence, the linear Hermite interpolation polynomials I 1 h u = I 1 h g of u and g, which are defined with respect to x * , coincide. With this we conclude again by applying standard interpolation error estimates for Hermite interpolation. (iv) Finally we consider the case where the boundary element τ 2 does not touch the active part act but joins at least one common node x 2 with a boundary element τ 1 , meas(τ 1 ∩ act ) > 0. In this case we have where the first part again is a standard interpolation error estimate. Let η(x) be a sufficiently smooth function satisfying follows. As in the third case we further obtain, by using |η| ≤ 1, For the remaining term we first consider the case σ ≥ 1. By using standard interpolation error estimates we then have In the two-dimensional case n = 2 and for σ ∈ ( 1 2 , 1) we first have and it remains to consider the Sobolev-Slobodeckii norm, i.e. by using the triangle inequality and |η| ≤ 1 we first obtain In fact, for x ∈ τ 2 and for σ ∈ ( 1 2 , 1) we find and therefore follows, since |hη (ξ )| is bounded, and by using the standard interpolation error estimate for the first summand. When combining this estimate with (3.3) we finally conclude Joining all four cases we have shown the desired error estimate.
As in the standard case of a variational equation, and by using an inverse estimate we are now able to give an error estimate in the energy norm.

Lemma 3.3 Let u ∈ K be the unique solution of the variational inequality (1.2).
Let u * h ∈ K h be as constructed in (3.2). Assume u, g ∈ H σ pw ( ) ∩ C( ) for some σ ∈ ( n− 1 2 , 2]. Let the mesh be globally quasi-uniform. Then there holds the error estimate Proof Let P h u ∈ S 1 h ( ) be the Galerkin projection of u defined as the unique solution of the variational problem Using standard techniques, see, e.g. [28], we conclude the error estimates Hence, by using the inverse inequality, and all previous error estimates, we conclude Note that the maximal value of σ as used in the error estimate (3.5) is determined by the choice of the basis functions, i.e. σ ≤ 2 when using piecewise linears, and by the regularity of the solution u ∈ H σ pw ( ). The latter may be obtained either from the properties of the underlying physical problem, or from the mapping properties of the involved operator A, see Theorem 2.2.

Error estimates in L 2 ( ): the Nitsche trick
In this section we present the main result of this paper, the Aubin-Nitsche trick to derive an error estimate in L 2 ( ) for the approximate solution of the variational inequality (1.2). Although the basic ideas of the proof follow the considerations for finite element approximations of variational inequalities in H 1 ( ), see [24,30], the proof as given here requires several considerations which are different.

The adjoint problem
To obtain error estimates in lower Sobolev norms we first assume that A : H 1/2+s ( ) → H −1/2+s ( ) is bijective and bounded for some s ∈ (0, 1 2 ]. For this s we consider a variational inequality which is adjoint with respect to the variational inequality (1.2), and we define where B 1/2−s : H 1/2−s ( ) → H −1/2+s ( ) is the associated Riesz operator. As for the solution u ∈ K of the primal variational inequality (1.2) we first state some regularity result for the solution z ∈ G of the variational inequality (4.1). For this we rewrite the variational inequality as a saddle point problem, i.e. by using the Lagrange multiplier μ ∈ H −1/2 ( act h ), μ ≤ 0 on act h , and α ∈ R + , we introduce the Lagrange functional The Karush-Kuhn-Tucker conditions which are related to the Lagrange functional (4.2) then read to find (z; μ, α) ∈ H 1/2 ( ) × H −1/2 ( ) × R such that

Hence we find from (4.3)
Az and it remains to consider two cases: (i) For α = 0 we find that z ∈ H 1/2 ( \ act,z h ) is the unique solution of the operator equation Az = B 1/2−s (u − u h ) and the assertion follows from the assumptions on A, i.e.
For s ∈ (0, 1 2 ] we finally conclude the assertion. (ii) For α > 0 we conclude λ, z = 0 and it remains to consider the subspace Again the assertion follows from the mapping properties of A and B 1/2−s .

Quasi-interpolation
In order to prove error estimates in lower order Sobolev spaces, for z ∈ G we need to construct a suitable approximation z h ∈ S 1 h ( ) which allows an error estimate in negative Sobolev norms, and which retains the inequality constraints on act h . For this we consider a quasi-interpolation, see also [6,26]. Since the discrete active set act h is given as the union of boundary elements τ , we define a dual boundary mesh as follows, see Fig. 2.
If x k ∈ act h is an interior node of the discrete active set, see for example x 4 and x 5 , the dual element ω k is defined by the midpoints of the adjacent primal elements, in the three-dimensional case we consider the midpoints of the adjacent edges in addition. If x k is a boundary node of the discrete active set, e.g. x 3 , the dual element ω k is the related part of the primal element in act h . For all other nodes we define the dual elements accordingly, where the missing parts of the boundary nodes, e.g., x 3 , are added to the dual element of a related node of the inactive set, see, e.g., ω 2 . For a globally quasi-uniform boundary element mesh h with mesh size h we conclude that all dual elements are of the same mesh size h.
With respect to the dual boundary element mesh we define the piecewise constant L 2 projection satisfying the local error estimate, for z ∈ H 1/2+s ( ), s ∈ (0, 1 2 and by applying the standard Aubin-Nitsche trick For x ∈ ω k ⊂ act h we have z ≤ 0 and therefore z 0 h ≤ 0 follows. Now we are in the position to define the quasi-interpolation z h ∈ S 1 h ( ), Lemma 4.2 Let z ∈ H 1/2+s ( ), s ∈ (0, 1 2 ], and let z h ∈ S 1 h ( ) be the quasiinterpolation as given in (4.9). Then there hold the error estimates z − z h H 1/2 ( ) ≤ ch s |z| H 1/2+s ( ) (4.10) and for x ∈ ω k we then conclude and therefore For what follows, let us recall some basic estimates from the numerical analysis of finite and boundary element methods, see, e.g., [1,28]. Let τ be the reference element to describe all boundary elements τ , in particular for n = 2 we have τ = (0, 1), while for n = 3 we have For x ∈ τ we obtain the local parametrisation x = J (η) and for the measure of a boundary element we conclude

For a boundary element basis function
as well as Hence we conclude the estimate ϕ j L 2 (ω k ) ≤ c h (n−1)/2 , and by using we obtain from the error estimate (4.7) For suppϕ j ∩ ω k = ∅ and x ∈ ω k we have and by using standard techniques we obtain With this we finally conclude the local error estimate and by summation over all elements of the dual mesh we obtain the global error estimate Note that in the three-dimensional case we use that all boundary elements are assumed to be shape regular, so that the number of terms in the local error estimate is bounded. Now the error estimate (4.10) in H 1/2 ( ) follows by using the standard H 1/2 ( ) projection on S 1 h ( ), and the inverse inequality, see, e.g. [28], and the proof of Lemma 3.3.
The error estimate (4.11) in H −1 ( ) follows as the error estimate in L 2 ( ), instead of (4.7) we now use (4.8), and With this we finally conclude the local error estimate and by localizing the dual norm in H −1 ( ) we obtain (4.11).

Error estimates in lower Sobolev norms
Now we are in a position to present error estimates in lower Sobolev norms for the solution u h of the discrete variational inequality (2.3). The main result as given in Theorem 4.5 is based on the following two estimates.

Lemma 4.3 Let u ∈ K and u h ∈ K h be the unique solutions of the variational inequalities (1.2) and (2.3), respectively. Let z ∈ G be the unique solution of the adjoint variational inequality (4.1), and let A :
h ( ) be an appropriate approximation of z ∈ G. Then there holds the error estimate as well as, by using λ ≤ 0 and u h ≤ g h on , Hence we can chose v = z +u −u h +g h −g ∈ G as a test function in the variational inequality (4.1) to obtain Now the assertion follows from the boundedness of A : H 1/2+s ( ) → H −1/2+s ( ), and from the regularity result (4.6).

Lemma 4.4
Let u ∈ K, u h ∈ K h , and z ∈ G be the unique solutions of the variational inequalities (1.2), (2.3), and (4.1), respectively. Let z h ∈ S 1 h ( ) be as constructed in (4.9). Then there holds the estimate (4.13) Proof By using λ k = 0 for u k < g k and λ k ≤ 0, z k ≤ 0 for u k = g k we first have Hence we obtain Now we state the main theorem of this paper. H −1/2+s ( ) be bounded for some s ∈ (0, 1 2 ], and let z h ∈ S 1 h ( ) be given as in (4.9). Assume λ ∈ H σ −1 ( ) for some σ ∈ ( n− 1 2 , 2]. Then there holds the error estimate (4.14) Proof Combing the error estimates (4.12) and (4.13) we first conclude By using the error estimates (4.10) and (4.11), an interpolation argument, as well as the regularity estimate (4.6) we obtain the assertion.
Together with the energy error estimate (3.5) we now obtain the final result of this subsection. In particular for s = 1 2 , i.e. A : H 1 ( ) → L 2 ( ), and when assuming u, g ∈ H 2 pw ( ) and λ ∈ H 1 ( ), we obtain the L 2 ( ) error estimate i.e., as for the standard approximation property of piecewise linear polynomials we can expect a quadratic order of convergence when measuring the error in the L 2 norm.

Semi-smooth Newton method
For the solution of the discrete variational inequality (2.3) we consider the discrete complementarity conditions (2.5), i.e.
which are equivalent to Hence we have to solve the system of (non)linear equations where the nonlinear equations F 2 (u, λ) have to be considered compenentwise. Since is not differentiable in x = 0, we introduce the slant derivative and the application of a semi-smooth Newton method reads where the application of G has to be understood componentwise. The Newton method requires the solution of the linear system

and the first line gives
The second equation gives for all k = 1, . . . , M we then conclude This is just the active set strategy, see, e.g., [5,14,15,17,19].

Screen problem
As a first example we consider the variational inequality to find u ∈ K such that is satisfied for all v ∈ K where = (0, 1 2 ), f ≡ 1, and In (6.1), D : H 1/2 ( ) → H −1/2 ( ) is the hypersingular boundary integral operator defined as Note that the variational inequality (6.1) is related to two-dimensional screen and crack problems, e.g., [29].
Since the variational inequality (6.1) perfectly fits into the framework of the present paper, all theoretical results are valid. In particular, D : H 1/2 ( ) → H −1/2 ( ) is a bounded and H 1/2 ( )-elliptic operator, see [7,23]. Moreover, D : Therefore we can expect almost linear convergence. For a numerical experiment we consider a uniform decomposition of the interval = (0, 1 2 ) into N = 2 L+1 boundary elements. In addition to the variational inequality (6.1) we also consider the solution   Table 1. The numerical results show linear convergence for both problems, as expected. Also the number of semi-smooth Newton iterations indicates super-linear convergence. In Fig. 3 we plot the solutions of the unconstrained and of the constrained screen problem.

Signorini problem
As a second example we consider the Signorini problem for the Laplacian, see, e.g., [27], where the Steklov-Poincaré operator S : H 1/2 ( ) → H −1/2 ( ) can by described, by using boundary integral operators as, see, e.g., [28], Note that V is the single layer integral operator, K is the double layer integral operator and K its adjoint, and D is the hypersingular boundary integral operator. By introducing the solution of the Signorini boundary value problem is equivalent to find the solution u ∈ K of the variational inequality Note that the Steklov-Poincaré operator S is H 1/2 ( N ∪ S )-elliptic which ensures unique solvability of the variational inequality. Assuming sufficient regularity of the given data one can not expect more than u ∈ H 5/2−ε ( ) for the solution of the Signorini boundary value problem, i.e. not more than u | ∈ H 2−ε ( ). In the case of a Lipschitz boundary we find from the mapping properties of all boundary integral operators [7] that S : H 1 ( ) → L 2 ( ). Hence we can apply Corollary 4.6 for s = 1 2 and σ = 2 − ε to conclude the error estimate, for all sufficient small ε > 0, i.e. we can expect almost quadratic convergence when using piecewise linear boundary elements. In [27], a proof of the energy error estimate in H 1/2 ( ) is given for the particular case of a smooth boundary of a bounded domain in two dimensions, and numerical examples are given. For numerical results in the case of contact problems in linear elasticity, see for example [9], and for optimal Dirichlet boundary control problems [25].