Publicly Available Published by De Gruyter February 15, 2023

A Symmetric Interior Penalty Method for an Elliptic Distributed Optimal Control Problem with Pointwise State Constraints

Susanne C. Brenner , Joscha Gedicke and Li-Yeng Sung

From the journal Computational Methods in Applied Mathematics

https://doi.org/10.1515/cmam-2022-0148

Abstract

We construct a symmetric interior penalty method for an elliptic distributed optimal control problem with pointwise state constraints on general polygonal domains. The resulting discrete problems are quadratic programs with simple box constraints that can be solved efficiently by a primal-dual active set algorithm. Both theoretical analysis and corroborating numerical results are presented.

Keywords: Elliptic Distributed Optimal Control Problems; Pointwise State Constraints; Symmetric Interior Penalty Methods

MSC 2010: 49J40; 49M41; 65K15; 65N30

1 Introduction

Let Ω ⊂ R 2 be a bounded polygonal domain, y d ∈ L 2 ⁢ ( Ω ) , g ∈ H 4 ⁢ ( Ω ) , and let 𝛽 be a positive constant. The optimal control problem (cf. [31]) is to find

(1.1) ( y ¯ , u ¯ ) = argmin ( y , u ) ∈ K g 1 2 ⁢ [ ∥ y − y d ∥ L 2 ⁢ ( Ω ) 2 + β ⁢ ∥ u ∥ L 2 ⁢ ( Ω ) 2 ] ,

where ( y , u ) ∈ H 1 ⁢ ( Ω ) × L 2 ⁢ ( Ω ) belongs to K g if and only if

(1.2) ∫ Ω ∇ y ⋅ ∇ z ⁢ d ⁢ x = ∫ Ω u ⁢ z ⁢ d x for all ⁢ z ∈ H 0 1 ⁢ ( Ω ) ,

(1.3) y = g on ⁢ ∂ Ω ,

(1.4) y ≤ ψ a.e. in ⁢ Ω .

We assume that the function

(1.5) ψ ⁢ belongs to ⁢ W 3 , p ⁢ ( Ω ) ⁢ for ⁢ p > 2 ⁢ and ⁢ ψ > g ⁢ on ⁢ ∂ Ω .

Remark 1.1

Throughout this paper, we will follow the standard notation for differential operators, function spaces and norms that can be found for example in [1, 36, 23].

Observe that (1.2)–(1.3) is equivalent to y ∈ g + E ̊ ⁢ ( Δ ; L 2 ⁢ ( Ω ) ) , where

E ̊ ⁢ ( Δ ; L 2 ⁢ ( Ω ) ) = { z ∈ H 0 1 ⁢ ( Ω ) : Δ ⁢ z ∈ L 2 ⁢ ( Ω ) }

and Δ ⁢ z is understood in the sense of distributions. It is well known (cf. [46, 37]) that

(1.6) E ̊ ⁢ ( Δ ; L 2 ⁢ ( Ω ) ) ⊂ H 1 + α ⁢ ( Ω ) for some α ∈ ( 1 / 2 , 1 ]

and

(1.7) ∥ y ∥ H 1 + α ⁢ ( Ω ) ≤ C Ω ⁢ ∥ Δ ⁢ y ∥ L 2 ⁢ ( Ω ) for all ⁢ y ∈ E ̊ ⁢ ( Δ ; L 2 ⁢ ( Ω ) ) .

Hence functions in E ̊ ⁢ ( Δ ; L 2 ⁢ ( Ω ) ) are continuous by the Sobolev inequality (cf. [1]).

We can reformulate the optimal control problem (1.1)–(1.4) as the following minimization problem that only involves 𝑦: find

(1.8) y ¯ = argmin y ∈ K g 1 2 ⁢ [ ∥ y − y d ∥ L 2 ⁢ ( Ω ) 2 + β ⁢ ∥ Δ ⁢ y ∥ L 2 ⁢ ( Ω ) 2 ] ,

where

(1.9) K g = { y ∈ g + E ̊ ⁢ ( Δ ; L 2 ⁢ ( Ω ) ) : y ≤ ψ ⁢ in ⁢ Ω } .

Our goal is to solve the optimal control problem (1.1)–(1.4) by a symmetric interior penalty (SIP) method (cf. [61, 3]) that is based on the reformulation (1.8)–(1.9). This reformulation was discussed in [56], and the first numerical scheme based on this idea appeared in [53], where the analysis was carried out under certain ad hoc assumptions on the free boundary from [8]. These assumptions were later removed in the new convergence analysis in [24] by exploiting the regularity results in [43, 44, 30] for fourth-order elliptic variational inequalities. Various finite element methods based on this new approach have appeared in [28, 19, 21, 20, 29, 26, 22, 27].

Comparing with the more traditional approach in [52, 54, 49, 35, 32, 55] that is based on reducing the optimal control problem (1.1)–(1.4) to a problem that only involves the control, a distinct feature of the new approach is that the convergence of the state can also be established in the L ∞ norm. Another useful feature is that the discrete problems are quadratic programs with simple box constraints where the system matrices are available and consequently they can be solved efficiently by many optimization algorithms. Moreover, general polygonal/polyhedral domains can also be handled by the new approach (cf. [20, 26, 22]). These features are also shared by the SIP method in this paper. Compared to the P 1 finite element methods in [26, 22], the inverse of the block diagonal mass matrix of the SIP method of any order can be evaluated exactly, and hence the system matrix is available without mass lumping (cf. Remark 2.4 below).

The rest of the paper is organized as follows. We recall relevant results of the continuous problem and set up the discrete problem in Section 2. Technical tools for the error analysis are collected in Section 3, followed by an abstract error estimate in Section 4 and concrete error estimates in Section 5. Numerical results are presented in Section 6, and we end with some concluding remarks in Section 7. A heuristic justification of some of the numerical results is given in Appendix A.

Throughout the paper, we use 𝐶 (with or without subscripts) to denote a generic positive constant independent of the mesh size ℎ. We also use A ≲ B to represent the statement that A ≤ ( constant ) ⁢ B , where the positive constant is independent of ℎ.

2 Continuous and Discrete Problems

In this section, we recall relevant results of the continuous problem and construct the SIP method.

2.1 The Continuous Problem

It follows from the classical theory of calculus of variations (cf. [40, 51]) that (1.8)–(1.9) has a unique solution y ¯ ∈ g + E ̊ ⁢ ( Δ ; L 2 ⁢ ( Ω ) ) characterized by the variational inequality

(2.1) ∫ Ω ( y ¯ − y d ) ⁢ ( y − y ¯ ) ⁢ d x + β ⁢ ∫ Ω ( Δ ⁢ y ¯ ) ⁢ Δ ⁢ ( y − y ¯ ) ⁢ d x ≥ 0 for all ⁢ y ∈ g + E ̊ ⁢ ( Δ ; L 2 ⁢ ( Ω ) ) .

Remark 2.1

In the case where Ω is convex, the space E ̊ ⁢ ( Δ ; L 2 ⁢ ( Ω ) ) coincides with H 2 ⁢ ( Ω ) ∩ H 0 1 ⁢ ( Ω ) and (2.1) is a fourth-order variational inequality.

The variational inequality (2.1) is equivalent to the following generalized Karush–Kuhn–Tucker conditions:

(2.2) ∫ Ω ( y ¯ − y d ) ⁢ z ⁢ d x + β ⁢ ∫ Ω ( Δ ⁢ y ¯ ) ⁢ ( Δ ⁢ z ) ⁢ d x = ∫ Ω z ⁢ d μ for all ⁢ z ∈ E ̊ ⁢ ( Δ ; L 2 ⁢ ( Ω ) ) ,

where

(2.3) μ ⁢ is a nonpositive finite Borel measure

and

(2.4) ∫ Ω ( y ¯ − ψ ) ⁢ d μ = 0 .

Details for the derivation of (2.2)–(2.4) and for the regularity results below can be found in [24, 26] and the references therein.

Remark 2.2

It follows from the complementarity condition (2.4) that 𝜇 is supported on the active set 𝒜 for the constraint (1.4) given by A = { x ∈ Ω ¯ : y ¯ ⁢ ( x ) = ψ ⁢ ( x ) } . Note that 𝒜 is a compact subset of Ω by the assumption that ψ > g on ∂ Ω .

We have

(2.5) u ¯ = − Δ ⁢ y ¯ ∈ H 0 1 ⁢ ( Ω )

and

(2.6) y ¯ ∈ H loc 3 ⁢ ( Ω ) ∩ W loc 2 , ∞ ⁢ ( Ω ) ∩ H 1 + α ⁢ ( Ω ) ,

where α ∈ ( 1 / 2 , 1 ] is the index of elliptic regularity as in (1.6).

Remark 2.3

In the case where Ω is convex, we can replace 𝛼 in (2.6) with some α ~ ∈ ( 1 , 2 ] by exploiting (2.5).

Finally, the Lagrange multiplier 𝜇 belongs to H − 1 ⁢ ( Ω ) , i.e.,

(2.7) ∫ Ω z ⁢ d μ ≤ C ⁢ | z | H 1 ⁢ ( Ω ) for all ⁢ z ∈ H 1 ⁢ ( Ω ) .

2.2 The Discrete Problem

Let T h be a triangulation of Ω, let 𝑘 be a positive integer and

Y h = { y ∈ L 2 ⁢ ( Ω ) : y T = y | T ∈ P k ⁢ ( T ) ⁢ for all ⁢ T ∈ T h }

the discontinuous finite element space of degree at most 𝑘. The set of the edges of T h is denoted by E h , E h b ( ⊂ E h ) is the set of the boundary edges and | e | is the diameter of an edge 𝑒.

The discrete Laplace operator Δ h : ( g + E ̊ ⁢ ( Δ ; L 2 ⁢ ( Ω ) ) ) + Y h → Y h is defined by

(2.8) ∫ Ω ( Δ h ⁢ ζ ) ⁢ z h ⁢ d x = − a h ⁢ ( ζ , z h ) for all ⁢ z h ∈ Y h ,

where, following the convention for jumps and averages in [4],

a h ( y , z ) = ∑ T ∈ T h ∫ T ∇ y ⋅ ∇ z d x − ∑ e ∈ E h ∫ e ( { { ∇ y } } ⋅ ⟦ z ⟧ + { { ∇ z } } ⋅ ⟦ y ⟧ ) d s + ∑ e ∈ E h σ | e | ∫ e ⟦ y ⟧ ⋅ ⟦ z ⟧ d s

is the bilinear form for the SIP method with a sufficiently large penalty parameter 𝜎.

The SIP method is consistent (cf. [57]) in the sense that

(2.9) a h ⁢ ( ζ , z h ) = ∫ Ω ( − Δ ⁢ ζ ) ⁢ z h ⁢ d x + ∑ e ∈ E h b σ | e | ⁢ ∫ e g ⁢ z h ⁢ d s − ∑ e ∈ E h b ∫ e g ⁢ ( ∇ z h ⋅ n ) ⁢ d s

for all ζ ∈ g + E ̊ ⁢ ( Δ ; L 2 ⁢ ( Ω ) ) and z h ∈ Y h .

Let g h ∈ Y h be defined by

(2.10) ∫ Ω g h ⁢ z h ⁢ d x = ∑ e ∈ E h b ∫ e g ⁢ ( ∇ z h ⋅ n ) ⁢ d s − ∑ e ∈ E h b σ | e | ⁢ ∫ e g ⁢ z h ⁢ d s for all ⁢ z h ∈ Y h .

Then (2.8), (2.9) and (2.10) imply

(2.11) g h + Δ h ⁢ y = Q h ⁢ ( Δ ⁢ y ) for all ⁢ y ∈ g + E ̊ ⁢ ( Δ ; L 2 ⁢ ( Ω ) ) ,

where Q h is the orthogonal projection from L 2 ⁢ ( Ω ) onto Y h .

The discrete problem is to find

(2.12) y ¯ h = argmin y h ∈ K h 1 2 ⁢ [ ∥ y h − y d ∥ L 2 ⁢ ( Ω ) 2 + β ⁢ ∥ g h + Δ h ⁢ y h ∥ L 2 ⁢ ( Ω ) 2 ] ,

where

(2.13) K h = { y ∈ Y h : y T ⁢ ( p ) ≤ ψ ⁢ ( p ) ⁢ for all ⁢ p ∈ V T ⁢ and all ⁢ T ∈ T h } .

Here V T is the set of the three vertices of 𝑇.

The unique solution y ¯ h ∈ K h of (2.12)–(2.13) is characterized by the discrete variational inequality

(2.14) ∫ Ω ( y ¯ h − y d ) ⁢ ( y h − y ¯ h ) ⁢ d x + β ⁢ ∫ Ω ( g h + Δ h ⁢ y ¯ h ) ⁢ Δ h ⁢ ( y h − y ¯ h ) ⁢ d x ≥ 0 for all ⁢ y h ∈ K h .

Note that we can express the constraints in (2.13) concisely as

(2.15) I T ⁢ y T ≤ I T ⁢ ( ψ | T ) for all ⁢ T ∈ T h ,

where I T is the nodal interpolation operator for the Lagrange P 1 finite element on 𝑇.

Remark 2.4

The discrete problem is a quadratic program with simple box constraints. Let M h be the mass matrix that represents the L 2 inner product on Y h , and let A h be the stiffness matrix that represents the bilinear form a h ⁢ ( ⋅ , ⋅ ) on Y h . The system matrix for the quadratic program is then given by M h + β ⁢ A h ⁢ M h − 1 ⁢ A h , and M h − 1 is available because M h is a block diagonal matrix. Therefore, the discrete problem can be solved efficiently by the primal-dual active set algorithm in [9, 48].

For simplicity, we will focus on the analysis of the discrete problem for the case where g = 0 . The extension to the case of general 𝑔 is straightforward.

3 Some Technical Tools

We collect the results for some finite element tools in this section.

3.1 Mesh-Dependent Norms

Let 𝐷 be a subdomain of Ω. We will denote by T h ⁢ ( D ) the collection of elements in T h that have a nonempty intersection with 𝐷, i.e.,

(3.1) T h ⁢ ( D ) = { T ∈ T h : T ∩ D ≠ ∅ } .

(3.2) | | | z | | | h , D = ( ∑ T ∈ T h ⁢ ( D ) [ | z | H 1 ⁢ ( T ) 2 + ∑ e ⊂ ∂ T ( σ − 1 | e | ∥ { { ∇ z } } ⋅ n e ∥ L 2 ⁢ ( e ) 2 + σ | e | − 1 ∥ ⟦ z ⟧ ∥ L 2 ⁢ ( e ) 2 ) ] ) 1 2 ,

where n e is a unit vector normal to 𝑒, and | | | ⋅ | | | h , Ω is simply denoted by | | | ⋅ | | | h .

Note that

(3.3) a h ⁢ ( y , z ) ≤ | | | y | | | h ⁢ | | | z | | | h for all ⁢ y , z ∈ E ̊ ⁢ ( Δ ; L 2 ⁢ ( Ω ) ) + Y h

by the Cauchy–Schwarz inequality, and

(3.4) | | | y h | | | h 2 ≤ C ⁢ a h ⁢ ( y h , y h ) for all ⁢ y h ∈ Y h

provided that 𝜎 is sufficiently large (cf. [57]).

There is also a bound for | | | y | | | h .

Lemma 3.1

We have

(3.5) | | | y | | | h ≤ C ⁢ ( | y | H 1 ⁢ ( Ω ) + h α ⁢ ∥ Δ ⁢ y ∥ L 2 ⁢ ( Ω ) ) for all ⁢ y ∈ E ̊ ⁢ ( Δ ; L 2 ⁢ ( Ω ) ) .

Proof

Let y ∈ E ̊ ⁢ ( Δ ; L 2 ⁢ ( Ω ) ) . It follows from (3.2) that

(3.6) | | | y | | | h 2 = | y | H 1 ⁢ ( Ω ) 2 + ∑ T ∈ T h ∑ e ⊂ ∂ T σ − 1 ⁢ | e | ⁢ ∥ { { ∇ y } } ⋅ n e ∥ L 2 ⁢ ( e ) 2 ,

and we have a trace inequality with scaling,

(3.7) | e | ⁢ ∥ { { ∇ y } } ⋅ n e ∥ L 2 ⁢ ( e ) 2 ≤ C ⁢ ∑ T ∈ T e ( | y | H 1 ⁢ ( T ) 2 + h 2 ⁢ α ⁢ | y | H 1 + α ⁢ ( T ) 2 ) ,

where T e is the set of the triangles in T h that share 𝑒 as a common edge.

Estimate (3.5) follows from (1.7), (3.6) and (3.7). ∎

3.2 Triangulations

We will consider both quasi-uniform triangulations and triangulations that are graded around the reentrant corners (cf. [5, 6, 45, 2, 12]).

Let 𝐷 be a subdomain such that D ⋐ Ω , i.e., the closure of 𝐷 is a compact subset of Ω. Note that T h is quasi-uniform around 𝐷 for both types of triangulations. Therefore, we have

(3.8) | | | y h | | | h , D ≤ C ⁢ h − 1 ⁢ ∥ y h ∥ L 2 ⁢ ( Ω ) for all ⁢ y h ∈ Y h

by applying standard inverse estimates (cf. [36, 23]) to (3.2).

Moreover, we have

(3.9) ∥ y h ∥ L ∞ ⁢ ( Ω ) ≤ C ⁢ ( 1 + | ln ⁡ h | ) 1 2 ⁢ | | | y h | | | h for all ⁢ y h ∈ Y h

by the discrete Sobolev inequality in [17].

Remark 3.2

Estimate (3.9) was established in [17] for quasi-uniform triangulations. But the proof in [17] is also valid for meshes graded around the reentrant corners since the discrete Sobolev inequality for conforming Lagrange elements holds for such meshes (cf. [59, Lemma 6.4] and [23, Lemma 4.9.2]).

3.3 The Interpolation Operator Π h

Let Π T be the nodal interpolation operator for the P k Lagrange finite element on the triangle 𝑇. The following estimates (cf. [36, 23]) are standard:

(3.11) ∥ ζ − Π T ⁢ ζ ∥ L ∞ ⁢ ( T ) + h T ⁢ | ζ − Π T ⁢ ζ | W 1 , ∞ ⁢ ( T ) ≤ C ⁢ h T 2 ⁢ | ζ | W 2 , ∞ ⁢ ( T ) for all ⁢ ζ ∈ W 2 , ∞ ⁢ ( T ) .

Let V h = Y h ∩ H 0 1 ⁢ ( Ω ) be the conforming P k Lagrange finite element subspace of H 0 1 ⁢ ( Ω ) associated with T h . The operator Π h : E ̊ ⁢ ( Δ ; L 2 ⁢ ( Ω ) ) → V h is defined by

( Π h ⁢ ζ ) | T = Π T ⁢ ( ζ | T ) for all ⁢ T ∈ T h .

It follows from (3.10) (cf. [57, 18]) that

(3.12) ∥ ζ − Π h ⁢ ζ ∥ L 2 ⁢ ( Ω ) + h ⁢ ∥ ζ − Π h ⁢ ζ ∥ L ∞ ⁢ ( Ω ) + h ⁢ | | | ζ − Π h ⁢ ζ | | | h ≤ C ⁢ h 1 + τ ⁢ ∥ Δ ⁢ ζ ∥ L 2 ⁢ ( Ω )

for all ζ ∈ E ̊ ⁢ ( Δ ; L 2 ⁢ ( Ω ) ) , where

(3.13) τ = { α if ⁢ T h ⁢ is quasi-uniform , 1 if ⁢ T h ⁢ is graded around the reentrant corners ,

and 𝛼 is the index of elliptic regularity as in (1.6).

In the case where ζ ∈ H 2 ⁢ ( Ω ) ∩ H 0 1 ⁢ ( Ω ) , we have

(3.14) ∥ ζ − Π h ⁢ ζ ∥ L 2 ⁢ ( Ω ) + h ⁢ ∥ ζ − Π h ⁢ ζ ∥ L ∞ ⁢ ( Ω ) + h ⁢ | | | ζ − Π h ⁢ ζ | | | h ≤ C ⁢ h 2 ⁢ | ζ | H 2 ⁢ ( Ω )

for both quasi-uniform and graded meshes. In particular, for a subdomain D ⋐ Ω , interior elliptic regularity (cf. [42]) and (3.14) imply

(3.15) ∥ ζ − Π h ⁢ ζ ∥ L 2 ⁢ ( D ) + h ⁢ ∥ ζ − Π h ⁢ ζ ∥ L ∞ ⁢ ( D ) + h ⁢ | | | ζ − Π h ⁢ ζ | | | h , D ≤ C ⁢ h 2 ⁢ ∥ Δ ⁢ ζ ∥ L 2 ⁢ ( Ω )

for all ζ ∈ E ̊ ⁢ ( Δ ; L 2 ⁢ ( Ω ) ) .

The following result provides an estimate of Δ h ∘ Π h on smooth functions.

Lemma 3.3

Let 𝜙 be a C ∞ function with compact support in Ω. There exists a positive constant 𝐶 independent of ℎ such that ∥ Δ h ⁢ ( Π h ⁢ ϕ ) ∥ L 2 ⁢ ( Ω ) ≤ C ⁢ ∥ ϕ ∥ H 2 ⁢ ( Ω ) .

Proof

Let D ⋐ Ω be an open neighborhood of the support of 𝜙, and let z h ∈ Y h be arbitrary. It follows from (2.8), (3.3), (3.8) and (3.14) that

∫ Ω [ Δ h ⁢ ( ϕ − Π h ⁢ ϕ ) ] ⁢ z h ⁢ d ⁢ x = − a h ⁢ ( ϕ − Π h ⁢ ϕ , z h ) ≤ | | | ϕ − Π h ⁢ ϕ | | | h ⁢ | | | z h | | | h , D ≤ C ⁢ h ⁢ | ϕ | H 2 ⁢ ( Ω ) ⁢ | | | z h | | | h , D ≤ C ⁢ | ϕ | H 2 ⁢ ( Ω ) ⁢ ∥ z h ∥ L 2 ⁢ ( Ω ) ,

which implies ∥ Δ h ⁢ ( ϕ − Π h ⁢ ϕ ) ∥ L 2 ⁢ ( Ω ) ≤ C ⁢ | ϕ | H 2 ⁢ ( Ω ) . Consequently, we have, in view of (2.11) (with g = 0 = g h ),

∥ Δ h ⁢ ( Π h ⁢ ϕ ) ∥ L 2 ⁢ ( Ω ) ≤ ∥ Δ h ⁢ ( ϕ − Π h ⁢ ϕ ) ∥ L 2 ⁢ ( Ω ) + ∥ Δ h ⁢ ϕ ∥ L 2 ⁢ ( Ω ) ≤ C ⁢ ∥ ϕ ∥ H 2 ⁢ ( Ω ) . ∎

3.4 The Ritz Operator R h

The operator R h : E ̊ ⁢ ( Δ ; L 2 ⁢ ( Ω ) ) → Y h is defined by

(3.16) a h ⁢ ( R h ⁢ ζ , z h ) = a h ⁢ ( ζ , z h ) for all ⁢ z h ∈ Y h ,

and we have the following well-known error estimates for the SIP method (cf. [57, 18]):

(3.18) ∥ ζ − R h ⁢ ζ ∥ L 2 ⁢ ( Ω ) ≤ C ⁢ h 2 ⁢ τ ⁢ ∥ Δ ⁢ ζ ∥ L 2 ⁢ ( Ω ) for all ⁢ ζ ∈ E ̊ ⁢ ( Δ ; L 2 ⁢ ( Ω ) ) .

Note that (2.8), (2.11) (with g = 0 = g h ) and (3.16) imply

(3.19) Δ h ⁢ ( R h ⁢ ζ ) = Δ h ⁢ ζ = Q h ⁢ ( Δ ⁢ ζ ) for all ⁢ ζ ∈ E ̊ ⁢ ( Δ ; L 2 ⁢ ( Ω ) ) .

3.5 Interior Estimates

Let D 1 and D 2 be subdomains of Ω such that D 1 ⋐ D 2 ⋐ Ω . We have an interior energy norm error estimate for the SIP method (cf. [34, 25])

and also an interior maximum norm error estimate (cf. [25])

that is valid for all ζ ∈ E ̊ ⁢ ( Δ ; L 2 ⁢ ( Ω ) ) , where the norm ∥ ⋅ ∥ W h 1 , ∞ ⁢ ( D ) is given by

(3.22) ∥ w ∥ W h 1 , ∞ ⁢ ( D ) = max T ∈ T h ⁢ ( D ) [ ∥ ∇ w ∥ L ∞ ⁢ ( T ) + max e ⊂ ∂ T ( ∥ { { ∇ w } } ⋅ n e ∥ L ∞ ⁢ ( e ) + | e | − 1 ∥ ⟦ w ⟧ ∥ L ∞ ⁢ ( e ) ) ]

for all w ∈ E ̊ ⁢ ( Δ ; L 2 ⁢ ( Ω ) ) + Y h .

3.6 The Connection Operator C h

Recall V h = Y h ∩ H 0 1 ⁢ ( Ω ) is the conforming P k Lagrange finite element subspace of H 0 1 ⁢ ( Ω ) associated with T h . We can construct an operator C h : Y h → V h by averaging,

(3.23) ( C h ⁢ y h ) ⁢ ( p ) = 1 | T h ⁢ ( p ) | ⁢ ∑ T ∈ T h ⁢ ( p ) ( y h | T ) ⁢ ( p )

for any node 𝑝 of the P k Lagrange finite element space interior to Ω, where T h ⁢ ( p ) is the set of the triangles in T h that share the node 𝑝.

It follows from (2.13) and (3.23) that

(3.24) C h ⁢ v h = v h for all ⁢ v h ∈ V h ,

(3.25) C h ⁢ y h ∈ K h for all ⁢ y h ∈ K h ,

and for any subdomain 𝐷 of Ω, we have (cf. [13, 14, 16])

(3.26) h − 2 ∥ y h − C h y h ∥ L 2 ⁢ ( D ) 2 + ∑ T ∈ T h ⁢ ( D ) | y h − C h y h | H 1 ⁢ ( T ) 2 ≤ C ∑ T ∈ T h * ⁢ ( D ) ∑ e ⊂ ∂ T | e | − 1 ∥ ⟦ y h ⟧ ∥ L 2 ⁢ ( e ) 2 ,

where

(3.27) T ∈ T h belongs to T h * ⁢ ( D ) if and only if S T ∩ D ≠ ∅ .

Here S T (the star of 𝑇) is the union of all the triangles in T h that share a common vertex with 𝑇.

Moreover, it follows from (3.23) and a scaling argument that

(3.28) ∥ C h ⁢ y h ∥ L ∞ ⁢ ( T ) ≤ C ⁢ ∥ y h ∥ L ∞ ⁢ ( S T ) for all ⁢ T ∈ T h .

3.7 The Smoothing Operator E h

The operator E h : Y h → E ̊ ⁢ ( Δ ; L 2 ⁢ ( Ω ) ) is defined by

(3.29) Δ ⁢ ( E h ⁢ y h ) = Δ h ⁢ y h .

It follows from (2.9) (with g = 0 ), (3.16) and (3.29) that

(3.30) y h = R h ⁢ ( E h ⁢ y h ) .

Consequently, we have, by (3.17), (3.18) and (3.29),

(3.32) ∥ E h ⁢ y h − y h ∥ L 2 ⁢ ( Ω ) ≤ C ⁢ h 2 ⁢ τ ⁢ ∥ Δ h ⁢ y h ∥ L 2 ⁢ ( Ω ) for all ⁢ y h ∈ Y h .

Let G ⁢ ( A ) ⋐ Ω be an open neighborhood of the active set 𝒜. Putting (3.15), (3.18), (3.20) and (3.29) together, we find

Note that (3.33) also holds if G ⁢ ( A ) is replaced by a subdomain 𝐷 such that G ⁢ ( A ) ⋐ D ⋐ Ω . It then follows from (3.2) that

(3.34) ∑ T ∈ T h * ⁢ ( G ⁢ ( A ) ) ∑ e ⊂ ∂ T | e | − 1 ∥ ⟦ y h ⟧ ∥ L 2 ⁢ ( e ) 2 ≤ C h 2 ∥ Δ h y h ∥ L 2 ⁢ ( Ω ) 2 for all y h ∈ Y h .

3.8 The Operator I h

The operator I h : E ̊ ⁢ ( Δ ; L 2 ⁢ ( Ω ) ) → V h is defined by ( I h ⁢ ζ ) | T = I T ⁢ ( ζ | T ) for all T ∈ T h , where I T is the nodal interpolation operator for the Lagrange P 1 finite element on 𝑇. It follows from the case k = 1 in Section 3.3 that all the estimates for Π h in Section 3.3 can be applied to I h .

In addition, we have an obvious estimate

(3.35) ∥ I h ⁢ ζ ∥ L ∞ ⁢ ( T ) ≤ ∥ ζ ∥ L ∞ ⁢ ( T ) for all ⁢ ζ ⁢ continuous on ⁢ T ¯ ,

and, by standard inverse and interpolation estimates (cf. [36, 23]),

(3.36) | q − I h ⁢ q | H 1 ⁢ ( T ) ≤ C ⁢ h T − 1 ⁢ | q − I h ⁢ q | L 2 ⁢ ( T ) ≤ C ⁢ | q | H 1 ⁢ ( T ) for all ⁢ q ∈ P k ⁢ ( T ) ,

where the positive constant 𝐶 only depends on the shape regularity of 𝑇 and 𝑘.

3.9 Estimates for y ¯

It follows from (2.5) and (3.19) and a standard interpolation error estimate (cf. [36, 23]) that

(3.37) ∥ Δ h ⁢ ( R h ⁢ y ¯ ) − Δ ⁢ y ¯ ∥ L 2 ⁢ ( Ω ) = ∥ Q h ⁢ ( Δ ⁢ y ¯ ) − Δ ⁢ y ¯ ∥ L 2 ⁢ ( Ω ) ≤ C ⁢ h ⁢ | Δ ⁢ y ¯ | H 1 ⁢ ( Ω ) .

Let D ⋐ Ω . We have, by (2.6), (3.1), (3.11), (3.22) and (3.27),

(3.38) ∥ y ¯ − Π h y ¯ ∥ L ∞ ⁢ ( D ) ≤ C h 2 max T ∈ T h ⁢ ( D ) | y ¯ | W 2 , ∞ ⁢ ( T ) ,

(3.39) ∥ y ¯ − Π h y ¯ ∥ W h 1 , ∞ ⁢ ( D ) ≤ C h max T ∈ T h * ⁢ ( D ) | y ¯ | W 2 , ∞ ⁢ ( T ) .

Combining (3.12), (3.18), (3.21), (3.38) and (3.39), we find

(3.40) ∥ y ¯ − R h ⁢ y ¯ ∥ L ∞ ⁢ ( D ) ≤ C ⁢ ( ( 1 + | ln ⁡ h | ) ⁢ h 2 + h 2 ⁢ τ )

for any subdomain D ⋐ Ω .

The following lemma provides a simple global error estimate in the maximum norm.

Lemma 3.4

We have

(3.41) ∥ y ¯ − R h ⁢ y ¯ ∥ L ∞ ⁢ ( Ω ) ≤ C ⁢ ( 1 + | ln ⁡ h | ) 1 2 ⁢ h τ .

Proof

According to (3.2) and (3.9),

(3.42) ∥ y ¯ − R h ⁢ y ¯ ∥ L ∞ ⁢ ( Ω ) ≤ ∥ y ¯ − Π h ⁢ y ¯ ∥ L ∞ ⁢ ( Ω ) + ∥ Π h ⁢ y ¯ − R h ⁢ y ¯ ∥ L ∞ ⁢ ( Ω ) ≲ ∥ y ¯ − Π h ⁢ y ¯ ∥ L ∞ ⁢ ( Ω ) + ( 1 + | ln ⁡ h | ) 1 2 ⁢ | | | Π h ⁢ y ¯ − R h ⁢ y ¯ | | | h ,

which together with (3.12) and (3.17) implies (3.41). ∎

Remark 3.5

More sophisticated maximum norm error estimates for discontinuous Galerkin methods under stronger assumptions can be found in [34, 33, 47].

4 An Abstract Error Estimate

We will measure the error in terms of the mesh-dependent norm ∥ ⋅ ∥ h defined by

(4.1) ∥ y ∥ h 2 = ∥ y ∥ L 2 ⁢ ( Ω ) 2 + β ⁢ ∥ Δ h ⁢ y ∥ L 2 ⁢ ( Ω ) 2 .

From here on, we also use ( ⋅ , ⋅ ) to denote the inner product on L 2 ⁢ ( Ω ) .

Our goal is to establish the abstract error estimate in the following theorem, where y ¯ ∈ E ̊ ⁢ ( Δ ; L 2 ⁢ ( Ω ) ) is the solution of (2.1) (with g = 0 ) and y ¯ h is the solution of (2.14) (with g h = 0 ).

Remark 4.1

For simplicity, we have absorbed various norms of 𝑦 involved in the error analysis into the generic constant 𝐶 that appears in this section and Section 5.

Theorem 4.2

There exists a positive constant 𝐶 independent of ℎ such that

(4.2) ∥ y ¯ − y ¯ h ∥ h ≤ C ⁢ ( h + inf y h ∈ K h [ ∥ y ¯ − y h ∥ h + ∥ I h ⁢ ( y ¯ − C h ⁢ y h ) ∥ L ∞ ⁢ ( A ) 1 2 ] ) .

Proof

Let y h ∈ K be arbitrary. It follows from (2.14), (4.1) and the Cauchy–Schwarz inequality that

(4.3) ∥ y h − y ¯ h ∥ h 2 = ( y h − y ¯ , y h − y ¯ h ) + β ⁢ ( Δ h ⁢ ( y h − y ¯ ) , Δ h ⁢ ( y h − y ¯ h ) ) + ( y ¯ − y d , y h − y ¯ h ) + β ⁢ ( Δ h ⁢ y ¯ , Δ h ⁢ ( y h − y ¯ h ) ) − ( y ¯ h − y d , y h − y ¯ h ) − β ⁢ ( Δ ⁢ y ¯ h , Δ ⁢ ( y h − y ¯ h ) ) ≤ ∥ y h − y ¯ ∥ h ⁢ ∥ y h − y ¯ h ∥ h + ( y ¯ − y d , y h − y ¯ h ) + β ⁢ ( Δ h ⁢ y ¯ , Δ h ⁢ ( y h − y ¯ h ) ) .

Since E h ⁢ ( y h − y ¯ h ) ∈ E ̊ ⁢ ( Δ ; L 2 ⁢ ( Ω ) ) , we can use (2.2) and (3.29) to write

(4.4) ( y ¯ − y d , y h − y ¯ h ) + β ⁢ ( Δ h ⁢ y ¯ , Δ h ⁢ ( y h − y ¯ h ) ) = ( y ¯ − y d , ( y h − y ¯ h ) − E h ⁢ ( y h − y ¯ h ) ) + ( y ¯ − y d , E h ⁢ ( y h − y ¯ h ) ) + β ⁢ ( Δ ⁢ y ¯ , Δ ⁢ ( E h ⁢ ( y h − y ¯ h ) ) ) = ( y ¯ − y d , ( y h − y ¯ h ) − E h ⁢ ( y h − y ¯ h ) ) + ∫ Ω E h ⁢ ( y h − y ¯ h ) ⁢ d ⁢ μ ,

and we have, by (3.32),

(4.5) ( y ¯ − y d , ( y h − y ¯ h ) − E h ⁢ ( y h − y ¯ h ) ) ≤ C ⁢ h 2 ⁢ τ ⁢ ∥ Δ h ⁢ ( y h − y ¯ h ) ∥ L 2 ⁢ ( Ω ) .

We can rewrite the last term on the right-hand side of (4.4) as

(4.6) ∫ Ω E h ⁢ ( y h − y ¯ h ) ⁢ d μ = ∫ Ω [ E h ⁢ ( y h − y ¯ h ) − C h ⁢ ( y h − y ¯ h ) ] ⁢ d μ + ∫ Ω I h ⁢ ( ψ − C h ⁢ y ¯ h ) ⁢ d μ + ∫ Ω [ I h ⁢ C h ⁢ ( y ¯ h − y h ) − C h ⁢ ( y ¯ h − y h ) ] ⁢ d μ + ∫ Ω I h ⁢ ( y ¯ − ψ ) ⁢ d μ + ∫ Ω I h ⁢ ( C h ⁢ y h − y ¯ ) ⁢ d μ = T 1 + T 2 + T 3 + T 4 + T 5 .

The estimates for T 2 , T 4 and T 5 are straightforward. In view of (2.3), (2.15) and (3.25), we immediately have

(4.7) T 2 = ∫ Ω I h ⁢ ( ψ − C h ⁢ y ¯ h ) ⁢ d μ ≤ 0 .

Using (1.5), (2.3), (2.4), (2.6) and (3.11) (applied to I T ), we obtain

(4.8) T 4 = ∫ Ω [ I h ⁢ ( y ¯ − ψ ) − ( y ¯ − ψ ) ] ⁢ d μ ≤ ∥ I h ⁢ ( y ¯ − ψ ) − ( y ¯ − ψ ) ∥ L ∞ ⁢ ( A ) ⁢ | μ ⁢ ( Ω ) | ≤ C ⁢ h 2 ,

(4.9) T 5 = ∫ Ω I h ⁢ ( C h ⁢ y h − y ¯ ) ⁢ d μ ≤ ∥ I h ⁢ ( y ¯ − C h ⁢ y h ) ∥ L ∞ ⁢ ( A ) ⁢ | μ ⁢ ( Ω ) | .

Next we turn to T 1 . Let G ⁢ ( A ) ⋐ Ω be an open neighborhood of the active set 𝒜. We have, by (2.7), (3.1), (3.2) and (3.26),

and hence

(4.10) T 1 ≤ C ⁢ h ⁢ ∥ Δ h ⁢ ( y h − y ¯ h ) ∥ L 2 ⁢ ( Ω )

by (3.33) and (3.34).

Finally, we consider T 3 . Let w h = I h ⁢ E h ⁢ ( y ¯ h − y h ) . We have, by (2.7), (3.1), (3.2), (3.12) (applied to I h ), (3.26), (3.33), (3.34) and (3.36),

(4.11) T 3 = ∫ Ω [ I h ⁢ C h ⁢ ( y ¯ h − y h ) − C h ⁢ ( y ¯ h − y h ) ] ⁢ d μ ≲ | I h ⁢ C h ⁢ ( y ¯ h − y h ) − C h ⁢ ( y ¯ h − y h ) | H 1 ⁢ ( G ⁢ ( A ) ) = | I h ⁢ [ C h ⁢ ( y ¯ h − y h ) − w h ] − [ C h ⁢ ( y ¯ h − y h ) − w h ] | H 1 ⁢ ( G ⁢ ( A ) ) ≲ ( ∑ T ∈ T h ⁢ ( G ⁢ ( A ) ) | C h ⁢ ( y ¯ h − y h ) − w h | H 1 ⁢ ( T ) 2 ) 1 2 ≲ ( ∑ T ∈ T h ⁢ ( G ⁢ ( A ) ) | C h ⁢ ( y ¯ h − y h ) − ( y ¯ h − y h ) | H 1 ⁢ ( T ) 2 ) 1 2 + ( ∑ T ∈ T h ⁢ ( G ⁢ ( A ) ) | ( y ¯ h − y h ) − E h ⁢ ( y ¯ h − y h ) | H 1 ⁢ ( T ) 2 ) 1 2 + ( ∑ T ∈ T h ⁢ ( G ⁢ ( A ) ) | E h ⁢ ( y ¯ h − y h ) − I h ⁢ E h ⁢ ( y ¯ h − y h ) | H 1 ⁢ ( T ) 2 ) 1 2 ≲ h ⁢ ∥ Δ h ⁢ ( y h − y ¯ h ) ∥ L 2 ⁢ ( Ω ) .

Putting (4.6)–(4.11) together, we find

∫ Ω E h ⁢ ( y h − y ¯ h ) ⁢ d μ ≤ C ⁢ ( h 2 + h ⁢ ∥ Δ h ⁢ ( y h − y ¯ h ) ∥ L 2 ⁢ ( Ω ) + ∥ I h ⁢ ( y ¯ − C h ⁢ y h ) ∥ L ∞ ⁢ ( A ) ) ,

which together with the inequality of arithmetic and geometric means, (4.1) and (4.3)–(4.5) implies

(4.12) ∥ y h − y ¯ h ∥ h ≤ C ⁢ ( h + ∥ y h − y ¯ ∥ h + ∥ I h ⁢ ( y ¯ − C h ⁢ y h ) ∥ L ∞ ⁢ ( A ) 1 2 ) for all ⁢ y h ∈ K h .

Estimate (4.2) follows from (4.12) and the triangle inequality. ∎

An Improved Abstract Error Estimate

Estimate (4.2), which is established under assumption (1.5), implies that ∥ y ¯ − y ¯ h ∥ h is at most O ⁢ ( h ) in general. However, it can be improved under the following additional regularity assumptions:

(4.13) ψ ∈ H 4 ⁢ ( Ω ) ,

(4.14) A ⁢ has a smooth boundary ,

(4.15) y ¯ ⁢ belongs to ⁢ H 4 ⁢ ( D ∖ A ) ,

for any D ⋐ Ω that is an open neighborhood of 𝒜. Note that it follows from the interior elliptic regularity for the biharmonic operator and (2.4) that it suffices to assume (4.15) for just one such 𝐷.

Relations (2.2), (2.4) and integration by parts imply

(4.16) ∫ Ω z d μ = ∫ A ( y ¯ − y d ) z d x + β ( ∫ ∂ A ⟦ ∂ ( Δ y ¯ ) / ∂ n ⟧ z d s + ∫ A ( Δ 2 y ¯ ) z d x )

for all z ∈ E ̊ ⁢ ( Δ ; L 2 ⁢ ( Ω ) ) , where 𝑛 is the outer unit normal on ∂ A and ⟦ ∂ ( Δ ⁢ y ¯ ) / ∂ n ⟧ equals ∂ ( Δ ⁢ y ¯ ) / ∂ n from the outside of 𝒜 minus ∂ ( Δ ⁢ y ¯ ) / ∂ n from the inside of 𝒜. Consequently, we have, by (4.16) and the trace theorem,

(4.17) | ∫ Ω z ⁢ d μ | ≤ C ϵ ⁢ ∥ z ∥ H 1 2 + ϵ ⁢ ( G ⁢ ( A ) ) for all ⁢ z ∈ E ̊ ⁢ ( Δ ; L 2 ⁢ ( Ω ) ) ,

for any ϵ > 0 .

Moreover, assumptions (4.13)–(4.15) also imply the estimate

(4.18) | y ¯ ⁢ ( x ) − ψ ⁢ ( x ) | ≤ C ϵ ⁢ d 3 − ϵ ,

which holds for any ϵ > 0 and any x ∈ D whose distance to 𝒜 is less than or equal to 𝑑 (cf. [20, Lemma 5.5]).

For simplicity, we assume τ = 1 in (3.13) in the discussion below.

In view of (4.18), we can improve estimate (4.8) to

(4.19) T 4 ≤ C ϵ ⁢ h 3 − ϵ for all ϵ > 0

provided that we choose G ⁢ ( A ) so that

(4.20) dist ⁢ ( Ω ∖ G ⁢ ( A ) , A ) ≤ h .

For the term T 1 , observe that we have the estimate

(4.21) ∥ E h ⁢ ( y h − y ¯ h ) − C h ⁢ ( y h − y ¯ h ) ∥ L 2 ⁢ ( G ⁢ ( A ) ) ≤ ∥ E h ⁢ ( y h − y ¯ h ) − ( y h − y ¯ h ) ∥ L 2 ⁢ ( G ⁢ ( A ) ) + ∥ ( y h − y ¯ h ) − C h ⁢ ( y h − y ¯ h ) ∥ L 2 ⁢ ( G ⁢ ( A ) ) ≤ C ⁢ h 2 ⁢ ∥ Δ h ⁢ ( y h − y ¯ h ) ∥ L 2 ⁢ ( Ω )

by (3.26), (3.32) and (3.34). Consequently, it follows from (4.10), (4.17), (4.21) and interpolation between Sobolev spaces (cf. [58]) that

(4.22) T 1 = ∫ Ω [ E h ⁢ ( y h − y ¯ h ) − C h ⁢ ( y h − y ¯ h ) ] ⁢ d ⁢ μ ≤ C ϵ ⁢ ∥ E h ⁢ ( y h − y ¯ h ) − C h ⁢ ( y h − y ¯ h ) ∥ H 1 2 + ϵ ⁢ ( G ⁢ ( A ) ) ≤ C ϵ ⁢ h 3 2 − ϵ ⁢ ∥ Δ h ⁢ ( y h − y ¯ h ) ∥ L 2 ⁢ ( Ω ) .

Note that we can assume G ⁢ ( A ) to be a smooth domain so that the interpolation between Sobolev spaces on G ⁢ ( A ) can be handled without difficulty.

Similarly, using (4.11) (where w h = I h ⁢ E h ⁢ ( y ¯ h − y h ) ) and a standard interpolation error estimate, we find

∥ I h ⁢ C h ⁢ ( y ¯ h − y h ) − C h ⁢ ( y ¯ h − y h ) ∥ L 2 ⁢ ( G ⁢ ( A ) ) ≤ ( ∑ T ∈ T h ⁢ ( G ⁢ ( A ) ) ∥ I h ⁢ [ C h ⁢ ( y ¯ h − y h ) − w h ] − [ C h ⁢ ( y ¯ h − y h ) − w h ] ∥ L 2 ⁢ ( T ) 2 ) 1 2 ≤ C ⁢ h ⁢ ( ∑ T ∈ T h ⁢ ( G ⁢ ( A ) ) | C h ⁢ ( y ¯ h − y h ) − w h | H 1 ⁢ ( T ) 2 ) 1 2 ≤ C ⁢ h 2 ⁢ ∥ Δ h ⁢ ( y h − y ¯ h ) ∥ L 2 ⁢ ( Ω ) ,

which together with (4.11), (4.17) and interpolation between Sobolev spaces yields

(4.23) T 3 = ∫ Ω [ I h ⁢ C h ⁢ ( y ¯ h − y h ) − C h ⁢ ( y ¯ h − y h ) ] ⁢ d ⁢ μ ≤ C ϵ ⁢ ∥ I h ⁢ C h ⁢ ( y ¯ h − y h ) − C h ⁢ ( y ¯ h − y h ) ∥ H 1 2 + ϵ ⁢ ( G ⁢ ( A ) ) ≤ C ϵ ⁢ h 3 2 − ϵ ⁢ ∥ Δ h ⁢ ( y h − y ¯ h ) ∥ L 2 ⁢ ( Ω ) .

Putting (4.1), (4.3), (4.4), (4.5) (with τ = 1 ), (4.6), (4.7), (4.9), (4.19), (4.22) and (4.23) together, we arrive at the improved abstract error estimate

(4.24) ∥ y ¯ − y ¯ h ∥ h ≤ C ϵ ⁢ ( h 3 2 − ϵ + inf y h ∈ K h [ ∥ y ¯ − y h ∥ h + ∥ I h ⁢ ( y ¯ − C h ⁢ y h ) ∥ L ∞ ⁢ ( A ) 1 2 ] ) .

5 Concrete Error Estimates

The key to derive concrete error estimates is to bound the infimum in (4.2) by constructing a function y h * ∈ K h with the desired properties.

Lemma 5.1

There exists y h * ∈ K h for ℎ sufficiently small such that

(5.1) ∥ y ¯ − y h * ∥ h + ∥ I h ⁢ ( y ¯ − C h ⁢ y h * ) ∥ L ∞ ⁢ ( A ) 1 2 ≤ C ⁢ ( ( 1 + | ln ⁡ h | ) 1 2 ⁢ h + h τ ) .

Proof

Let G ⁢ ( A ) ⋐ Ω be an open neighborhood of the active set 𝒜. According to (3.40), we have

(5.2) ϵ h = ∥ R h ⁢ y ¯ − y ¯ ∥ L ∞ ⁢ ( G ⁢ ( A ) ) ≤ C ⁢ ( ( 1 + | ln ⁡ h | ) ⁢ h 2 + h 2 ⁢ τ ) .

Let 𝜙 be a nonnegative C ∞ function with compact support in Ω such that ϕ = 1 on G ⁢ ( A ) , and let y h * ∈ Y h be defined by

(5.3) y h * = R h ⁢ y ¯ − ϵ h ⁢ Π h ⁢ ϕ .

First we show that y h * belongs to K h for sufficiently small ℎ. Indeed, we have ψ − y ¯ ≥ δ > 0 on Ω ∖ G ⁢ ( A ) and hence

y h * ≤ R h ⁢ y ¯ = y ¯ + ( R h ⁢ y ¯ − y ¯ ) ≤ ψ − δ + ( R h ⁢ y ¯ − y ¯ ) on Ω ∖ G ⁢ ( A ) ,

which, in view of (3.41), implies

y h * ⁢ ( p ) ≤ ψ ⁢ ( p ) for any vertex p of T h that belongs to Ω ∖ G ⁢ ( A )

provided that ℎ is sufficiently small. On the other hand, for any vertex 𝑝 of T h that belongs to G ⁢ ( A ) , we have, by (5.2) and (5.3),

y h * ⁢ ( p ) = ( R h ⁢ y ¯ ) ⁢ ( p ) − ϵ h = y ¯ ⁢ ( p ) + ( R h ⁢ y ¯ − y ¯ ) ⁢ ( p ) − ϵ h ≤ y ¯ ⁢ ( p ) ≤ ψ ⁢ ( p ) .

Next we estimate the two terms that appear on the left-hand side of (5.1). According to Lemma 3.3, (3.18), (3.19), (5.2) and (5.3), we have

(5.4) ∥ y ¯ − y h * ∥ L 2 ⁢ ( Ω ) ≤ ∥ y ¯ − R h ⁢ y ¯ ∥ L 2 ⁢ ( Ω ) + ϵ h ⁢ ∥ Π h ⁢ ϕ ∥ L 2 ⁢ ( Ω ) ≤ C ⁢ ϵ h ,

(5.5) ∥ Δ h ⁢ ( y ¯ − y h * ) ∥ L 2 ⁢ ( Ω ) = ϵ h ⁢ ∥ Δ h ⁢ ( Π h ⁢ ϕ ) ∥ L 2 ⁢ ( Ω ) ≤ C ⁢ ϵ h ,

and therefore, in view of (4.1) and (5.2),

(5.6) ∥ y ¯ − y h * ∥ h ≤ C ⁢ ( ( 1 + | ln ⁡ h | ) ⁢ h 2 + h 2 ⁢ τ ) .

For the second term, we find, by (3.24), (3.27), (3.28) and (3.35),

(5.7) ∥ I h ⁢ ( y ¯ − C h ⁢ y h * ) ∥ L ∞ ⁢ ( A ) ≤ max T ∈ T h ⁢ ( G ⁢ ( A ) ) ∥ y ¯ − C h y h * ∥ L ∞ ⁢ ( T ) ≤ max T ∈ T h ⁢ ( G ⁢ ( A ) ) ⁡ ( ∥ y ¯ − Π h ⁢ y ¯ ∥ L ∞ ⁢ ( T ) + ∥ C h ⁢ ( Π h ⁢ y ¯ − R h ⁢ y ¯ ) ∥ L ∞ ⁢ ( T ) + ϵ h ⁢ ∥ Π h ⁢ ϕ ∥ L ∞ ⁢ ( T ) ) ≲ max T ∈ T h * ⁢ ( G ⁢ ( A ) ) ⁡ ( ∥ y ¯ − Π h ⁢ y ¯ ∥ L ∞ ⁢ ( T ) + ∥ Π h ⁢ y ¯ − R h ⁢ y ¯ ∥ L ∞ ⁢ ( T ) ) + ϵ h ≲ max T ∈ T h * ⁢ ( G ⁢ ( A ) ) ⁡ ( ∥ y ¯ − Π h ⁢ y ¯ ∥ L ∞ ⁢ ( T ) + ∥ y ¯ − R h ⁢ y ¯ ∥ L ∞ ⁢ ( T ) ) + ϵ h ,

and hence, in view of (3.38), (3.40) and (5.2),

(5.8) ∥ I h ⁢ ( y ¯ − C h ⁢ y h * ) ∥ L ∞ ⁢ ( A ) 1 2 ≤ C ⁢ ( ( 1 + | ln ⁡ h | ) 1 2 ⁢ h + h τ ) .

Estimate (5.1) follows from (5.6)–(5.8). ∎

We can now establish several concrete error estimates.

Theorem 5.2

Let y ¯ h ∈ K h be the solution of (2.12)–(2.13) and u ¯ h = − Δ h ⁢ y ¯ h . There exists a positive constant 𝐶 independent of ℎ such that

(5.9) ∥ u ¯ − u ¯ h ∥ L 2 ⁢ ( Ω ) + ∥ y ¯ − y ¯ h ∥ L 2 ⁢ ( Ω ) ≤ C ⁢ ( ( 1 + | ln ⁡ h | ) 1 2 ⁢ h + h τ ) .

Proof

It follows from (4.1), Theorem 4.2 and Lemma 5.1 that

(5.10) ∥ y ¯ − y ¯ h ∥ L 2 ⁢ ( Ω ) + ∥ Δ h ⁢ ( y ¯ − y ¯ h ) ∥ L 2 ⁢ ( Ω ) ≤ C ⁢ ( ( 1 + | ln ⁡ h | ) 1 2 ⁢ h + h τ ) ,

which together with (2.11) (where g = 0 = g h ) and (3.37) yields

∥ u ¯ − u ¯ h ∥ L 2 ⁢ ( Ω ) ≤ ∥ Δ ⁢ y ¯ − Δ h ⁢ y ¯ ∥ L 2 ⁢ ( Ω ) + ∥ Δ h ⁢ ( y ¯ − y ¯ h ) ∥ L 2 ⁢ ( Ω )

(5.11) = ∥ Δ ⁢ y ¯ − Q h ⁢ ( Δ ⁢ y ¯ ) ∥ L 2 ⁢ ( Ω ) + ∥ Δ h ⁢ ( y ¯ − y ¯ h ) ∥ L 2 ⁢ ( Ω )

≤ C ⁢ ( ( 1 + | ln ⁡ h | ) 1 2 ⁢ h + h τ ) . ∎

The following lemmas will enable us to establish estimates for y ¯ − y ¯ h in other norms.

Lemma 5.3

There exists a positive constant 𝐶 independent of ℎ such that

| | | y h | | | h ≤ C ⁢ ∥ Δ h ⁢ y h ∥ L 2 ⁢ ( Ω ) for all ⁢ y h ∈ Y h .

Proof

In view of (3.3), (3.4), (3.16), (3.30), we have

| | | y h | | | h 2 ≲ a h ⁢ ( y h , y h ) = a h ⁢ ( E h ⁢ y h , y h ) ≤ | | | E h ⁢ y h | | | h ⁢ | | | y h | | | h

and hence

| | | y h | | | h ≲ | | | E h ⁢ y h | | | h ≲ ∥ Δ ⁢ ( E h ⁢ y h ) ∥ L 2 ⁢ ( Ω ) = ∥ Δ h ⁢ y h ∥ L 2 ⁢ ( Ω )

by (1.7), Lemma 3.1 and (3.29). ∎

Lemma 5.4

There exists a positive constant 𝐶 independent of ℎ such that

(5.12) ∥ y h ∥ L ∞ ⁢ ( Ω ) ≤ C ⁢ ∥ Δ h ⁢ y h ∥ L 2 ⁢ ( Ω ) for all ⁢ y h ∈ Y h .

Proof

It follows from (3.9), (3.12), (3.29) and (3.31) that

(5.13) ∥ y h − E h ⁢ y h ∥ L ∞ ⁢ ( Ω ) ≤ ∥ y h − Π h ⁢ ( E h ⁢ y h ) ∥ L ∞ ⁢ ( Ω ) + ∥ Π h ⁢ ( E h ⁢ y h ) − E h ⁢ y h ∥ L ∞ ⁢ ( Ω ) ≲ ( 1 + | ln ⁡ h | ) 1 2 ⁢ | | | y h − Π h ⁢ ( E h ⁢ y h ) | | | h + h τ ⁢ ∥ Δ ⁢ ( E h ⁢ y h ) ∥ L 2 ⁢ ( Ω ) ≲ ( 1 + | ln ⁡ h | ) 1 2 ⁢ ( | | | y h − E h ⁢ y h | | | h + | | | E h ⁢ y h − Π h ⁢ ( E h ⁢ y h ) | | | h ) + h τ ⁢ ∥ Δ ⁢ ( E h ⁢ y h ) ∥ L 2 ⁢ ( Ω ) ≲ ( 1 + | ln ⁡ h | ) 1 2 ⁢ h τ ⁢ ∥ Δ ⁢ ( E h ⁢ y h ) ∥ L 2 ⁢ ( Ω ) = ( 1 + | ln ⁡ h | ) 1 2 ⁢ h τ ⁢ ∥ Δ h ⁢ y h ∥ L 2 ⁢ ( Ω ) .

Moreover, we have, by (1.7), (3.29) and the Sobolev inequality,

(5.14) ∥ E h ⁢ y h ∥ L ∞ ⁢ ( Ω ) ≲ ∥ E h ⁢ y h ∥ H 1 + α ⁢ ( Ω ) ≲ ∥ Δ h ⁢ ( E h ⁢ y h ) ∥ L 2 ⁢ ( Ω ) = ∥ Δ h ⁢ y h ∥ L 2 ⁢ ( Ω ) .

Estimate (5.12) follows from (5.13) and (5.14). ∎

Theorem 5.5

There exists a positive constant 𝐶 independent of ℎ such that

(5.15) | | | y ¯ − y ¯ h | | | h ≤ C ⁢ ( ( 1 + | ln ⁡ h | ) 1 2 ⁢ h + h τ ) ,

(5.16) ∥ y ¯ − y ¯ h ∥ L ∞ ⁢ ( Ω ) ≤ C ⁢ ( ( 1 + | ln ⁡ h | ) 1 2 ⁢ h + h τ ) .

Proof

In view of (3.19), (5.10) and Lemma 5.3, we have

| | | R h ⁢ y ¯ − y ¯ h | | | h ≲ ∥ Δ h ⁢ ( R h ⁢ y ¯ − y ¯ h ) ∥ L 2 ⁢ ( Ω ) = ∥ Δ h ⁢ ( y ¯ − y ¯ h ) ∥ L 2 ⁢ ( Ω ) ≲ ( 1 + | ln ⁡ h | ) 1 2 ⁢ h + h τ ,

and therefore

| | | y ¯ − y ¯ h | | | h ≤ | | | y ¯ − R h ⁢ y ¯ | | | h + | | | R h ⁢ y ¯ − y ¯ h | | | h ≲ ( 1 + | ln ⁡ h | ) 1 2 ⁢ h + h τ

by (3.17).

According to (3.19), (5.10) and Lemma 5.4, we have

∥ R h ⁢ y ¯ − y ¯ h ∥ L ∞ ⁢ ( Ω ) ≲ ∥ Δ h ⁢ ( R h ⁢ y ¯ − y ¯ h ) ∥ L 2 ⁢ ( Ω ) = ∥ Δ h ⁢ ( y ¯ − y ¯ h ) ∥ L 2 ⁢ ( Ω ) ≲ ( 1 + | ln ⁡ h | ) 1 2 ⁢ h + h τ ,

which together with (3.41) and the triangle inequality implies (5.16). ∎

Improved Error Estimates

We can improve estimate (5.1) under additional regularity assumptions. For the discussion below, we assume (4.13), (4.14) and strengthen (4.15) to

(5.17) y ¯ ∈ H 4 ⁢ ( Ω ∖ A ) .

Combining y ¯ ∈ H loc 3 ⁢ ( Ω ) with (4.13), (4.14) and (5.17), we see that y ¯ belongs to C 2 ⁢ ( Ω ¯ ) and

(5.18) y ¯ ∈ H 7 2 − ϵ ⁢ ( Ω ) for all ϵ > 0 .

For simplicity, we assume τ = 1 in (3.13) and k = 2 in the discussion below.

It follows from (5.18) that

(5.19) | | | y ¯ − Π h ⁢ y ¯ | | | h ≤ C ⁢ h 2 ⁢ | y ¯ | H 3 ⁢ ( Ω ) ,

and estimates (3.17) and (3.18) can then be improved to

(5.20) ∥ y ¯ − R h ⁢ y ¯ ∥ L 2 ⁢ ( Ω ) + h ⁢ | | | y ¯ − R h ⁢ y ¯ | | | h ≤ C ⁢ h 3 ⁢ | y ¯ | H 3 ⁢ ( Ω ) .

Assumption (4.13) and the Sobolev inequality imply that ψ ∈ W 3 , p ⁢ ( Ω ) for any p < ∞ and hence

(5.21) ∥ ψ − Π T ⁢ ψ ∥ L ∞ ⁢ ( T ) ≤ C ϵ ⁢ h T 3 − ϵ for all ⁢ T ∈ T h , ϵ > 0

by the Bramble–Hilbert lemma (cf. [11, 39]) and scaling. We can combine estimate (4.18) with (5.21) to obtain

(5.22) ∥ y ¯ − Π h ⁢ y ¯ ∥ L ∞ ⁢ ( G ⁢ ( A ) ) ≤ ∥ ( y ¯ − ψ ) − Π h ⁢ ( y ¯ − ψ ) ∥ L ∞ ⁢ ( G ⁢ ( A ) ) + ∥ ψ − Π h ⁢ ψ ∥ L ∞ ⁢ ( G ⁢ ( A ) ) ≤ C ϵ ⁢ h 3 − ϵ

provided that G ⁢ ( A ) satisfies (4.20). Similarly, the estimate

(5.23) ∥ y ¯ − Π h ⁢ y ¯ ∥ W h 1 , ∞ ⁢ ( G ⁢ ( A ) ) ≤ C ϵ ⁢ h 2 − ϵ

follows from (3.22), (4.18) and (5.21). Combining (3.21) and (5.20)–(5.23), we have the estimate

(5.24) ϵ h = ∥ R h ⁢ y ¯ − y ¯ ∥ L ∞ ⁢ ( G ⁢ ( A ) ) ≤ C ϵ ⁢ h 3 − ϵ

that improves (5.2). Consequently, estimate (5.8) is improved to

(5.25) ∥ I h ⁢ ( y ¯ − C h ⁢ y h * ) ∥ L ∞ ⁢ ( A ) 1 2 ≤ C 3 ⁢ h 3 2 − ϵ

by (5.7), (5.22) and (5.24). Putting (5.4), (5.5), (5.24) and (5.25) together, we obtain the estimate

(5.26) ∥ y ¯ − y h * ∥ h + ∥ I h ⁢ ( y ¯ − C h ⁢ y h * ) ∥ L ∞ ⁢ ( A ) 1 2 ≤ C ϵ ⁢ h 3 2 − ϵ

that improves (5.1). Hence we have

∥ y ¯ − y ¯ h ∥ h ≤ C ϵ ⁢ h 3 2 − ϵ

by (4.24) and (5.26), and estimate (5.10) becomes

∥ y ¯ − y ¯ h ∥ L 2 ⁢ ( Ω ) + ∥ Δ h ⁢ ( y ¯ − y ¯ h ) ∥ L 2 ⁢ ( Ω ) ≤ C ϵ ⁢ h 3 2 − ϵ ,

and thus, in view of (5.18), estimate (5.11) becomes

∥ u ¯ − u ¯ h ∥ L 2 ⁢ ( Ω ) ≤ ∥ Δ ⁢ y ¯ − Δ h ⁢ y ¯ ∥ L 2 ⁢ ( Ω ) + ∥ Δ h ⁢ ( y ¯ − y ¯ h ) ∥ L 2 ⁢ ( Ω ) ≤ ∥ Δ ⁢ y ¯ − Q h ⁢ ( Δ ⁢ y ¯ ) ∥ L 2 ⁢ ( Ω ) + ∥ Δ h ⁢ ( y ¯ − y ¯ h ) ∥ L 2 ⁢ ( Ω ) ≤ C ϵ ⁢ h 3 2 − ϵ .

Therefore, we have the estimate

(5.27) ∥ u ¯ − u ¯ h ∥ L 2 ⁢ ( Ω ) + ∥ y ¯ − y ¯ h ∥ L 2 ⁢ ( Ω ) ≤ C ϵ ⁢ h 3 2 − ϵ

that improves (5.9). Similarly estimates (5.15) and (5.16) can be improved to

| | | y ¯ − y ¯ h | | | h ≤ C ϵ ⁢ h 3 2 − ϵ and ∥ y ¯ − y ¯ h ∥ L ∞ ⁢ ( Ω ) ≤ C ϵ ⁢ h 3 2 − ϵ

through (3.42), (5.19) and (5.20).

6 Numerical Results

We have performed numerical experiments on the SIP method for k = 1 and k = 2 , where the penalty parameter σ = 3 ⁢ k ⁢ ( k + 1 ) (cf. [41]). The mesh parameter for the 𝑗-th level refinement is h j = 8 / 2 j . The discrete problems are solved by a primal-dual active set algorithm (cf. [7, 9, 48, 50]).

We report the errors of the optimal state in ∥ ⋅ ∥ L 2 ⁢ ( Ω ) and | | | ⋅ | | | h , and the errors of the optimal control in ∥ ⋅ ∥ L 2 ⁢ ( Ω ) . We also report the approximation of ∥ y ¯ − y ¯ h ∥ L ∞ ⁢ ( Ω ) by ∥ Π h ⁢ y ¯ − y ¯ h ∥ ∞ , which is the ∞-norm of the finite element coefficient vector representing Π h ⁢ y ¯ − y ¯ h .

The approximate solution y ¯ h j and the operator Π h j are denoted by y ¯ j and Π j , respectively, in all the tables.

Remark 6.1

We have also tested the SIP method based on the P 3 finite element. The results are similar to the ones for the P 2 element for the examples in Section 6.1 and Section 6.2.

6.1 Square Example

Let Ω = ( − 4 , 4 ) 2 . We consider the optimal control problem (1.1) with ψ ⁢ ( x ) = | x | 2 − 1 , β = 1 , g = 0 and

y d ⁢ ( x ) = { Δ 2 ⁢ y ¯ ⁢ ( x ) + y ¯ ⁢ ( x ) if ⁢ | x | > 1 , Δ 2 ⁢ y ¯ ⁢ ( x ) + y ¯ ⁢ ( x ) + 2 if ⁢ | x | ≤ 1 ,

where the optimal state y ¯ is given by

y ¯ ⁢ ( x ) = { | x | 2 − 1 if ⁢ | x | ≤ 1 , v ⁢ ( | x | ) + [ 1 − ϕ ⁢ ( | x | ) ] ⁢ w ⁢ ( x ) if ⁢ 1 ≤ | x | ≤ 3 , w ⁢ ( x ) if ⁢ | x | ≥ 3 ,

with

v ⁢ ( r ) = ( r 2 − 1 ) ⁢ ( 1 − r − 1 2 ) 4 + 1 4 ⁢ ( r − 1 ) 2 ⁢ ( r − 3 ) 4 , ϕ ⁢ ( r ) = [ 1 + 4 ⁢ ( r − 1 2 ) + 10 ⁢ ( r − 1 2 ) 2 + 20 ⁢ ( r − 1 2 ) 3 ] ⁢ ( 1 − r − 1 2 ) 4 , w ⁢ ( x ) = 2 ⁢ sin ⁡ ( π 8 ⁢ ( x 1 + 4 ) ) ⁢ sin ⁡ ( π 8 ⁢ ( x 2 + 4 ) ) .

This example from [26] is designed such that 𝒜 is the disc defined by | x | ≤ 1 , and y ¯ and u ¯ have zero traces on the boundary. The numerical results for k = 1 and k = 2 are reported in Table 1 and Table 2, respectively.

Table 1

Example on square ( k = 1 )

𝑗	∥ y ¯ − y ¯ j ∥ L 2 ⁢ ( Ω )	Order	\| \| \| y ¯ − y ¯ j \| \| \| h	Order	∥ Π j ⁢ y ¯ − y ¯ j ∥ ∞	Order	∥ u ¯ − u ¯ j ∥ L 2 ⁢ ( Ω )	Order
0	4.5810 × 10 1	—	5.3417 × 10 1	—	9.2024 × 10 0	—	1.8989 × 10 1	–
1	3.5379 × 10 1	0.37	3.6566 × 10 1	0.55	9.5624 × 10 0	−0.06	4.5579 × 10 1	−1.26
2	3.0432 × 10 0	3.54	4.0927 × 10 0	3.16	1.1012 × 10 0	3.12	7.5447 × 10 0	2.59
3	1.3713 × 10 0	1.15	2.7151 × 10 0	0.59	4.1099 × 10 − 1	1.42	6.1330 × 10 0	0.30
4	1.9044 × 10 − 1	2.85	1.1673 × 10 0	1.22	8.0634 × 10 − 2	2.35	2.9217 × 10 0	1.07
5	4.8280 × 10 − 2	1.98	5.5406 × 10 − 1	1.08	2.5790 × 10 − 2	1.64	1.0700 × 10 0	1.45
6	1.7648 × 10 − 2	1.45	2.7205 × 10 − 1	1.03	7.9764 × 10 − 3	1.69	3.6576 × 10 − 1	1.55
7	1.2534 × 10 − 2	0.49	1.3587 × 10 − 1	1.00	2.9755 × 10 − 3	1.42	1.6622 × 10 − 1	1.14
8	2.4529 × 10 − 3	2.35	6.7554 × 10 − 2	1.01	6.5540 × 10 − 4	2.18	5.8690 × 10 − 2	1.50
9	5.1696 × 10 − 4	2.25	3.3747 × 10 − 2	1.00	1.3571 × 10 − 4	2.27	2.0999 × 10 − 2	1.48

Table 2

Example on square ( k = 2 )

𝑗	∥ y ¯ − y ¯ j ∥ L 2 ⁢ ( Ω )	Order	\| \| \| y ¯ − y ¯ j \| \| \| h	Order	∥ Π j ⁢ y ¯ − y ¯ j ∥ ∞	Order	∥ u ¯ − u ¯ j ∥ L 2 ⁢ ( Ω )	Order
0	3.7032 × 10 0	—	1.4658 × 10 1	—	1.1510 × 10 0	—	9.1462 × 10 0	–
1	8.3726 × 10 0	−1.18	7.3305 × 10 0	1.00	1.7694 × 10 0	−0.62	9.2196 × 10 0	−0.01
2	3.2350 × 10 0	1.37	5.3145 × 10 0	0.46	1.3622 × 10 0	0.38	1.4794 × 10 1	−0.68
3	8.3648 × 10 − 1	1.95	1.5522 × 10 0	1.78	2.6935 × 10 − 1	2.34	6.0855 × 10 0	1.28
4	1.2604 × 10 − 1	2.73	3.3384 × 10 − 1	2.22	5.7337 × 10 − 2	2.23	2.1357 × 10 0	1.51
5	2.8159 × 10 − 2	2.16	7.5093 × 10 − 2	2.15	6.8294 × 10 − 3	3.07	8.1279 × 10 − 1	1.39
6	7.2141 × 10 − 3	1.96	1.8852 × 10 − 2	1.99	2.2932 × 10 − 3	1.57	2.5827 × 10 − 1	1.65
7	1.2074 × 10 − 3	2.58	4.1633 × 10 − 3	2.18	4.1299 × 10 − 4	2.47	9.9289 × 10 − 2	1.38
8	4.1976 × 10 − 4	1.52	1.1245 × 10 − 3	1.89	1.3026 × 10 − 4	1.66	3.3727 × 10 − 2	1.56

Note that the additional regularity assumptions (4.13), (4.14) and (5.17) are satisfied for this example. In view of (5.18), we have

(6.1) u ¯ = − Δ ⁢ y ¯ ∈ H 3 2 − ϵ ⁢ ( Ω ) for all ϵ > 0 ,

and also the Hölder regularity

(6.2) y ¯ ∈ C 2 , 1 2 − ϵ ⁢ ( Ω ¯ ) for all ϵ > 0

by the Sobolev Embedding Theorem (cf. [60, Section 2.8]).

The following are the best possible results allowed by the regularity (5.18), (6.1) and (6.2):

(6.3) ∥ y ¯ − y ¯ h ∥ L 2 ⁢ ( Ω ) ≤ { C ⁢ h 2 , k = 1 , C ⁢ h 3 , k = 2 ,

(6.4) | | | y ¯ − y ¯ h | | | h ≤ { C ⁢ h , k = 1 , C ⁢ h 2 , k = 2 ,

(6.5) ∥ y ¯ − y ¯ h ∥ L ∞ ⁢ ( Ω ) ≤ { C ⁢ h 2 , k = 1 , C ϵ ⁢ h 5 2 − ϵ , k = 2 ,

(6.6) ∥ u ¯ − u ¯ h ∥ L 2 ⁢ ( Ω ) ≤ { C ϵ ⁢ h 3 2 − ϵ , k = 1 , C ϵ ⁢ h 3 2 − ϵ , k = 2 .

The numerical results for | | | y ¯ − y ¯ h | | | h and ∥ u ¯ − u ¯ h ∥ L 2 ⁢ ( Ω ) in Table 1 and Table 2 match exactly (6.4) and (6.6). The orders of convergence for ∥ y ¯ − y ¯ h ∥ L 2 ⁢ ( Ω ) and ∥ Π h ⁢ y ¯ − y ¯ h ∥ L ∞ ⁢ ( Ω ) are less stable, possibly due to the fact that the free boundary is a circle that is not captured exactly by the mesh, and that the quality of the representation of the circle by the meshes varies from refinement level to refinement level. Nevertheless, the numerical results for ∥ Π h ⁢ y ¯ − y ¯ h ∥ L ∞ ⁢ ( Ω ) in Table 1 and in Table 2 (up to the seventh level refinement) indicate (6.5) asymptotically. The order of convergence for ∥ y ¯ − y ¯ h ∥ L 2 ⁢ ( Ω ) in Table 1 also agrees with (6.3). Finally, the order of convergence for ∥ y ¯ − y ¯ h ∥ L 2 ⁢ ( Ω ) in Table 2 appears to be at most 5 / 2 , which is less than the optimal order in (6.3) for k = 2 .

Note that the O ⁢ ( h ) convergence for | | | y ¯ − y ¯ h | | | h in Table 1 agrees with estimate (5.15) (for τ = 1 ), and the O ⁢ ( h 3 2 ) convergence for ∥ u ¯ − u ¯ h ∥ L 2 ⁢ ( Ω ) in Table 2 agrees with the estimate in (5.27). While these two estimates are the only ones we can establish within our theoretical framework that match the numerical results exactly, we have provided a heuristic justification of the other numerical results in Appendix A.

Pictures for the optimal state, the optimal control and the active set computed by the P 1 SIP method at level 7 are depicted in Figure 1.

Figure 1
State, control and active set computed by a uniform mesh (

k
=
1

k=1

,

j
=
7

j=7

)

Figure 1

State, control and active set computed by a uniform mesh ( k = 1 , j = 7 )

6.2 L-Shaped Domain Example

We extend the previous example to the L-shaped domain Ω = ( − 8 , 8 ) 2 ∖ ( [ 0 , 8 ] × [ − 8 , 0 ] ) by shifting the active set to the center of the upper left quadrant (cf. Figure 2) and by adding the corner singular function

ψ s ⁢ ( x ) = r 2 3 ⁢ sin ⁡ ( 2 3 ⁢ θ ) ,

where ( r , θ ) are the polar coordinates of 𝑥.

More precisely, we take β = 1 , g ∈ H 4 ⁢ ( Ω ) such that g = 4 ⁢ ψ s on ∂ Ω , x * = ( − 4 , 4 ) (the center of the upper left quadrant),

ψ ⁢ ( x ) = | x − x * | 2 − 1 + 4 ⁢ ψ s ⁢ ( x ) ,

and

y d ⁢ ( x ) = { Δ 2 ⁢ y ¯ ⁢ ( x ) + y ¯ ⁢ ( x ) if ⁢ | x − x * | > 1 , Δ 2 ⁢ y ¯ ⁢ ( x ) + y ¯ ⁢ ( x ) + 2 if ⁢ | x − x * | ≤ 1 ,

where

y ¯ ⁢ ( x ) = 4 ⁢ ψ s ⁢ ( x ) + { | x − x * | 2 − 1 if ⁢ | x − x * | ≤ 1 , v ⁢ ( | x − x * | ) + [ 1 − ϕ ⁢ ( | x − x * | ) ] ⁢ w ⁢ ( x − x * ) if ⁢ 1 ≤ | x − x * | ≤ 3 , w ⁢ ( x − x * ) if ⁢ | x − x * | ≥ 3 ,

with 𝑣, 𝜙, 𝑤 as in the previous example.

The numerical results for k = 1 and k = 2 on uniform and graded meshes are displayed in Tables 3 and 4 and Tables 5 and 6, respectively.

Due to the singularity at the reentrant corner, we observe O ⁢ ( h 2 3 ) convergence for | | | y ¯ − y ¯ h | | | h on uniform meshes in Table 3 and Table 5, which agrees with estimate (5.15) for τ = 2 3 − ϵ . For graded mesh refinement with grading parameter 0.6 for k = 1 and 0.3 for k = 2 , the convergence of | | | y ¯ − y ¯ h | | | h is improved to O ⁢ ( h ) and O ⁢ ( h 2 ) , respectively, as in the square domain example of Section 6.1. The performance of other approximations on graded meshes is also similar to the ones in Section 6.1.

Table 3

Example on L-shaped domain for uniform meshes ( k = 1 )

𝑗	∥ y ¯ − y ¯ j ∥ L 2 ⁢ ( Ω )	Order	\| \| \| y ¯ − y ¯ j \| \| \| h	Order	∥ Π j ⁢ y ¯ − y ¯ j ∥ ∞	Order	∥ u ¯ − u ¯ j ∥ L 2 ⁢ ( Ω )	Order
0	3.9920 × 10 1	—	4.4007 × 10 1	—	9.4995 × 10 0	—	2.2870 × 10 1	–
1	3.9904 × 10 1	0.00	3.3403 × 10 1	0.40	1.0482 × 10 1	−0.14	4.0431 × 10 1	−0.82
2	3.4537 × 10 0	3.53	5.7522 × 10 0	2.54	1.4600 × 10 0	2.84	7.4429 × 10 0	2.44
3	1.4574 × 10 0	1.24	3.6802 × 10 0	0.64	9.7520 × 10 − 1	0.58	6.1192 × 10 0	0.28
4	2.1943 × 10 − 1	2.73	1.9689 × 10 0	0.90	6.1792 × 10 − 1	0.66	2.9231 × 10 0	1.07
5	5.7923 × 10 − 2	1.92	1.1369 × 10 0	0.79	3.8999 × 10 − 1	0.66	1.0723 × 10 0	1.45
6	1.8490 × 10 − 2	1.65	6.7892 × 10 − 1	0.74	2.4632 × 10 − 1	0.66	3.6672 × 10 − 1	1.55
7	1.4113 × 10 − 2	0.39	4.1360 × 10 − 1	0.72	1.5534 × 10 − 1	0.67	1.6617 × 10 − 1	1.14
8	2.7271 × 10 − 3	2.37	2.5482 × 10 − 1	0.70	9.7860 × 10 − 2	0.67	5.8672 × 10 − 2	1.50

Table 4

Example on L-shaped domain for graded meshes with grading μ = 0.6 ( k = 1 )

𝑗	∥ y ¯ − y ¯ j ∥ L 2 ⁢ ( Ω )	Order	\| \| \| y ¯ − y ¯ j \| \| \| h	Order	∥ Π j ⁢ y ¯ − y ¯ j ∥ ∞	Order	∥ u ¯ − u ¯ j ∥ L 2 ⁢ ( Ω )	Order
0	3.9920 × 10 1	—	4.4007 × 10 1	—	9.4995 × 10 0	—	2.2870 × 10 1	–
1	2.6192 × 10 1	0.61	2.8876 × 10 1	0.61	1.0147 × 10 1	−0.10	3.3432 × 10 1	−0.55
2	2.6055 × 10 0	3.33	4.4817 × 10 0	2.69	1.1320 × 10 0	3.16	7.5177 × 10 0	2.15
3	1.5013 × 10 0	0.80	2.9781 × 10 0	0.59	4.1680 × 10 − 1	1.44	6.1210 × 10 0	0.30
4	2.0284 × 10 − 1	2.89	1.3167 × 10 0	1.18	1.5541 × 10 − 1	1.42	2.9203 × 10 0	1.07
5	5.7017 × 10 − 2	1.83	6.3409 × 10 − 1	1.05	7.3336 × 10 − 2	1.08	1.0716 × 10 0	1.45
6	2.0081 × 10 − 2	1.51	3.1756 × 10 − 1	1.00	3.8830 × 10 − 2	0.92	3.6665 × 10 − 1	1.55
7	1.1602 × 10 − 2	0.79	1.5707 × 10 − 1	1.02	1.5405 × 10 − 2	1.33	1.6612 × 10 − 1	1.14
8	2.5718 × 10 − 3	2.17	7.8324 × 10 − 2	1.00	7.2773 × 10 − 3	1.08	5.8662 × 10 − 2	1.50

Table 5

Example on L-shaped domain for uniform meshes ( k = 2 )

𝑗	∥ y ¯ − y ¯ j ∥ L 2 ⁢ ( Ω )	Order	\| \| \| y ¯ − y ¯ j \| \| \| h	Order	∥ Π j ⁢ y ¯ − y ¯ j ∥ ∞	Order	∥ u ¯ − u ¯ j ∥ L 2 ⁢ ( Ω )	Order
0	4.2702 × 10 0	—	1.1243 × 10 1	—	2.2648 × 10 0	—	8.6780 × 10 0	–
1	8.2827 × 10 0	−0.96	7.3681 × 10 0	0.61	1.7615 × 10 0	0.36	9.1700 × 10 0	−0.08
2	3.3045 × 10 0	1.33	5.5664 × 10 0	0.40	1.3632 × 10 0	0.37	1.4789 × 10 1	−0.69
3	8.3352 × 10 − 1	1.99	1.8441 × 10 0	1.59	4.3771 × 10 − 1	1.64	6.0881 × 10 0	1.28
4	1.3230 × 10 − 1	2.66	7.0735 × 10 − 1	1.38	2.7683 × 10 − 1	0.66	2.1354 × 10 0	1.51
5	3.2286 × 10 − 2	2.03	3.9996 × 10 − 1	0.82	1.7469 × 10 − 1	0.66	8.1272 × 10 − 1	1.39
6	7.3717 × 10 − 3	2.13	2.4830 × 10 − 1	0.69	1.1000 × 10 − 1	0.67	2.5828 × 10 − 1	1.65
7	1.6831 × 10 − 3	2.13	1.5603 × 10 − 1	0.67	6.9299 × 10 − 2	0.67	9.9293 × 10 − 2	1.38

Table 6

Example on L-shaped domain for graded meshes with grading μ = 0.3 ( k = 2 )

𝑗	∥ y ¯ − y ¯ j ∥ L 2 ⁢ ( Ω )	Order	\| \| \| y ¯ − y ¯ j \| \| \| h	Order	∥ Π j ⁢ y ¯ − y ¯ j ∥ ∞	Order	∥ u ¯ − u ¯ j ∥ L 2 ⁢ ( Ω )	Order
0	4.2702 × 10 0	—	1.1243 × 10 1	—	2.2648 × 10 0	—	8.6780 × 10 0	–
1	1.2868 × 10 1	−1.59	1.7278 × 10 1	−0.62	5.9386 × 10 0	−1.39	2.4659 × 10 1	−1.51
2	3.3202 × 10 0	1.95	5.3727 × 10 0	1.69	1.4139 × 10 0	2.07	1.4917 × 10 1	0.73
3	8.3547 × 10 − 1	1.99	1.5466 × 10 0	1.80	2.6935 × 10 − 1	2.39	6.0817 × 10 0	1.29
4	1.2014 × 10 − 1	2.80	3.2962 × 10 − 1	2.23	5.7276 × 10 − 2	2.23	2.1342 × 10 0	1.51
5	2.7127 × 10 − 2	2.15	7.4216 × 10 − 2	2.15	6.8314 × 10 − 3	3.07	8.1287 × 10 − 1	1.39
6	8.0864 × 10 − 3	1.75	1.8834 × 10 − 2	1.98	2.3745 × 10 − 3	1.52	2.5826 × 10 − 1	1.65
7	1.4089 × 10 − 3	2.52	4.0875 × 10 − 3	2.20	4.8728 × 10 − 4	2.28	9.9273 × 10 − 2	1.38

The pictures of the optimal state, the optimal control and the active set computed by the P 1 SIP method on the uniform mesh at level 6 are depicted in Figure 2. We have also included the picture of the active set computed on the graded mesh, which shows a better approximation.

Figure 2
State, control and active set computed by a uniform mesh and the active set computed by a graded mesh for the L-shaped domain example (

k
=
1

k=1

,

j
=
6

j=6

)

Figure 2

State, control and active set computed by a uniform mesh and the active set computed by a graded mesh for the L-shaped domain example ( k = 1 , j = 6 )

6.3 Cube Example

We have also tested the SIP method on a three-dimensional example, which is an extension of the example in Section 6.1.

Let Ω = ( − 4 , 4 ) 3 . We consider the optimal control problem with ψ ⁢ ( x ) = | x | 2 − 1 , β = 1 , g = 0 and

y d ⁢ ( x ) = { Δ 2 ⁢ y ¯ ⁢ ( x ) + y ¯ ⁢ ( x ) if ⁢ | x | > 1 , Δ 2 ⁢ y ¯ ⁢ ( x ) + y ¯ ⁢ ( x ) + 2 if ⁢ | x | ≤ 1 ,

where the optimal state y ¯ is given by

y ¯ ⁢ ( x ) = { | x | 2 − 1 if ⁢ | x | ≤ 1 , v ⁢ ( | x | ) + [ 1 − ϕ ⁢ ( | x | ) ] ⁢ w ⁢ ( x ) if ⁢ 1 ≤ | x | ≤ 3 , w ⁢ ( x ) if ⁢ | x | ≥ 3 ,

with radial functions 𝑣, 𝜙 as in the two-dimensional example, and

w ⁢ ( x ) = 2 ⁢ sin ⁡ ( π 8 ⁢ ( x 1 + 4 ) ) ⁢ sin ⁡ ( π 8 ⁢ ( x 2 + 4 ) ) ⁢ sin ⁡ ( π 8 ⁢ ( x 3 + 4 ) ) .

The numerical results are reported in Table 7 for k = 1 and Table 8 for k = 2 . The performance of the SIP method is similar to that in Section 6.1.

Table 7

Example on cube ( k = 1 )

𝑗	∥ y ¯ − y ¯ j ∥ L 2 ⁢ ( Ω )	Order	\| \| \| y ¯ − y ¯ j \| \| \| h	Order	∥ Π j ⁢ y ¯ − y ¯ j ∥ ∞	Order	∥ u ¯ − u ¯ j ∥ L 2 ⁢ ( Ω )	Order
0	1.6410 × 10 1	—	2.3815 × 10 1	—	9.7347 × 10 − 1	—	2.3152 × 10 1	–
1	7.5940 × 10 0	1.11	1.4455 × 10 1	0.72	1.1095 × 10 0	−0.19	2.0836 × 10 1	0.15
2	2.7130 × 10 0	1.48	4.8940 × 10 0	1.56	5.5711 × 10 − 1	0.99	9.4985 × 10 0	1.13
3	1.3542 × 10 0	1.00	3.3239 × 10 0	0.56	4.2026 × 10 − 1	0.41	7.6705 × 10 0	0.31
4	9.5533 × 10 − 2	3.83	1.5610 × 10 0	1.09	8.2658 × 10 − 2	2.35	3.4858 × 10 0	1.14
5	4.5772 × 10 − 2	1.06	7.7556 × 10 − 1	1.01	1.8483 × 10 − 2	2.16	1.3494 × 10 0	1.37

Table 8

Example on cube ( k = 2 )

𝑗	∥ y ¯ − y ¯ j ∥ L 2 ⁢ ( Ω )	Order	\| \| \| y ¯ − y ¯ j \| \| \| h	Order	∥ Π j ⁢ y ¯ − y ¯ j ∥ ∞	Order	∥ u ¯ − u ¯ j ∥ L 2 ⁢ ( Ω )	Order
0	2.6455 × 10 1	—	4.8283 × 10 1	—	2.0914 × 10 0	—	2.8380 × 10 1	–
1	3.0710 × 10 0	3.11	3.9681 × 10 0	3.60	2.7922 × 10 − 1	2.90	7.0809 × 10 0	2.00
2	2.3105 × 10 0	0.41	2.9753 × 10 0	0.42	1.5391 × 10 0	−2.46	9.9390 × 10 0	−0.49
3	2.6989 × 10 − 1	3.10	6.1120 × 10 − 1	2.28	6.7021 × 10 − 2	4.52	3.1649 × 10 0	1.65
4	1.0293 × 10 − 2	4.71	1.6018 × 10 − 1	1.93	2.1690 × 10 − 2	1.63	9.7816 × 10 − 1	1.69
5	4.2952 × 10 − 3	1.26	4.5806 × 10 − 2	1.81	3.9147 × 10 − 3	2.47	3.8867 × 10 − 1	1.33

7 Concluding Remarks

We have developed a SIP method for an elliptic distributed optimal control problem with pointwise state constraints on general polygonal domains. The resulting discrete problems are quadratic programs with box constraints that can be solved efficiently by a primal-dual active set method.

By using graded meshes, we can show O ⁢ ( h ) convergence (up to a log factor) for the optimal control in the L 2 norm, and for the optimal state in the L 2 norm, energy norm and L ∞ norm under general assumptions on the data. Better convergence can also be established under additional regularity assumptions.

We have tested the SIP method based on the P 1 and P 2 discontinuous finite element spaces. The numerical results are better than our theoretical results. A heuristic justification of the numerical results is provided in the appendix where an optimal control problem without the pointwise state constraints, but with the correct assumption on the regularity of the data, is considered.

We expect the interior maximum norm estimate (3.21) can be extended to the other discontinuous Galerkin methods in [4] along the lines taken in [47], in which case the analysis in this paper can also be applied to the discretizations of (2.1) based on these methods.

Even though our theory is purely for two-dimensional domains, numerical results indicate that the SIP method has similar performance in three dimensions. Therefore, the discontinuous Galerkin approach may provide higher-order methods for the optimal control problem on general polyhedral domains where the discrete problems are quadratic programs with box constraints that can be solved efficiently by a primal-dual active set algorithm. A rigorous analysis of the SIP method in three dimensions will require an interior maximum norm error estimate for discontinuous Galerkin methods in three dimensions, which is still absent from the literature.

Funding source: National Science Foundation

Award Identifier / Grant number: DMS-19-13035

Award Identifier / Grant number: DMS-22-08404

Funding statement: The work of the first and third authors was supported in part by the National Science Foundation under Grant No. DMS-19-13035 and Grant No. DMS-22-08404.

A An Optimal Control Problem Without State Constraints

In order to give a heuristic justification of the numerical results observed in Section 6.1, we consider the elliptic optimal control problem (1.1) without the state constraint (1.4) on a square Ω. The optimal control problem is equivalent to the following problem after dropping the bars over 𝑦 and 𝑢 (cf. [52]): find ( u , y ) ∈ H 0 1 ⁢ ( Ω ) × H 0 1 ⁢ ( Ω ) such that

(A.1) ∫ Ω y ⁢ v ⁢ d x + β ⁢ ∫ Ω ∇ u ⋅ ∇ v ⁢ d ⁢ x = ∫ Ω y d ⁢ v ⁢ d x for all ⁢ v ∈ H 0 1 ⁢ ( Ω ) ,

(A.2) ∫ Ω ∇ y ⋅ ∇ z ⁢ d ⁢ x − ∫ Ω u ⁢ z ⁢ d x = 0 for all ⁢ z ∈ H 0 1 ⁢ ( Ω ) .

To match the regularity of the Lagrange multiplier in (4.17), we replace y d ∈ L 2 ⁢ ( Ω ) by a linear functional 𝐹 that belongs to H − 1 2 − ϵ ⁢ ( Ω ) for any ϵ > 0 and rewrite (A.1)–(A.2) concisely as

(A.3) B ⁢ ( ( u , y ) , ( v , z ) ) = F ⁢ ( v ) for all ⁢ ( v , z ) ∈ H 0 1 ⁢ ( Ω ) × H 0 1 ⁢ ( Ω ) ,

where

B ⁢ ( ( u , y ) , ( v , z ) ) = β ⁢ ∫ Ω ∇ u ⋅ ∇ v ⁢ d ⁢ x + ∫ Ω y ⁢ v ⁢ d x − ∫ Ω u ⁢ z ⁢ d x + ∫ Ω ∇ y ⋅ ∇ z ⁢ d ⁢ x .

Remark A.1

Equation (2.2) is identical to (A.3) with

F ⁢ ( v ) = ∫ Ω y d ⁢ v ⁢ d x + ∫ Ω v ⁢ d μ ,

where the Lagrange multiplier 𝜇 satisfies (4.17).

The bilinear form B ⁢ ( ⋅ , ⋅ ) is clearly bounded on H 0 1 ⁢ ( Ω ) × H 0 1 ⁢ ( Ω ) . It is also coercive because

B ⁢ ( ( v , z ) , ( v , z ) ) = β ⁢ | v | H 1 ⁢ ( Ω ) 2 + | z | H 1 ⁢ ( Ω ) 2 for all ⁢ ( y , z ) ∈ H 0 1 ⁢ ( Ω ) × H 0 1 ⁢ ( Ω ) .

Since Ω is convex, elliptic regularity and (A.1) (with y d replaced by 𝐹) imply (cf. [37])

(A.4) u ∈ H 3 2 − ϵ ⁢ ( Ω ) for any ϵ > 0 ,

and then, since Ω is a square and u ∈ H 3 2 − ϵ ⁢ ( Ω ) ∩ H 0 1 ⁢ ( Ω ) , elliptic regularity and (A.2) imply (cf. [38, 15])

(A.5) y ∈ H 7 2 − ϵ ⁢ ( Ω ) for any ϵ > 0 .

The SIP method for (A.3) is to find ( u h , y h ) ∈ Y h × Y h such that

(A.6) ∫ Ω y h ⁢ v ⁢ d x + β ⁢ a h ⁢ ( u h , v ) = F ⁢ ( v ) for all ⁢ v ∈ Y h ,

(A.7) a h ⁢ ( y h , z ) − ∫ Ω u h ⁢ z ⁢ d x = 0 for all ⁢ z ∈ Y h .

Remark A.2

The discussion below can also be carried out for the SIP method defined by (A.6) and (A.7) that involves additional technicalities. For simplicity, we consider instead conforming finite element methods for (A.3).

Recall V h = Y h ∩ H 0 1 ⁢ ( Ω ) is the P k ( k = 1 , 2 ) Lagrange finite element space associated with a uniform triangulation T h of Ω. The discrete problem is to find ( u h , y h ) ∈ V h × V h such that

(A.8) ∫ Ω y h ⁢ v ⁢ d x + β ⁢ ∫ Ω ∇ u h ⋅ ∇ v ⁢ d ⁢ x = F ⁢ ( v ) for all ⁢ v ∈ V h ,

(A.9) ∫ Ω ∇ y h ⋅ ∇ z ⁢ d ⁢ x − ∫ Ω u h ⁢ z ⁢ d x = 0 for all ⁢ z ∈ V h ,

or equivalently,

B ⁢ ( ( u h , y h ) , ( v , z ) ) = F ⁢ ( v ) for all ⁢ ( v , z ) ∈ V h × V h .

Comparing (A.1)–(A.2) (where y d is replaced by 𝐹) with (A.8)–(A.9), we have the Galerkin relations

(A.10) ∫ Ω ( y − y h ) ⁢ v ⁢ d x + β ⁢ ∫ Ω ∇ ( u − u h ) ⋅ ∇ v ⁢ d ⁢ x = 0 for all ⁢ v ∈ V h ,

(A.11) ∫ Ω ∇ ( y − y h ) ⋅ ∇ z ⁢ d ⁢ x − ∫ Ω ( u − u h ) ⁢ z ⁢ d x = 0 for all ⁢ z ∈ V h .

Combing the boundedness and coercivity of B ⁢ ( ⋅ , ⋅ ) , the regularity (A.4)–(A.5), the Galerkin relations (A.10)–(A.11) and standard interpolation error estimates, we obtain the error estimate

(A.12) | u − u h | H 1 ⁢ ( Ω ) + | y − y h | H 1 ⁢ ( Ω ) ≤ C ⁢ inf ( v , z ) ∈ V h × V h [ | u − v | H 1 ⁢ ( Ω ) + | y − z | H 1 ⁢ ( Ω ) ] ≤ C ϵ ⁢ h 1 2 − ϵ ,

which then implies

(A.13) ∥ u − u h ∥ L 2 ⁢ ( Ω ) + ∥ y − y h ∥ L 2 ⁢ ( Ω ) ≤ C ϵ ⁢ h 3 2 − ϵ

by a standard duality argument.

In the case where k = 1 , we can improve the estimate for | y − y h | H 1 ⁢ ( Ω ) in (A.12) through relation (A.11) with z = Π h ⁢ y − y h , which yields

(A.14) ∫ Ω ∇ ( y − y h ) ⋅ ∇ ( y − y h ) ⁡ d ⁢ x + ∫ Ω ∇ ( y − y h ) ⋅ ∇ ( Π h ⁢ y − y ) ⁡ d ⁢ x = ∫ Ω ( u − u h ) ⁢ ( y − y h ) ⁢ d x + ∫ Ω ( u − u h ) ⁢ ( Π h ⁢ y − y ) ⁢ d x .

It follows from (A.5), (A.13), (A.14) and standard interpolation error estimates that

| y − y h | H 1 ⁢ ( Ω ) 2 ≤ C 1 ⁢ h ⁢ | y − y h | H 1 ⁢ ( Ω ) + C 2 , ϵ 1 ⁢ h 3 − ϵ 1 + C 3 , ϵ 2 ⁢ h 5 2 − ϵ 2 ,

which implies (with ϵ 1 = 1 and ϵ 2 = 1 2 )

(A.15) | y − y h | H 1 ⁢ ( Ω ) ≤ C ⁢ h .

Remark A.3

The estimate in (A.13) for ∥ u − u h ∥ L 2 ⁢ ( Ω ) and estimate (A.15) for | y − y h | H 1 ⁢ ( Ω ) match the results observed in Table 1.

In the case where k = 2 , we need to use a duality argument that is based on the equivalence of (A.3) with the following fourth-order boundary value problem (cf. (2.2) and Remark A.1):

(A.16) β ⁢ ∫ Ω ( Δ ⁢ y ) ⁢ ( Δ ⁢ z ) ⁢ d x + ∫ Ω y ⁢ z ⁢ d x = F ⁢ ( z ) for all ⁢ z ∈ H 2 ⁢ ( Ω ) ∩ H 0 1 ⁢ ( Ω ) .

Note that (A.3) can be interpreted as a mixed formulation of (A.16).

Let ϕ ∈ H 2 ⁢ ( Ω ) ∩ H 0 1 ⁢ ( Ω ) be defined by

(A.17) β ⁢ ∫ Ω ( Δ ⁢ v ) ⁢ ( Δ ⁢ ϕ ) ⁢ d x + ∫ Ω v ⁢ ϕ ⁢ d x = ∫ Ω v ⁢ ( y − y h ) ⁢ d x for all ⁢ v ∈ H 2 ⁢ ( Ω ) ∩ H 0 1 ⁢ ( Ω ) .

Since Ω is a square, we have

(A.18) ∥ ϕ ∥ H 3 − ϵ ⁢ ( Ω ) ≤ C ϵ ⁢ ∥ y − y h ∥ L 2 ⁢ ( Ω )

by the elliptic regularity result in [10] for the biharmonic equation with the boundary conditions of simply supported plates.

Let y ~ h ∈ H 0 1 ⁢ ( Ω ) be defined by

(A.19) ∫ Ω ∇ y ~ h ⋅ ∇ z ⁢ d ⁢ x = ∫ Ω u h ⁢ z ⁢ d x for all ⁢ z ∈ H 0 1 ⁢ ( Ω ) .

Then we have

(A.20) ∥ y ~ h ∥ H 3 − ϵ ⁢ ( Ω ) ≤ C ϵ ⁢ ∥ u h ∥ H 1 ⁢ ( Ω )

by the elliptic regularity for a square (cf. [37]).

In view of (A.9) and (A.19), y h ∈ V h is precisely the approximation of y ~ h generated by the standard P 2 Lagrange finite element method. Consequently, we have

| y ~ h − y h | H 1 ⁢ ( Ω ) ≤ C ϵ ⁢ h 2 − ϵ ⁢ | y ~ h | H 3 − ϵ ⁢ ( Ω ) ≤ C ϵ ⁢ h 2 − ϵ

by (A.20) and the fact that ∥ u h ∥ H 1 ⁢ ( Ω ) is uniformly bounded by (A.12), and then the estimate

(A.21) ∥ y ~ h − y h ∥ L 2 ⁢ ( Ω ) ≤ C ⁢ h ⁢ | y ~ h − y h | H 1 ⁢ ( Ω ) ≤ C ϵ ⁢ h 3 − ϵ

follows from a standard duality argument.

Now we take v = y − y ~ h in (A.17) and then use (A.2) and (A.19) to obtain the relation

β ⁢ ∫ Ω ∇ ( u − u h ) ⁢ ∇ ϕ ⁢ d ⁢ x + ∫ Ω ( y − y h ) ⁢ ϕ ⁢ d x + ∫ Ω ( y h − y ~ h ) ⁢ ϕ ⁢ d x = ∫ Ω ( y − y h ) 2 ⁢ d x + ∫ Ω ( y h − y ~ h ) ⁢ ( y − y h ) ⁢ d x ,

which implies, through (A.10),

∥ y − y h ∥ L 2 ⁢ ( Ω ) 2 = β ⁢ ∫ Ω ∇ ( u − u h ) ⋅ ∇ ( ϕ − Π h ⁢ ϕ ) ⁡ d ⁢ x + ∫ Ω ( y − y h ) ⁢ ( ϕ − Π h ⁢ ϕ ) ⁢ d x + ∫ Ω ( y h − y ~ h ) ⁢ ϕ ⁢ d x − ∫ Ω ( y h − y ~ h ) ⁢ ( y − y h ) ⁢ d x ≤ β ⁢ | u − u h | H 1 ⁢ ( Ω ) ⁢ | ϕ − Π h ⁢ ϕ | H 1 ⁢ ( Ω ) + ∥ y − y h ∥ L 2 ⁢ ( Ω ) ⁢ ∥ ϕ − Π h ⁢ ϕ ∥ L 2 ⁢ ( Ω ) + ∥ y h − y ~ h ∥ L 2 ⁢ ( Ω ) ⁢ ∥ ϕ ∥ L 2 ⁢ ( Ω ) + ∥ y h − y ~ h ∥ L 2 ⁢ ( Ω ) ⁢ ∥ y − y h ∥ L 2 ⁢ ( Ω ) .

In view of (A.13), (A.18), (A.21) and standard interpolation error estimates, we conclude that

∥ y − y h ∥ L 2 ⁢ ( Ω ) 2 ≤ C 1 , ϵ 1 ⁢ h 5 2 − ϵ 1 ⁢ ∥ y − y h ∥ L 2 ⁢ ( Ω ) + C 2 , ϵ 2 ⁢ h 3 − ϵ 2 ⁢ ∥ y − y h ∥ L 2 ⁢ ( Ω ) ,

and hence

(A.22) ∥ y − y h ∥ L 2 ⁢ ( Ω ) ≤ C ϵ ⁢ h 5 2 − ϵ

by the inequality of arithmetic and geometric means, which improves the estimate in (A.13).

Finally, we return to relation (A.14) and use (A.13), (A.22) to obtain the estimate

| y − y h | H 1 ⁢ ( Ω ) 2 ≤ C 1 ⁢ h 2 ⁢ | y − y h | H 1 ⁢ ( Ω ) + C 2 , ϵ 1 ⁢ h 4 − ϵ 1 + C 3 , ϵ 2 ⁢ h 9 2 − ϵ 2 ,

and hence

(A.23) | y − y h | H 1 ⁢ ( Ω ) ≤ C ϵ ⁢ h 2 − ϵ .

Remark A.4

The estimate in (A.13) for ∥ u − u h ∥ L 2 ⁢ ( Ω ) , the estimate in (A.22) for ∥ y − y h ∥ L 2 ⁢ ( Ω ) and the estimate in (A.23) for | y − y h | H 1 ⁢ ( Ω ) match the results observed in Table 2.

References

[1] R. A. Adams and J. J. F. Fournier, Sobolev Spaces, 2nd ed., Pure Appl. Math. (Amsterdam) 140, Elsevier/Academic Press, Amsterdam, 2003. Search in Google Scholar

[2] T. Apel, A.-M. Sändig and J. R. Whiteman, Graded mesh refinement and error estimates for finite element solutions of elliptic boundary value problems in non-smooth domains, Math. Methods Appl. Sci. 19 (1996), no. 1, 63–85. 10.1002/(SICI)1099-1476(19960110)19:1<63::AID-MMA764>3.0.CO;2-SSearch in Google Scholar

[3] D. N. Arnold, An interior penalty finite element method with discontinuous elements, SIAM J. Numer. Anal. 19 (1982), no. 4, 742–760. 10.1137/0719052Search in Google Scholar

[4] D. N. Arnold, F. Brezzi, B. Cockburn and L. D. Marini, Unified analysis of discontinuous Galerkin methods for elliptic problems, SIAM J. Numer. Anal. 39 (2001/02), no. 5, 1749–1779. 10.1137/S0036142901384162Search in Google Scholar

[5] I. Babuška, Finite element method for domains with corners, Computing 6 (1970), 264–273. 10.1007/BF02238811Search in Google Scholar

[6] I. Babuška, R. B. Kellogg and J. Pitkäranta, Direct and inverse error estimates for finite elements with mesh refinements, Numer. Math. 33 (1979), no. 4, 447–471. 10.1007/BF01399326Search in Google Scholar

[7] M. Bergounioux, K. Ito and K. Kunisch, Primal-dual strategy for constrained optimal control problems, SIAM J. Control Optim. 37 (1999), no. 4, 1176–1194. 10.1137/S0363012997328609Search in Google Scholar

[8] M. Bergounioux and K. Kunisch, Augmented Lagrangian techniques for elliptic state constrained optimal control problems, SIAM J. Control Optim. 35 (1997), no. 5, 1524–1543. 10.1137/S036301299529330XSearch in Google Scholar

[9] M. Bergounioux and K. Kunisch, Primal-dual strategy for state-constrained optimal control problems, Comput. Optim. Appl. 22 (2002), no. 2, 193–224. Search in Google Scholar

[10] H. Blum and R. Rannacher, On the boundary value problem of the biharmonic operator on domains with angular corners, Math. Methods Appl. Sci. 2 (1980), no. 4, 556–581. 10.1002/mma.1670020416Search in Google Scholar

[11] J. H. Bramble and S. R. Hilbert, Estimation of linear functionals on Sobolev spaces with application to Fourier transforms and spline interpolation, SIAM J. Numer. Anal. 7 (1970), 112–124. 10.1137/0707006Search in Google Scholar

[12] J. J. Brannick, H. Li and L. T. Zikatanov, Uniform convergence of the multigrid 𝑉-cycle on graded meshes for corner singularities, Numer. Linear Algebra Appl. 15 (2008), no. 2–3, 291–306. 10.1002/nla.574Search in Google Scholar

[13] S. C. Brenner, Two-level additive Schwarz preconditioners for nonconforming finite elements, Domain Decomposition Methods in Scientific and Engineering Computing, Contemp. Math. 180, American Mathematical Society, Providence (1994), 9–14. 10.1090/conm/180/01951Search in Google Scholar

[14] S. C. Brenner, Two-level additive Schwarz preconditioners for nonconforming finite element methods, Math. Comp. 65 (1996), no. 215, 897–921. 10.1090/S0025-5718-96-00746-6Search in Google Scholar

[15] S. C. Brenner, Overcoming corner singularities using multigrid methods, SIAM J. Numer. Anal. 35 (1998), no. 5, 1883–1892. 10.1137/S0036142996308022Search in Google Scholar

[16] S. C. Brenner, Poincaré–Friedrichs inequalities for piecewise H 1 functions, SIAM J. Numer. Anal. 41 (2003), no. 1, 306–324. 10.1137/S0036142902401311Search in Google Scholar

[17] S. C. Brenner, Discrete Sobolev and Poincaré inequalities for piecewise polynomial functions, Electron. Trans. Numer. Anal. 18 (2004), 42–48. Search in Google Scholar

[18] S. C. Brenner, J. Cui and L.-Y. Sung, Multigrid methods for the symmetric interior penalty method on graded meshes, Numer. Linear Algebra Appl. 16 (2009), no. 6, 481–501. 10.1002/nla.630Search in Google Scholar

[19] S. C. Brenner, C. B. Davis and L.-Y. Sung, A partition of unity method for a class of fourth order elliptic variational inequalities, Comput. Methods Appl. Mech. Engrg. 276 (2014), 612–626. 10.1016/j.cma.2014.04.004Search in Google Scholar

[20] S. C. Brenner, J. Gedicke and L.-Y. Sung, C 0 interior penalty methods for an elliptic distributed optimal control problem on nonconvex polygonal domains with pointwise state constraints, SIAM J. Numer. Anal. 56 (2018), no. 3, 1758–1785. 10.1137/17M1140649Search in Google Scholar

[21] S. C. Brenner, M. Oh, S. Pollock, K. Porwal, M. Schedensack and N. S. Sharma, A C 0 interior penalty method for elliptic distributed optimal control problems in three dimensions with pointwise state constraints, Topics in Numerical Partial Differential Equations and Scientific Computing, IMA Vol. Math. Appl. 160, Springer, New York (2016), 1–22. 10.1007/978-1-4939-6399-7_1Search in Google Scholar

[22] S. C. Brenner, M. Oh and L.-Y. Sung, P 1 finite element methods for an elliptic state-constrained distributed optimal control problem with Neumann boundary conditions, Results Appl. Math. 8 (2020), Paper No. 100090. 10.1016/j.rinam.2019.100090Search in Google Scholar

[23] S. C. Brenner and L. R. Scott, The Mathematical Theory of Finite Element Methods, 3rd ed., Texts Appl. Math. 15, Springer, New York, 2008. 10.1007/978-0-387-75934-0Search in Google Scholar

[24] S. C. Brenner and L.-Y. Sung, A new convergence analysis of finite element methods for elliptic distributed optimal control problems with pointwise state constraints, SIAM J. Control Optim. 55 (2017), no. 4, 2289–2304. 10.1137/16M1088090Search in Google Scholar

[25] S. C. Brenner and L.-Y. Sung, An interior maximum norm error estimate for the symmetric interior penalty method on planar polygonal domains, Comput. Methods Appl. Math. 23 (2023), no. 1, 49–63. 10.1515/cmam-2022-0127Search in Google Scholar

[26] S. C. Brenner, L.-Y. Sung and J. Gedicke, P 1 finite element methods for an elliptic optimal control problem with pointwise state constraints, IMA J. Numer. Anal. 40 (2020), no. 1, 1–28. 10.1093/imanum/dry071Search in Google Scholar

[27] S. C. Brenner, L.-Y. Sung and Z. Tan, A C 1 virtual element method for an elliptic distributed optimal control problem with pointwise state constraints, Math. Models Methods Appl. Sci. 31 (2021), no. 14, 2887–2906. 10.1142/S0218202521500640Search in Google Scholar

[28] S. C. Brenner, L.-Y. Sung and Y. Zhang, A quadratic C 0 interior penalty method for an elliptic optimal control problem with state constraints, Recent Developments in Discontinuous Galerkin Finite Element Methods for Partial Differential Equations, IMA Vol. Math. Appl. 157, Springer, Cham (2014), 97–132. 10.1007/978-3-319-01818-8_4Search in Google Scholar

[29] S. C. Brenner, L.-Y. Sung and Y. Zhang, C 0 interior penalty methods for an elliptic state-constrained optimal control problem with Neumann boundary condition, J. Comput. Appl. Math. 350 (2019), 212–232. 10.1016/j.cam.2018.10.015Search in Google Scholar

[30] L. A. Caffarelli and A. Friedman, The obstacle problem for the biharmonic operator, Ann. Sc. Norm. Super. Pisa Cl. Sci. (4) 6 (1979), no. 1, 151–184. Search in Google Scholar

[31] E. Casas, Control of an elliptic problem with pointwise state constraints, SIAM J. Control Optim. 24 (1986), no. 6, 1309–1318. 10.1137/0324078Search in Google Scholar

[32] E. Casas, M. Mateos and B. Vexler, New regularity results and improved error estimates for optimal control problems with state constraints, ESAIM Control Optim. Calc. Var. 20 (2014), no. 3, 803–822. 10.1051/cocv/2013084Search in Google Scholar

[33] H. Chen, Pointwise error estimates of the local discontinuous Galerkin method for a second order elliptic problem, Math. Comp. 74 (2005), no. 251, 1097–1116. 10.1090/S0025-5718-04-01700-4Search in Google Scholar

[34] Z. Chen and H. Chen, Pointwise error estimates of discontinuous Galerkin methods with penalty for second-order elliptic problems, SIAM J. Numer. Anal. 42 (2004), no. 3, 1146–1166. 10.1137/S0036142903421527Search in Google Scholar

[35] S. Cherednichenko and A. Rösch, Error estimates for the discretization of elliptic control problems with pointwise control and state constraints, Comput. Optim. Appl. 44 (2009), no. 1, 27–55. 10.1007/s10589-008-9186-5Search in Google Scholar

[36] P. G. Ciarlet, The Finite Element Method for Elliptic Problems, Stud. Math. Appl. 4, North-Holland, Amsterdam, 1978. 10.1115/1.3424474Search in Google Scholar

[37] M. Dauge, Elliptic Boundary Value Problems on Corner Domains, Lecture Notes in Math. 1341, Springer, Berlin, 1988. 10.1007/BFb0086682Search in Google Scholar

[38] M. Dauge, S. Nicaise, M. Bourlard and J. M.-S. Lubuma, Coefficients des singularités pour des problèmes aux limites elliptiques sur un domaine à points coniques. II. Quelques opérateurs particuliers, RAIRO Modél. Math. Anal. Numér. 24 (1990), no. 3, 343–367. 10.1051/m2an/1990240303431Search in Google Scholar

[39] T. Dupont and R. Scott, Polynomial approximation of functions in Sobolev spaces, Math. Comp. 34 (1980), no. 150, 441–463. 10.1090/S0025-5718-1980-0559195-7Search in Google Scholar

[40] I. Ekeland and R. Témam, Convex Analysis and Variational Problems, Class. Appl. Math. 28, Society for Industrial and Applied Mathematics, Philadelphia, 1999. 10.1137/1.9781611971088Search in Google Scholar

[41] Y. Epshteyn and B. Rivière, Estimation of penalty parameters for symmetric interior penalty Galerkin methods, J. Comput. Appl. Math. 206 (2007), no. 2, 843–872. 10.1016/j.cam.2006.08.029Search in Google Scholar

[42] L. C. Evans, Partial Differential Equations, 2nd ed., Grad. Stud. Math. 19, American Mathematical Society, Providence, 2010. Search in Google Scholar

[43] J. Frehse, Zum Differenzierbarkeitsproblem bei Variationsungleichungen höherer Ordnung, Abh. Math. Semin. Univ. Hamburg 36 (1971), 140–149. 10.1007/BF02995917Search in Google Scholar

[44] J. Frehse, On the regularity of the solution of the biharmonic variational inequality, Manuscripta Math. 9 (1973), 91–103. 10.1007/BF01320669Search in Google Scholar

[45] R. Fritzsch and P. Oswald, Zur optimalen Gitterwahl bei Finite-Elemente-Approximationen, Wiss. Z. Tech. Univ. Dresden 37 (1988), no. 3, 155–158. Search in Google Scholar

[46] P. Grisvard, Elliptic Problems in Nonsmooth Domains, Monogr. Stud. Math. 24, Pitman, Boston, 1985. Search in Google Scholar

[47] J. Guzmán, Pointwise error estimates for discontinuous Galerkin methods with lifting operators for elliptic problems, Math. Comp. 75 (2006), no. 255, 1067–1085. 10.1090/S0025-5718-06-01823-0Search in Google Scholar

[48] M. Hintermüller, K. Ito and K. Kunisch, The primal-dual active set strategy as a semismooth Newton method, SIAM J. Optim. 13 (2002), no. 3, 865–888. 10.1137/S1052623401383558Search in Google Scholar

[49] M. Hinze, R. Pinnau, M. Ulbrich and S. Ulbrich, Optimization with PDE Constraints, Math. Model. Theory Appl. 23, Springer, New York, 2009. Search in Google Scholar

[50] K. Ito and K. Kunisch, Lagrange Multiplier Approach to Variational Problems and Applications, Adv. Des. Control 15, Society for Industrial and Applied Mathematics, Philadelphia, 2008. 10.1137/1.9780898718614Search in Google Scholar

[51] D. Kinderlehrer and G. Stampacchia, An Introduction to Variational Inequalities and Their Applications, Class. Appl. Math. 31, Society for Industrial and Applied Mathematics, Philadelphia, 2000. 10.1137/1.9780898719451Search in Google Scholar

[52] J.-L. Lions, Optimal Control of Systems Governed by Partial Differential Equations, Grundlehren Math. Wiss. 170, Springer, New York, 1971. 10.1007/978-3-642-65024-6Search in Google Scholar

[53] W. Liu, W. Gong and N. Yan, A new finite element approximation of a state-constrained optimal control problem, J. Comput. Math. 27 (2009), no. 1, 97–114. Search in Google Scholar

[54] C. Meyer, Error estimates for the finite-element approximation of an elliptic control problem with pointwise state and control constraints, Control Cybernet. 37 (2008), no. 1, 51–83. Search in Google Scholar

[55] I. Neitzel, J. Pfefferer and A. Rösch, Finite element discretization of state-constrained elliptic optimal control problems with semilinear state equation, SIAM J. Control Optim. 53 (2015), no. 2, 874–904. 10.1137/140960645Search in Google Scholar

[56] M. Pierre and J. Sokoł owski, Differentiability of projection and applications, Control of Partial Differential Equations and Applications, Lect. Notes Pure Appl. Math. 174, Dekker, New York (1996), 231–240. Search in Google Scholar

[57] B. Rivière, Discontinuous Galerkin Methods for Solving Elliptic and Parabolic Equations, Front. Appl. Math. 35, Society for Industrial and Applied Mathematics, Philadelphia, 2008. 10.1137/1.9780898717440Search in Google Scholar

[58] L. Tartar, An Introduction to Sobolev Spaces and Interpolation Spaces, Lect. Notes Unione Mat. Ital. 3, Springer, Berlin, 2007. Search in Google Scholar

[59] V. Thomée, Galerkin Finite Element Methods for Parabolic Problems, 2nd ed., Springer Ser. Comput. Math. 25, Springer, Berlin, 2006. Search in Google Scholar

[60] H. Triebel, Interpolation Theory, Function Spaces, Differential Operators, North-Holland Math. Libr. 18, North-Holland, Amsterdam, 1978. Search in Google Scholar

[61] M. F. Wheeler, An elliptic collocation-finite element method with interior penalties, SIAM J. Numer. Anal. 15 (1978), no. 1, 152–161. 10.1137/0715010Search in Google Scholar

Received: 2022-07-21

Accepted: 2022-09-26

Published Online: 2023-02-15

Published in Print: 2023-07-01

A Symmetric Interior Penalty Method for an Elliptic Distributed Optimal Control Problem with Pointwise State Constraints

Abstract

1 Introduction

2 Continuous and Discrete Problems

2.1 The Continuous Problem

2.2 The Discrete Problem

3 Some Technical Tools

3.1 Mesh-Dependent Norms

Proof

3.2 Triangulations

3.3 The Interpolation Operator Π h

Proof

3.4 The Ritz Operator R h

3.5 Interior Estimates

3.6 The Connection Operator C h

3.7 The Smoothing Operator E h

3.8 The Operator I h

3.9 Estimates for y ¯

Proof

4 An Abstract Error Estimate

Proof

An Improved Abstract Error Estimate

5 Concrete Error Estimates

Proof

Proof

Proof

Proof

Proof

Improved Error Estimates

6 Numerical Results

6.1 Square Example

6.2 L-Shaped Domain Example

6.3 Cube Example

7 Concluding Remarks

A An Optimal Control Problem Without State Constraints

References

Journal and Issue

Articles in the same Issue