Convergence analysis of Galerkin finite element approximations to shape gradients in eigenvalue optimization1footnote 11footnote 1This work was supported in part by Science and Technology Commission of Shanghai Municipality (No. 18dz2271000) and the National Natural Science Foundation of China under grants 11201153 and 11571115.

# Convergence analysis of Galerkin finite element approximations to shape gradients in eigenvalue optimization111This work was supported in part by Science and Technology Commission of Shanghai Municipality (No. 18dz2271000) and the National Natural Science Foundation of China under grants 11201153 and 11571115.

Shengfeng Zhu Department of Data Mathematics & Shanghai Key Laboratory of Pure Mathematics and Mathematical Practice, School of Mathematical Sciences, East China Normal University, Shanghai 200241, China Xianliang Hu School of Mathematical Sciences, Zhejiang University, Hangzhou 310027, Zhejiang, China
###### Abstract

Numerical computation of shape gradients from Eulerian derivatives is essential to wildly used gradient type methods in shape optimization. Boundary type Eulerian derivatives are popularly used in literature. The volume type Eulerian derivatives hold more generally, but are rarely noticed and used numerically. We investigate thoroughly the accuracy of Galerkin finite element approximations of the two type shape gradients for optimization of elliptic eigenvalues. Under certain regularity assumptions on domains, we show a priori error estimates for the two approximate shape gradients. The convergence analysis shows that the volume integral formula converges faster and generally offers better accuracy. Numerical experiments verify theoretical results for the Dirichlet case. For the Neumann case, however, the boundary formulation surprisingly converges as fast as the volume one. Numerical results are presented.

###### keywords:
Shape optimization, shape gradient, eigenvalue problem, error estimate, finite element, multiple eigenvalue

## 1 Introduction

With the development of computer sciences, shape optimization has become important and promising in many fields of engineering (cf Defour (); HN (); HP (); Pironneau (); Sok1 ()), e.g., structural mechanics Bendsoe (), acoustics Osher2 (), computational fluid dynamics (see e.g., JHH (); Pironneau (); QR ()), etc. Analytic approaches can be applied to study existence and regularity BB (); Defour (). For shape design of complex systems in practice, numerical methods such as gradient-type algorithms are usually employed instead to seek “approximately” optimal shapes. One can adopt two strategies: optimize-then-discretize HN (); Sok1 () and discretize-then-optimize Bendsoe (). They are not equivalent at certain circumstances. For the former, the so-called Eulerian derivative and corresponding shape gradient of a shape functional are usually required by shape sensitivity analysis for sensitivity calculation of shape functionals with respect to domain variations (see e.g., Defour (); Sok1 ()). In 1907, Hadamard computed the Eulerian derivative of the first eigenvalue of a clamped plate with smooth boundary Hadamard (). Later, a structure theorem was developed by Zolésio for general shape functionals on -domains (). By the structure theorem, the Eulerian derivative can be expressed as a boundary integral. Due to the attractive concise appearance, this type Eulerian derivative has caused much attention both in shape optimization theory Defour (); Sok1 () and most existing numerical shape gradient algorithms (see e.g. Pironneau ()). But this type Eulerian derivative actually fails to hold when the boundary does not satisfy the required smoothness. Another more general type Eulerian derivative expressed as a domain integral then can be used instead Defour (). The two type Eulerian derivatives are equivalent through integration by parts if the boundary is regular enough.

For shape gradient computations, numerical approaches such as finite elements, finite differences (see e.g. ZhuJCP ()) and boundary element methods AF () are used to solve the state and possible adjoint constraints. Galerkin finite element method (cf. Brenner ()) is popular for discretizations of PDEs in shape optimization (see e.g. HN (); LiuZhu (); WHZ (); ZhuANM (); ZHW ()). This method based on domain triangulation is flexible to shape representation and shape changes in shape optimization.

The accuracy of finite element approximations of shape gradients seems to be essential for implementation of numerical optimization algorithms. Delfour and Zolésio (Remark 2.3 on pp. 531 Defour ()) pointed that “the boundary integral expression is not suitable since the finite element solution does not have the appropriate smoothness under which the boundary integral formula is obtained”. Pironneau et al. (pp 210 MP ()) presented convergence analysis for consistent approximations of boundary shape gradients in linear elliptic problems. Bergren Berggren () remarked that “the sensitivity information-directional derivatives of objective functions and constraints needs to be very accurately computed in order for the optimization algorithms to fully converge”. The use of domain expressions of Eulerian derivatives seems to be more promising. Recently, Hiptmair et al. Hiptmair () first showed that the Galerkin finite element approximation of shape gradient in the volume integral type converges faster and is more accurate than that in the boundary integral for linear elliptic problems. The volume formulations of shape gradients are derived and used for numerical shape optimization algorithms in magnetic induction tomography HLY () and parabolic diffusion problems SSW (), respectively. Shape gradients are popular in boundary form when combined with the level set method for shape optimization (see e.g. Allaire (); Burger (); LDZKL (); LiuZhu (); Osher2 ()). Recently, the volume type Eulerian derivatives were also incorporated numerically into the level set method for shape and topology optimization Burman (); Laurain (). In Kiniger (), an eigenvalue shape optimization problem was transformed to be an optimal control problem and a priori error estimates were obtained after finite element discretizations.

In this paper, we prove convergence for Galerkin finite element approximations of shape gradients in eigenvalue optimization. The motivations arise from the following aspects. First, eigenvalue problems in optimal shape design have fundamental importance for science and engineering, especially in structural mechanics (see. e.g., Allaire2 (); Bendsoe (); Defour (); HenriBook (); Rousselet83 (); Sok1 (); ZhuJCP ()). Second, finite element approximations to shape gradients of eigenvalues in boundary formulations are wildly used for numerical algorithms in eigenvalue optimization (see e.g., Allaire2 (); AF (); Defour (); HenriBook (); Osher2 (); Oudet (); Sok1 () and references therein). Third, numerical evidence shows that the potential advantages for using the new volume shape gradients in algorithms for shape optimization of Laplace eigenvalue problems ZhuJOTA (). To the best of our knowledge, there is no literature reported on convergence analysis of the approximate shape gradients in boundary/volume formulation for eigenvalue problems. For both Dirichlet and Neumann Laplace eigenvalue problems AF (); BB (); HenriBook (); KS (); LNS (); Oudet (), we prove convergence of the finite element approximations to shape gradients in both boundary and volume formulations. A priori error estimates are presented first in an infinite-dimensional operator norm. Numerical results are presented for verifying convergence of approximate shape gradients.

The rest of the paper is organized as follows. Laplace eigenvalue problems with Dirichlet and Neumann boundary conditions are presented in Section 2. Shape calculus is briefly introduced to give the Eulerian derivatives of eigenvalues in the forms of boundary and volume integrals. In Section 3, under certain regularity assumptions on domains, we present a priori convergence analysis of the finite element approximations to shape gradients in boundary and volume formulations. In Section 4, numerical results are presented for verifying convergence of approximate shape gradients as well as effectiveness of shape gradient algorithms in shape optimization. Brief conclusions are drawn in Section 5.

## 2 Problem formulation

Let be a bounded domain in with Lipschitz continuous boundary . We consider the Laplace eigenvalue problem:

 ⎧⎨⎩−Δu=λuin Ωu=0or∂u∂n=0on ∂Ω, (1)

where is the Laplacian. The homogeneous Dirichlet or Neumann boundary condition physically corresponds to a vibrating planar membrane () being fixed or free, respectively.

Let us first introduce briefly some notations on Sobolev spaces Adams (). For , the Banach space consists of measurable functions such that the associated norm . For each integer , the Banach space is equipped with the norm . In particular, when , we write instead of and instead of . Notice that is indeed a Hilbert space with respect to the scalar product with being the usual inner product. Denote by the closure of with respect to the norm . We write instead of when . For readability, we use the same notations for Sobolev norms of vector-valued and scalar functions.

The variational formulation of (1) is to find such that Babuska ()

 (∇u,∇v)=λ(u,v)∀v∈V, (2)

where () for the Dirichlet (Neumann) boundary condition. Due to the positiveness, self-adjointness and the compactness of the inverse of negative Laplacian operator, there exists a sequence of eigenpairs as solutions of (2) with eigenvalues

 0<λ1≤λ2≤⋯↗+∞ (3)

and corresponding eigenfunctions , which can be normalized with

 (ui,uj)=δij, (4)

where is the Kronecker delta. Note that for the case of Neumann boundary condition, . A typical optimization problem consists of minimizing some eigenvalue of (3) subject to a volume constraint AF (); BB (); HenriBook (); KS (); Oudet ().

### 2.1 A priori error estimates for Laplace eigenvalue problems

We consider the standard Ritz-Galerkin finite element method Brenner () for discretization and approximation of the variational formulation (2) Boffi (). For the shape gradient deformation algorithm we shall present, the domain here at each iteration is naturally assumed to be a polygon/polyhedron, which can be triangulated exactly with no geometric error introduced.

###### Remark 2.1.

For Dirichlet Laplace eigenvalue problems on planar polygonals (see Remark 4.2 Babuska ()), we have the following regularity result similar as linear elliptic problems Grisvard (). Let () be the maximal interior angle of the vertices of . Then, we have the for a Laplace eigenfunction such that

 u∈Hr(Ω)∩H10(Ω)with1

i.e., can be for any small . If is a convex polygon, then . For , more delicated discussions required to be made for the Poisson equation and we assume that () for simplicity.

Consider a family of triangulations satisfying that , where the mesh size with for any . Let be a family of finite-dimensional subspaces of . For the linear Lagrange elements, in the Dirichlet case with denoting the set of piecewise linear polynomials on and in the Neumann case. Denote with . Throughout, we shall denote by a general constant, which may differ at different occurrences and may depend on the mesh aspect ratio and the shape of , but it is always independent of eigenfunction and the mesh size . We assume that the mesh family is regular so that the following approximation property holds Brenner ():

 infvh∈Vh(∥u−vh∥L2(Ω)+h∥∇u−∇vh∥L2(Ω))≤Chr|u|Hr(Ω)∀u∈Hr(Ω). (5)

Suppose moreover that the mesh is quasi-uniform, i.e.,

 minK∈ThhK≥Ch∀h>0,

based on which the inverse inequality holds (see e.g. Theorem 4.5.11 Brenner ()). The weak formulation for conforming finite element approximation of the problem (2) reads: find and such that

 (∇uh,∇vh)=λh(uh,vh)∀vh∈Vh, (6)

For (6), there exist a finite sequence of eigenvalues

 0<λ1,h≤λ2,h≤⋯≤λN,h,N=dim Vh,

and corresponding eigenvectors

 u1,h,u2,h⋯uN,h,

which can be assumed to satisfy

 (ui,h,uj,h)=δij. (7)

For , we suppose that is the lowest index of the th distinct eigenvalue of (2) with being its multiplicity. More precisely, suppose that

 λki−1+li−1−1=λki−1<λki=λki+1=⋯=λki+li−1<λki+li=λki+1.

We have the following a priori error estimates on approximate eigenvalues and eigenfunctions.

###### Lemma 2.2.

Assume that is a polygon/polyhedron and are quasi-uniform. Let and be an eigenpair of (2) and (6), respectively with (). Then,

 λki+j−1≤λki+j−1,h≤λki+j−1+Ch2s|uki+j−1|H1+s(Ω),

and can be chosen so that (7) holds and

 ∥uki+j−1−uki+j−1,h∥L2(Ω)+h∥∇uki+j−1−∇uki+j−1,h∥L2(Ω)≤Ch1+s|uki+j−1|H1+s(Ω),
 ∥∇uki+j−1−∇uki+j−1,h∥H−s(Ω)≤Ch2s|uki+j−1|H1+s(Ω),

where with .

###### Proof.

By combining (3.25) KO () and Theorem 5.1 BO1989 (), we obtain the first two a priori error estimates above. To prove the last inequality, we first have

 ∥∇u−∇uh∥H−1(Ω) =sup0≠v∈H10(Ω)d⟨∇u−∇uh,v⟩∥v∥H1(Ω) (8) =sup0≠v∈H10(Ω)d(u−uh,divv)∥v∥H1(Ω) ≤∥u−uh∥L2(Ω)sup0≠v∈H10(Ω)d∥divv∥L2(Ω)∥v∥H1(Ω) ≤C∥u−uh∥L2(Ω) ≤Ch2|u|H2(Ω)for u∈H2(Ω),

where we have omitted the subscript for notational simplicity. On the other hand,

 ∥∇u−∇uh∥L2(Ω)≤C|u|H1(Ω)for u∈H1(Ω). (9)

Using the operator interpolation theorem (Lemma 22.3 Tartar ()) to (8) and (9), we have

 ∥∇u−∇uh∥H−s(Ω)≤Ch2s|u|H1+s(Ω)for u∈H1+s(Ω),0≤s≤1.

Note that this Lemma has included results for the special case of simple eigenvalues, i.e., . In following, we omit the index number for a specific eigenvalue and eigenfunction for simplicity. Let and be eigenpairs of (2) and (6), respectively.

Let us first define the Ritz projection such that

 (∇Phu,∇vh)=(∇u,∇vh)∀vh∈Vh. (10)
###### Lemma 2.3.

Let assumptions in Lemma 2.2 hold with . Then,

 ∥∇(Phu−uh)∥L2(Ω)≤Ch2|u|H2(Ω).
###### Proof.

We take in(2), (6) and (10). Then, we have

 (∇(Phu−uh),∇(Phu−uh)) =(λu−λhuh,Phu−uh) =(λ(u−uh)+(λ−λh)uh,Phu−uh).

Then, by the Cauchy-Schwarz inequality and triangle inequality,

 ∥∇(Phu−uh)∥2L2(Ω) ≤(λ∥u−uh∥L2(Ω)+|λ−λh|∥uh∥L2(Ω))∥Phu−uh∥L2(Ω) =(λ∥u−uh∥L2(Ω)+|λ−λh|)∥Phu−uh∥L2(Ω) ≤C(λ∥u−uh∥L2(Ω)+|λ−λh|)∥∇(Phu−uh)∥L2(Ω),

where the Poincaré inequality is used in the last inequality. Therefore,

 ∥∇(Phu−uh)∥L2(Ω) ≤C(λ∥u−uh∥L2(Ω)+|λ−λh|) (11) ≤C(λCh2+Ch2)|u|H2(Ω) ≤Ch2|u|H2(Ω)

using Lemma 2.2. ∎

###### Lemma 2.4.

Assume that . Then,

 ∥∇u−∇uh∥L4(Ω)≤Ch|u|W2,4(Ω). (12)
###### Proof.

By triangle inequality, we have

 ∥∇u−∇uh∥L4(Ω)≤∥∇u−∇Phu∥L4(Ω)+∥∇Phu−∇uh∥L4(Ω). (13)

By (8.5.4) on pp. 230 Brenner () and the approximation property (4.4.28) on pp. 110 Brenner (),

 ∥∇u−∇Phu∥L4(Ω) ≤Cinfv∈Vh∥∇u−∇v∥L4(Ω) (14) ≤Ch|u|W2,4(Ω),

By inverse inequality and Lemma 2.3, we have

 ∥∇Phu−∇uh∥L4(Ω) ≤Ch−d4∥∇Phu−∇uh∥L2(Ω) (15) ≤Ch2−d4|u|H2(Ω).

A combination of (13), (14) and (15) allows the conclusion to hold. ∎

### 2.2 Shape sensitivity analysis

As a tool in shape optimization, shape calculus/shape sensitivity analysis can be performed by the velocity (speed) method Defour (); Sok1 () and the perturbation of identity method Pironneau (). The two approaches are equivalent in the sense of first-order expansion with respect to domain perturbations. We recall basic shape calculus using the speed method (Sec. 2.9, pp. 54 and pp. 98 of Sok1 ()).

For a variable with , we introduce a velocity field with being the space of continuously differentiable transformations of . Then, we define a family of transformations with . For with , it satisfies the following flow system

 dxdt(t,X)=V(t,x(t,X)),  x(0,X)=X. (16)

For some domain , a shape functional depending on the shape is denoted by with . Denote in the following for simplicity.

###### Definition 2.5.

The Eulerian derivative of at in the direction is defined by

 dJ(Ω;V):=limt↘0J(Ωt)−J(Ω)t (17)

if the limit exists Defour ().

###### Definition 2.6.

The shape functional is called shape differentiable at if (i) there exist Eulerian derivatives for all directions ;
(ii) the map is linear and continuous from to .

We remark that non-differentiable cases occur when, e.g., the Eulerian derivative exists but the mapping is nonlinear. Such cases occur for the shape functionals of multiple eigenvalues.

###### Definition 2.7.

The material derivative in some Sobolev space of a state variable in a direction is denoted as

 ˙u(Ω;V):=limt↘0u(Ωt)∘Tt(V)(Ω)−u(Ω)t (18)

if the limit exists.

Material derivatives on the boundary can be defined analogously. When taking into account the strong (or weak) convergence in for the limit, the material derivative will be more specified as a strong (or weak) version.

The structure theorem [10, Corollary 1, pp. 480] states that the boundary Eulerian derivative of shape functional depends only on the normal part of the velocity on the boundary when certain smoothness of the boundary is satisfied. The volume formulation of Eulerian derivative actually holds with less smoothness requirement on boundary Laurain () and offers more accuracy Hiptmair (). By the speed method, we can obtain the Eulerian derivatives of simple as well as multiple eigenvalues for both Dirichlet and Neumann boundary conditions. Denote . The Eulerian derivative of an eigenvalue (depending on ) in the direction is defined to be

 dλ(Ω;V):=limt↘0λ(Ωt)−λ(Ω)t. (19)

For a simple eigenvalue, let be an eigenpair of the problem (2). Then, is shape differentiable and

 dλ(Ω;V)=∫Ω[−2∇u⋅DV∇u+divV(|∇u|2−λu2)]dx, (20)

where denotes the Jacobian of . If, furthermore, is convex or if it is of class , then the boundary Eulerian derivative of Dirichlet eigenvalue

 dλ(Ω;V)=−∫∂Ω(∂u∂n)2Vnds. (21)

If is of class for the Neumann case, then the boundary Eulerian derivative

 dλ(Ω;V)=∫∂Ω(|∇Γu|2−λu2)Vnds, (22)

 ∇Γu:=∇u−∂u∂nn.

(See in Appendix the formal derivations of (20)-(22) for self-containedness of the paper).

For sufficiently regular, the expression (20) corresponds to those appearing concise in HenriBook ():

 dλ(Ω;V)= −∫Ωdiv(|∇u|2V)dxforDirichlet (23) (or ∫Ωdiv((|∇u|2−λu2)V)dxfor  Neumann).

However, it is not appropriate to use (23) for discretization. The usual Lagrange finite element discretization of (23) fails to hold since is not continuously differentiable. We thus consider (20) for discretization.

For the multiple eigenvalue case, we simplify as and let be its eigenfunctions satisfying (4). Then, is no longer shape differentiable. Two strategies can be considered: the sub-differential and directional derivatives HenriBook (); Rousselet83 (); Sok1 (). We adopt the latter to follow closely the derivations of directional derivatives for the Dirichlet case as in Theorem 2.5.8 HenriBook () or Rousselet83 (). Then we can have the following results for both boundary conditions:

###### Proposition 2.8.

Assume that is a multiple eigenvalue of order for (1) with Dirichlet boundary condition. Let be an -orthonormal basis of the eigenspace associated with , then the Eulerian derivative is one of the eigenvalues of the matrix with the entry

 mi,j=∫Ω[−(DV+DVT)∇ui⋅∇uj+divV(∇ui⋅∇uj−λuiuj)]dx. (24)

If, furthermore, is convex or if it is of class , then

 mi,j=−∫∂Ω∂ui∂n∂uj∂nVnds. (25)

In the case of Neumann boundary condition,

 mi,j=∫Ω[−(DV+DVT)∇ui⋅∇uj+divV(∇ui⋅∇uj−λuiuj)]dx. (26)

If is of class , then

 mi,j=∫∂Ω(∇Γui⋅∇Γuj−λuiuj)Vnds. (27)

with .

## 3 A priori error estimates of approximate shape gradients in Eulerian derivatives

With the Galerkin finite element method for discretizations of the Laplace eigenvalue problem, we compute the approximate Eulerian derivatives and resulting shape gradients. We will analyze the convergence rates with a priori error estimates in an infinite-dimensional operator norm. For simplicity, we will only discuss the Dirichlet case. The results below however can be similarly extended to the Neumann case. We will consider both cases of simple and multiple eigenvalues. We first discuss the case of simple eigenvalues. In last section, is now simplified as . In order to differentiate notations for the boundary and volume type Eulearian derivatives, we denote (20) and (21) by and , respectively. The finite element approximations of (20) and (21) then read respectively:

 λ′(Ω;V)Ω,h:=∫Ω[−2∇uh⋅DV∇uh+divV(|∇uh|2−λhu2h)]dx (28)

and

 λ′(Ω;V)∂Ω,h:=−∫∂Ω(∂uh∂n)2Vnds. (29)

In the continuous setting, if is . With being discretized by finite elements, we have .

For the case of multiple eigenvalues, denote the matrices (resp. eigenvalues) (resp. ) and (resp. ) corresponding to (24) and (25), respectively. The approximations of (24) and (25) are

 mhi,j(Ω;V)Ω:=∫Ω[−(DV+DVT)∇ui,h⋅∇uj,h+divV(∇ui,h⋅∇uj,h−λui,huj,h)]dx (30)

and

 mhi,j(Ω;V)∂Ω:=−∫∂Ω∂ui,h∂n∂uj,h∂nVnds (31)

with . The corresponding matrices (resp. eigenvalues) are denoted by (resp. ) and (resp. ), respectively.

For each simple/multiple eigenvalue case, a priori error estimates are presented for two type (volume and boundary) finite element approximations of Eulerian derivatives and corresponding shape gradients. We first consider the case of simple eigenvalues and then the multiple case.

###### Remark 3.1.

In most cases for shape gradient algorithms, the domain is polyhedral and thus no geometric errors are introduced after triangulations. We will not consider the geometric approximation errors when performing convergence analysis of approximate Eulerian derivatives. When the domain is smooth, e.g., , for boundary formulas of Eulerian derivatives to hold, we assume that the geometric errors can be negligible by using isoparametric finite elements or fine meshes on boundaries.

### 3.1 Simple eigenvalue case

For the continuous formulas (20)-(21), we present convergence analysis of the approximate Eulerian derivatives with the volume integral (28) and boundary integral (29), respectively. For the volume type, we have

###### Theorem 3.2.

Let assumptions in Lemma 2.2 hold. Let be a single eigenpair of (2) and be its Galerkin Lagrange finite element approximation in (6). Then,

 |λ′(Ω;V)Ω−λ′(Ω;V)Ω,h|≤Ch2s|u|H1+s(Ω)|V|W1,∞(Ω),0

If and , we further have

 |λ′(Ω;V)Ω−λ′(Ω;V)Ω,h|≤Ch2|u|W2,4(Ω)|V|H2(Ω). (33)
###### Proof.

First, we have by (20), (28) and the triangle inequality

 |λ′(Ω;V)−λ′(Ω;V)Ω,h|≤∣∣∣∫Ω(2∇u⋅DV∇u−2∇uh⋅DV∇uh)dx∣∣∣ (34) +∣∣∣∫ΩdivV(|∇u|2−|∇uh|2)dx∣∣∣+∣∣∣∫ΩdivV(λ2u2−λ2hu2h)dx∣∣∣.

For the first term in R.H.S. of (34),

 ∣∣∣∫Ω(2∇u⋅DV∇u−2∇uh⋅DV∇uh)dx∣∣∣ (35) =∣∣∣∫Ω2(∇u−∇uh)⋅(DV+DVT)∇udx−∫Ω2(∇u−∇uh)⋅DV(∇u−∇uh)dx∣∣∣ ≤∣∣∣∫Ω2(∇u−∇uh)⋅(DV+DVT)∇udx∣∣∣+∣∣∣∫Ω2(∇u−∇uh)⋅DV(∇u−∇uh)dx∣∣∣ ≤2∥∇u−∇uh∥H−s(Ω)∥(DV+DVT)∇u∥Hs(Ω)+2∥DV∥L∞(Ω)∥∇u−∇uh∥2L2(Ω) ≤Ch2s|u|2H1+s(Ω)|V|W1,∞(Ω)+Ch2s|u|2H1+s(Ω)|V|W1,∞(Ω),

where Lemma 2.2 is used in the last inequality. For the second term in R.H.S. of (34), we have analogously

 ∣∣∣∫ΩdivV(|∇u|2−|∇uh|2)dx∣∣∣ (36) =∣∣∣2∫ΩdivV∇u⋅(∇u−∇uh)dx+∫ΩdivV|∇u−∇uh|2dx∣∣∣ ≤∣∣∣2∫ΩdivV∇u⋅(∇u−∇uh)dx∣∣∣+∣∣∣∫ΩdivV|∇u−∇uh|2dx∣∣∣ ≤2∥divV∥L∞(Ω)∥∇u∥Hs(Ω)∥∇u−∇uh∥H−s(Ω)+∥divV∥L∞(Ω)∥∇u−∇uh∥2L2(Ω) ≤Ch2s|u|2H1+s(Ω)|V|W1,∞(Ω)+Ch2s|u|2H1+s(Ω)|V|W1,∞(Ω).

For the third term in R.H.S. of (34), simple estimations yield

 ∣∣∣∫ΩdivV(λ2u2−λ2hu2h)dx∣∣∣ (37) =∣∣∣∫ΩdivV[2λu(λu−λhuh)−(λu−λhuh)2]dx∣∣∣ ≤∥divV∥L∞(Ω)(∫Ω|2λu(λu−λhuh)|dx+∫Ω|λu−λhuh|2dx) ≤C|V|W1,∞(Ω)(2λ∥u∥L2(Ω)∥λu−λhuh∥L2(Ω)+∥λu−λhuh∥2L2(Ω)),

where Cauchy-Schwarz inequality is used. In (37), and

 ∥λu−λhuh∥L2(Ω) =∥λ(u−uh)+(λ−λh)uh∥L2(Ω) (38) ≤λ∥u−uh∥L2(Ω)+|λ−λh|∥uh∥L2(Ω) ≤Cλh1+s|u|H1+s(Ω)+Ch2s|u|H1+s(Ω)

with and Lemma 2.2 being used. Therefore,

 ∣∣∣∫ΩdivV(λ2u2−λ2hu2h)dx∣∣∣ ≤C|V|W1,∞(Ω)(2λ(Cλh1+s+Ch2s)|u|H1+s(Ω) (39) +(Cλh1+s+Ch2s)2|u|2H1+s(Ω)) ≤Cλ2h2s|V|W1,∞(Ω)|u|H1+s(Ω).

Substituting (35), (36) and (39) into (34) allows (32) to hold.

Now we prove (33) when and . Each term on R.H.S. of (34) is estimated differently from above due to different regularities on and . By Hölder’s inequality, Lemma 2.2, Lemma 2.4, and Sobolev embedding theorem, we have

 ∣∣∣∫Ω(2∇u⋅DV∇u−2∇uh⋅DV∇uh)dx∣∣∣ (40) ≤2∥∇u−∇uh∥H−1(Ω)∥(DV+DVT)∇u∥H1(Ω) +2∥DV∥L4(Ω)∥∇u−∇uh∥L4(Ω)∥∇u−∇uh∥L2(Ω) ≤Ch2|u|H2(Ω)(|V|H2(Ω)|u|W1,∞(Ω)+|V|W1,4(Ω)|u|W2,4(Ω)) +Ch2|u|H2(Ω)|u|W2,4(Ω)|V|W1,4(Ω) ≤Ch2|u|H2(Ω)|V|H2(Ω)|u|W2,4(Ω)+Ch2|u|H2(Ω)|u|W2,4(Ω)|V|H2(Ω)</