Solving ill-posed control problems by stabilized finite element methods: an alternative to Tikhonov regularization

This is the published version of a paper published in Inverse Problems.

Citation for the original published paper (version of record):

Burman, E., Hansbo, P., Larson, M G. (2018)

Solving ill-posed control problems by stabilized finite element methods: an alternative to Tikhonov regularization

Inverse Problems, 34(3): 035004

https://doi.org/10.1088/1361-6420/aaa32b

Access to the published version may require subscription.

N.B. When citing this work, cite the original published paper.

Permanent link to this version:

http://urn.kb.se/resolve?urn=urn:nbn:se:umu:diva-144929

PAPER • OPEN ACCESS

Solving ill-posed control problems by stabilized finite element methods:

an alternative to Tikhonov regularization

To cite this article: Erik Burman et al 2018 Inverse Problems 34 035004

View the article online for updates and enhancements.

Inverse Problems

Solving ill-posed control problems by stabilized finite element methods:

an alternative to Tikhonov regularization

Erik Burman

, Peter Hansbo

and Mats G Larson

1 Department of Mathematics, University College London, London, UK, WC1E 6BT, United Kingdom

2 Department of Mechanical Engineering, Jönköping University, SE-55111 Jönköping, Sweden

3 Department of Mathematics and Mathematical Statistics, Umeå University, SE, 901 87 Umeå, Sweden

E-mail: e.burman@ucl.ac.uk

Received 23 September 2016, revised 16 November 2017 Accepted for publication 20 December 2017

Published 31 January 2018 Abstract

Tikhonov regularization is one of the most commonly used methods for the regularization of ill-posed problems. In the setting of finite element solutions of elliptic partial differential control problems, Tikhonov regularization amounts to adding suitably weighted least squares terms of the control variable, or derivatives thereof, to the Lagrangian determining the optimality system. In this note we show that the stabilization methods for discretely ill- posed problems developed in the setting of convection-dominated convection – diffusion problems, can be highly suitable for stabilizing optimal control problems, and that Tikhonov regularization will lead to less accurate discrete solutions. We consider some inverse problems for Poisson ’s equation as an illustration and derive new error estimates both for the reconstruction of the solution from the measured data and reconstruction of the source term from the measured data. These estimates include both the effect of the discretization error and error in the measurements.

Keywords: optimal control problem, data assimilation, source identification, finite elements, regularization

(Some figures may appear in colour only in the online journal)

E Burman et al

Solving ill-posed control problems by stabilized finite element methods: an alternative to Tikhonov regularization

Printed in the UK 035004

INPEEY

© 2018 IOP Publishing Ltd 34

Inverse Problems

IP

1361-6420

10.1088/1361-6420/aaa32b

Paper

3

1

36

Inverse Problems

Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

2018

4 Author to whom any correspondence should be addressed.

Inverse Problems 34 (2018) 035004 (36pp) https://doi.org/10.1088/1361-6420/aaa32b

1. Introduction

In this note we propose an alternative to the classical Tikhonov regularization approach in the finite element approximations of optimal control problems governed by elliptic partial differ- ential equations. We shall, following [5], consider problems of the type

J(u, q) → min!, A(u) = f + B(q),

(1) where J is a cost functional, A is an elliptic differential operator for the state variable u, and B is an impact operator for the control variable q. Introducing the costate variable λ, this prob- lem can be formulated as finding saddle points for the Lagrangian functional

L( u, q, λ) := J(u, q) + (λ, A(u) − f − B(q)),

(2) where ( ·, ·) denotes the L

inner product, determined by the system

 



 

dd1

L(u + 

v, q, λ)|

= 0,

dd2

L(u, q + 

r, λ)|

= 0,

dd3

L(u, q, λ + 

µ) |

= 0.

(3)

In a finite element setting, the continuous states, controls and costates ( u, q, λ) ∈ V × Q × V are replaced by their discrete counterparts (u

, q

, λ

) ∈ V

× Q

× V

, where V

and Q

are finite dimensional counterparts of the appropriate Hilbert spaces V and Q, respectively.

Typically, the cost functional measures some distance between the discrete state and a known or sampled function u

over a subdomain _{M ⊆ Ω} , where Ω ⊂ R

, and d = 2,3 is the polyhedral (polygonal) domain of computation.

J(u, q) := 1

2 u − u

(4) which may not lead to a well-posed problem. A classical regularization method due to Tikhonov, see [30], is to add a stabilizing functional n(q, q) ,

J(u, q) := 1

2 u − u

+ n(q, q),

(5) where typically

n(q, q) := α

q − q

+ α

∇(q − q

) 

,

(6) where α

and α

are regularization parameters and q

is the background state, or first guess state. The role of the background state is to diminish the inconsistent character of the Tikhonov regularization and implies additional a priori knowledge of the system beyond the samples u

. In this note we will assume that no such additional a priori data is at hand, q

= 0, and that there is no physical justification for the addition of the term n(q, q) —or that the parameters ^α

and α

given by the application are too small to provide sufficient stabilization of the system for computational purposes. The actual choice of the parameters α

, α

typically depends on the level of noise in the data —see for example [ 32, 36] —and can be determined iteratively.

Work has also been done using the Morozov discrepancy principle to choose both regulariza-

tion parameters and discretization space [1, 3, 28]. More generally, in the computation the

mesh size is assumed to be small enough for the regularization to stabilize the numerical

scheme as well. The rule of thumb is then that the discretization error should be of the order

of the noise level [40, 43]. None of the above references consider analysis of the discretization error and perturbations in a unified setting.

It is well known that an inverse problem based on a linear operator with finite dimensional range is always well posed in the sense of Nashed [39], since the discrete forward operator has a closed range. If the linear system comes from the finite element discretization of an ill-posed PDE, then even if the system happens to be square with non-zero eigenvalues, it will typically be increasingly ill-conditioned as the dimension of the discretization space increases. This ill-conditioning is expressed as non-uniformity in the norms determining the stability proper- ties of the system: as the dimension increases, the stability bounds degenerate. This in its turn may have a detrimental effect on the accuracy of a numerical approximation. Indeed, for a computational PDE this phenomenon is not restricted to ill-posed problems, but arises typi- cally in situations where coercivity arguments fail, see for instance [21, section 5]. The aim of the present paper is to consider the discretizations of PDE problems that are ill-posed in a finite element framework, and provide regularization methods that allow for an analysis that considers both the discretization error and perturbation errors in a unified framework. More precisely, we are interested in problems that have some conditional stability property (see for instance [33, definition 4.3]). The key point is that the existence of a solution implies unique- ness and that a conditional stability estimate is available. The approach that we will follow is to eliminate the Tikhonov regularization on the continuous level and instead regularize the discrete formulation, hence making the regularization part of the computational method in the form of a weakly consistent stabilization. The regularization/stabilization parameter will then be linked to the mesh size and perturbations in the data through the error estimates. The ter- minology stabilization versus regularization is slightly ambiguous, but in classical numerical analysis the method of modified equations [27] provides a link between these concepts. The accuracy that is achievable depends on the regularity of the exact solution and the size of the perturbations in the data, but is optimal with respect to the order of the method and the stabil- ity of the problem. Note in particular that in this case the optimal parameter is independent of the stability of the quantity of interest. To make the comparison with the standard regulariza- tion technique easier we also add the analysis outline of such methods, relating the regulariza- tion parameter to the mesh-size. In contrast with our proposed method, the optimal choice of parameter for the Tikhonov regularization requires detailed knowledge of the stability proper- ties of the quantity of interest of the continuous problem. Even if this parameter is optimized with respect to the stability of the target quantity, the theoretical convergence order of the finite element method is always superior to that of Tikhonov regularization. On the other hand stabilizations such as those proposed herein change the matrix structure of the Galerkin method by increasing the bandwidth. This raises questions on how to solve the linear systems efficiently, which are more challenging compared to classical Tikhonov regularization.

Finally, let us point out that the use of conditional stability estimates has already been used for the choice of regularizing parameters [16, 34, 35], adaptive approaches for the choice of the parameter were explored in [26, 32], and iterative methods to find an optimal parameter proposed in [32]. The idea of using conditional stability estimates for analysis of the recon- struction method was used, for instance, in [17] for reconstructions using the Landweber itera- tion and in [29] in the context of variational regularization.

1.1. Model problems

We will discuss two model problems below, both based on a Poisson-type elliptic problem. In

the first, data is given in a subset of the bulk instead of on the boundary, and in the second the

solution is known in the domain and we wish to reconstruct the source term. The first problem can be associated with an ill-posed boundary control problem or a data assimilation problem, and the second is related to an ill-posed distributed control problem.

1.1.1. Model problem one: unique continuation.

Given a subset _{M ⊂ Ω} that is open, non- empty and distinct, the data u

∈ H

( M) and source term _{f ∈ L}

(Ω) , we wish to find u ∈ H

(Ω) such that

u = u

in M

(7) and u is a weak solution to

−∆u = f in Ω. (8)

This is not possible for all u

and f, but we will assume that the data is such that a solution exists. We will refer to this problem as the unique continuation problem. The ill-posedness of this problem follows by the following argument. Consider the operator

T : H

(Ω) → H

(Ω) × H

( M), Tu = (∆u, u|

).

The best approximate solution to (7) and (8) is given by the pseudoinverse T

( f , u

) when- ever this is well defined. It is well-known, see [20, proposition 2.4], that T

is continuous if and only if the range R(T) is closed. In this case T

is also well-defined on the whole

H

(Ω) × H

( M) . It makes sense to say that the minimization problem min Tu − ( f , u

) 

is well posed if and only if T

is continuous.

However, R(T) is not closed. Indeed, to get a contradiction suppose that this is the case.

We know that T is injective. Then T : H ^˜

(Ω) → R(T) , T = T ^˜ is bijective. As R(T) is closed, it is a Hilbert space, and T ^˜

: R(T) → H

(Ω) is continuous by the open mapping theorem.

But this implies that

u

= ˜T

Tu

 CTu

= C(∆u

+ u

).

Finally by choosing u as the solution of the classical Hadamard counterexample for the Cauchy problem we obtain a contradiction. We can for example take u

= sin(nx) sinh(ny) , so that ∆u

= 0 and assume that _M is centered around the origin; then C cannot be fixed independently of n.

It is known that in case the solution exists, a conditional stability estimate holds. The esti- mate can be quantified in a three-sphere inequality. Let _M include a ball B

(x

) with radius r

, centered at x

. We can then show the stability of the solution u in B

(x

) , with r

> r

under the a priori assumption that B

(x

) ⊂ B

(x

) ⊂ Ω and _{u ∈ H}

(Ω) is a weak solution to (8).

Lemma 1 (Three-sphere inequality). Assume that _{u : Ω → R} is a weak solution of (8) with _{f ∈ H}

(Ω) such that _f

 ε for some ε > 0 . For every r

, r

, r

, r such that 0 < r

< r

< r

< r and for every x

∈ Ω such that dist(x

, ∂Ω) > r there holds

u

 C 

u

+ ε 

· 

u

+ ε 

(9)

where B

, i = 1,2,3 are balls centered at x

, C > 0 and τ, 0 < τ < 1 only depend on the ge-

ometry of Ω, and the ratios r

/r

and r

/r

.

Proof. For a proof in the non-homogeneous case see Allessandrini et al [2, theorem 1.10].

□ Note that similar results hold in the H

-semi-norm —see for example [ 7] —and the results below can be extended to this case, under suitable modifications of the scheme. The sensitivity on the perturbations in the data is, however, expected to increase, since in this case the H

- semi-norm of the perturbations will come into play in the estimates.

We now cast the problem in the form of a constrained, regularized minimization problem:

find _{u ∈ H}

(Ω) minimizing 1

2 u − u

+ α

2 ∇ u

(10) subject to

−∆u = f , in Ω. (11)

Introducing the Lagrange multiplier λ ∈ H

(Ω) , we have the optimality system (see [30, 31]):

find (u

, λ

) ∈ H

(Ω) × H

(Ω) such that

M

u

v dΩ + 

Ω

∇λ

· ∇v dΩ + α

Ω

∇u

· ∇v dΩ =

M

u

v dΩ ∀v ∈ H

(Ω), (12)

Ω

∇u

· ∇µ dΩ =

Ω

f µ dΩ ∀µ ∈ H

(Ω).

(13) For α > 0 this system is well-posed by the Babuska –Lax–Milgram lemma [ 4], but even if u

is such that a solution to (11) exists with u = u

in _M , the solution will in general not be a solution to the optimality system (12) and (13). This is only satisfied by u = lim

u

. Here u denotes the function satisfying (11) and such that _u|

= u

. The regularized solution how- ever must fit the data to order O(α

) . This follows by testing equation (12) with u

− u and using that _u|

= u

leading to

2u

− u

+ α ∇u

+ α ∇(u

− u)

= α ∇u

. By (13) there holds ∆u

= ∆u independent of α so that _∆u

= ∆u

. Collecting the above bounds we deduce that for all α > 0

u

− u

+ α

∇u

 2α

∇u

and ∆u

= ∆u

. (14) It also follows from this relation that _∇u

is bounded in the limit. Observe that in the special case when _{u ∈ H}

(Ω) satisfies (11) with _u|

= u

, then u

= u for all α > 0 . This is easily verified by choosing u

= u and λ

= −αu in (12) and (13).

In the case where α = 0 we have the optimality system: find ( u, λ) ∈ H

(Ω) × H

(Ω) such that

M

u v dΩ + 

Ω

∇λ · ∇v dΩ =

M

u

v dΩ ∀v ∈ H

(Ω),

(15)

Ω

∇u · ∇µ dΩ =

Ω

f µ dΩ ∀µ ∈ H

(Ω).

(16)

Clearly the solution to (7) and (8) solves (15) and (16) with λ = 0 . But the unregularized minimization problem is ill-posed, similar to (7) and (8).

Below we will assume that u

∈ H

( M) is the unperturbed measurement for which the unique solution exists and consider a numerical method for the approximation of u given some perturbed data, ˜ u

:= u|

+ δu , with δ u ∈ L

( M) ; hence ˜ u

∈ L

( M) , which is why we minimize the L

-norm, _{u − u}

.

Remark 1 (Relation to boundary control problems). We can also consider the prob- lem of finding a function _{q : ∂Ω → R} minimizing

1

2 u(q) − u

+ α

2 ∇ u(q)

(17) subject to

−∆u = f , in Ω, u = q on ∂Ω.

(18) Clearly if q is found, which minimizes (17), the above discussed unique continuation problem is also solved and vice versa. Indeed, by solving the unique continuation problem, we find the optimal q by taking the trace of u. In the boundary control case it may be convenient to make the regularizing least squares term in (17) act directly on q. Finite element methods for this problem have been proposed in [15, 18, 24, 25, 38, 44], for instance, typically with a regular- izing term on the form _q

, but there appears to be no work in the case of the vanishing regularization, using the conditional stability of the limit problem, which is the topic of the present paper.

1.2. Model problem two: source reconstruction

Our second example considers the case where the data is available in the whole domain, M ≡ Ω , but the source term is unknown and must be reconstructed. The challenge here is that the application of the Laplacian is unstable. This case will be referred to as source reconstruc- tion below, but is also related to a distributed control problem. We consider the elementary problem: minimize

1

2 u − u

+ α

2  q

(19) subject to

−∆u = q in Ω, u = 0 on ∂Ω. (20)

Here, u

is known data and Ω is a convex polygonal (polyhedral) subset of _R

, d = 2,3, with outward-pointing normal n. We assume that we wish to solve (19) and (20) in the situation where u

is a measurement on a system that is of the form (20). This means that if no perturba- tions are present in the data and measurements are available in every point of Ω, the minimizer for α = 0 is u = u

∈ H

(Ω) ∩ H

(Ω) and an associated q = −∆u ∈ L

(Ω) exists so that (20) is satisfied. Also assume that the regularity bound

u

 Cq

(21)

holds. Below we will consider the problem of reconstructing q from u

, using a stabilized finite element method, and provide an analysis accounting for both the discretization error and the error due to errors in the measured data u

.

Introducing the Lagrange multiplier λ ∈ H

(Ω) , we have the optimality system (see [30, 31]): find ( u, q, λ) ∈ H

(Ω) × L

(Ω) × H

(Ω) such that

Ω

u v dΩ + 

Ω

∇λ · ∇v dΩ =

Ω

u

v dΩ ∀v ∈ H

(Ω),

(22)

α

Ω

q r dΩ + 

Ω

λ r dΩ = 0 ∀r ∈ L

(Ω),

(23)

Ω

∇u · ∇µ dΩ =

Ω

q µ dΩ ∀µ ∈ H

(Ω).

(24) We note here that the trace of λ is zero on the boundary and that this introduces an artificial boundary condition on q through the regularization term. Since q is equal to λ almost every- where we also observe that the regularization in L

actually imposes an H

-regularity on the source term. These artefacts carry over to the approximate solution when the regularized sys- tem is discretized. Once again we are interested here in the optimality system of (19) and (20) when α = 0 : find ( u, q, λ) ∈ H

(Ω) × L

(Ω) × H

(Ω) such that

Ω

u v dΩ + 

Ω

∇λ · ∇v dΩ =

Ω

u

v dΩ ∀v ∈ H

(Ω),

(25)

Ω

λ r dΩ = 0 ∀r ∈ L

(Ω),

(26)

Ω

∇u · ∇µ dΩ =

Ω

q µ dΩ ∀µ ∈ H

(Ω).

(27) Here we also assume that only the perturbed data ˜ u

= u|

+ δu with δ u ∈ L

( M) is avail- able for the reconstruction, and we will estimate quantitatively the effect of δu on the bounds below.

Remark 2 (Relation to distributed control problems). The related distributed con- trol problem is that of finding a function _{q : Ω → R} minimizing

1

2 u(q) − u

+ α

2  q

(28) subject to

−∆u = q, in Ω, u = 0 on ∂Ω.

(29) Here we consider the simplest context where the data u

is known over the whole domain Ω.

Finite element methods for the distributed control problems have been considered in [6, 22,

23, 37], for instance. Often some additional pointwise constraints on the control variable q are

introduced. To the best of our knowledge no works consider the situation of vanishing regu-

larization, which is at the interface between optimal control and inverse problems.

2. Discretization of ill-posed problems

Since the equations to which we wish to compute solutions are ill-posed, both the well-posed- ness of the forward and the adjoint problems may be compromised, and a naive discretization of the problem can not be expected to be successful. Indeed it is known, for instance, that the finite element equivalent of the unique extension of harmonic functions, which is at the basis of the stability of lemma 1, does not hold. To rectify the situation we draw on the experiences of the stabilization of well-posed, but numerically unstable problems, such as indefinite prob- lems, convection –diffusion equations or Stokes equations [ 8, 12, 13] and propose to stabilize the discrete ill-posed problems by penalizing the fluctuations of the solution. This results in a regularization that has no obvious interpretation on the continuous level, but ensures that the discrete system is invertible. The stability of the continuous problem may then be applied to the error of the computational approximation, resulting in error estimates.

2.1. Derivation of the discrete model

Let _{T

}

denote a family of shape regular and quasi-uniform tesselations of Ω in non-overlap- ping simplices, such that for any two different simplices K, K

∈ T

, _{K ∩ K}

consists of either the empty set, a common face/edge or a common vertex. The outwardly pointing normal of a simplex K will be denoted n

and the diameter h

. The global mesh parameter of _T

is defined by h := max

h

. We denote the set of interior element faces F in _T

by _F

. To each face we associate its diameter h

and a normal n

, whose orientation is arbitrary but fixed. We define the standard finite element space of continuous piecewise affine functions on _T

V

:= {v

∈ C

(Ω) : v

|

∈ P

( K), ∀K ∈ T

},

where P

(K) denotes the set of affine functions on K. We also define V

:= V

∩ H

(Ω) . The following bilinear forms defined on V

× V

will be useful for the formulation of our finite element methods,

m

(v

, w

) := 

X

v

w

dΩ, for X ⊆ Ω,

(30)

s

(v

, w

) := 

F∈FI

γ

F

h

[[ ∇v

· n

]][[ ∇w

· n

]] ds

(31) where [[y

]] |

:= lim

(y

( x − n

) − y

(x + n

)) denotes the jump of the quantity y

over the face F, with the normal n

, and γ

∈ R

denotes a dimensionless parameter independent of h. Finally

a(v

, w

) := 

Ω

∇v

· ∇w

dΩ.

(32)

2.2. Discretization of the unique continuation problem

If we write the minimization problem (10) and (11) on the discrete space V

, we may replace

the term _∇u

with the form s

(u

, u

) defined in (31). We then want to find u

∈ V

minimizing

1

2 u

− u

+ 1

2 s

(u

, u

)

(33) under the constraint

a(u

, v

) = m

( f , v

) , ∀v

∈ V

.

Observe that the constraint equation is underdetermined since the trial space is larger than the test space.

Remark 3. In the discrete setting we may take both the trial and the test space to be V

, but in this case, for consistency, we must add a boundary term to the definition (32) of a leading to the modified form

a(v

, w

) := 

Ω

∇v

· ∇w

dΩ − 

∂Ω

∇v

· nw

ds.

This term does not make sense on the continuous level, but can be included in the analysis below after minor modifications. For conciseness we leave the details to the reader.

Considering the optimality of the above discrete minimization problem we obtain the finite element formulation: find u

, λ

∈ V

× V

such that

m

(u

, v

) + s

(u

, v

) + a(v

, λ

) = m

(˜ u

, v

) ∀v

∈ V

,

(34) a(u

, µ

) = m

(˜ f , µ

) ∀µ

∈ V

,

(35) where ˜f := f + δf and ˜ u

:= u

+ δu

, with δ f ∈ H

(Ω) and δu

∈ L

( M) denoting mea- surement errors in the source term and data.

This may then be written on the compact form, find u

, λ

∈ V

, with _V

:= V

× V

such that

A

[(u

, λ

), (v

, µ

)] = m

(˜ f , µ

) + m

(˜ u

, v

) , ∀v

, µ

∈ V

, with

A

[(u

, λ

), (v

, µ

)] := m

(u

, v

) + s

(u

, v

) + a(v

, λ

) + a(u

, µ

).

Remark 4. An analogy can be made with the stabilized finite element discretizations for singularly perturbed convection –diffusion equations of the form

u + a · ∇u − ε∆u = f ,

(36) where _{a ∈ [W}

(Ω)]

and ε ∈ R

. When the P éclet number ( Pe(h) = (|a|h)/ε , where h is the mesh size), is large, the H

-stability of the scheme _ε

_∇u

 ε

 f 

obtained by the standard energy estimate is insufficient to ensure the stability of the discretization and spurious oscillations appear that compromise the accuracy. Reference [12] proposed adding a stabilizing term on the form s

( ·, ·) to the standard FEM formulation, which was shown to improve the estimates in the spirit of other known methods for this problem (see [41]). The key to the enhanced stability is that the addition of the stabilizing term makes the following stability estimate hold

u

+ h

a · ∇u

+ ε

∇u

 C f 

.

The stabilization improves the robustness in the limit ε → 0 , and in the limit _{h → 0} it limits the rate of blow-up of _{a · ∇u}

.

If we compare it with the unique continuation problem, the unsymmetric operator is the (discrete) Laplacian ∆

u

∈ V

defined by

m

(∆

u

, v

) = a(u

, v

) ∀v

∈ V

. Taking v

= h

∆

u

we obtain

h∆

u

= m

(∆

u

, h

∆

u

) = a(u

, h

∆

u

).

One may then easily show, using integration by parts and the techniques introduced below (anticipating the inequalities (41)), that

a(u

, h

∆

u

)  

F∈FI

F

h

|[[∇u

· n

]] |h

|h

∆

u

|  C|u

|

h∆

u

.

It follows that

h∆

u

 C|u

|

,

showing that the addition of the penalty to the gradient fluctuation enhances the stability on the discrete level, similar to the case of convection –diffusion equations. This estimate, how- ever, does not make sense for continuous equations.

2.3. Discretization of the source reconstruction problem

If we write the minimization problem (19) and (20) on the discrete space V

, we may replace the term _q

by the form s

(q

, q

) defined in (31). We then want to find u

∈ V

minimizing

1

2 u

− u

+ 1

2 s

(u

, u

) + 1

2 s

(q

, q

)

(37) under the constraint

a(u

, v

) = m

(q

, v

) , ∀v

∈ V

.

Observe that we also need to stabilize the solution u

. This is due to the large kernel of the operator s

, which makes the natural stability of the constraint equation insufficient. Below we will also discuss other more conventional regularizations in our framework. Considering the optimality of the above discrete minimization problem we obtain the finite element for- mulation: find u

, q

, λ

∈ V

× V

× V

such that

m

(u

, v

) + s

(u

, v

) + a(v

, λ

) = m

(˜ u

, v

) ∀v

∈ V

,

(38) where we have used the perturbed data ˜ u

= u

+ δu , where δ u ∈ L

(Ω) ,

m

(λ

, w

) − s

(q

, w

) = 0 ∀w

∈ V

(39) and

a(u

, µ

) = m

(q

, µ

) ∀µ

∈ V

(40)

with bilinear forms given by (30) –( 32) above. Below we will distinguish the stabilization parameters of s

( ·, ·) and s

( ·, ·) and denote them by γ

and γ

respectively.

This may then be written on the compact form: find u

, q

, λ

∈ V

with V

:= V

× V

× V

, such that

A

[(u

, q

, λ

), (v

, w

, µ

)] = m

(˜ u

, v

) , ∀v

, r

, µ

∈ V

, with

A

[(u

, q

, λ

), (v

, w

, µ

)] := m

(u

, v

) + s

(u

, v

)

+ a(v

, λ

) + s

(q

, w

) − m

(λ

, w

)

− a(u

, µ

) + m

(q

, µ

).

3. Preliminary technical results

First we define the semi-norms associated with the stabilization operator s

( ·, ·) ,

|x

|

= s

(x

, x

)

, ∀x

∈ V

.

We recall the following well-known inverse and trace inequalities (see for instance [19, section 1.4.3])

v

 C

(h

v

+ h

∇v

) , ∀v ∈ H

(K), h

∇v

· n

 C

∇v

, ∀v

∈ P

(K), h

∇v

+ h

v

 C

v

, ∀v

∈ P

(K).

(41)

As an immediate consequence of (41) we have the following stabilities for some C

> 0 depending only on the mesh geometry:

|x (42)

|

 C

h

x

, |x

|

 C

h

∇x

. This follows from the definition of s

( ·, ·) and the inverse inequalities (41)

|x

|

:= 

F∈FI

h

[[ ∇x

· n

]] 

 C 

K∈Th

h

∇x

 C 

K∈Th

h

x

. Let i

: H

(Ω) → V

denote the Scott –Zhang interpolant and ^π

: L

→ V

and π

: L

→ V

denote the L

-projections on the respective finite element spaces. The following error estimate is known to hold both for i

and π

:

u − i (43)

u

+ h∇(u − i

u)

 Ch

|u|

, t = 1, 2.

To prove the stability of our formulations below we need to show that for any function v

∈ V

the L

-norm is equivalent to _π

v

+ |v

|

. We prove the result in this technical lemma.

Lemma 2. There exists C

> 0 such that for all v

∈ V

there holds hv

 C

( v

+ |v

|

).

(44) There exists c

, c

> 0 such that for any function v

∈ V

c

( π

v

+ |v

|

)  v

 c

( π

v

+ |v

|

).

(45)

Proof. The discrete Poincar é type inequality ( 44) may be proved using a compactness ar- gument similar to that of [14]. To lessen the technical detail, we use here an approach with continuous Poincar é inequalities and discrete interpolation instead. Let I

: ∇V

→ [V

]

be a quasi-interpolation operator [11, section 5] such that

∇v (46)

− I

∇v

 C|v

|

, I

∇v

 C∇v

where ∆

:= ∪

. The following Poincar é inequality is well known (see [ 21, lemma B.63]). If f : H

(Ω) → R is a linear functional that is non-zero for constant functions then

u

 C

( | f (u)| + ∇u

) , ∀u ∈ H

(Ω).

For instance, we may take f (u) = 

M

u dΩ  Cu

.

As an immediate consequence we have the bound

v (47)

 C(v

+ ∇v

).

Now let _M

_{⊂ M} be the set of interior triangles of _M , M

:= {K ∈ T

: ∆

⊂ M}.

It then follows by the stability of _I

that _I

∇v

 C∇v

. Adding and sub- tracting _I

∇v in the second term on the right-hand side of (47) gives

v

 C(v

+ ∇v

− I

∇v

+ I

∇v

)

 C(v

+ |v

|

+ I

∇v

) (48) where we used (46) in the second inequality. For the third term on the right-hand side of (48) we once again use Poincar é’s inequality, the stability of I

and the inverse inequality (41) to conclude that

I

∇v

 C(I

∇v

+ ∇I

∇v

)

 C(h

v

+ ∇I

∇v

).

Using the fact that _∇v

is constant on each element we may write

∇I

∇v

= 

K∈Th

∇(∇v

− I

∇v

) 

.

Then using the inverse inequality (41) of each term in the sum, and the left relation of (46) it follows that

∇I

∇v

 Ch

∇v

− I

∇v

 Ch

|v

|

. Collecting these bounds we have shown that

I

∇v

 Ch

( v

+ |v

|

),