Numerical Analysis for the Pure Neumann Control Problem Using the Gradient Discretisation Method

Jérome Droniou; Neela Nataraj; Devika Shylaja

doi:10.1515/cmam-2017-0054

Publicly Available Published by De Gruyter December 5, 2017

Numerical Analysis for the Pure Neumann Control Problem Using the Gradient Discretisation Method

Jérome Droniou , Neela Nataraj and Devika Shylaja

From the journal Computational Methods in Applied Mathematics

https://doi.org/10.1515/cmam-2017-0054

Abstract

The article discusses the gradient discretisation method (GDM) for distributed optimal control problems governed by diffusion equation with pure Neumann boundary condition. Using the GDM framework enables to develop an analysis that directly applies to a wide range of numerical schemes, from conforming and non-conforming finite elements, to mixed finite elements, to finite volumes and mimetic finite differences methods. Optimal order error estimates for state, adjoint and control variables for low-order schemes are derived under standard regularity assumptions. A novel projection relation between the optimal control and the adjoint variable allows the proof of a super-convergence result for post-processed control. Numerical experiments performed using a modified active set strategy algorithm for conforming, non-conforming and mimetic finite difference methods confirm the theoretical rates of convergence.

Keywords: Elliptic Equations; Optimal Control; Neumann Boundary Conditions; Error Estimates,Super-Convergence; GDM; Finite Elements; Mimetic Finite Differences

MSC 2010: 49J20; 49M25; 65N15; 65N30

1 Introduction

We consider the following distributed optimal control problem governed by the diffusion equation with Neumann boundary condition:

(1.1)

(1.1a)minu∈𝒰ad⁡J⁢(y,u)⁢ ⁢subject to

(1.1b)-div⁢(A⁢∇⁡y)=u+f ⁢in ⁢Ω,

(1.1c)A⁢∇⁡y⋅𝐧Ω=0 ⁢on ∂⁡Ω,⨍Ωy⁢(𝒙)⁢d𝒙=0.

Here, Ω⊂ℝd (d≤3) is a bounded domain with boundary ∂⁡Ω, 𝐧Ω is the outer unit normal to Ω and

⨍Ωy⁢(𝒙)⁢d𝒙:=1|Ω|⁢∫Ωy⁢(𝒙)⁢d𝒙

denotes the average value of the function y over Ω (here and in the rest of the paper, |E| is the Lebesgue measure of a set E⊂ℝd). The cost functional, dependent on the state variable y and the control variable u, is given by

(1.2)J⁢(y,u):=12⁢∥y-y¯d∥L2⁢(Ω)2+α2⁢∥u∥L2⁢(Ω)2

with α>0. The desired state variable y¯d∈L2⁢(Ω) is chosen to satisfy ⨍Ωy¯d⁢(𝒙)⁢d𝒙=0. The source term f∈L2⁢(Ω) also satisfies the zero average condition ⨍Ωf⁢(𝒙)⁢d𝒙=0. The diffusion matrix A:Ω→ℳd⁢(ℝ) is a measurable, bounded and uniformly elliptic matrix-valued function such that A⁢(𝒙) is symmetric for a.e. 𝒙∈Ω. Finally, the admissible set of controls 𝒰ad is the non-empty convex set defined by

(1.3)𝒰ad={u∈L2⁢(Ω):a≤u≤b⁢ and ⁢⨍Ωu⁢(𝒙)⁢d𝒙=0},

where a and b are constants in [-∞,+∞] with a<0<b (this condition is necessary to ensure that 𝒰ad is not empty or reduced to {0}).

In this paper, we discuss the discretisation of control and state by the gradient discretisation method (GDM). The GDM is a generic framework for the convergence analysis of numerical methods (finite elements, mixed finite elements, finite volume, mimetic finite difference methods, etc.) for linear and non-linear elliptic and parabolic diffusion equations (including degenerate equations), the Navier–Stokes equations, variational inequalities, Darcy flows in fractured media, etc. See for example [16, 17, 18, 1, 22], and the monograph [15] for a complete presentation of the GDM. The contributions of this article are summarised now. Here, we

establish basic error estimates that provide 𝒪⁢(h) convergence rate for all the three variables (control, state and adjoint) for low-order schemes under standard regularity assumptions for the pure Neumann problem, without reaction term. Given that the optimal control u¯ is approximated by piecewise constant functions, the convergence rates are optimal.
prove super-convergence result for post-processed optimal controls, state and adjoint variables.
establish a projection relation (see Lemma 3.10) between control and adjoint variables. This relation, which is non-standard since it has to account for the zero average constraints, is the key to prove the super-convergence result for all three variables.
design a modified active set strategy algorithm for GDM that is adapted to this non-standard projection relation.
discuss some numerical results that confirm the theoretical rates of convergence for conforming, non-conforming and mimetic finite difference methods.

The literature contains several contributions to numerical analysis for second-order distributed optimal control problems governed by diffusion equation with Dirichlet and Neumann boundary conditions (BCs) (we refer to [20, 30] and the references therein). For Dirichlet BCs, the super-convergence of post-processed controls for conforming finite element (FE) methods has been investigated in [30]. Recently, this result was extended to the GDM in [20]. Carrying out this analysis in the context of the GDM means that it readily applies to various schemes, including non-conforming ℙ1 finite elements and hybrid mimetic mixed schemes (HMM), which contains mixed-hybrid mimetic finite differences; and for these schemes, the analysis of [20] provides novel estimates. We refer to [8, 9, 10, 11, 12, 13] for an analysis of non-linear elliptic optimal control problems. Several works cover optimal control for second-order Neumann boundary value problems, albeit with an additional (linear or non-linear) reaction term which makes the state equation naturally well-posed, without zero average constraint, see [2, 3, 11, 26, 29]. To the best of our knowledge, the numerical analysis of pure Neumann control problems, without reaction term and thus with the integral constraint, is open even for finite element methods. Our results therefore seem to be new and, being established in the GDM framework, cover a range of numerical methods, including conforming Galerkin methods, non-conforming finite elements, and mimetic finite differences. Although done on the simple problem (1.1), the analysis uncovers some properties, such as the specific relation formula between the adjoint and control variables and a modified active set algorithm used to compute the solution of the numerical scheme. These have a wider application potential in optimisation problems involving an integral constraint. For example, in the model of [14, 31] describing the miscible displacement of one fluid by another in a porous medium, the pressure is subjected to an elliptic equation with homogeneous Neumann boundary conditions. In this equation, the source terms model the injection and production wells, and are typically the only quantities that engineers can adjust (to some extent). Hence, considering these source terms as controls of the pressure may lead to optimal control problems as in (1.1), with the exact same boundary conditions and integral constraints on the control terms.

The paper is organised as follows. Section 2 recalls the GDM for elliptic problems with Neumann BCs and the properties needed to prove its convergence. Some classical examples of GDM are presented in Section 2.2. Section 3 deals with the GDM for the optimal control problem (1.1). The basic error estimates and super-convergence results are presented in Sections 3.2 and 3.3. The first super-convergence result provides a nearly quadratic convergence rate for a post-processed control. Under an L∞ stability assumption of the GDM, which stems for most schemes from the quasi-uniformity of the mesh, the second super-convergence theorem establishes a full quadratic super-convergence rate. Discussions on post-processed controls and the projection relation between control and proper adjoint are presented in Section 3.4. The active set strategy is an algorithm to solve the non-linear Karush–Kuhn–Tucker (KKT) formulation of the optimal control problem [32]. Section 4.1 presents a modification of this algorithm that accounts for the zero average constraint on the control. This modified active set algorithm also automatically selects the proper discrete adjoint whose projection provides the discrete control variable. In Section 4.2, we present the results of some numerical experiments. The paper ends with an appendix, Section A, where we prove the results stated in Sections 3.2 and 3.3.

Before concluding this introduction, we discuss the optimality conditions for (1.1). For a given u∈𝒰ad, there exists a unique weak solution y⁢(u)∈H⋆1⁢(Ω):={w∈H1⁢(Ω):⨍Ωw⁢(𝒙)⁢d𝒙=0} of (1.1b)–(1.1c). That is, for u∈𝒰ad, there exists a unique y⁢(u)∈H⋆1⁢(Ω) such that for all w∈H⋆1⁢(Ω),

a⁢(y,w)=∫Ωu⁢w⁢d𝒙,

where

a⁢(ϕ,ψ)=∫ΩA⁢∇⁡ϕ⋅∇⁡ψ⁢d⁢𝒙 for all ϕ,ψ∈H1⁢(Ω).

The term y⁢(u) is the state associated with the control u.

In the following, the norm and scalar product in L2⁢(Ω) (or L2⁢(Ω)d for vector-valued functions) are denoted by ∥⋅∥ and (⋅,⋅). The convex control problem (1.1) has a unique solution (y¯,u¯)∈H⋆1⁢(Ω)×𝒰ad and there exists a co-state p¯∈H1⁢(Ω) such that the triplet (y¯,p¯,u¯)∈H⋆1⁢(Ω)×H1⁢(Ω)×𝒰ad satisfies the KKT optimality conditions [28, Chapter 2]:

(1.4)

(1.4a)a⁢(y¯,w)=(u¯+f,w) ⁢for all ⁢w∈H⋆1⁢(Ω),

(1.4b)a⁢(z,p¯)=(y¯-y¯d,z) ⁢for all ⁢z∈H1⁢(Ω),

(1.4c)(p¯+α⁢u¯,v-u¯)≥0 ⁢for all ⁢v∈𝒰ad.

Remark 1.1 (Choice of the Adjoint State).

Several co-states satisfy the optimality conditions (1.4), as p¯ is only determined up to an additive constant by (1.4). The same will be true for the discrete co-state, solution to a discrete version of these KKT equations. Establishing error estimates require the continuous and discrete co-states to have the same average. The usual choice is to fix this average as zero. However, for the control problem with pure Neumann conditions, this is not the best choice. Indeed, as seen in Lemma 3.10, establishing a proper relation between the control and co-state requires a certain zero average of a non-linear function of this co-state. A more efficient approach, that we will adopt, to fix the proper co-states is thus the following:

Design an algorithm (the modified active set algorithm of Section 4.1) that computes a discrete co-state with the proper condition, so that the discrete control can be easily obtained in terms of this discrete co-state.
Fix the average of the continuous co-state p¯ to be the same as the average of the discrete co-state obtained above.

As we will see, an algebraic relation between this p¯ and the continuous control u¯ can still be written, upon selecting a proper (but non-explicit) translation of p¯.

Remark 1.2 (Zero Average Constraint on the Source Term and Desired State).

(i) If we consider (1.1) without the constraint ⨍Ωf⁢(𝒙)⁢d𝒙=0 on the source term, the set of admissible controls needs to be modified into

𝒰ad={u∈L2⁢(Ω):a≤u≤b⁢ and ⁢⨍Ω(u+f)⁢d𝒙=0}.

In this case, a simple transformation can bring us back to the case of a source term with zero average. Rewrite the state equation (1.1b) as -div⁢(A⁢∇⁡y)=u⋆+f⋆ with u⋆=u+⨍Ωf⁢d𝒙 and f⋆=f-⨍Ωf⁢d𝒙. Then ⨍Ωf⋆⁢d𝒙=0 and u⋆∈𝒰ad⋆, where

𝒰ad⋆={u⋆∈L2⁢(Ω):a⋆≤u⋆≤b⋆⁢ and ⁢⨍Ωu⋆⁢d𝒙=0}

with

a⋆=a+⨍Ωf⁢d𝒙 and b⋆=b+⨍Ωf⁢d𝒙.

(ii) If the desired state y¯d∈L2⁢(Ω) is such that ⨍Ωy¯d(𝒙)d𝒙=:m≠0, then it is natural to select states y in (1.1) with the same average m (since the average of these states can be freely fixed, and the choice made in (1.1c) is arbitrary). This ensures the best possible approximation of the desired state y¯d. In that case, working with y-m and y¯d-m instead of y and y¯d brings back to the original formulation (1.1) with a desired state y¯d-m having a zero average.

Remark 1.3 (Non-homogeneous BCs).

The study of the second-order distributed control problem (1.1) with non-homogeneous boundary conditions A⁢∇⁡y⋅𝐧Ω=g on ∂⁡Ω (with g∈L2⁢(∂⁡Ω)) follows in a similar way. In this case, the source terms and boundary condition are supposed to satisfy the compatibility condition

∫Ωf⁢d𝒙+∫∂⁡Ωg⁢ds⁢(𝒙)=0.

The controls are still taken in 𝒰ad defined by (1.3) and the KKT optimality condition is [28]: Seek a function (y¯,p¯,u¯)∈H⋆1⁢(Ω)×H1⁢(Ω)×𝒰ad such that

a⁢(y¯,w)=(u¯+f,w)+(g,γ⁢(w))∂for all ⁢w∈H⋆1⁢(Ω),

a⁢(z,p¯)=(y¯-y¯d,z)for all ⁢z∈H1⁢(Ω),

(p¯+α⁢u¯,v-u¯)≥0for all ⁢v∈𝒰ad,

where γ:H1⁢(Ω)→L2⁢(∂⁡Ω) is the trace operator and (⋅,⋅)∂ is the inner product in L2⁢(∂⁡Ω).

2 GDM for Elliptic PDE with Neumann BCs

The GDM consists in writing numerical schemes, called gradient schemes (GS), by replacing the continuous space and operators by discrete ones in the weak formulation of the problem [15, 16, 21]. These discrete space and operators are given by a gradient discretisation (GD).

2.1 Gradient Discretisation and Gradient Scheme

A notion of gradient discretisation for Neumann BCs is given in [15, Definition 3.1]. The following extends this definition by demanding the existence of the element 1𝒟 and is always satisfied in practical applications. This existence ensures that the zero average condition can be put in the discretisation space or in the bilinear form as for the continuous formulation, see Remark 2.2.

Definition 2.1 (Gradient Discretisation for Neumann BCs).

A gradient discretisation (GD) for homogeneous Neumann BCs is given by 𝒟=(X𝒟,Π𝒟,∇𝒟) such that:

X𝒟 is a finite-dimensional vector space on ℝ.
Π𝒟:X𝒟→L2⁢(Ω) and ∇𝒟:X𝒟→L2⁢(Ω)d are linear mappings.
The quantity
(2.1)∥w∥𝒟2:=∥∇𝒟w∥2+|⨍ΩΠ𝒟w(𝒙)d𝒙|2
is a norm on X𝒟.
There exists 1𝒟∈X𝒟 such that Π𝒟⁢1𝒟=1 on Ω and ∇𝒟⁡1𝒟=0 on Ω.

The flexibility of the GDM analysis framework comes from the wide possible range of choices for the gradient discretisation (X𝒟,Π𝒟,∇𝒟). Each of these choices correspond to a particular numerical scheme (see Section 2.2 for a few examples). The space X𝒟 represents the degrees of freedom (unknowns) of the method; a vector in X𝒟 gathers values for such unknowns. The operators Π𝒟 and ∇𝒟 reconstruct, from a set of values of these unknowns, a scalar (respectively, vector) function on the entire set Ω. The scalar function is suppose to play the role of the solution/test functions itself in the weak formulation of the PDE; the vector-function, reconstructed “gradient”, is used in lieu of the gradients of these solution/test functions. Performing these substitutions in the weak formulation leads to a finite-dimensional system of equations (on the unknowns), which is dubbed the gradient scheme corresponding to the gradient discretisation 𝒟.

If F∈L2⁢(Ω) is such that ⨍ΩF⁢(𝒙)⁢d𝒙=0, the weak formulation of the Neumann boundary value problem

(2.2){-div⁢(A⁢∇⁡ψ)=Fin Ω,A⁢∇⁡ψ⋅𝐧Ω=0on ∂⁡Ω,

is given as follows:

(2.3)Find ψ∈H⋆1⁢(Ω) such that, for all w∈H⋆1⁢(Ω), ⁢a⁢(ψ,w)=∫ΩF⁢w⁢d𝒙.

As explained above, a gradient scheme for (2.2) is then obtained from a GD 𝒟 by writing the weak formulation (2.3) with the continuous spaces, functions and gradients replaced with their discrete counterparts:

(2.4)Find ψ𝒟∈X𝒟,⋆ such that, for all w𝒟∈X𝒟,⋆, ⁢a𝒟⁢(ψ𝒟,w𝒟)=∫ΩF⁢Π𝒟⁢w𝒟⁢d𝒙,

where

a𝒟⁢(ϕ𝒟,z𝒟)=∫ΩA⁢∇𝒟⁡ϕ𝒟⋅∇𝒟⁡z𝒟⁢d⁢𝒙 for all ϕ𝒟,z𝒟∈X𝒟,

and

X𝒟,⋆={w𝒟∈X𝒟:⨍ΩΠ𝒟⁢w𝒟⁢d𝒙=0}.

Remark 2.2.

As for the continuous formulation (2.3), using the element 1𝒟∈X𝒟 actually enables us to consider in (2.4) test functions w𝒟 in X𝒟, rather than just X𝒟,⋆. The simplest technique to achieve this is to use a quadratic penalty method [24, Chapter 11]. For any ρ>0, (2.4) can be shown equivalent to

(2.5)

Find ψ𝒟∈X𝒟 such that, for all w𝒟∈X𝒟,

a𝒟⁢(ψ𝒟,w𝒟)+ρ⁢(⨍ΩΠ𝒟⁢ψ𝒟⁢d𝒙)⁢(⨍ΩΠ𝒟⁢w𝒟⁢d𝒙)=∫ΩF⁢Π𝒟⁢w𝒟⁢d𝒙.

Indeed, considering w𝒟=1𝒟 in (2.5) shows that the solution to this problem belongs to X𝒟,⋆ and is therefore a solution to (2.4). The converse is straightforward.

2.2 Examples of Gradient Discretisations

We briefly present here a few examples based on known numerical methods. We refer to [15] for a detailed analysis of these methods, and more examples of gradient discretisations. As demonstrated by these examples, the GDM cover a wide range of different numerical methods. This means that the analysis carried out in the GDM framework for the control problem in (1.1) readily applies to all these methods. In particular, this makes the control problem accessible to numerical schemes not usually considered but relevant to diffusion models, such as schemes applicable on generic meshes (not just triangular/quadrangular meshes) as encountered for example in reservoir engineering applications.

Let us consider a mesh 𝒯 of Ω. A precise definition can be found for example in [15, Definition 7.2] but, for our purpose, the intuitive understanding of 𝒯 as a partition of Ω in polygonal/polyhedral sets is sufficient.

Conforming ℙ1 Finite Elements.

The simplest gradient discretisation is perhaps obtained by considering conforming ℙ1 finite elements. The mesh is made of triangles (in 2D) or tetrahedra (in 3D), with no hanging nodes. Each v𝒟∈X𝒟 is a vector of values at the vertices of the mesh (the standard unknowns of conforming ℙ1 finite elements), Π𝒟⁢v𝒟 is the continuous piecewise linear function on the mesh which takes these values at the vertices, and ∇𝒟⁡v𝒟=∇⁡(Π𝒟⁢v𝒟). Then (2.4) is the standard ℙ1 finite element scheme for (2.2).

Non-conforming ℙ1 Finite Elements.

As above, the mesh is made of conforming triangles or tetrahedra. Each v𝒟∈X𝒟 is a vector of values at the centers of mass of the edges/faces, Π𝒟⁢v𝒟 is the piecewise linear function on the mesh which takes these values at these centers of mass, and ∇𝒟⁡v𝒟=∇𝒯⁡(Π𝒟⁢v𝒟) is the broken gradient of Π𝒟⁢v𝒟. In that case, (2.4) gives the non-conforming ℙ1 finite element approximation of (2.2).

Mass-Lumped Non-conforming ℙ1 Finite Elements.

Still considering a conforming triangular/tetrahedral mesh, the space X𝒟 and gradient reconstruction ∇𝒟 are identical to those of the non-conforming ℙ1 finite elements described above, but the function reconstruction is modified to be piecewise constant. For each edge/face σ of the mesh, we consider the diamond Dσ around σ constructed from the edge/face and the one or two cell centers on each side (see Figure 1). Then, for v=(vσ)σ∈X𝒟, the reconstructed function Π𝒟⁢v is the piecewise constant function on the diamonds, equal to vσ on Dσ for all edge/face σ.

$Figure 1 Diamonds for the definition of the mass-lumped non-conforming ℙ1{\mathbb{P}_{1}} finite elements.$

Figure 1

Diamonds for the definition of the mass-lumped non-conforming ℙ1 finite elements.

Hybrid Mixed Mimetic Method (HMM).

We consider a generic mesh 𝒯 (not necessarily triangular/tetrahedral) with one point 𝒙K chosen in each cell K∈𝒯 such that K is strictly star-shaped with respect to 𝒙K; see Figure 2 for some notations. A vector v∈X𝒟 is made of cell (vK)K and face (vσ)σ values, and the operator Π𝒟 reconstructs a piecewise constant function from the cell values: for any cell K, Π𝒟⁢v=vK on K. The gradient reconstruction ∇𝒟 is built in two pieces: a consistent gradient ∇¯K constant over the cell and stabilisation terms constant over the half-diamonds DK,σ (and akin to the remainders of first-order Taylor expansions between the cell and face values). For any cell K and any face σ of K, we set

∇𝒟⁡v=∇¯K⁢v+ddK,σ⁢(vσ-vK-∇¯K⁢v⋅(𝒙¯σ-𝒙K))⁢𝐧K,σ on K,

where dK,σ is the orthogonal distance between 𝒙K and σ, 𝒙¯σ is the center of mass of σ, 𝐧K,σ is the outer normal to K on σ and, denoting by ℰK the set of faces of K,

∇¯K⁢v=1|K|⁢∑σ∈ℰK|σ|⁢vσ⁢𝐧K,σ.

Once used in the gradient scheme (2.4), this HMM gradient discretisation gives rise to a numerical method that can be applied on any mesh (including with hanging nodes, non-convex cells, etc.). This scheme can also be re-interpreted as a finite volume method [15, Section 13.3].

Figure 2

Notations for the construction of the HMM gradient discretisation.

2.3 Error Estimates for the GDM for the Neumann Problem

The accuracy of a gradient scheme (2.4) is measured by three quantities. The first one is a discrete Poincaré constant C𝒟, which ensures the coercivity of the method:

(2.6)C𝒟:=maxw∈X𝒟∖{0}⁡∥Π𝒟⁢w∥∥w∥𝒟.

The second quantity is the interpolation error S𝒟, which measures what is called, in the GDM framework, the GD-consistency of 𝒟:

(2.7)S𝒟⁢(φ)=minw∈X𝒟⁡(∥Π𝒟⁢w-φ∥+∥∇𝒟⁡w-∇⁡φ∥) for all ⁢φ∈H1⁢(Ω).

Finally, we measure the limit-conformity of a GD by defining

(2.8)

W𝒟(𝝋)=maxw∈X𝒟∖{0}1∥w∥𝒟|∫Ω(Π𝒟wdiv(𝝋)+∇𝒟w⋅𝝋)d𝒙| for all 𝝋∈Hdiv0(Ω),

where H0div⁢(Ω)={𝝋∈L2⁢(Ω)d:div⁢(𝝋)∈L2⁢(Ω),γ𝐧⁢(𝝋)=0} with γ𝐧 being the normal trace of 𝝋 on ∂⁡Ω.

By using these quantities, an error estimate can be established for GS. We refer to [15] for a proof of the following theorem. Here and in the rest of the paper,

(2.9)X≲Y⁢ means that ⁢X≤C⁢Y⁢ for some C depending only on Ω, A and an upper bound of C𝒟.

Theorem 2.3 (Error Estimate for the GDM).

Let D be a GD in the sense of Definition 2.1, let ψ be the solution to (2.3), and let ψD be the solution to (2.4). Then

(2.10)∥Π𝒟⁢ψ𝒟-ψ∥+∥∇𝒟⁡ψ𝒟-∇⁡ψ∥≲WS𝒟⁢(ψ),

where

(2.11)WS𝒟⁢(ψ)=W𝒟⁢(A⁢∇⁡ψ)+S𝒟⁢(ψ).

Remark 2.4 (Rates of Convergence).

For all classical low-order methods based on meshes (such as ℙ1 conforming and non-conforming finite element methods, finite volume methods, etc.), if A is Lipschitz continuous and ψ∈H2⁢(Ω), then 𝒪⁢(h) estimates can be obtained for W𝒟⁢(A⁢∇⁡ψ) and S𝒟⁢(ψ) (see [15]). Theorem 2.3 then gives a linear rate of convergence for these methods.

Remark 2.5.

Note that Theorem 2.3 also holds if we replace the zero average condition on ψ and Π𝒟⁢ψ𝒟 with

⨍ΩΠ𝒟⁢ψ𝒟⁢d𝒙=⨍Ωψ⁢d𝒙.

In this case, estimate (2.10) can be obtained by considering the translation of ψ𝒟 and ψ. Set ψ~𝒟=ψ𝒟-c⁢1𝒟 and ψ~=ψ-c⁢1, where

c=⨍ΩΠ𝒟⁢ψ𝒟⁢d𝒙=⨍Ωψ⁢d𝒙

and 1 is the constant function. Using Definition 2.1, we find Π𝒟⁢ψ~𝒟=Π𝒟⁢ψ𝒟-c, ∇𝒟⁡ψ~𝒟=∇𝒟⁡ψ𝒟 and ∇⁡ψ~=∇⁡ψ. This gives

⨍ΩΠ𝒟⁢ψ~𝒟⁢d𝒙=⨍Ωψ~⁢d𝒙=0.

Applying Theorem 2.3,

∥Π𝒟⁢ψ~𝒟-ψ~∥+∥∇𝒟⁡ψ~𝒟-∇⁡ψ~∥≲WS𝒟⁢(ψ~),

which implies

∥Π𝒟⁢ψ𝒟-ψ∥+∥∇𝒟⁡ψ𝒟-∇⁡ψ∥≲WS𝒟⁢(ψ~)=WS𝒟⁢(ψ).

The following stability result, useful to our analysis, is straightforward.

Proposition 2.6 (Stability of the GDM).

Let a¯ be a coercivity constant of A. If ψD is the solution to the gradient scheme (2.4), then

∥∇𝒟⁡ψ𝒟∥≤C𝒟a¯⁢∥F∥ 𝑎𝑛𝑑 ∥Π𝒟⁢ψ𝒟∥≤C𝒟2a¯⁢∥F∥.

Proof.

Take w𝒟=ψ𝒟 in (2.4) and use the definition of C𝒟 to write

a¯⁢∥∇𝒟⁡ψ𝒟∥2≤∥F∥⁢∥Π𝒟⁢ψ𝒟∥≤C𝒟⁢∥F∥⁢∥ψ𝒟∥𝒟.

Since ⨍ΩΠ𝒟⁢ψ𝒟⁢d𝒙=0, recalling the definition (2.1) of ∥⋅∥𝒟 shows that ∥ψ𝒟∥𝒟=∥∇𝒟⁡ψ𝒟∥ and the proof of first estimate is complete. The second estimate follows from the definition of C𝒟. ∎

3 GDM for the Control Problem and Main Results

This section starts with a description of GDM for the optimal control problem and is followed by the basic error estimates and super-convergence results in Sections 3.2 and 3.3. The super-convergence results for a post-processed control are presented. A discussion on post-processed controls and the projection relation between control and proper adjoint are presented in Section 3.4.

3.1 GDM for the Optimal Control Problem

Let 𝒟 be a GD as in Definition 2.1. The space 𝒰h is defined as the space of piecewise constant functions on a mesh 𝒯 of Ω. The space 𝒰ad,h=𝒰ad∩𝒰h is a finite-dimensional subset of 𝒰ad. A gradient scheme for (1.4) consists in seeking (y¯𝒟,p¯𝒟,u¯h)∈X𝒟,⋆×X𝒟×𝒰ad,h such that

(3.1)

(3.1a)a𝒟⁢(y¯𝒟,w𝒟)=(u¯h+f,Π𝒟⁢w𝒟) ⁢for all ⁢w𝒟∈X𝒟,⋆,

(3.1b)a𝒟⁢(z𝒟,p¯𝒟)=(Π𝒟⁢y¯𝒟-y¯d,Π𝒟⁢z𝒟) ⁢for all ⁢z𝒟∈X𝒟,

(3.1c)(Π𝒟⁢p¯𝒟+α⁢u¯h,vh-u¯h)≥0 ⁢for all ⁢vh∈𝒰ad,h.

As in the continuous KKT conditions (1.4), these equations do not define p¯𝒟 uniquely. One possible constraint that fixes a unique p¯𝒟 is described in Lemma 3.10. This particular choice ensures a simple projection relation between p¯𝒟 and u¯h.

Two projection operators play a major role throughout the paper: the orthogonal projection on piecewise constant functions on 𝒯, namely 𝓟𝒯:L1⁢(Ω)→𝒰h and the cut-off function P[a,b]:ℝ→[a,b]. They are defined as

(3.2)(𝓟𝒯⁢v)|K:=⨍Kv⁢(𝒙)⁢d𝒙 ⁢for all ⁢v∈L1⁢(Ω)⁢ and all ⁢K∈𝒯,

(3.3)P[a,b]⁢(s):=min⁡(b,max⁡(a,s)) ⁢for all ⁢s∈ℝ.

3.2 Basic Error Estimate for the GDM for the Control Problem

To state the error estimates, we define the projection error Eh by

(3.4)Eh⁢(W)=∥W-𝓟𝒯⁢W∥ for all ⁢W∈L2⁢(Ω).

The proofs of the following basic error estimates are provided in Section A. They follow by adapting the corresponding proofs in [20] to account for the pure Neumann boundary conditions and integral constraints.

Theorem 3.1 (Control Estimate).

Let D be a GD, let (y¯,p¯,u¯) be a solution to (1.4) and let (y¯D,p¯D,u¯h) be a solution to (3.1) such that ⨍ΩΠD⁢p¯D⁢d𝐱=⨍Ωp¯⁢d𝐱. Then, recalling (2.9), (2.11) and (3.4), there exists a constant C depends only on α such that

(3.5)

∥u¯-u¯h∥≲C⁢(Eh⁢(p¯)+Eh⁢(u¯)+WS𝒟⁢(p¯)+WS𝒟⁢(y¯)).

Proposition 3.2 (State and Adjoint Error Estimates).

Let D be a GD, let (y¯,p¯,u¯) be a solution to (1.4) and let (y¯D,p¯D,u¯h) be a solution to (3.1). Assume that ⨍ΩΠD⁢p¯D⁢d𝐱=⨍Ωp¯⁢d𝐱. Then the following error estimates hold:

(3.6)∥Π𝒟⁢y¯𝒟-y¯∥+∥∇𝒟⁡y¯𝒟-∇⁡y¯∥≲∥u¯-u¯h∥+WS𝒟⁢(y¯),

(3.7)∥Π𝒟⁢p¯𝒟-p¯∥+∥∇𝒟⁡p¯𝒟-∇⁡p¯∥≲∥u¯-u¯h∥+WS𝒟⁢(y¯)+WS𝒟⁢(p¯).

Remark 3.3 (Rates of Convergence for the Control Problem).

As in Remark 2.4, if A is Lipschitz continuous and (y¯,p¯,u¯)∈H2⁢(Ω)2×H1⁢(Ω), then (3.5), (3.6) and (3.7) give linear rates of convergence for all classical first-order methods.

3.3 Super-Convergence for Post-Processed Controls

In this subsection, we define the post-processed continuous and discrete controls (see (3.11)) and state the super-convergence results. The proofs are presented in Section A.

We make here the following assumptions.

Assumption A1 (Interpolation Operator).

For each w∈H2⁢(Ω), there exists w𝒯∈L2⁢(Ω) such that:

If w∈H2⁢(Ω) solves -div⁢(A⁢∇⁡w)=g∈H1⁢(Ω), and w𝒟 is the solution to the corresponding GS with
⨍ΩΠ𝒟⁢w𝒟⁢d𝒙=⨍Ωw⁢d𝒙,
then
(3.8)∥Π𝒟⁢w𝒟-w𝒯∥≲h2⁢∥g∥H1⁢(Ω).
For any w∈H2⁢(Ω), it holds
(3.9)
|(w-w𝒯,Π𝒟⁢v𝒟)|≲h2⁢∥w∥H2⁢(Ω)⁢∥Π𝒟⁢v𝒟∥ for all ⁢v𝒟∈X𝒟
and
(3.10)∥𝓟𝒯⁢(w-w𝒯)∥≲h2⁢∥w∥H2⁢(Ω).

Assumption A2.

The estimate ∥Π𝒟⁢v𝒟-𝓟𝒯⁢(Π𝒟⁢v𝒟)∥≲h⁢∥∇𝒟⁡v𝒟∥ holds for any v𝒟∈X𝒟.

Assumption A3 (Discrete Sobolev Imbedding).

For all v𝒟∈X𝒟, it holds

∥Π𝒟⁢v𝒟∥L2*⁢(Ω)≲∥v𝒟∥𝒟,

where 2* is a Sobolev exponent of 2, that is, 2*∈[2,∞) if d=2, and 2*=6 if d=3.

Let

𝒯2={K∈𝒯:u¯=a a.e. on K, or u¯=b a.e. on K, or a<u¯<b a.e. on K},

and 𝒯1=𝒯∖𝒯2. That is, 𝒯1 is the set of cells where u¯ crosses at least one constraint a or b. For i=1,2, we let Ωi,𝒯=int⁢(⋃K∈𝒯iK¯). The space W1,∞⁢(𝒯1) is the usual broken Sobolev space, endowed with its broken norm. Our last assumption is:

Assumption A4.

We have |Ω1,𝒯|≲h and u¯|Ω1,𝒯∈W1,∞⁢(𝒯1).

Possible choices of the mapping w↦w𝒯, depending on the considered gradient discretisation (that is, the considered numerical method), is discussed in Remark 3.4 below. Note that Assumptions A1–A4 are similar to that in [20] with X𝒟,0 substituted by X𝒟, and an additional average condition in Assumption A1. We also refer to [20] for a detailed discussion on Assumptions A1–A4.

Assuming p¯∈H2⁢(Ω) (see Theorem 3.5) and letting p¯𝒯 be defined as in Assumption A1, the post-processed continuous and discrete controls are given by

(3.11)

u~⁢(𝒙)=P[a,b]⁢(-1α⁢p¯𝒯⁢(𝒙)) and u~h⁢(𝒙)=P[a,b]⁢(-1α⁢Π𝒟⁢p¯𝒟⁢(𝒙)).

For a detailed discussion on the post-processed controls, we refer the reader to Section 3.4.

Remark 3.4 (Choice of wT in Assumption A1).

When Π𝒟 is a piecewise linear reconstruction, the super-convergence result (3.8) usually holds with w𝒯=w. This is for example well known for conforming and non-conforming ℙ1 FE method. When Π𝒟⁢v𝒟 is piecewise constant on 𝒯 for all v𝒟∈X𝒟, the super-convergence (3.8) requires to project the exact solution on piecewise constant functions on the mesh. This is usually done by setting w𝒯⁢(𝒙)=⨍Kw⁢(𝒙)⁢d𝒙 for all 𝒙∈K and all K∈𝒯 (or, equivalently at order 𝒪⁢(h2), w𝒯⁢(𝒙)=w⁢(𝒙¯K) with 𝒙¯K the center of mass of K). With this choice, the super-convergence result is known, e.g., for mixed/hybrid and nodal mimetic finite difference schemes (see [4, 19]). As a consequence:

For FE methods, u~=P[a,b]⁢(-α-1⁢p¯).
For mimetic finite difference methods, u~|K=P[a,b]⁢(-α-1⁢p¯⁢(𝒙¯K)) for all K∈𝒯.

For K∈𝒯 and 𝒙¯K the center of mass of K, let ρK=max⁡{r>0:B⁢(𝒙¯K,r)⊂K} be the maximal radius of balls centred at 𝒙¯K and included in K. Fixing η>0 such that

η≥diam⁢(K)ρK for all ⁢K∈𝒯,

we use the following extension of the notation (2.9):

X≲ηY⁢ means that ⁢X≤C⁢Y⁢ for some C depending only on Ω, A, an upper bound of C𝒟, and η.

Theorem 3.5 (Super-Convergence for Post-Processed Controls I).

Let D be a GD and let T be a mesh. Assume that

Assumptions A1–A4 hold,
y¯ and p¯ belong to H2⁢(Ω),
y¯d and f belong to H1⁢(Ω),

and let u~, u~h be the post-processed controls defined by (3.11), where p¯ and p¯D are chosen such that

⨍ΩΠ𝒟⁢p¯𝒟⁢d𝒙=⨍Ωp¯⁢d𝒙.

Then there exists C depending only on α in (1.2) such that

∥u~-u~h∥≲ηC⁢h2-12*⁢∥u¯∥W1,∞⁢(𝒯1)+C⁢h2⁢ℱ⁢(y¯d,f,y¯,p¯)

and

ℱ⁢(y¯d,f,y¯,p¯)=∥y¯d∥H1⁢(Ω)+∥f∥H1⁢(Ω)+∥y¯∥H2⁢(Ω)+∥p¯∥H2⁢(Ω).

Theorem 3.6 (Super-Convergence for Post-Processed Controls II).

Let the assumptions and notations of Theorem 3.5 hold, except Assumption A3 which is replaced by

(3.12)there exists δ>0 such that, for any F∈L2⁢(Ω), the solution ψD to (2.4) satisfies ∥ΠD⁢ψD∥L∞⁢(Ω)≤δ⁢∥F∥.

Then there exists C depending only on α and δ such that

∥u~-u~h∥≲ηC⁢h2⁢[∥u¯∥W1,∞⁢(𝒯1)+ℱ⁢(y¯d,f,y¯,p¯)].

Remark 3.7.

For most methods, assumption (3.12) is satisfied if the mesh is quasi-uniform (see [23] for conforming and non-conforming ℙ1 finite element method, and [20, Theorem 7.1] for HMM methods).

Corollary 3.8 (Super-Convergence for the State and Adjoint Variables).

Let the assumptions of Theorem 3.5 be fulfilled. Then the following error estimates hold, with C depending only on α:

(3.13)∥y¯𝒯-Π𝒟⁢y¯𝒟∥≲ηC⁢hr⁢∥u¯∥W1,∞⁢(𝒯1)+C⁢h2⁢ℱ⁢(y¯d,f,y¯,p¯),

(3.14)∥p¯𝒯-Π𝒟⁢p¯𝒟∥≲ηC⁢hr⁢∥u¯∥W1,∞⁢(𝒯1)+C⁢h2⁢ℱ⁢(y¯d,f,y¯,p¯),

where y¯T and p¯T are defined as in Assumption A1, and r=2-12*.

Under the assumptions of Theorem 3.6, (3.13) and (3.14) hold with r=2 and C depending only α and δ.

3.4 Discussion on Post-Processed Controls

In this subsection, we present a detailed analysis of post-processed controls given by (3.11). This analysis is performed under the assumptions of Section 3.3, and by also assuming that WS𝒟⁢(φ)≲h for all φ∈H2⁢(Ω) (see Remark 2.4). We begin by stating and proving two lemmas which discuss projection relations between control and adjoint variables for the pure Neumann problem, both at the continuous level and at the discrete level. We then show that the post-processed controls remain O⁢(h) close to their corresponding original controls, see (3.20) and (3.24). Hence, the super-convergence result makes sense: since u¯h is piecewise constant, it is impossible to expect more than O⁢(h) approximation on the controls; but by “moving” these controls by a specific O⁢(h), we obtain computable post-processed controls that enjoy an O⁢(h2) convergence result.

Lemma 3.9.

Let -∞≤a<0<b≤∞ and ϕ∈L1⁢(Ω). Define Γ:R→R by

Γ⁢(c)=∫ΩP[a,b]⁢(ϕ+c)⁢d𝒙,

where P[a,b] is given by (3.3). Set m=a-ess sup⁢(ϕ)∈[-∞,+∞) and M=b-ess inf⁢(ϕ)∈(-∞,+∞]. Then we have the following:

Γ is Lipschitz continuous.
limc→m⁡Γ⁢(c)=a⁢|Ω|, limc→M⁡Γ⁢(c)=b⁢|Ω|, and there is c⋆∈(m,M) such that Γ⁢(c⋆)=0.
If ϕ∈H1⁢(Ω), then for any compact interval Q in (m,M), there exists ρQ>0 such that if c,c′∈Q with c<c′, then
(3.15)Γ⁢(c′)-Γ⁢(c)≥ρQ⁢(c′-c).
As a consequence, the real number c⋆ in item 2 is unique.

Proof.

Item 1 is obvious since P[a,b] is Lipschitz continuous.

Let us now analyse the limits in item 2. Let (cn) be a sequence in ℝ such that cn→M as n→∞. By the definition of M, this implies P[a,b]⁢(ϕ+cn)→b a.e. on Ω. Let (cn) be bounded below by R and note that ϕ+R∈L1⁢(Ω). Moreover, for s∈ℝ, a<0<b implies P[a,b]⁢(s)≥min⁡(s,0) so

P[a,b]⁢(ϕ+cn)≥min⁡(ϕ+cn,0)≥min⁡(ϕ+R,0)∈L1⁢(Ω).

By Fatou’s lemma,

∫Ωb⁢d𝒙≤lim infn→∞⁡∫ΩP[a,b]⁢(ϕ⁢(𝒙)+cn)⁢d𝒙,

which gives b⁢|Ω|≤lim infn→∞⁡Γ⁢(cn). Since Γ⁢(cn)≤b⁢|Ω| (because P[a,b]⁢(s)≤b), we infer that

limn→∞⁡Γ⁢(cn)=b⁢|Ω|,

and thus that limc→M⁡Γ⁢(c)=b⁢|Ω|. In a similar way, we deduce that limc→m⁡Γ⁢(c)=a⁢|Ω|.

The existence of c⋆ such that Γ⁢(c⋆)=0 then follows from the Intermediate Value Theorem and

limc→m⁡Γ⁢(c)=a⁢|Ω|⁢<0⁢<b|⁢Ω|=limc→M⁡Γ⁢(c).

We now assume that ϕ∈H1⁢(Ω) and we turn to item 3. For a.e. c∈ℝ, Γ′⁢(c)=∫Ω𝟙(a,b)⁢(ϕ⁢(𝒙)+c)⁢d𝒙, where 𝟙(a,b) is the characteristic function of (a,b). Define Θ⁢(c)=∫Ω𝟙(a,b)⁢(ϕ⁢(𝒙)+c)⁢d𝒙 for all c∈ℝ. We claim that

Θ is lower semi-continuous,
Θ⁢(c)>0 for all c∈(m,M).

To prove that Θ is lower semi-continuous, let cn→c as n→∞. Since 𝟙(a,b) is lower semi-continuous on ℝ, we have, for all 𝒙∈Ω,

𝟙(a,b)⁢(ϕ⁢(𝒙)+c)≤lim infn→∞⁡𝟙(a,b)⁢(ϕ⁢(𝒙)+cn).

Applying Fatou’s lemma,

Θ⁢(c)≤lim infn→∞⁡∫Ω𝟙(a,b)⁢(ϕ⁢(𝒙)+cn)⁢d𝒙=lim infn→∞⁡Θ⁢(cn).

It follows that Θ is lower semi-continuous. We now show that Θ>0 on the interval (m,M). Let c∈(m,M). Then I=(a-c,b-c)∩(ess⁢inf⁡ϕ,ess⁢sup⁡ϕ) is an interval of positive length, since a-c<ess⁢sup⁡ϕ and b-c>ess⁢inf⁡ϕ. The set WI,c={𝒙:ϕ⁢(𝒙)∈I} has a non-zero measure because ϕ∈H1⁢(Ω) and Ω is connected. To see this, let α<β be the endpoints of I and assume that ϕ∈H1⁢(Ω) takes some values less than α on a non-null set, some values greater than β on a non-null set, but that WI,c is a null set. Then P[α,β]⁢(ϕ)∈H1⁢(Ω) exactly takes the values α and β (outside a set of zero measure). Hence ∇⁡P[α,β]⁢(ϕ)=𝟙[α,β]⁢(ϕ)⁢∇⁡ϕ=0 and P[α,β]⁢(ϕ) should be constant, since Ω is connected, which is a contradiction. Thus, WI,c has a non-zero measure. Since {𝒙:ϕ⁢(𝒙)+c∈(a,b)}⊇WI,c, this gives Θ⁢(c)≥|WI,c|>0.

Coming back to item 3, let Q be a compact interval in (m,M). We know that Θ>0 on Q and Θ is lower semi-continuous. Hence Θ reaches its minimum on Q and infQ⁡Θ=Θ⁢(c0)>0 for some c0∈Q. Since Γ′=Θ a.e., we have Γ′≥infQ⁡Θ a.e. on Q and, Γ being Lipschitz and [c,c′]⊂Q, we infer

Γ⁢(c′)-Γ⁢(c)=∫cc′Γ′⁢(s)⁢ds≥[infQ⁡Θ]⁢(c′-c),

which establishes (3.15). The uniqueness of c⋆ such that Γ⁢(c⋆)=0 follows from this inequality, which shows that Γ is strictly increasing on (m,M). ∎

Lemma 3.10 (Projection Formulas for the Controls).

If p¯∈H1⁢(Ω) is a co-state and c¯∈R is such that

⨍ΩP[a,b]⁢(-1α⁢p¯⁢(𝒙)+c¯)⁢d𝒙=0,

then the continuous optimal control u¯ in (1.4) can be expressed in terms of the projection formula

(3.16)u¯⁢(𝒙)=P[a,b]⁢(-1α⁢p¯⁢(𝒙)+c¯).

If D is a GD and p¯D is chosen such that

(3.17)⨍ΩP[a,b]⁢(𝓟𝒯⁢(-1α⁢Π𝒟⁢p¯𝒟))⁢d𝒙=0,

then the discrete optimal control in (3.1) is given by

(3.18)u¯h⁢(𝒙)=P[a,b]⁢(𝓟𝒯⁢(-1α⁢Π𝒟⁢p¯𝒟⁢(𝒙))).

Proof.

Set p~=p¯-α⁢c¯. Clearly, p~∈H1⁢(Ω). From the optimality condition for the control problem (1.4c), we deduce that

(p~+α⁢u¯,v-u¯)≥0 for all ⁢v∈𝒰ad,

since ∫Ωu¯⁢d𝒙=∫Ωv⁢d𝒙=0. Set U=P[a,b]⁢(-α-1⁢p~), i.e.,

U={aon ⁢Ω+={𝒙∈Ω:p~⁢(𝒙)+α⁢U⁢(𝒙)>0},-α-1⁢p~on ⁢Ω0={𝒙∈Ω:p~⁢(𝒙)+α⁢U⁢(𝒙)=0},bon ⁢Ω-={𝒙∈Ω:p~⁢(𝒙)+α⁢U⁢(𝒙)<0}.

It is then straightforward to see that U∈𝒰ad, i.e., U∈[a,b] and ⨍ΩU⁢(𝒙)⁢𝑑𝒙=0 (by the choice of c¯). Then, using the definitions of Ω+, Ω0 and Ω-, since v≥a=U on Ω+ and v≤b=U on Ω-,

(p~+α⁢U,v-U)=∫Ω+(p~+α⁢U)⁢(v-U)⁢d𝒙+∫Ω0(p~+α⁢U)⁢(v-U)⁢d𝒙+∫Ω-(p~+α⁢U)⁢(v-U)⁢d𝒙≥0.

Recall that the optimality condition is nothing but a characterisation of the L2⁢(Ω) orthogonal projection of -α-1⁢p~ on 𝒰ad and, as such, defines a unique element u¯ of 𝒰ad. We just proved that U=P[a,b]⁢(-α-1⁢p~) satisfies this optimality condition, which shows that it is equal to u¯. The proof of (3.16) is complete.

The second relation follows in a similar way by noticing that, since controls are piecewise-constants on 𝒯, (3.1c) is equivalent to (𝓟𝒯⁢(Π𝒟⁢p¯𝒟)+α⁢u¯h,vh-u¯h)≥0 for all vh∈𝒰ad,h. We also notice that, by the definition of 𝓟𝒯 and assumption (3.17), P[a,b]⁢(𝓟𝒯⁢(-1α⁢Π𝒟⁢p¯𝒟))∈𝒰ad,h. ∎

Remark 3.11.

There is at least one adjoint p¯𝒟 such that (3.17) is satisfied: start from any adjoint p¯𝒟0 and, by applying Lemma 3.9 (item 2) to ϕ=𝓟𝒯⁢(-α-1⁢Π𝒟⁢p¯𝒟0) and by noticing that ϕ+c⋆=𝓟𝒯⁢(-α-1⁢Π𝒟⁢p¯𝒟0)+c⋆, find c⋆ such that p¯𝒟=p¯𝒟0-α⁢c⋆⁢1𝒟 satisfies (3.17).

Since the discrete co-state p¯𝒟 is a computable quantity, its average is easier to fix than the average of the non-computable p¯. Hence, the projection relation (3.18) is the most natural choice to express the discrete control u¯h in terms of the discrete adjoint variable. This is the choice made in the modified active set strategy presented in Section 4.1. Once this choice is made, since p¯ must have the same average as Π𝒟⁢p¯𝒟 for u~ defined in (3.11) to satisfy super-convergence estimates, it is clear that P[a,b]⁢(-1α⁢p¯) will not have a zero average in general. Hence, if we want to express the continuous control in terms of p¯, we need to offset this p¯ by the correct c¯, as stated in the lemma.

Lemma 3.12 (Stability of the Discrete States).

Let D be a GD, let (y¯,p¯,u¯) be a solution to (1.4) and let (y¯D,p¯D,u¯h) be a solution to (3.1). Assume that ⨍ΩΠD⁢p¯D⁢d𝐱=⨍Ωp¯⁢d𝐱. Then

(3.19)∥Π𝒟⁢y¯𝒟∥+∥∇𝒟⁡y¯𝒟∥+∥Π𝒟⁢p¯𝒟∥+∥∇𝒟⁡p¯𝒟∥≲1.

Proof.

Let ϕ∈H0div⁢(Ω). Taking w=0 in (2.7), we get S𝒟⁢(ϕ)≤∥ϕ∥+∥∇⁡ϕ∥. By the Cauchy–Schwarz inequality, using (2.6) and (2.1), for w∈X𝒟,

∫Ω(Π𝒟⁢w⁢div⁢(ϕ)+∇𝒟⁡w⋅ϕ)⁢d𝒙≤∥Π𝒟⁢w∥⁢∥div⁢(ϕ)∥+∥∇𝒟⁡w∥⁢∥ϕ∥≤C𝒟⁢∥w∥𝒟⁢(∥div⁢(ϕ)∥+∥ϕ∥).

With (2.8), this implies WS𝒟⁢(ϕ)≤C𝒟⁢(∥div⁢(ϕ)∥+∥ϕ∥).

Therefore, for A⁢∇⁡ψ∈H0div⁢(Ω), recalling the definition (2.11) of WS𝒟, we have

WS𝒟⁢(ψ)≲∥ψ∥H1⁢(Ω)+∥div⁢(A⁢∇⁡ψ)∥+∥A⁢∇⁡ψ∥≲∥ψ∥H1⁢(Ω)+∥div⁢(A⁢∇⁡ψ)∥≲1.

By using Proposition 3.2, Theorem 3.1, this shows that

∥∇𝒟⁡y¯𝒟∥+∥Π𝒟⁢y¯𝒟∥≲∥Π𝒟⁢y¯𝒟-y¯∥+∥∇𝒟⁡y¯𝒟-∇⁡y¯∥≲1.

The result for the adjoint variable can be derived similarly and hence (3.19) follows. ∎

In the rest of this section, we establish 𝒪⁢(h) estimates between the controls u¯, u¯h and their post-processed versions u~, u~h. These estimates justify that the post-processed controls are indeed meaningful quantities.

By using (3.11), (3.18), the Lipschitz continuity of P[a,b], Assumption A2 and Lemma 3.12, we have the following estimate between u¯h and u~h:

(3.20)∥u~h-u¯h∥≤α-1⁢∥Π𝒟⁢p¯𝒟-𝓟𝒯⁢(Π𝒟⁢p¯𝒟)∥≲α-1⁢h⁢∥∇𝒟⁡p¯𝒟∥≲α-1⁢h.

Let us now turn to estimating u¯-u~. The co-state p¯∈H1⁢(Ω) in (1.4) is still taken such that

⨍Ωp¯⁢d𝒙=⨍ΩΠ𝒟⁢p¯𝒟⁢d𝒙.

From Lemma 3.9, it follows that there exists a unique constant c¯∈(m,M) such that ∫ΩP[a,b]⁢(-1α⁢p¯+c¯)⁢d𝒙=0, where m and M are defined as in Lemma 3.9. Using Lemma 3.10 and recalling (3.11),

(3.21)u¯⁢(𝒙)=P[a,b]⁢(-1α⁢p¯⁢(𝒙)+c¯) and u~⁢(𝒙)=P[a,b]⁢(-1α⁢p¯𝒯⁢(𝒙)).

Starting from (3.21) and using the Lipschitz continuity of P[a,b], the assumption WS𝒟⁢(φ)≲h, Corollary 3.8 and Proposition 3.2, we get a constant C depending only on α, f, y¯d, p¯, y¯ and u¯ such that

∥u¯-u~∥≤α-1⁢∥p¯𝒯-p¯+α⁢c¯∥≲α-1⁢∥p¯𝒯-p¯∥+|c¯|

≲α-1⁢(∥p¯𝒯-Π𝒟⁢p¯𝒟∥+∥Π𝒟⁢p¯𝒟-p¯∥)+|c¯|

(3.22)≲ηC⁢h+|c¯|.

To estimate the last term in (3.22), recall the definition of Γ⁢(c) from Lemma 3.9 for ϕ=α-1⁢p¯:

Γ⁢(c)=∫ΩP[a,b]⁢(-1α⁢p¯+c)⁢d𝒙.

By the choice of c¯, Γ⁢(c¯)=0. From Lemma 3.10, the choice of p¯𝒟 shows that

Γ⁢(0)=∫ΩP[a,b]⁢(-1α⁢p¯)⁢d𝒙=∫Ω(P[a,b]⁢(-1α⁢p¯)-P[a,b]⁢(𝓟𝒯⁢(-1α⁢Π𝒟⁢p¯𝒟)))⁢d𝒙.

Let q𝒟 be the solution to (3.1b) with source term y¯-y¯d (that is, the solution to the GS for equation (1.4b) satisfied by p¯) such that ⨍ΩΠ𝒟⁢q𝒟⁢d𝒙=⨍ΩΠ𝒟⁢p¯𝒟⁢d𝒙=⨍Ωp¯⁢d𝒙. Using the Lipschitz continuity of P[a,b], the Cauchy–Schwarz inequality, the triangle inequality, Remark 2.5, Proposition 2.6, Assumption A2, Theorem 2.3 and Lemma 3.12, we obtain

|Γ(0)|≤∫Ω|P[a,b](-1αp¯)-P[a,b](𝓟𝒯(-1αΠ𝒟p¯𝒟))|d𝒙

≲α-1⁢∥p¯-𝓟𝒯⁢(Π𝒟⁢p¯𝒟)∥

≲α-1⁢∥p¯-Π𝒟⁢q𝒟∥+α-1⁢∥Π𝒟⁢q𝒟-Π𝒟⁢p¯𝒟∥+α-1⁢∥Π𝒟⁢p¯𝒟-𝓟𝒯⁢(Π𝒟⁢p¯𝒟)∥

≲α-1⁢(WS𝒟⁢(p¯)+α-1⁢∥y¯-Π𝒟⁢y¯𝒟∥+h⁢∥∇𝒟⁡p¯𝒟∥)

(3.23)≲α-1⁢(WS𝒟⁢(p¯)+WS𝒟⁢(y¯)+h⁢∥∇𝒟⁡p¯𝒟∥)≲α-1⁢h.

Let m,M be as in Lemma 3.9 for ϕ=α-1⁢p¯. Relation (3.23) shows that a⁢|Ω|⁢<Γ⁢(0)⁢<b|⁢Ω| if h is small enough; hence, in this case, 0∈(m,M). There is therefore a compact interval Q in (m,M) depending only on p¯ such that 0 and c¯ belong to Q. Without loss of generality, we can assume that c¯≥0. A use of Lemma 3.9 leads to Γ⁢(c¯)-Γ⁢(0)≥ρQ⁢c¯, where ρQ>0. This implies 0≤c¯≲α-1⁢hρQ, using (3.23) and the fact that Γ⁢(c¯)=0. Combining this with (3.22), we infer that

(3.24)∥u~-u¯∥≲η(C+α-1ρQ)⁢h.

4 Numerical Experiments

In this section, we first present the modified active set strategy. This is followed by results of numerical experiments for conforming, non-conforming and mimetic finite difference methods.

4.1 A Modified Active Set Strategy

The interest of choosing an adjoint given by (3.17) is highlighted in Lemma 3.10: we have the projection relation (3.18) between the discrete control and adjoint. Such a relation is at the core of the (standard) active set algorithm [32]. For a detailed analysis of this method, we refer the reader to [5, 6, 27]. Here, we propose a modified active set algorithm that enforces the proper zero average condition, and thus the proper relation between discrete adjoint and control.

We first notice that, when selecting the p¯𝒟 such that (3.17) holds, the KKT optimality conditions (3.1) can be rewritten as: Seek (y¯𝒟,p¯𝒟,u¯h)∈X𝒟×X𝒟×𝒰ad,h such that

(4.1)

(4.1a)a𝒟⁢(y¯𝒟,w𝒟)+ρ⁢(⨍ΩΠ𝒟⁢y¯𝒟⁢d𝒙)⁢(⨍ΩΠ𝒟⁢w𝒟⁢d𝒙)=(u¯h+f,Π𝒟⁢w𝒟),w𝒟∈X𝒟,

(4.1b)a𝒟⁢(z𝒟,p¯𝒟)+ρ⁢(⨍ΩP[a,b]⁢[𝓟𝒯⁢(-α-1⁢Π𝒟⁢p¯𝒟)]⁢d𝒙)⁢(⨍ΩΠ𝒟⁢z𝒟⁢d𝒙)=(Π𝒟⁢y¯𝒟-y¯d,Π𝒟⁢z𝒟),z𝒟∈X𝒟,

(4.1c)(Π𝒟⁢p¯𝒟+α⁢u¯h,vh-u¯h)≥0,vh∈𝒰ad,h,

where ρ>0 is constant.

Set μ¯h=-(α-1⁢Π𝒟⁢p¯𝒟+u¯h). As the original active set strategy [32], the modified active set strategy is an iterative algorithm. As initial guesses, two arbitrary functions, uh0,μh0 are chosen. In the nth step of the algorithm, we define the set of active and inactive restrictions by

Aa,hn⁢(𝒙)={𝒙:uhn-1⁢(𝒙)+μhn-1⁢(𝒙)<a},Ab,hn⁢(𝒙)={𝒙:uhn-1⁢(𝒙)+μhn-1⁢(𝒙)>b},Ihn=Ω∖(Aa,hn∪Ab,hn).

max⁡(∥uhn-uhn-1∥L∞⁢(Ω)∥uhn-1∥L∞⁢(Ω),∥Π𝒟⁢p𝒟n-Π𝒟⁢p𝒟n-1∥L∞⁢(Ω)∥Π𝒟⁢p𝒟n-1∥L∞⁢(Ω))≤10-10,

then we terminate the algorithm. In this case, we notice that the relative L∞ difference between Π𝒟⁢y𝒟n-1 and Π𝒟⁢y𝒟n is less than 10-6 for all examples in Section 4.2. Else we find y𝒟n,p𝒟n and uhn solution to the system

(4.2)

(4.2a)a𝒟⁢(y𝒟n,w𝒟)+ρ⁢(⨍ΩΠ𝒟⁢y𝒟n⁢d𝒙)⁢(⨍ΩΠ𝒟⁢w𝒟⁢d𝒙)=(uhn+f,Π𝒟⁢w𝒟),w𝒟∈X𝒟,

(4.2b)a𝒟⁢(z𝒟,p𝒟n)+ρ⁢(⨍ΩP[a,b]⁢[𝓟𝒯⁢(-α-1⁢Π𝒟⁢p𝒟n)]⁢d𝒙)⁢(⨍ΩΠ𝒟⁢z𝒟⁢d𝒙)=(Π𝒟⁢y𝒟n-y¯d,Π𝒟⁢z𝒟),z𝒟∈X𝒟,

(4.2c)uhn={aon ⁢Aa,hn,𝓟𝒯⁢(-α-1⁢Π𝒟⁢p𝒟n)on ⁢Ihn,bon ⁢Ab,hn.

The above algorithm consists of non-linear equations. It can however be approximated by a linearized system in the following way, thus leading to our final modified active set algorithm. Instead of solving (4.2), we solve

(4.3)

(4.3a)a𝒟⁢(y𝒟n,w𝒟)+ρ⁢(⨍ΩΠ𝒟⁢y𝒟n⁢d𝒙)⁢(⨍ΩΠ𝒟⁢w𝒟⁢d𝒙)=(uhn+f,Π𝒟⁢w𝒟),w𝒟∈X𝒟,

(4.3b)a𝒟⁢(z𝒟,p𝒟n)+ρ⁢(⨍ΩΠ𝒟⁢p𝒟n⁢d𝒙)⁢(⨍ΩΠ𝒟⁢z𝒟⁢d𝒙)=(Π𝒟⁢y𝒟n-y¯d,Π𝒟⁢z𝒟)+ρ⁢S𝒟n-1,z𝒟∈X𝒟,

(4.3c)uhn={aon ⁢Aa,hn,𝓟𝒯⁢(-α-1⁢Π𝒟⁢p𝒟n)on ⁢Ihn,bon ⁢Ab,hn,

where

S𝒟n-1=(⨍ΩΠ𝒟⁢z𝒟⁢d𝒙)⁢(⨍Ω{Π𝒟⁢p𝒟n-1-P[a,b]⁢[𝓟𝒯⁢(-α-1⁢Π𝒟⁢p𝒟n-1)]}⁢d𝒙).

Note that (4.3c) can be re-written in the following more commonly used form:

uhn+(1-𝟙a,hn-𝟙b,hn)⁢α-1⁢Π𝒟⁢p𝒟n=𝟙a,hn⁢a+𝟙b,hn⁢b,

where 𝟙a,hn and 𝟙b,hn denote the characteristic functions of the sets Aa,hn and Ab,hn, respectively.

Remark 4.1.

The convergence analysis of the proposed algorithm is a plan for future study. However, if (Π𝒟⁢y𝒟n,Π𝒟⁢p𝒟n) converges weakly to (Π𝒟⁢y¯𝒟,Π𝒟⁢p¯𝒟) in H1⁢(Ω) and uhn converges to u¯h in L2⁢(Ω), then the solution to (4.3) converges to the solution of (4.1) as n→∞.

4.2 Examples

In this subsection, we illustrate examples for the numerical solution of (1.1). We use three specific schemes for the state and adjoint variables: conforming finite element (FE) method, non-conforming finite element (ncℙ1FE) method and hybrid mimetic mixed (HMM) method. All three are GDMs with gradient discretisations with bounds on C𝒟, order h estimate on WS𝒟, and satisfying Assumptions A1–A4, and (3.12) on quasi-uniform meshes; see [15, 20] and Remarks 2.4 and 3.4.

The control variable is discretised using piecewise constant functions on the corresponding meshes. The discrete solution is computed using the modified active set algorithm mentioned in Section 4.1 with zero as an initial guess for both u and μ. Here, Ua and Ya denote the average values of the computed control u¯h and the reconstructed state solution Π𝒟⁢y¯𝒟 respectively. Let ni denote the number of iterations required for the convergence of the modified active set algorithm, and fa denote the numerical average of the source term f calculated using the same quadrature rule as in the implementation of the schemes, i.e.,

fa=1|Ω|⁢∑K∈𝒯|K|⁢f⁢(𝒙¯K),

where 𝒙¯K denotes the center of mass of the cell K. This numerical average enables us to evaluate the quality of the quadrature rule for each mesh; in particular, since f has a zero average, any quantity of the order of fa can be considered to be equal to zero, up to quadrature error. The relative errors are denoted by

𝖾𝗋𝗋𝒟⁢(y¯):=∥Π𝒟⁢y¯𝒟-y¯𝒯∥∥y¯∥,𝖾𝗋𝗋𝒟⁢(∇⁡y¯):=∥∇𝒟⁡y¯𝒟-∇⁡y¯∥∥∇⁡y¯∥,

𝖾𝗋𝗋𝒟⁢(p¯):=∥Π𝒟⁢p¯𝒟-p¯𝒯∥∥p¯∥,𝖾𝗋𝗋𝒟⁢(∇⁡p¯):=∥∇𝒟⁡p¯𝒟-∇⁡p¯∥∥∇⁡p¯∥,

𝖾𝗋𝗋⁢(u¯):=∥u¯h-u¯∥∥u¯∥,𝖾𝗋𝗋⁢(u~):=∥u~h-u~∥∥u¯∥.

The data in the optimal control problem (1.1) are chosen as follows:

y¯=2⁢cos⁡(π⁢x)⁢cos⁡(π⁢y),p¯=2⁢cos⁡(π⁢x)⁢cos⁡(π⁢y),α=1,𝒰ad=[a,b],u¯=P[a,b]⁢(-p¯+c¯),

where c¯ is chosen to ensure that ⨍Ωu¯⁢d𝒙=0. The matrix-valued function is given by A=Id unless otherwise specified. The source term f and the desired state y¯d are then computed using

f=-Δ⁢y¯-u¯,y¯d=y¯+Δ⁢p¯.

4.2.1 Example 1: Ω=(0,1)2, ρ=10-4, a=-1, b=1

We here consider the computational domain Ω=(0,1)2. We have p¯⁢(x,y)=-p¯⁢(1-x,y) and, since P[-1,1] is odd, P[-1,1]⁢(-p¯)⁢(1-x,y)=-P[-1,1]⁢(-p¯)⁢(x,y). Integrating this relation over Ω shows that P[-1,1]⁢(-p¯) has a zero average and thus, by Lemma 3.10, that c¯=0. We thus see that u¯=P[-1,1]⁢(-p¯).

Conforming FEM.

The discrete solution is computed on a family of uniform grids with mesh sizes h=12i, i=2,…,6. Due to the symmetry of the mesh and of the solution, approximate solutions are also symmetric and thus have zero average at an order compatible with the stopping criterion in the active set algorithm (the discrete solutions of (3.1) are only approximated by this algorithm), see Table 1. As also seen in this table, the number of iterations of the modified active set algorithm remains very small, and independent on the mesh size. The error estimates and the convergence rates of the control, the post-processed control, the state and the adjoint variables are presented in Table 2. The numerical results corroborate Theorem 3.1, Theorem 3.6 and Corollary 3.8.

Table 1

Example 1: Conforming FEM.

h	Ua	fa	Ya	ni
0.250000	0.002752×10-13	0.20699×10-14	-0.008396×10-13	2
0.125000	-0.008049×10-13	0.20912×10-14	-0.004684×10-13	3
0.062500	-0.001370×10-13	0.20548×10-14	0.010486×10-13	3
0.031250	-0.032432×10-13	0.21299×10-14	0.050725×10-13	3
0.015625	-0.917129×10-13	0.20367×10-14	-0.495753×10-13	3

Table 2

Convergence results. Example 1: Conforming FEM.

h	𝖾𝗋𝗋𝒟⁢(y¯)	Order	𝖾𝗋𝗋𝒟⁢(∇⁡y¯)	Order	𝖾𝗋𝗋𝒟⁢(p¯)	Order
0.250000	0.325104	–	0.293424	–	0.333213	–
0.125000	0.086450	1.9110	0.129922	1.1753	0.089153	1.9021
0.062500	0.022176	1.9628	0.064767	1.0043	0.022967	1.9567
0.031250	0.005591	1.9879	0.032578	0.9914	0.005798	1.9860
0.015625	0.001402	1.9960	0.016337	0.9958	0.001453	1.9960

h	𝖾𝗋𝗋𝒟⁢(∇⁡p¯)	Order	𝖾𝗋𝗋⁢(u¯)	Order	𝖾𝗋𝗋⁢(u~)	Order
0.250000	0.300144	–	0.464300	–	0.222006	–
0.125000	0.131070	1.1953	0.254036	0.8700	0.065430	1.7626
0.062500	0.064930	1.0134	0.126358	1.0075	0.016668	1.9728
0.031250	0.032599	0.9941	0.063453	0.9938	0.004226	1.9797
0.015625	0.016339	0.9965	0.031778	0.9977	0.001047	2.0136

Non-conforming FEM.

For comparison, we compute the solutions of the ncℙ1 finite element method on the same grids. As for conforming FE, the symmetry of the problem ensures that the approximation solutions have a zero average at an order dictated by the stopping criterion used in the active set algorithm. The results in Tables 3 and 4 are similar to those obtained with the conforming FE.

Table 3

Example 1: ncℙ1FEM.

h	Ua	fa	Ya	ni
0.250000	-0.003919×10-13	0.206991×10-14	0.000518×10-12	3
0.125000	-0.067706×10-13	0.209121×10-14	-0.003856×10-12	3
0.062500	0.030900×10-13	0.205478×10-14	0.017217×10-12	3
0.031250	-0.075427×10-13	0.212989×10-14	-0.041154×10-12	3
0.015625	-0.187208×10-13	0.203674×10-14	0.933499×10-12	3

Table 4

Convergence results. Example 1: ncℙ1FEM.

h	𝖾𝗋𝗋𝒟⁢(y¯)	Order	𝖾𝗋𝗋𝒟⁢(∇⁡y¯)	Order	𝖾𝗋𝗋𝒟⁢(p¯)	Order
0.250000	0.148286	–	0.409750	–	0.146306	–
0.125000	0.033274	2.1559	0.189599	1.1118	0.033499	2.1268
0.062500	0.008134	2.0324	0.093105	1.0260	0.008122	2.0443
0.031250	0.002023	2.0077	0.046348	1.0064	0.002025	2.0036
0.015625	0.000505	2.0019	0.023148	1.0016	0.000505	2.0041

h	𝖾𝗋𝗋𝒟⁢(∇⁡p¯)	Order	𝖾𝗋𝗋⁢(u¯)	Order	𝖾𝗋𝗋⁢(u~)	Order
0.250000	0.408120	–	0.473176	–	0.284795	–
0.125000	0.189770	1.1047	0.250457	0.9178	0.071206	1.9999
0.062500	0.093102	1.0274	0.126078	0.9902	0.017716	2.0069
0.031250	0.046349	1.0063	0.063407	0.9916	0.004440	1.9965
0.015625	0.023149	1.0016	0.031770	0.9970	0.001109	2.0007

Table 5

Example 1: HMM.

h	Ua	fa	Ya	ni
0.250000	-0.016326	0.016324	-0.017271	4
0.125000	-0.005300	0.005300	-0.004968	4
0.062500	-0.001503	0.001503	-0.001277	3
0.031250	-0.000352	0.000352	-0.000321	3

Table 6

Convergence results. Example 1: HMM.

h	𝖾𝗋𝗋𝒟⁢(y¯)	Order	𝖾𝗋𝗋𝒟⁢(∇⁡y¯)	Order	𝖾𝗋𝗋𝒟⁢(p¯)	Order
0.250000	0.025586	–	0.143963	–	0.033104	–
0.125000	0.006764	1.9194	0.070970	1.0204	0.010044	1.7207
0.062500	0.001709	1.9847	0.035358	1.0052	0.002443	2.0397
0.031250	0.000429	1.9958	0.017663	1.0013	0.000619	1.9811

h	𝖾𝗋𝗋𝒟⁢(∇⁡p¯)	Order	𝖾𝗋𝗋⁢(u¯)	Order	𝖾𝗋𝗋⁢(u~)	Order
0.250000	0.144012	–	0.214573	–	0.034890	–
0.125000	0.070972	1.0209	0.109352	0.9725	0.009603	1.8613
0.062500	0.035359	1.0052	0.055045	0.9903	0.002403	1.9989
0.031250	0.017663	1.0013	0.027551	0.9985	0.000605	1.9893

HMM Scheme.

This scheme was tested on a series of regular triangular meshes from [25] where the points 𝒫 (see [17, Definition 2.21]) are located at the center of mass of the cells. These meshes are no longer symmetric and thus the symmetry of the approximate solution is lost. Zero averages are thus obtained up to quadrature error, see Table 5. It has been proved in [7, 19] that the state and adjoint equations enjoy a super-convergence property in L2 norm for such a sequence of meshes; hence, as expected from Theorem 3.6, so does the scheme for the entire control problem after post-processing of the control. The errors in the energy norm and the L2 norm, together with their orders of convergence, are presented in Table 6.

For all three methods (conforming ℙ1 FE, ncℙ1 FE and HMM), the theoretical rates of convergence are confirmed by the numerical outputs. Without post-processing, an 𝒪⁢(h) convergence rate is obtained on the controls, which validates Theorem 3.1. With post-processing of the controls, the order of convergence of Theorem 3.6 is recovered. We also notice that the super-convergence on the state and adjoint stated in Corollary 3.8 is confirmed, provided that the exact state and adjoint are properly projected (usage of the functions y¯𝒯 and p¯𝒯 in 𝖾𝗋𝗋𝒟⁢(y¯) and 𝖾𝗋𝗋𝒟⁢(p¯)).

Remark 4.2.

As seen in Table 5, the modified active set algorithm converges in very few iterations if ρ=10-4. We however found that, if ρ=1, the modified active set algorithm no longer converges. Further work will investigate in more depth the convergence analysis of the modified active set algorithm, to understand better its dependency with respect to ρ.

4.2.2 Example 2: A=100⁢Id, Ω=(0,1)2, ρ=10-2, a=-1, b=1

In this subsection, we present some numerical results for the control problem defined on the unit square domain Ω=(0,1)2 and A=100⁢Id. As explained in Example 1, a=-1 and b=1 imply c¯=0.

Table 7

Example 2: Conforming FEM.

h	Ua	fa	Ya	ni
0.250000	0.002280×10-11	0.209361×10-12	0.009117×10-11	2
0.125000	0.018065×10-11	0.209375×10-12	0.024496×10-11	3
0.062500	0.027564×10-11	0.209400×10-12	-0.012778×10-11	3
0.031250	-0.103755×10-11	0.209420×10-12	-0.028850×10-11	3
0.015625	-0.168277×10-11	0.209297×10-12	0.624160×10-11	3

Table 8

Convergence results. Example 2: Conforming FEM.

h	𝖾𝗋𝗋𝒟⁢(y¯)	Order	𝖾𝗋𝗋𝒟⁢(∇⁡y¯)	Order	𝖾𝗋𝗋𝒟⁢(p¯)	Order
0.250000	0.328223	–	0.296014	–	0.328304	–
0.125000	0.087182	1.9126	0.130232	1.1846	0.087209	1.9125
0.062500	0.022409	1.9600	0.064814	1.0067	0.022417	1.9599
0.031250	0.005653	1.9870	0.032584	0.9921	0.005655	1.9870
0.015625	0.001417	1.9963	0.016338	0.9960	0.001417	1.9963

h	𝖾𝗋𝗋𝒟⁢(∇⁡p¯)	Order	𝖾𝗋𝗋⁢(u¯)	Order	𝖾𝗋𝗋⁢(u~)	Order
0.250000	0.296080	–	0.463836	–	0.218080	–
0.125000	0.130243	1.1848	0.253857	0.8696	0.064390	1.7600
0.062500	0.064816	1.0068	0.126333	1.0068	0.016358	1.9768
0.031250	0.032584	0.9922	0.063449	0.9936	0.004145	1.9807
0.015625	0.016338	0.9960	0.031778	0.9976	0.001026	2.0139

Conforming FEM.

We provide in Table 7 the details of active set algorithm for the conforming finite element method. As expected, the symmetries of the problem provide approximate solutions with a nearly perfect average. For such grids, we obtain super-convergence result for the post-processed control. The errors between the true and computed solutions are computed for different mesh sizes and presented in Table 8. They still follow the expected theoretical rates, and the number of iterations of the active set algorithm remain small.

Non-conforming FEM.

The results, presented in Tables 9 and 10, are similar to those for the conforming FE scheme.

Table 9

Example 2: ncℙ1FEM.

h	Ua	fa	Ya	ni
0.250000	-0.000977×10-10	0.209361×10-12	0.000739×10-10	3
0.125000	-0.006518×10-10	0.209375×10-12	-0.000960×10-10	3
0.062500	-0.004320×10-10	0.209400×10-12	0.002330×10-10	3
0.031250	-0.007236×10-10	0.209420×10-12	0.054029×10-10	3
0.015625	0.346321×10-10	0.209297×10-12	0.276914×10-10	3

Table 10

Convergence results. Example 2: ncℙ1FEM.

h	𝖾𝗋𝗋𝒟⁢(y¯)	Order	𝖾𝗋𝗋𝒟⁢(∇⁡y¯)	Order	𝖾𝗋𝗋𝒟⁢(p¯)	Order
0.250000	0.148286	–	0.409750	–	0.148262	–
0.125000	0.033263	2.1564	0.189599	1.1118	0.033265	2.1561
0.062500	0.008131	2.0324	0.093105	1.0260	0.008131	2.0325
0.031250	0.002022	2.0077	0.046348	1.0064	0.002022	2.0076
0.015625	0.000505	2.0019	0.023148	1.0016	0.000505	2.0019

h	𝖾𝗋𝗋𝒟⁢(∇⁡p¯)	Order	𝖾𝗋𝗋⁢(u¯)	Order	𝖾𝗋𝗋⁢(u~)	Order
0.250000	0.409732	–	0.473136	–	0.286537	–
0.125000	0.189600	1.1117	0.250390	0.9181	0.071004	2.0128
0.062500	0.093105	1.0260	0.126079	0.9899	0.017719	2.0026
0.031250	0.046348	1.0064	0.063407	0.9916	0.004439	1.9970
0.015625	0.023148	1.0016	0.031770	0.9970	0.001109	2.0004

HMM Scheme.

The results are presented in Tables 11 and 12. They are qualitatively similar to those for Example 1. As mentioned before, the algorithm is not convergent for ρ=1.

Table 11

Example 2: HMM.

h	Ua	fa	Ya	ni
0.250000	-1.000000	1.817180	81.718002	–
0.125000	-0.528617	0.523335	-0.528289	6
0.062500	-0.136036	0.134678	-0.135812	5
0.031250	-0.034208	0.033866	-0.034178	5

Table 12

Convergence results. Example 2: HMM.

h	𝖾𝗋𝗋𝒟⁢(y¯)	Order	𝖾𝗋𝗋𝒟⁢(∇⁡y¯)	Order	𝖾𝗋𝗋𝒟⁢(p¯)	Order
0.250000	81.717083	–	0.143996	–	392673.3	–
0.125000	0.528572	7.2724	0.070971	1.0207	0.876152	18.7737
0.062500	0.135884	1.9597	0.035359	1.0052	0.216641	2.0159
0.031250	0.034196	1.9905	0.017663	1.0013	0.054238	1.9979

h	𝖾𝗋𝗋𝒟⁢(∇⁡p¯)	Order	𝖾𝗋𝗋⁢(u¯)	Order	𝖾𝗋𝗋⁢(u~)	Order
0.250000	0.143987	–	1.686504	–	1.691346	–
0.125000	0.070970	1.0207	0.878849	0.9403	0.874778	0.9512
0.062500	0.035358	1.0052	0.237129	1.8899	0.230873	1.9218
0.031250	0.017663	1.0013	0.064673	1.8744	0.058612	1.9778

4.2.3 Example 3: Ω=(0,1)2, ρ=10-4, a=-0.5, b=1

In this case, since P[a,b] is no longer odd, P[a,b]⁢(-p¯) no longer has a zero average and, to compute 𝖾𝗋𝗋𝒟⁢(u¯), we need to find c¯ such that

⨍ΩP[a,b]⁢(-p¯+c¯)⁢d𝒙=0.

This c¯ can be found by a bisection method, by computing the averages on a very thin mesh and bisecting until we find a proper c¯. Using a mesh of size h=0.00195, we find c≈-0.24596797.

Conforming FEM.

The numerical results obtained using conforming finite element method are shown in Tables 13 and 14, respectively. Since there is a loss of symmetry, the approximate solutions have zero averages only up to quadrature error (compare Ua and fa in Table 13). Here, we observed that the modified active set algorithm converges only when ρ≤10-1. When it does, though, the number of iterations remain very small. As in Examples 1 and 2, the theoretical rates of convergence are confirmed by these numerical outputs.

Table 13

Example 3: Conforming FEM.

h	Ua	fa	Ya	ni
0.250000	0.0020160	-0.0020160	0.201602×10-6	4
0.125000	0.0055595	-0.0055595	0.555952×10-6	4
0.062500	-0.0004794	0.0004795	-0.047944×10-6	4
0.031250	0.0001470	-0.0001470	0.014705×10-6	5
0.015625	-0.0000136	0.0000136	-0.001362×10-6	5

Table 14

Convergence results. Example 3: Conforming FEM.

h	𝖾𝗋𝗋𝒟⁢(y¯)	Order	𝖾𝗋𝗋𝒟⁢(∇⁡y¯)	Order	𝖾𝗋𝗋𝒟⁢(p¯)	Order
0.250000	0.325266	–	0.293567	–	0.346894	–
0.125000	0.086733	1.9070	0.130041	1.1747	0.097046	1.8378
0.062500	0.022291	1.9601	0.064790	1.0051	0.025081	1.9521
0.031250	0.005624	1.9868	0.032581	0.9917	0.006219	2.0117
0.015625	0.001410	1.9963	0.016337	0.9959	0.001569	1.9865

h	𝖾𝗋𝗋𝒟⁢(∇⁡p¯)	Order	𝖾𝗋𝗋⁢(u¯)	Order	𝖾𝗋𝗋⁢(u~)	Order
0.250000	0.300149	–	0.466701	–	0.234197	–
0.125000	0.131075	1.1953	0.268982	0.7950	0.064265	1.8656
0.062500	0.064931	1.0134	0.138258	0.9602	0.016053	2.0012
0.031250	0.032599	0.9941	0.069620	0.9898	0.003996	2.0064
0.015625	0.016339	0.9965	0.034944	0.9945	0.001002	1.9950

Non-conforming FEM.

The results are similar to those obtained with the conforming FE method (see Tables 15 and 16).

Table 15

Example 3: ncℙ1FEM.

h	Ua	fa	Ya	ni
0.250000	0.002016	-0.002016	0.0201601×10-5	4
0.125000	0.005560	-0.005559	-0.1301803×10-5	5
0.062500	-0.000480	0.000479	-0.0017424×10-5	5
0.031250	0.000147	-0.000147	0.0011093×10-5	5
0.015625	-0.000014	0.000014	-0.0001436×10-5	5

Table 16

Convergence results. Example 3: ncℙ1FEM.

h	𝖾𝗋𝗋𝒟⁢(y¯)	Order	𝖾𝗋𝗋𝒟⁢(∇⁡y¯)	Order	𝖾𝗋𝗋𝒟⁢(p¯)	Order
0.250000	0.148286	–	0.409750	–	0.141781	–
0.125000	0.033270	2.1561	0.189600	1.1118	0.032889	2.1080
0.062500	0.008133	2.0324	0.093105	1.0260	0.008041	2.0322
0.031250	0.002022	2.0077	0.046348	1.0064	0.001994	2.0118
0.015625	0.000505	2.0019	0.023148	1.0016	0.000498	2.0015

h	𝖾𝗋𝗋𝒟⁢(∇⁡p¯)	Order	𝖾𝗋𝗋⁢(u¯)	Order	𝖾𝗋𝗋⁢(u~)	Order
0.250000	0.408120	–	0.494425	–	0.269888	–
0.125000	0.189770	1.1047	0.269866	0.8735	0.080165	1.7513
0.062500	0.093102	1.0274	0.138223	0.9652	0.019967	2.0054
0.031250	0.046349	1.0063	0.069625	0.9893	0.005091	1.9715
0.015625	0.023149	1.0016	0.034941	0.9947	0.001283	1.9883

HMM Scheme.

Tables 17 and 18 show that the HMM scheme behave similarly to the FE schemes. Note that, here too, the convergence of the modified active set algorithm is only observed if ρ≤10-1.

Table 17

Example 3: HMM.

h	Ua	fa	Ya	ni
0.250000	-0.019043	0.019041	-0.017271	5
0.125000	-0.005459	0.005459	-0.004968	5
0.062500	-0.001300	0.001300	-0.001277	5
0.031250	-0.000331	0.000331	-0.000321	5

Table 18

Convergence results. Example 3: HMM.

h	𝖾𝗋𝗋𝒟⁢(y¯)	Order	𝖾𝗋𝗋𝒟⁢(∇⁡y¯)	Order	𝖾𝗋𝗋𝒟⁢(p¯)	Order
0.250000	0.026037	–	0.144014	–	0.055044	–
0.125000	0.006841	1.9284	0.070972	1.0209	0.013361	2.0425
0.062500	0.001728	1.9853	0.035359	1.0052	0.003342	1.9995
0.031250	0.000433	1.9956	0.017663	1.0013	0.000843	1.9869

h	𝖾𝗋𝗋𝒟⁢(∇⁡p¯)	Order	𝖾𝗋𝗋⁢(u¯)	Order	𝖾𝗋𝗋⁢(u~)	Order
0.250000	0.144013	–	0.237647	–	0.050184	–
0.125000	0.070972	1.0209	0.120112	0.9844	0.013482	1.8962
0.062500	0.035359	1.0052	0.061226	0.9722	0.003468	1.9586
0.031250	0.017663	1.0013	0.030583	1.0014	0.000872	1.9916

Funding source: Australian Research Council

Award Identifier / Grant number: DP170100605

Funding source: Department of Science and Technology, Ministry of Science and Technology

Award Identifier / Grant number: SR/S4/MS/808/12

Funding statement: The first author acknowledges the funding support from the Australian Government through the Australian Research Council’s Discovery Projects funding scheme (project number DP170100605). The second and third authors acknowledge the funding support from the DST project SR/S4/MS/808/12.

A Appendix

The proofs of error estimates for control, state and adjoint variables are obtained by modifying the proofs of the corresponding results in [20]. For the sake of completeness and readability, we provide here detailed proofs, highlighting in chosen places where modifications are required due to the pure Neumann boundary conditions (which mostly amount to making sure that certain averages have been properly fixed).

Proof of Theorem 3.1.

Define the following auxiliary discrete problem: Seek

(y𝒟⁢(u¯),p𝒟⁢(u¯))∈X𝒟,⋆×X𝒟

such that

(A.1)

(A.1a)a𝒟⁢(y𝒟⁢(u¯),w𝒟)=(f+u¯,Π𝒟⁢w𝒟)for all ⁢w𝒟∈X𝒟,⋆,

(A.1b)a𝒟⁢(z𝒟,p𝒟⁢(u¯))=(y¯-y¯d,Π𝒟⁢z𝒟)for all ⁢z𝒟∈X𝒟,

where the co-state p𝒟⁢(u¯) is chosen such that

⨍ΩΠ𝒟⁢p𝒟⁢(u¯)⁢d𝒙=⨍ΩΠ𝒟⁢p¯𝒟⁢d𝒙=⨍Ωp¯⁢d𝒙.

For Neumann boundary conditions, this particular choice is essential as it ensures that p𝒟⁢(u¯)-p¯𝒟∈X𝒟,⋆ can be used as a test function w𝒟 in (3.1a) and (A.1a). Recalling that 𝓟𝒯 is the orthogonal projection on piecewise constant functions on 𝒯, we obtain 𝓟𝒯⁢(𝒰ad)⊂𝒰h. Also, for u∈𝒰ad and K∈𝒯, 𝓟𝒯⁢u|K=⨍Ku⁢d𝒙∈[a,b] and, using (1.3),

⨍Ω𝓟𝒯⁢u⁢d𝒙=∑K∈𝒯⨍K𝓟𝒯⁢u⁢d𝒙=∑K∈𝒯⨍Ku⁢d𝒙=⨍Ωu⁢d𝒙=0.

Hence, 𝓟𝒯⁢(𝒰ad)⊂𝒰ad,h.

Set P𝒟,α⁢(u¯)=α-1⁢Π𝒟⁢p𝒟⁢(u¯), P¯𝒟,α=α-1⁢Π𝒟⁢p¯𝒟 and P¯α=α-1⁢p¯. Since u¯h∈𝒰ad,h⊂𝒰ad and 𝓟𝒯⁢u¯∈𝒰ad,h, from the optimality conditions ((1.4c) and (3.1c)),

-α⁢(P¯α+u¯,u¯-u¯h)≥0,

α⁢(P¯𝒟,α+u¯h,u¯-u¯h)≥α⁢(P¯𝒟,α+u¯h,u¯-𝓟𝒯⁢u¯).

Adding these two inequalities yields

α⁢∥u¯-u¯h∥2≤-α⁢(P¯𝒟,α+u¯h,u¯-𝓟𝒯⁢(u¯))+α⁢(P¯𝒟,α-P¯α,u¯-u¯h)

(A.2)=-α⁢(P¯𝒟,α+u¯h,u¯-𝓟𝒯⁢u¯)+α⁢(P¯𝒟,α-P𝒟,α⁢(u¯),u¯-u¯h)-α⁢(P¯α-P𝒟,α⁢(u¯),u¯-u¯h).

By the orthogonality property of 𝓟𝒯, we have (u¯h,u¯-𝓟𝒯⁢u¯)=0 and (𝓟𝒯⁢P¯α,u¯-𝓟𝒯⁢u¯)=0. Therefore, the first term in the right-hand side of (A.2) can be re-cast as

-α⁢(P¯𝒟,α+u¯h,u¯-𝓟𝒯⁢u¯)=-α⁢(P¯α-𝓟𝒯⁢P¯α,u¯-𝓟𝒯⁢u¯)+α⁢(P¯α-P𝒟,α⁢(u¯),u¯-𝓟𝒯⁢u¯)

(A.3)+α⁢(P𝒟,α⁢(u¯)-P¯𝒟,α,u¯-𝓟𝒯⁢u¯).

By the Cauchy–Schwarz inequality, the first term on the right-hand side of (A.3) is estimated as

(A.4)-α⁢(P¯α-𝓟𝒯⁢P¯α,u¯-𝓟𝒯⁢u¯)≤Eh⁢(p¯)⁢Eh⁢(u¯).

Equation (A.1b) shows that p𝒟⁢(u¯) is the solution of the GS corresponding to the adjoint problem (1.4b), whose solution is p¯. Therefore, using the fact that ⨍ΩΠ𝒟⁢p𝒟⁢(u¯)⁢d𝒙=⨍Ωp¯⁢d𝒙 (note that the specific relation between the continuous and discrete co-states is essential here), by Theorem 2.3,

(A.5)∥P¯α-P𝒟,α⁢(u¯)∥=α-1⁢∥p¯-Π𝒟⁢p𝒟⁢(u¯)∥≲α-1⁢WS𝒟⁢(p¯).

Hence, using the Cauchy–Schwarz inequality,

(A.6)α⁢(P¯α-P𝒟,α⁢(u¯),u¯-𝓟𝒯⁢u¯)≲WS𝒟⁢(p¯)⁢Eh⁢(u¯).

Using the definitions of C𝒟, ∥⋅∥𝒟 and the fact that p𝒟⁢(u¯)-p¯𝒟∈X𝒟,⋆, we find that

(A.7)∥Π𝒟⁢p𝒟⁢(u¯)-Π𝒟⁢p¯𝒟∥2≲∥p𝒟⁢(u¯)-p¯𝒟∥𝒟2=∥∇𝒟⁡p𝒟⁢(u¯)-∇𝒟⁡p¯𝒟∥2.

By writing the difference of (A.1b) and (3.1b), we see that p𝒟⁢(u¯)-p¯𝒟 is the solution to the GS (2.4) with source term F=y¯-Π𝒟⁢y¯𝒟. i.e., for all z𝒟∈X𝒟,

a𝒟⁢(z𝒟,p𝒟⁢(u¯)-p¯𝒟)=(y¯-Π𝒟⁢y¯𝒟,Π𝒟⁢z𝒟).

Choose z𝒟=p𝒟⁢(u¯)-p¯𝒟 in the above equality and use it in (A.7) to obtain

∥Π𝒟⁢p𝒟⁢(u¯)-Π𝒟⁢p¯𝒟∥2≲∥∇𝒟⁡p𝒟⁢(u¯)-∇𝒟⁡p¯𝒟∥2≲∥y¯-Π𝒟⁢y¯𝒟∥⁢∥Π𝒟⁢p𝒟⁢(u¯)-Π𝒟⁢p¯𝒟∥.

As a consequence,

∥P𝒟,α⁢(u¯)-P¯𝒟,α∥=α-1⁢∥Π𝒟⁢p𝒟⁢(u¯)-Π𝒟⁢p¯𝒟∥≲α-1⁢∥y¯-Π𝒟⁢y¯𝒟∥

≲α-1⁢∥y¯-Π𝒟⁢y𝒟⁢(u¯)∥+α-1⁢∥Π𝒟⁢y𝒟⁢(u¯)-Π𝒟⁢y¯𝒟∥.

Use Theorem 2.3 with ψ=y¯ to bound the first term in the above expression. This along with an application of Young’s inequality yields an estimate for the last term in (A.3) as

(A.8)α⁢(P𝒟,α⁢(u¯)-P¯𝒟,α,u¯-𝓟𝒯⁢u¯)≤C1⁢Eh⁢(u¯)⁢WS𝒟⁢(y¯)+C1⁢Eh⁢(u¯)2+14⁢∥Π𝒟⁢y𝒟⁢(u¯)-Π𝒟⁢y¯𝒟∥2,

where C1 depends only on Ω, A and an upper bound of C𝒟. Plugging (A.4), (A.6) and (A.8) into (A.3) yields

-α⁢(P¯𝒟,α+u¯h,u¯-𝓟𝒯⁢u¯)≤Eh⁢(p¯)⁢Eh⁢(u¯)+C2⁢Eh⁢(u¯)⁢WS𝒟⁢(p¯)+C1⁢Eh⁢(u¯)⁢WS𝒟⁢(y¯)

(A.9)+C1⁢Eh⁢(u¯)2+14⁢∥Π𝒟⁢y𝒟⁢(u¯)-Π𝒟⁢y¯𝒟∥2,

where C2 is the hidden constant in (A.6). Let us turn to the second term on the right-hand side of (A.2). From (3.1b) and (A.1b), for all z𝒟∈X𝒟,

(A.10)a𝒟⁢(z𝒟,p¯𝒟-p𝒟⁢(u¯))=(Π𝒟⁢y¯𝒟-y¯,Π𝒟⁢z𝒟).

Also, from (3.1a) and (A.1a), for all w𝒟∈X𝒟,⋆,

(A.11)a𝒟⁢(y¯𝒟-y𝒟⁢(u¯),w𝒟)=(u¯h-u¯,Π𝒟⁢w𝒟).

Choose z𝒟=y¯𝒟-y𝒟⁢(u¯)∈X𝒟 in (A.10), w𝒟=p¯𝒟-p𝒟⁢(u¯)∈X𝒟,⋆ in (A.11), use the symmetry of a𝒟⁢(⋅,⋅), Theorem 2.3 with ψ=y¯ and Young’s inequality to obtain

α⁢(P¯𝒟,α-P𝒟,α⁢(u¯),u¯-u¯h)=-(Π𝒟⁢y¯𝒟-y¯,Π𝒟⁢y¯𝒟-Π𝒟⁢y𝒟⁢(u¯))

=(y¯-Π𝒟⁢y𝒟⁢(u¯),Π𝒟⁢y¯𝒟-Π𝒟⁢y𝒟⁢(u¯))-∥Π𝒟⁢y¯𝒟-Π𝒟⁢y𝒟⁢(u¯)∥2

≲WS𝒟⁢(y¯)⁢∥Π𝒟⁢y¯𝒟-Π𝒟⁢y𝒟⁢(u¯)∥-∥Π𝒟⁢y¯𝒟-Π𝒟⁢y𝒟⁢(u¯)∥2

(A.12)≤C3⁢WS𝒟⁢(y¯)2+14⁢∥Π𝒟⁢y¯𝒟-Π𝒟⁢y𝒟⁢(u¯)∥2-∥Π𝒟⁢y¯𝒟-Π𝒟⁢y𝒟⁢(u¯)∥2,

where C3 only depends on Ω, A and an upper bound of C𝒟. The last term on the right-hand side of (A.2) can be estimated using (A.5) and Young’s inequality:

(A.13)-α⁢(P¯α-P𝒟,α⁢(u¯),u¯-u¯h)≤α2⁢∥u¯-u¯h∥2+C4⁢WS𝒟⁢(p¯)2,

where C4 only depends on Ω, A, α and an upper bound of C𝒟. Substitute (A.9), (A.12) and (A.13) into (A.2), apply the Young’s inequality and ∑iai2≤∑iai to complete the proof. ∎

Proof of Proposition 3.2.

An application of triangle inequality yields

∥Π𝒟⁢y¯𝒟-y¯∥+∥∇𝒟⁡y¯𝒟-∇⁡y¯∥≤∥Π𝒟⁢y¯𝒟-Π𝒟⁢y𝒟⁢(u¯)∥+∥∇𝒟⁡y¯𝒟-∇𝒟⁡y𝒟⁢(u¯)∥

(A.14)+∥Π𝒟⁢y𝒟⁢(u¯)-y¯∥+∥∇𝒟⁡y𝒟⁢(u¯)-∇⁡y¯∥.

Subtracting (3.1a) and (A.1a), and using the stability property of GS (Proposition 2.6), the first two terms on the right-hand side of the above inequality can be estimated as

∥Π𝒟⁢y¯𝒟-Π𝒟⁢y𝒟⁢(u¯)∥+∥∇𝒟⁡y¯𝒟-∇𝒟⁡y𝒟⁢(u¯)∥≲∥u¯-u¯h∥.

The last two terms on the right-hand side of the (A.14) are estimated using Theorem 2.3 as

∥Π𝒟⁢y𝒟⁢(u¯)-y¯∥+∥∇𝒟⁡y𝒟⁢(u¯)-∇⁡y¯∥≲WS𝒟⁢(y¯).

A combination of the above two results yields the error estimates (3.6) for the state variable. A use of ⨍ΩΠ𝒟⁢p𝒟⁢(u¯)⁢d𝒙=⨍ΩΠ𝒟⁢p¯𝒟⁢d𝒙 in Proposition 2.6 leads to the error estimates for the adjoint variable in a similar way. ∎

Proof of Theorem 3.5.

Consider the auxiliary problem defined by: For g∈L2⁢(Ω), let p𝒟*⁢(g)∈X𝒟 solve

(A.15)a𝒟⁢(z𝒟,p𝒟*⁢(g))=(Π𝒟⁢y𝒟⁢(g)-y¯d,Π𝒟⁢z𝒟) for all ⁢z𝒟∈X𝒟,

where y𝒟⁢(g) is given by (A.1a) with u¯ replaced by g. We fix p𝒟*⁢(g) by imposing ⨍ΩΠ𝒟⁢p𝒟*⁢(g)⁢d𝒙=⨍Ωp¯⁢d𝒙. This choice is dictated by the pure Neumann boundary condition and will be essential.

For K∈𝒯, let 𝒙¯K be the centroid (centre of mass) of K. A standard approximation property (see e.g. [19, Lemma A.7] with wK≡1) yields

(A.16)∀K∈𝒯,∀ϕ∈H2⁢(K),∥𝓟𝒯⁢ϕ-ϕ⁢(𝒙¯K)∥L2⁢(K)≲ηdiam⁢(K)2⁢∥ϕ∥H2⁢(K).

Define u^ and p^ a.e. on Ω by: For all K∈𝒯 and all 𝒙∈K, u^⁢(𝒙)=u¯⁢(𝒙¯K) and p^⁢(𝒙)=p¯⁢(𝒙¯K). From (3.11) and the Lipschitz continuity of P[a,b], we obtain

∥u~-u~h∥≤α-1⁢∥Π𝒟⁢p¯𝒟-p¯𝒯∥

≤α-1⁢∥p¯𝒯-Π𝒟⁢p𝒟*⁢(u¯)∥+α-1⁢∥Π𝒟⁢p𝒟*⁢(u¯)-Π𝒟⁢p𝒟*⁢(u^)∥

+α-1⁢∥Π𝒟⁢p𝒟*⁢(u^)-Π𝒟⁢p¯𝒟∥

(A.17)=:α-1T1+α-1T2+α-1T3.

Step 1: Estimate of T1. A use of triangle inequality, (1.4b), (A.1b) and Assumption A1 (i) leads to

T1≤∥p¯𝒯-Π𝒟⁢p𝒟⁢(u¯)∥+∥Π𝒟⁢p𝒟⁢(u¯)-Π𝒟⁢p𝒟*⁢(u¯)∥

(A.18)≲h2⁢∥y¯-y¯d∥H1⁢(Ω)+∥Π𝒟⁢p𝒟⁢(u¯)-Π𝒟⁢p𝒟*⁢(u¯)∥.

We now estimate the last term in this inequality. Use the definitions of C𝒟, ∥⋅∥𝒟 and the fact that

⨍ΩΠ𝒟⁢p𝒟⁢(u¯)⁢d𝒙=⨍ΩΠ𝒟⁢p𝒟*⁢(u¯)⁢d𝒙

to obtain

(A.19)∥Π𝒟⁢(p𝒟⁢(u¯)-p𝒟*⁢(u¯))∥2≲∥∇𝒟⁡(p𝒟⁢(u¯)-p𝒟*⁢(u¯))∥2.

Subtract (A.15) with g=u¯ from (A.1b), substitute z𝒟=p𝒟⁢(u¯)-p𝒟*⁢(u¯), use property (3.9) in Assumption A1 (ii) and the Cauchy–Schwarz inequality to obtain

∥∇𝒟⁡(p𝒟⁢(u¯)-p𝒟*⁢(u¯))∥2≲a𝒟⁢(p𝒟⁢(u¯)-p𝒟*⁢(u¯),p𝒟⁢(u¯)-p𝒟*⁢(u¯))

=(y¯-Π𝒟⁢y𝒟⁢(u¯),Π𝒟⁢(p𝒟⁢(u¯)-p𝒟*⁢(u¯)))

=(y¯-y¯𝒯,Π𝒟⁢(p𝒟⁢(u¯)-p𝒟*⁢(u¯)))+(y¯𝒯-Π𝒟⁢y𝒟⁢(u¯),Π𝒟⁢(p𝒟⁢(u¯)-p𝒟*⁢(u¯)))

≲h2⁢∥y¯∥H2⁢(Ω)⁢∥Π𝒟⁢(p𝒟⁢(u¯)-p𝒟*⁢(u¯))∥+∥y¯𝒯-Π𝒟⁢y𝒟⁢(u¯)∥⁢∥Π𝒟⁢(p𝒟⁢(u¯)-p𝒟*⁢(u¯))∥.

A use of (A.19) and Assumption A1 (i) leads to

∥Π𝒟⁢p𝒟⁢(u¯)-Π𝒟⁢p𝒟*⁢(u¯)∥≲h2⁢∥y¯∥H2⁢(Ω)+h2⁢∥f+u¯∥H1⁢(Ω).

Plugged into (A.18), this estimate yields

(A.20)T1≲h2⁢(∥y¯-y¯d∥H1⁢(Ω)+∥y¯∥H2⁢(Ω)+∥f+u¯∥H1⁢(Ω)).

Step 2: Estimate of T2. Subtract equations (A.15) satisfied by p𝒟*⁢(u¯) and p𝒟*⁢(u^) to obtain, for all z𝒟∈X𝒟,

(A.21)a𝒟⁢(z𝒟,p𝒟*⁢(u¯)-p𝒟*⁢(u^))=(Π𝒟⁢y𝒟⁢(u¯)-Π𝒟⁢y𝒟⁢(u^),Π𝒟⁢z𝒟).

Since p𝒟*⁢(u^)-p𝒟*⁢(u¯)∈X𝒟,⋆, a use of Proposition 2.6 in (A.21) yields

(A.22)T2=∥Π𝒟⁢p𝒟*⁢(u¯)-Π𝒟⁢p𝒟*⁢(u^)∥≲∥Π𝒟⁢y𝒟⁢(u¯)-Π𝒟⁢y𝒟⁢(u^)∥.

Choose z𝒟=y𝒟⁢(u¯)-y𝒟⁢(u^) in equation (A.21), subtract equations (A.1a) satisfied by y𝒟⁢(u¯) and y𝒟⁢(u^), since p𝒟*⁢(u¯)-p𝒟*⁢(u^)∈X𝒟,⋆, to obtain

∥Π𝒟⁢(y𝒟⁢(u¯)-y𝒟⁢(u^))∥2=a𝒟⁢(y𝒟⁢(u¯)-y𝒟⁢(u^),p𝒟*⁢(u¯)-p𝒟*⁢(u^))

=(u¯-u^,Π𝒟⁢p𝒟*⁢(u¯)-Π𝒟⁢p𝒟*⁢(u^)).

Set w𝒟=p𝒟*⁢(u¯)-p𝒟*⁢(u^), use the orthogonality of 𝓟𝒯, the Cauchy–Schwarz inequality and Assumption A2 to infer

∥Π𝒟⁢(y𝒟⁢(u¯)-y𝒟⁢(u^))∥2=(u¯-u^,Π𝒟⁢w𝒟)

=(u¯-𝓟𝒯⁢u¯,Π𝒟⁢w𝒟-𝓟𝒯⁢(Π𝒟⁢w𝒟))+(𝓟𝒯⁢u¯-u^,Π𝒟⁢w𝒟)

(A.23)≲ηh⁢∥u¯∥H1⁢(Ω)⁢h⁢∥∇𝒟⁡w𝒟∥+∫Ω1,𝒯(𝓟𝒯⁢u¯-u^)⁢Π𝒟⁢w𝒟⁢d𝒙+∫Ω2,𝒯(𝓟𝒯⁢u¯-u^)⁢Π𝒟⁢w𝒟⁢d𝒙.

Equation (A.21) and the stability of the GDM (Proposition 2.6) show that

(A.24)∥∇𝒟⁡w𝒟∥=∥∇𝒟⁡(p𝒟*⁢(u¯)-p𝒟*⁢(u^))∥≲∥Π𝒟⁢(y𝒟⁢(u¯)-y𝒟⁢(u^))∥.

A use of Hölder’s inequality, Assumptions A4 and A3, the fact that w𝒟∈X𝒟,⋆ and (A.24) yields an estimate for the second term on the right-hand side of (A.23) as follows:

∫Ω1,𝒯(𝓟𝒯⁢u¯-u^)⁢Π𝒟⁢w𝒟⁢d𝒙≤∥𝓟𝒯⁢u¯-u^∥L2⁢(Ω1,𝒯)⁢∥Π𝒟⁢w𝒟∥L2⁢(Ω1,𝒯)

≤h⁢∥u¯∥W1,∞⁢(𝒯1)⁢|Ω1,𝒯|12⁢∥Π𝒟⁢w𝒟∥L2*⁢(Ω)⁢|Ω1,𝒯|12-12*

≲h2-12*⁢∥u¯∥W1,∞⁢(𝒯1)⁢∥w𝒟∥𝒟

=h2-12*⁢∥u¯∥W1,∞⁢(𝒯1)⁢∥∇𝒟⁡w𝒟∥

(A.25)≲h2-12*⁢∥u¯∥W1,∞⁢(𝒯1)⁢∥Π𝒟⁢(y𝒟⁢(u¯)-y𝒟⁢(u^))∥.

Consider now the last term on the right-hand side of (A.23). For any K∈𝒯2, we have u¯=a on K, u¯=b on K, or, by (3.16), u¯=-α-1⁢p¯+c¯ on K. Hence, on K, 𝓟𝒯⁢u¯-u^=0 or 𝓟𝒯⁢u¯-u^=α-1⁢(p^-𝓟𝒯⁢p¯). This leads to |𝓟𝒯⁢u¯-u^|≤α-1⁢|p^-𝓟𝒯⁢p¯| on Ω2,𝒯. Use (A.16), the definition of C𝒟, the fact that w𝒟∈X𝒟,⋆ and (A.24) to obtain

∫Ω2,𝒯(𝓟𝒯⁢u¯-u^)⁢Π𝒟⁢w𝒟⁢d𝒙≤∥𝓟𝒯⁢u¯-u^∥L2⁢(Ω2,𝒯)⁢∥Π𝒟⁢w𝒟∥

≤α-1⁢∥𝓟𝒯⁢p¯-p^∥L2⁢(Ω2,𝒯)⁢∥Π𝒟⁢w𝒟∥

≲ηh2⁢α-1⁢∥p¯∥H2⁢(Ω2,𝒯)⁢∥∇𝒟⁡w𝒟∥

(A.26)≲ηh2⁢α-1⁢∥p¯∥H2⁢(Ω2,𝒯)⁢∥Π𝒟⁢(y𝒟⁢(u¯)-y𝒟⁢(u^))∥.

Plug (A.24), (A.25) and (A.26) into (A.23) and then into (A.22) to get

(A.27)T2≲ηh2-12*⁢∥u¯∥W1,∞⁢(𝒯1)+h2⁢(∥u¯∥H1⁢(Ω)+α-1⁢∥p¯∥H2⁢(Ω2,𝒯)).

Step 3: Estimate of T3. Subtract (3.1b) from (A.15) with g=u^ and (3.1a) from (A.1a) with u^ instead of u¯, we obtain for all z𝒟∈X𝒟 and w𝒟∈X𝒟,⋆,

(A.28)a𝒟⁢(z𝒟,p𝒟*⁢(u^)-p¯𝒟)=(Π𝒟⁢y𝒟⁢(u^)-Π𝒟⁢y¯𝒟,Π𝒟⁢z𝒟)

and

(A.29)a𝒟⁢(y𝒟⁢(u^)-y¯𝒟,w𝒟)=(u^-u¯h,Π𝒟⁢w𝒟).

Substitute z𝒟=p𝒟*⁢(u^)-p¯𝒟∈X𝒟,⋆ in (A.28), w𝒟=y𝒟⁢(u^)-y¯𝒟∈X𝒟,⋆ in (A.29) and use Proposition 2.6 to obtain

(A.30)T3=∥Π𝒟⁢p𝒟*⁢(u^)-Π𝒟⁢p¯𝒟∥≲∥Π𝒟⁢y𝒟⁢(u^)-Π𝒟⁢y¯𝒟∥≲∥u^-u¯h∥.

The optimality condition (1.4c) (see [30, Lemma 3.5]) yields for a.e. 𝒙∈Ω,

(p¯⁢(𝒙)+α⁢u¯⁢(𝒙))⁢(v⁢(𝒙)-u¯⁢(𝒙))≥0 for all ⁢v∈𝒰ad.

Since u¯, p¯ and u¯h are continuous at the centroid 𝒙¯K, we can choose 𝒙=𝒙¯K and v(𝒙¯K)=u¯h(𝒙¯K)(=u¯h on K). All the involved functions being constants over K, this gives

(p^+α⁢u^)⁢(u¯h-u^)≥0 on K, for all K∈𝒯.

Integrate over K and sum over K∈𝒯 to deduce

(p^+α⁢u^,u¯h-u^)≥0.

Choose vh=u^ in the discrete optimality condition (3.1c) to establish

(Π𝒟⁢p¯𝒟+α⁢u¯h,u^-u¯h)≥0.

Add the above two inequalities to obtain

(p^-Π𝒟⁢p¯𝒟+α⁢(u^-u¯h),u¯h-u^)≥0,

and thus

α⁢∥u^-u¯h∥2≤(p^-Π𝒟⁢p¯𝒟,u¯h-u^)

=(p^-p¯𝒯,u¯h-u^)+(p¯𝒯-Π𝒟⁢p𝒟*⁢(u^),u¯h-u^)+(Π𝒟⁢p𝒟*⁢(u^)-Π𝒟⁢p¯𝒟,u¯h-u^)

(A.31)=:M1+M2+M3.

Since u¯h-u^ is piecewise constant on 𝒯, the orthogonality property of 𝓟𝒯, (A.16) and (3.10) in Assumption A1 (ii) lead to

M1=(p^-𝓟𝒯⁢p¯𝒯,u¯h-u^)

=(p^-𝓟𝒯⁢p¯,u¯h-u^)+(𝓟𝒯⁢(p¯-p¯𝒯),u¯h-u^)

(A.32)≲ηh2⁢∥p¯∥H2⁢(Ω)⁢∥u¯h-u^∥.

By the Cauchy–Schwarz inequality, the triangle inequality and the notations in (A.17), we obtain

(A.33)M2≤∥p¯𝒯-Π𝒟⁢p𝒟*⁢(u^)∥⁢∥u¯h-u^∥≲(T1+T2)⁢∥u¯h-u^∥.

Subtract equations (3.1a) and (A.1a) (with u^ instead of u¯) satisfied by y¯𝒟 and y𝒟⁢(u^), choose w𝒟=p𝒟*⁢(u^)-p¯𝒟 and use equations (3.1b) and (A.15) on p¯𝒟 and p𝒟*⁢(u^) to arrive at

M3=(Π𝒟⁢(p𝒟*⁢(u^)-p¯𝒟),u¯h-u^)

=a𝒟⁢(y¯𝒟-y𝒟⁢(u^),p𝒟*⁢(u^)-p¯𝒟)

(A.34)=(Π𝒟(y𝒟(u^)-y¯𝒟),Π𝒟(y¯𝒟-y𝒟(u^))≤0.

A substitution of (A.32)–(A.34) (together with estimates (A.20) and (A.27) of T1 and T2) into (A.31) yields an estimate on ∥u¯h-u^∥ which, when plugged into (A.30), gives

(A.35)T3≲ηα-1h2-12*∥u¯∥W1,∞⁢(𝒯1)+α-1h2(∥y¯-y¯d∥H1⁢(Ω)+∥y¯∥H2⁢(Ω)+(1+α-1)∥p¯∥H2⁢(Ω)+∥f+u¯∥H1⁢(Ω)+∥u¯∥H1⁢(Ω)).

Step 4: Conclusion. A use of (1.2) and the fact that u¯ is optimal leads to

α2⁢∥u¯∥2≤J⁢(y¯,u¯)≤J⁢(y⁢(0),0)=12⁢∥y⁢(0)-y¯d∥2,

where y⁢(0) is the solution to the state equation (1.1b) with u=0. Hence,

(A.36)∥u¯∥≲α-1⁢(∥f∥+∥y¯d∥).

From (3.16) and (3.3), we have

∇⁡u¯=∇⁡P[a,b]⁢(-α-1⁢p¯+c¯)=𝟙(-α-1⁢p¯+c¯)∈[a,b]⁢∇⁡(-α-1⁢p¯+c¯).

Note that |∇⁡(-α-1⁢p¯+c¯)|=α-1⁢|∇⁡p¯|. Therefore,

(A.37)∥∇⁡u¯∥2=∫Ω|∇⁡u¯|2⁢d𝒙=∫Ω|𝟙(-α-1⁢p¯+c¯)∈[a,b]⁢∇⁡(-α-1⁢p¯+c¯)|2⁢d𝒙≲α-2⁢∥∇⁡p¯∥2.

Combine (A.36) and (A.37) to obtain

(A.38)∥u¯∥H1⁢(Ω)≲α-1⁢(∥f∥+∥y¯d∥)+α-1⁢∥∇⁡p¯∥.

Use (A.38) in (A.20), (A.27) and (A.35) and plug the resulting estimates into (A.17) to complete the proof. ∎

Proof of Theorem 3.6.

The proof of this theorem is identical to the proof of Theorem 3.5, except for the estimate of T2. This estimate is the only source of the 2-12* power (instead of 2), and the only place where we used Assumption A3, here replaced by (3.12). Recall Assumption A4 and use (3.12) in (A.21) satisfied by p𝒟*⁢(u¯)-p𝒟*⁢(u^) to write

∫Ω1,𝒯(𝓟𝒯⁢u¯-u^)⁢Π𝒟⁢w𝒟⁢d𝒙=∫Ω1,𝒯(𝓟𝒯⁢u¯-u^)⁢(Π𝒟⁢p𝒟*⁢(u¯)-Π𝒟⁢p𝒟*⁢(u^))⁢d𝒙

≲∥𝓟𝒯⁢u¯-u^∥L∞⁢(Ω1,𝒯)⁢∥Π𝒟⁢p𝒟*⁢(u¯)-Π𝒟⁢p𝒟*⁢(u^)∥L∞⁢(Ω1,𝒯)⁢|Ω1,𝒯|

≲h2⁢∥u¯∥W1,∞⁢(𝒯1)⁢∥Π𝒟⁢p𝒟*⁢(u¯)-Π𝒟⁢p𝒟*⁢(u^)∥L∞⁢(Ω)

≲h2⁢∥u¯∥W1,∞⁢(𝒯1)⁢δ⁢∥Π𝒟⁢y𝒟⁢(u¯)-Π𝒟⁢y𝒟⁢(u^)∥.

The rest of the proof follows from this estimate. ∎

Proof of Corollary 3.8.

The result for the state and adjoint variables can be derived exactly as in [20, Corollary 3.7]. ∎

References

[1] Y. Alnashri and J. Droniou, Gradient schemes for the Signorini and the obstacle problems, and application to hybrid mimetic mixed methods, Comput. Math. Appl. 72 (2016), 2788–2807. 10.1016/j.camwa.2016.10.004Search in Google Scholar

[2] T. Apel, J. Pfefferer and A. Rösch, Finite element error estimates for Neumann boundary control problems on graded meshes, Comput. Optim. Appl. 52 (2012), no. 1, 3–28. 10.1007/s10589-011-9427-xSearch in Google Scholar

[3] T. Apel, J. Pfefferer and A. Rösch, Finite element error estimates on the boundary with application to optimal control, Math. Comp. 84 (2015), no. 291, 33–70. 10.1090/S0025-5718-2014-02862-7Search in Google Scholar

[4] L. Beirão da Veiga, K. Lipnikov and G. Manzini, The Mimetic Finite Difference Method for Elliptic Problems, MS&A. Model. Simul. Appl. 11, Springer, Cham, 2014. 10.1007/978-3-319-02663-3Search in Google Scholar

[5] 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

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

[7] F. Brezzi, K. Lipnikov and V. Simoncini, A family of mimetic finite difference methods on polygonal and polyhedral meshes, Math. Models Methods Appl. Sci. 15 (2005), no. 10, 1533–1551. 10.1142/S0218202505000832Search in Google Scholar

[8] E. Casas and K. Kunisch, Optimal control of semilinear elliptic equations in measure spaces, SIAM J. Control Optim. 52 (2014), no. 1, 339–364. 10.1137/13092188XSearch in Google Scholar

[9] E. Casas, M. Mateos and J.-P. Raymond, Error estimates for the numerical approximation of a distributed control problem for the steady-state Navier–Stokes equations, SIAM J. Control Optim. 46 (2007), no. 3, 952–982. 10.1137/060649999Search in Google Scholar

[10] E. Casas, M. Mateos and J.-P. Raymond, Penalization of Dirichlet optimal control problems, ESAIM Control Optim. Calc. Var. 15 (2009), no. 4, 782–809. 10.1051/cocv:2008049Search in Google Scholar

[11] E. Casas, M. Mateos and F. Tröltzsch, Error estimates for the numerical approximation of boundary semilinear elliptic control problems, Comput. Optim. Appl. 31 (2005), no. 2, 193–219. 10.1007/s10589-005-2180-2Search in Google Scholar

[12] E. Casas and J.-P. Raymond, Error estimates for the numerical approximation of Dirichlet boundary control for semilinear elliptic equations, SIAM J. Control Optim. 45 (2006), no. 5, 1586–1611. 10.1137/050626600Search in Google Scholar

[13] E. Casas and F. Tröltzsch, A general theorem on error estimates with application to a quasilinear elliptic optimal control problem, Comput. Optim. Appl. 53 (2012), no. 1, 173–206. 10.1007/s10589-011-9453-8Search in Google Scholar

[14] G. Chavent and J. Jaffré, Mathematical Models and Finite Elements for Reservoir Simulation, Stud. Math. Appl. 17, North-Holland, Amsterdam, 1986. Search in Google Scholar

[15] J. Droniou, R. Eymard, T. Gallouët, C. Guichard and R. Herbin, The gradient discretisation method: A framework for the discretisation and numerical analysis of linear and nonlinear elliptic and parabolic problems, preprint (2017), https://hal.archives-ouvertes.fr/hal-01382358, Version 4. Search in Google Scholar

[16] J. Droniou, R. Eymard, T. Gallouët and R. Herbin, Gradient schemes: A generic framework for the discretisation of linear, nonlinear and nonlocal elliptic and parabolic equations, Math. Models Methods Appl. Sci. 23 (2013), no. 13, 2395–2432. 10.1142/S0218202513500358Search in Google Scholar

[17] J. Droniou, R. Eymard and R. Herbin, Gradient schemes: Generic tools for the numerical analysis of diffusion equations, ESAIM Math. Model. Numer. Anal. 50 (2016), no. 3, 749–781. 10.1051/m2an/2015079Search in Google Scholar

[18] J. Droniou, J. Hennicker and R. Masson, Numerical analysis of a two-phase flow discrete fracture model, preprint (2016),https://arxiv.org/abs/1612.07373. 10.1007/s00211-018-0994-ySearch in Google Scholar

[19] J. Droniou and N. Nataraj, Improved L2 estimate for gradient schemes, and super-convergence of HMM and TPFA finite volume methods, preprint (2017), https://arxiv.org/abs/1602.07359; IMA J. Numer. Anal. (2017), DOI 10.1093/imanum/drx028. Search in Google Scholar

[20] J. Droniou, N. Nataraj and D. Shylaja, The gradient discretisation method for optimal control problems, with super-convergence for non-conforming finite elements and mixed-hybrid mimetic finite differences, SIAM J. Control Optim. 55 (2017), no. 6, 3640–3672. 10.1137/17M1117768Search in Google Scholar

[21] R. Eymard, C. Guichard and R. Herbin, Small-stencil 3d schemes for diffusive flows in porous media, ESAIM Math. Model. Numer. Anal. 46 (2012), 265–290. 10.1051/m2an/2011040Search in Google Scholar

[22] R. Eymard, C. Guichard, R. Herbin and R. Masson, Gradient schemes for two-phase flow in heterogeneous porous media and Richards equation, ZAMM Z. Angew. Math. Mech. 94 (2014), no. 7–8, 560–585. 10.1002/zamm.201200206Search in Google Scholar

[23] L. Gastaldi and R. Nochetto, Optimal L∞-error estimates for nonconforming and mixed finite element methods of lowest order, Numer. Math. 50 (1987), no. 5, 587–611. 10.1007/BF01408578Search in Google Scholar

[24] M. S. Gockenbach, Understanding and Implementing the Finite Element Method, Society for Industrial and Applied Mathematics, Philadelphia, 2006. 10.1137/1.9780898717846Search in Google Scholar

[25] R. Herbin and F. Hubert, Benchmark on discretization schemes for anisotropic diffusion problems on general grids, Finite Volumes for Complex Applications V, ISTE, London (2008), 659–692. Search in Google Scholar

[26] K. Krumbiegel and J. Pfefferer, Superconvergence for Neumann boundary control problems governed by semilinear elliptic equations, Comput. Optim. Appl. 61 (2015), no. 2, 373–408. 10.1007/s10589-014-9718-0Search in Google Scholar

[27] K. Kunisch and A. Rösch, Primal-dual active set strategy for a general class of constrained optimal control problems, SIAM J. Optim. 13 (2002), no. 2, 321–334. 10.1137/S1052623499358008Search in Google Scholar

[28] 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

[29] M. Mateos and A. Rösch, On saturation effects in the Neumann boundary control of elliptic optimal control problems, Comput. Optim. Appl. 49 (2011), no. 2, 359–378. 10.1007/s10589-009-9299-5Search in Google Scholar

[30] C. Meyer and A. Rösch, Superconvergence properties of optimal control problems, SIAM J. Control Optim. 43 (2004), no. 3, 970–985. 10.1137/S0363012903431608Search in Google Scholar

[31] D. Peaceman, Improved treatment of dispersion in numerical calculation of multidimensional miscible displacement, Soc. Pet. Eng. J. 6 (1966), no. 3, 213–216. 10.2118/1362-PASearch in Google Scholar

[32] F. Tröltzsch, Optimal Control of Partial Differential Equations, Grad. Stud. Math. 112, American Mathematical Society, Providence, 2010. 10.1090/gsm/112Search in Google Scholar

Received: 2017-05-10

Revised: 2017-09-13

Accepted: 2017-11-09

Published Online: 2017-12-05

Published in Print: 2018-10-01

Numerical Analysis for the Pure Neumann Control Problem Using the Gradient Discretisation Method

Abstract

1 Introduction

Remark 1.1 (Choice of the Adjoint State).

Remark 1.2 (Zero Average Constraint on the Source Term and Desired State).

Remark 1.3 (Non-homogeneous BCs).

2 GDM for Elliptic PDE with Neumann BCs

2.1 Gradient Discretisation and Gradient Scheme

Definition 2.1 (Gradient Discretisation for Neumann BCs).

Remark 2.2.

2.2 Examples of Gradient Discretisations

Conforming ℙ1 Finite Elements.

Non-conforming ℙ1 Finite Elements.

Mass-Lumped Non-conforming ℙ1 Finite Elements.

Hybrid Mixed Mimetic Method (HMM).

2.3 Error Estimates for the GDM for the Neumann Problem

Theorem 2.3 (Error Estimate for the GDM).

Remark 2.4 (Rates of Convergence).

Remark 2.5.

Proposition 2.6 (Stability of the GDM).

Proof.

3 GDM for the Control Problem and Main Results

3.1 GDM for the Optimal Control Problem

3.2 Basic Error Estimate for the GDM for the Control Problem

Theorem 3.1 (Control Estimate).

Proposition 3.2 (State and Adjoint Error Estimates).

Remark 3.3 (Rates of Convergence for the Control Problem).

3.3 Super-Convergence for Post-Processed Controls

Assumption A1 (Interpolation Operator).

Assumption A2.

Assumption A3 (Discrete Sobolev Imbedding).

Assumption A4.

Remark 3.4 (Choice of wT in Assumption A1).

Theorem 3.5 (Super-Convergence for Post-Processed Controls I).

Theorem 3.6 (Super-Convergence for Post-Processed Controls II).

Remark 3.7.

Corollary 3.8 (Super-Convergence for the State and Adjoint Variables).

3.4 Discussion on Post-Processed Controls

Lemma 3.9.

Proof.

Lemma 3.10 (Projection Formulas for the Controls).

Proof.

Remark 3.11.

Lemma 3.12 (Stability of the Discrete States).

Proof.

4 Numerical Experiments

4.1 A Modified Active Set Strategy

Remark 4.1.

4.2 Examples

4.2.1 Example 1: Ω=(0,1)2, ρ=10-4, a=-1, b=1

Conforming FEM.

Non-conforming FEM.

HMM Scheme.

Remark 4.2.

4.2.2 Example 2: A=100⁢Id, Ω=(0,1)2, ρ=10-2, a=-1, b=1

Conforming FEM.

Non-conforming FEM.

HMM Scheme.

4.2.3 Example 3: Ω=(0,1)2, ρ=10-4, a=-0.5, b=1

Conforming FEM.

Non-conforming FEM.

HMM Scheme.

A Appendix

Proof of Theorem 3.1.

Proof of Proposition 3.2.

Proof of Theorem 3.5.

Proof of Theorem 3.6.

Proof of Corollary 3.8.

References

Journal and Issue

Articles in the same Issue