Finite element error estimates for normal derivatives on boundary concentrated meshes

# Finite element error estimates for normal derivatives on boundary concentrated meshes

Johannes Pfefferer Technical University of Munich, Department of Mathematics, Chair of Optimal Control, Boltzmannstraße 3, 85748 Garching by Munich, Germany. pfefferer@ma.tum.de    Max Winkler Chemnitz University of Technology, Faculty of Mathematics, Professorship Numerical Mathematics (Partial Differential Equations), Straße der Nationen 62, 09111 Chemnitz, Germany. max.winkler@mathematik.tu-chemnitz.de
###### Abstract

This paper is concerned with approximations and related discretization error estimates for the normal derivatives of solutions of linear elliptic partial differential equations. In order to illustrate the ideas, we consider the Poisson equation with homogeneous Dirichlet boundary conditions and use standard linear finite elements for its discretization. The underlying domain is assumed to be polygonal but not necessarily convex. Approximations of the normal derivatives are introduced in a standard way as well as in a variational sense. On general quasi-uniform meshes, one can show that these approximate normal derivatives possess a convergence rate close to one in as long as the singularities due to the corners are mild enough. Using boundary concentrated meshes, we show that the order of convergence can even be doubled in terms of the mesh parameter while increasing the complexity of the discrete problems only by a small factor. As an application, we use these results for the numerical analysis of Dirichlet boundary control problems, where the control variable corresponds to the normal derivative of some adjoint variable.

f

inite element error estimates, local mesh refinement, boundary concentrated meshes, Dirichlet boundary control, surface flux, normal derivatives

\minisec

AMS subject classification 35J05, 49J20, 65N15, 65N30

## 1 Introduction

The main purpose of this paper is to investigate convergence properties of two types of approximations to the normal derivative of the weak solution of the Poisson equation

 −Δu=finΩ,u=0on∂Ω,

posed in polygonal domains . The first approximation, denoted by , is defined in a classical way, whereas the second one, denoted by , is introduced in a variational sense. Both of them require the knowledge about discrete solutions to the Poisson equation. In this regard and also to illustrate the ideas, we choose standard linear finite elements for the discretization.

In the recent past, error estimates for the two different approximations have been established. In general, the quality of the estimates do not only depend on the regularity of the solution but also on the structure of the underlying computational meshes. We emphasize that in the present case of polygonal domains the regularity of the solution may additionally be lowered by the appearance of corner singularities even though the input datum is arbitrarily smooth. In the following, for a concise discussion of the results from literature, we assume that the regularity of the solution is only limited by the singular terms coming from the corners and not by rough input data.

On general quasi-uniform meshes, the classical and the discrete variational normal derivative converge in with the rate (up to logarithmic factors), provided that the largest interior angle in the domain is less than . For larger interior angles the convergence rate is reduced due to the corner singularities. More precisely, the convergence rate fulfills . The corresponding results for the classical approximation of the normal derivative have been discussed in [14, 20], whereas related results for the discrete variational normal derivative can be found in [4]. On certain superconvergence meshes, where, roughly speaking, neighboring elements almost form a parallelogram, the convergence rate for the discrete variational normal derivative can be improved to (again up to logarithmic factors) if the largest interior angle is less than . Otherwise, the convergence rate satisfies the condition from before, cf. [4]. The convergence rates for the different approximations of the normal derivatives are illustrated in Figure 1 depending on the largest interior angle and the structure of the underlying computational meshes.

As the quantities of interest live on the boundary, it might be promising to appropriately refine the mesh towards the boundary. In this regard, we consider a certain class of boundary concentrated meshes. These are isotropically refined towards the boundary such that the element diameter at the boundary is of order with denoting the maximal element diameter in the interior of the domain. As we will see, the number of elements corresponding to such meshes is of order and no longer of order as in case of quasi-uniform meshes. However, with this slight increase in the number of elements, it is possible to double the convergence rates of the two different approximations in terms of the maximal element diameter (compared to general quasi-uniform meshes). More precisely, the convergence rate in is two (up to logarithmic factors) as long as the largest interior angle is less than . For larger interior angles we obtain a rate fulfilling , see Figure 1 for an illustration.

Our proof of error estimates for the approximating normal derivatives heavily relies on finite element error estimates in weighted -norms. Thereby, the weight is a regularized distance function with respect to the boundary. In order to bound these finite element errors on graded meshes appropriately, one has to be able to handle the weights within the estimates. This requires to establish regularity results in weighted Sobolev spaces with the aforementioned regularized distance function as weight. Based on this, the weighted -errors can then be treated by an adapted duality argument employing a dyadic decomposition of the domain with respect to its boundary and local energy norm estimates on the subsets. These techniques are known for instance from maximum norm error estimates [2, 27, 28] or from finite element error estimates on the boundary for the Neumann problem [6, 29]. However, in all these references the weights, and hence the dyadic decomposition used in the proofs, are related to the corners of the polygonal domain.

As an application of the discrete variational normal derivative we consider Dirichlet boundary control problems with -regularization, where this type of normal derivative naturally arises in the discrete optimality system. In the last decade, Dirichlet boundary control problems have been under active research. We start with mentioning the contribution [10], where a control constrained problem subject to a semilinear elliptic equation is considered. There, a convergence order of is proved for the error of the controls in . This means a rate close to one is only possible if the largest interior angle is less than . However, non-convex domains are excluded in that reference. Later on, in [19], comparable results for the controls are provided in case of linear problems without control constraints. In addition, the authors of this reference show that the states exhibit better convergence properties. The proof relies on a duality argument and estimates for the controls in weaker norms than . We note that, to the best of our knowledge, such an argumentation is restricted to problems without control constraints. For a certain time, this was the state of the art. Nevertheless, numerical experiments indicated that the controls converge with an order close to one also for larger interior angles, and can even achieve a rate close to if the underlying meshes satisfy certain superconvergence properties. For smooth domains , where no corner singularities appear, these convergence rates for general and superconvergence meshes are shown in [11]. Therein, the domain is approximated by a sequence of polygons, on which the discrete approximations are posed. The first contributions dealing with quasi-optimal convergence rates for quasi-uniform and superconvergence meshes in polygonal domains are [3] and [4]. More precisely, in [3], accurate regularity results are derived for the solution of the optimal control problem. In [4], these are applied within the proofs of the error estimates for the control. The rates of convergence for the controls in the unconstrained case coincide with those from above for the discrete variational normal derivative , as such an error for the adjoint state is one of the dominating error contributions. In these references, the control constrained case is discussed is as well. While in convex domains the error estimates for the controls are similar to those in the unconstrained case (depending on the specific choice of the control bounds), in non-convex domains the convergence rates are considerably larger. This is due to a smoothing effect of the control bounds on the continuous solution. For a more detailed discussion, we refer to the introduction of [4]. In the present paper, we only consider the case without control constraints. However, we notice that the estimates can be extend to the control constrained case as well. In the unconstrained case, if we use boundary concentrated meshes, we obtain a rate of two (up to logarithmic factors) as long as the interior angles are less than . Otherwise, we get the reduced rate . This is quite natural as we have already observed that the error for the discrete variational normal derivative is the limiting term within the error estimates.

Finally, we notice that there is an alternative approach to the -regularization. Several articles, for instance [16, 23, 24], consider a regularization in the -norm instead. This guarantees a higher regularity of the solution. In numerical experiments it turns out that the convergence rate in case of a standard discretization on quasi-uniform meshes seems to be one order higher than for the -regularization in case of quasi-uniform meshes, this is (up to logarithmic factors) if , and if . To the best of our knowledge, the corresponding estimates in the literature show a lower rate of convergence. This is mainly due to the fact that either standard techniques are used to bound the error for the discrete variational normal derivative or lower regularity of the data is assumed. A proof of the convergence rates stated above will be subject of a forthcoming article.

The paper is organized as follows: In Section 2, we introduce the variational formulation to the Poisson equation and establish regularity results in weighted Sobolev spaces, where the weight is a regularized distance function with respect to the boundary. Moreover, we collect several regularity results in different weighted Sobolev spaces from the literature for the later error analysis. The discretization of the Poisson equation and the boundary concentrated meshes are introduced at the beginning of Section 3. Moreover, the discretization error estimates in weighted -norms are proven in this section. These are applied in Section 4 in order to derive the error estimates for the two different approximations to the normal derivative of the solution of the Poisson equation. In addition, numerical experiments are included in this section which underline the theoretical findings. The numerical analysis for Dirichlet boundary control problems is outlined in Section 5. Moreover, numerical examples are presented which exactly show the convergence rates from the theory.

In the following will denote a generic constant which is always independent of the mesh parameter . We will use the notation to indicate that and .

## 2 Weighted regularity for elliptic problems

Let us first introduce some notation which is used in this paper. We consider computational domains that are bounded by a polygon . The corner points are numerated counter-clockwise and are denoted by , . The interior angle at a corner point is denoted by . The index set collects all indices with , i.e., the indices corresponding to non-convex corners. The boundary edge having endpoints and () is denoted by , . The classical Sobolev spaces are denoted as usual by for , , and by in case of . The corresponding norms are denoted by and , respectively. By we denote the completion of functions with respect to . Moreover, we use the notation and for the norm and the inner product in . An analogous notation is used for the spaces defined on the boundary.

For we consider the Poisson equation in variational form:

 Find u∈H10(Ω):(∇u,∇v)L2(Ω)=(f,v)L2(Ω)∀v∈H10(Ω). (1)

We introduce the regularized distance function

 σ(x):=dI+ρ(x)withρ(x):=dist(x,Γ):=infy∈Γ|x−y|, (2)

and some number exactly specified later. In the following we investigate regularity results in weighted spaces containing as weight function.

###### Lemma 2.1.

There exists a constant independent of such that

 (i) ∥σ−1u∥L2(Ω)+∥∇u∥L2(Ω) ≤c∥σf∥L2(Ω), (ii) ∥σ−1/2u∥L2(Ω)+|lndI|∥σ1/2∇u∥L2(Ω) ≤c|lndI|2∥σ3/2f∥L2(Ω),%if Ω is convex.
###### Proof.

As illustrated in Figure 2 we associate to each edge , , the subsets

 ΩΓi:={x∈Ω:ρ(x)=dist(x,Γi)},

and to each non-convex corner with the subsets

 Ωcj:={x∈Ω:ρ(x)=|x−cj|}

such that

 ¯Ω=(d⋃i=1¯ΩΓi)∪⎛⎜⎝⋃j∈Cnon¯Ωcj⎞⎟⎠.

For each set , we introduce the local coordinates with an affine linear map . Here, is a rotation matrix and a translation vector chosen in such a way that and . Moreover, we introduce the bounds and such that can be parameterized in local coordinates by

 ΩΓi={(x,y)⊤=Fi(xi,yi)∈R2:xi∈(0,¯¯¯xi), yi∈(0,¯¯¯yi(xi))}. (3)

The regularized distance function satisfies for

 (4)

Moreover, we write and confirm that

 ∇ΓiuΓi(xi,yi)=B⊤i∇u(Fi(xi,yi))with ∇Γi=(∂/∂xi,∂/∂yi)⊤.

To describe the sets , we instead use polar coordinates and located at the corner such that . Then, we find a representation of the form

 Ωcj={(x,y)⊤=(rjcosφj,rjsinφj)⊤∈R2:φj∈(π2,ωj−π2), rj∈(0,¯rj(φj))}

with an appropriate function . Within we write . There is the relation

 ∇cjucj(rj,φj)=(cosφjsinφj−rjsinφjrjcosφj)∇u(rjcosφj,rjsinφj)

with . Moreover, the weight function possesses for the representation

 σ(x,y)=dI+rj(x,y).

The result follows from the weak formulation of (1) and the Cauchy-Schwarz inequality:

 ∥∇u∥2L2(Ω)=(∇u,∇u)L2(Ω)=(f,u)L2(Ω)≤∥σf∥L2(Ω)∥σ−1u∥L2(Ω). (5)

Once we have shown , the result is proven. For that purpose, we consider the sub-domains and separately. Integration by parts using the local coordinates and the fact that implies together with the Cauchy-Schwarz inequality

 12∥σ−1u∥2L2(ΩΓi)=12∫¯xi0∫¯yi(xi)0u2Γi(xi,yi)(dI+yi)2dyidxi =−12∫¯xi0u2Γi(xi,yi)dI+yi∣∣ ∣∣¯yi(xi)yi=0dx+∫¯xi0∫¯yi(xi)0uΓi(xi,yi)∂yiuΓi(xi,yi)dI+yidyidxi ≤∥σ−1u∥L2(ΩΓi)∥∇ΓiuΓi∥L2(ΩΓi)=∥σ−1u∥L2(ΩΓi)∥∇u∥L2(ΩΓi).

In case of the sub-domains , we first use the property . In a second step we enlarge the domain to a circular sector with radius and containing . Afterwards, we use the fact that in combination with a Poincaré type inequality on the enlarged domain. This leads to

 ∥σ−1 u∥2L2(Ωcj)≤∫ωj−π2π2∫¯rj(φj)0r−1jucj(rj,φi)2drjdφj≤∫ωj0∫^rj0r−1jucj(rj,φi)2drjdφj ≤c∫^rj0∫ωj0r−1j(∂φjucj(rj,φj))2dφjdrj≤c∫^rj0∫ωj0rj|∇u(rjcosφj,rjsinφj)|2dφjdrj ≤c∥∇u∥2L2(Ω).

Summation over all subsets , , and , , yields the desired estimate

 ∥σ−1u∥L2(Ω)≤c∥∇u∥L2(Ω). (6)

To show the second estimate , we apply the Leibniz rule:

 ∥σ1/2∇u∥2L2(Ω)=∫Ωσ∇u⋅∇u=∫Ω∇u⋅∇(σu)−∫Ωu∇u⋅∇σ. (7)

The variational formulation (1) with leads to

 ∫Ω∇u⋅∇(σu)=(f,σu)L2(Ω)≤c∥σ3/2f∥L2(Ω)∥σ−1/2u∥L2(Ω). (8)

For the second term on the right-hand side of (7) we get from (4)

 ∫Ωu∇u⋅∇σ=d∑i=1∫ΩΓiuΓi(xi,yi)∇ΓiuΓi(xi,yi)⋅(01)dxidyi=d∑i=1∫ΩΓi12∂yiuΓi(xi,yi)2dxidyi. (9)

Integration by parts and exploiting the fact that vanishes on yields

 ≤∫¯¯¯xi0∫¯yi(xi)012∂yiuΓi(xi,yi)2dyidxi=12∫¯¯¯xi0u2Γi(xi,¯¯¯yi(xi))dxi ≥c∗2∫¯¯¯xi0u2Γi(xi,¯¯¯yi(xi))√1+¯¯¯y′i(xi)2dxi=c∗2∥u∥2L2(∂ΩΓi∖Γ).

The constant

 c∗:=mini∈Cessinfxi∈(0,¯¯¯xi)1/√1+¯¯¯y′i(xi)2

depends solely on the geometry of . Insertion into (9) yields

 ∫Ωu∇u⋅∇σ≥c∗∥u∥2L2(~Γ),~Γ:=d⋃i,j=1∂Ωi∩∂Ωj. (10)

Combining the estimates (7), (8) and (10) leads with Young’s inequality to

 ∥σ1/2∇u∥2L2(Ω)+c∗∥u∥2L2(~Γ)≤c|lndI|2∥σ3/2f∥2L2(Ω)+ε|lndI|−2∥σ−1/2u∥2L2(Ω). (11)

It remains to appropriately bound the latter term in (11) to show a weighted -estimate for . The decomposition into the subsets , integration by parts and Young’s inequality yield

 ∥σ−1/2u∥2L2(Ω) =d∑i=1∫¯¯¯xi0∫¯yi(xi)01dI+yiuΓi(xi,yi)2dyidxi =d∑i=1(∫¯¯¯xi0ln(dI+¯¯¯yi(xi))uΓi(xi,¯¯¯yi(xi))2dxi −∫¯¯¯xi0∫¯yi(xi)0ln(dI+yi)2uΓi(xi,yi)∂yiuΓi(xi,yi)dyidxi) ≤c∗∗∥u∥2L2(~Γ)+^c|lndI|2∥σ1/2∇u∥2L2(Ω)+12∥σ−1/2u∥2L2(Ω), (12)

with

 ^c>0andc∗∗:=2maxi∈Cesssupxi∈(0,¯¯¯xi)ln(1+¯yi(xi))√1+¯¯¯y′i(xi)2.

The latter term in (2) may be kicked back to the left-hand side. Inserting (2) into (11) and choosing

 ε=14min{c∗c∗∗,1^c}

yields

 ∥σ1/2∇u∥2L2(Ω)+c∗∥u∥2L2(~Γ) ≤c|lndI|2∥σ3/2f∥2L2(Ω)+12(∥σ1/2∇u∥2L2(Ω)+c∗|lndI|−2∥u∥2L2(~Γ)). (13)

Due to , a kick-back argument leads to the desired estimate for the term . Using the estimates (2) and (2) we finally confirm that

 ∥σ−1/2u∥2L2(Ω)≤c(∥u∥2L2(~Γ)+|lndI|2∥σ1/2∇u∥2L2(Ω))≤c|lndI|4∥σ3/2f∥2L2(Ω).

For technical reasons we decompose our domain into dyadic subsets

 ΩJ:={x∈Ω:ρ(x)∈(dJ+1,dJ)},J=0,…,I, (14)

where given , the numbers fulfill , for , and . Without loss of generality we assume that , otherwise a simple scaling argument can be used to achieve this. By means of the subsets , we will be able to handle the weight function within the proofs. More precisely, we will especially use that

 infx∈ΩJσ(x)∼supx∈ΩJσ(x)∼dJ,J=0,…,I. (15)

In the next lemma we show a local regularity result on the subsets and, as a consequence of this, global a priori estimates for second derivatives in weighted norms. Due to pollution effects we have to take into account the patches defined by

 Ω′J:=ΩJ−1∪ΩJ∪ΩJ+1,Ω′′J:=Ω′J−1∪Ω′J∪Ω′J+1

with the obvious modifications for the cases , and .

###### Lemma 2.2.

Let be a polygonal domain. The solution of (1) satisfies

 (i) ∥∇2u∥L2(ΩJ) ≤c(d−1J∥∇u∥L2(Ω′J)+∥f∥L2(Ω′J)),J=0,…,I, provided that u∈H2(Ω′J). Moreover, if Ω is convex, the estimates (ii) ∥σ∇2u∥L2(Ω) ≤c∥σf∥L2(Ω), (iii) ∥σ3/2∇2u∥L2(Ω) ≤c|lndI|∥σ3/2f∥L2(Ω)

are fulfilled.

###### Proof.

We consider a covering of consisting of finitely many balls , , with radius and centers . This implies . In a first step, we appropriately bound on . For this purpose, we introduce a smooth cut-off function satisfying in and . In case of the balls may overlap the boundary. Then, and are the intersections of each ball with . It is possible to construct in such a way that . Next, let be defined by

 ¯u:={|BdJ/4(xi)|−1∫BdJ/4(xi)uif suppη⊂Ω,0otherwise.

The function is the weak solution of the boundary value problem

 −Δ(η(u−¯u))=−Δη(u−¯u)−2∇u⋅∇η+ηfin Ω,η(u−¯u)=0on Γ.

Note that the right-hand side belongs to . With standard regularity results [13] we conclude

 ∥∇2u∥L2(BdJ/8(xi)) ≤c∥∇2(η(u−¯u))∥L2(Ω)≤c∥Δ(η(u−¯u))∥L2(Ω) ≤c(d−2J∥u−¯u∥L2(BdJ/4(xi))+d−1J∥∇u∥L2(BdJ/4(xi))+∥f∥L2(BdJ/4(xi))).

An application of the Poincaré inequality allows to bound the first term on the right-hand side by the second one. Note that a careful choice of the midpoints according to [20, Lemma A.1] guarantees that in each point only a finite number of balls overlap. The number depends only on the spatial dimension. Hence, we get the desired estimate (i),

 ∥∇2u∥2L2(ΩJ) ≤N∑i=1∥∇2u∥2L2(BdJ/8(xi))≤cN∑i=1(d−2J∥∇u∥2L2(BdJ/4(xi))+∥f∥2L2(BdJ/4(xi))) ≤c(d−2J∥∇u∥2L2(Ω′J)+∥f∥2L2(Ω′J)).

The estimate (ii) follows from (i) taking into account the relation (15). From this we deduce

 ∥σ∇2u∥2L2(Ω) ≤cI∑J=0d2J∥∇2u∥2L2(ΩJ)≤cI∑J=0(∥∇u∥2L2(Ω′J)+d2J∥f∥2L2(Ω′J)) ≤c(∥∇u∥2L2(Ω)+∥σf∥2L2(Ω)).

With Lemma 2.1 we conclude the assertion. In the same way the estimate (iii) follows. ∎

As we consider computational domains having a polygonal boundary, we have to deal with singularities occurring in the vicinity of vertices of the domain as well. For an accurate description of these singularities, we exploit regularity results in weighted Sobolev spaces with weights related to the corners. We denote the distance functions to the corners by , . Moreover, we introduce the regions , and choose appropriately such that these domains do not intersect. Furthermore, we introduce the region . On each we define for , and the local norm

 ∥v∥pVk,pβj(ΩjR):=∑|α|≤k∥rβj+|α|−kjDαv∥pLp(ΩjR), for p∈[1,∞), ∥v∥Vk,∞βj(ΩjR):=max|α|≤k∥rβj+|α|−kjDαv∥L∞(ΩjR).

The weighted Sobolev space with weight vector is defined as the set of measurable functions with finite norm

 ∥v∥Vk,p→β(Ω):=∥v∥Wk,p(^ΩR/2)+d∑j=1∥v∥Vk,pβj(ΩjR). (16)

We will frequently use these norms on subdomains . In this case, the weight functions are still related to the corners of and not of .

Under certain assumptions on the input data, one can show that the solution of (1) belongs to these weighted Sobolev space provided that the weights are sufficiently large. The lower bounds for the weights depend on the singular exponents

 λj:=π/ωj, j∈C.

The following result is taken from [21, §1.3, Theorem 3.1] for , and [18, Theorem 2.6.1] for .

###### Lemma 2.3.

Let with , and satisfying for all . Then, the solution of (1) belongs to and satisfies

 ∥u∥V2,p→β(Ω)≤c∥f∥V0,p→β(Ω).
###### Remark 2.4.

According to [17, Theorem 7.1.1] (see also [25, Lemma 2.32]) there holds

 d∑j=1⎛⎝∑|α|≤1∣∣(Dαv)(cj)∣∣⎞⎠=0

if with for . Thus, if with for , then the normal derivative is equal to zero at each convex corner. At non-convex corners it has a pole in general, see also the discussions in [3].

In order to derive optimal error estimates, we need a similar result for the case , which is excluded in the previous lemma. However, taking regularity results in weighted Hölder spaces into account (see [18]) the assertion of Lemma 2.3 remains true when assuming slightly more regularity for the right-hand side. For the proof of the following result we refer to [28, Lemma 4.2].

###### Lemma 2.5.

Assume that with some . Let the weight vector be chosen such that for all . Then, the solution of (1) belongs to