Numerical analysis for the pure Neumann control problem using the gradient discretisation method

TL;DR

This paper presents a unified analysis of the gradient discretisation method for Neumann boundary control problems, deriving optimal error estimates and demonstrating super-convergence through numerical experiments across various numerical schemes.

Contribution

It introduces a general framework that applies to multiple numerical schemes for Neumann control problems, including new super-convergence results for post-processed control variables.

Findings

Optimal order error estimates for state, adjoint, and control variables.

Super-convergence of post-processed control variables.

Numerical experiments confirm theoretical convergence rates.

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, nonconforming and mimetic finite difference methods confirm the theoretical rates of convergence.

Figures

Tables

Table 1. Table 1. Example 1, conforming FEM

$h$	$U_a$	$f_a$	$Y_a$	ni
0.250000	0.002752$\times {10}^{-13}$	0.20699$\times {10}^{-14}$	-0.008396$\times {10}^{-13}$	2
0.125000	-0.008049$\times {10}^{-13}$	0.20912$\times {10}^{-14}$	-0.004684$\times {10}^{-13}$	3
0.062500	-0.001370$\times {10}^{-13}$	0.20548$\times {10}^{-14}$	0.010486$\times {10}^{-13}$	3
0.031250	-0.032432$\times {10}^{-13}$	0.21299$\times {10}^{-14}$	0.050725$\times {10}^{-13}$	3
0.015625	-0.917129$\times {10}^{-13}$	0.20367$\times {10}^{-14}$	-0.495753$\times {10}^{-13}$	3

Table 2. Table 2. Convergence results, Example 1, conforming FEM

$h$	${\mathrm{𝖾𝗋𝗋}}_{𝒟}\left(\overline{y}\right)$	Order	${\mathrm{𝖾𝗋𝗋}}_{𝒟}\left(\nabla \overline{y}\right)$	Order	${\mathrm{𝖾𝗋𝗋}}_{𝒟}\left(\overline{p}\right)$	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

Table 3. Table 3. Example 1, nc ℙ 1 subscript ℙ 1 \mathbb{P}_{1} FEM

$h$	$U_a$	$f_a$	$Y_a$	ni
0.250000	-0.003919$\times {10}^{-13}$	0.206991$\times {10}^{-14}$	0.000518$\times {10}^{-12}$	3
0.125000	-0.067706$\times {10}^{-13}$	0.209121$\times {10}^{-14}$	-0.003856$\times {10}^{-12}$	3
0.062500	0.030900$\times {10}^{-13}$	0.205478$\times {10}^{-14}$	0.017217$\times {10}^{-12}$	3
0.031250	-0.075427$\times {10}^{-13}$	0.212989$\times {10}^{-14}$	-0.041154$\times {10}^{-12}$	3
0.015625	-0.187208$\times {10}^{-13}$	0.203674$\times {10}^{-14}$	0.933499$\times {10}^{-12}$	3

Table 4. Table 4. Convergence results, Example 1, nc ℙ 1 subscript ℙ 1 \mathbb{P}_{1} FEM

$h$	${\mathrm{𝖾𝗋𝗋}}_{𝒟}\left(\overline{y}\right)$	Order	${\mathrm{𝖾𝗋𝗋}}_{𝒟}\left(\nabla \overline{y}\right)$	Order	${\mathrm{𝖾𝗋𝗋}}_{𝒟}\left(\overline{p}\right)$	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

Table 5. Table 5. Example 1, HMM

$h$	$U_a$	$f_a$	$Y_a$	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. Table 6. Convergence results, Example 1, HMM

$h$	${\mathrm{𝖾𝗋𝗋}}_{𝒟}\left(\overline{y}\right)$	Order	${\mathrm{𝖾𝗋𝗋}}_{𝒟}\left(\nabla \overline{y}\right)$	Order	${\mathrm{𝖾𝗋𝗋}}_{𝒟}\left(\overline{p}\right)$	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

Table 7. Table 7. Example 2, conforming FEM

$h$	$U_a$	$f_a$	$Y_a$	ni
0.250000	0.002280$\times {10}^{-11}$	0.209361$\times {10}^{-12}$	0.009117$\times {10}^{-11}$	2
0.125000	0.018065$\times {10}^{-11}$	0.209375$\times {10}^{-12}$	0.024496$\times {10}^{-11}$	3
0.062500	0.027564$\times {10}^{-11}$	0.209400$\times {10}^{-12}$	-0.012778$\times {10}^{-11}$	3
0.031250	-0.103755$\times {10}^{-11}$	0.209420$\times {10}^{-12}$	-0.028850$\times {10}^{-11}$	3
0.015625	-0.168277$\times {10}^{-11}$	0.209297$\times {10}^{-12}$	0.624160$\times {10}^{-11}$	3

Table 8. Table 8. Convergence results, Example 2, conforming FEM

$h$	${\mathrm{𝖾𝗋𝗋}}_{𝒟}\left(\overline{y}\right)$	Order	${\mathrm{𝖾𝗋𝗋}}_{𝒟}\left(\nabla \overline{y}\right)$	Order	${\mathrm{𝖾𝗋𝗋}}_{𝒟}\left(\overline{p}\right)$	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

Table 9. Table 9. Example 2, nc ℙ 1 subscript ℙ 1 \mathbb{P}_{1} FEM

$h$	$U_a$	$f_a$	$Y_a$	ni
0.250000	-0.000977$\times {10}^{-10}$	0.209361$\times {10}^{-12}$	0.000739 $\times {10}^{-10}$	3
0.125000	-0.006518$\times {10}^{-10}$	0.209375$\times {10}^{-12}$	-0.000960$\times {10}^{-10}$	3
0.062500	-0.004320$\times {10}^{-10}$	0.209400$\times {10}^{-12}$	0.002330$\times {10}^{-10}$	3
0.031250	-0.007236$\times {10}^{-10}$	0.209420 $\times {10}^{-12}$	0.054029 $\times {10}^{-10}$	3
0.015625	0.346321$\times {10}^{-10}$	0.209297$\times {10}^{-12}$	0.276914 $\times {10}^{-10}$	3

Table 10. Table 10. Convergence results, Example 2, nc ℙ 1 subscript ℙ 1 \mathbb{P}_{1} FEM

$h$	${\mathrm{𝖾𝗋𝗋}}_{𝒟}\left(\overline{y}\right)$	Order	${\mathrm{𝖾𝗋𝗋}}_{𝒟}\left(\nabla \overline{y}\right)$	Order	${\mathrm{𝖾𝗋𝗋}}_{𝒟}\left(\overline{p}\right)$	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

Table 11. Table 11. Example 2, HMM

$h$	$U_a$	$f_a$	$Y_a$	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. Table 12. Convergence results, Example 2, HMM

$h$	${\mathrm{𝖾𝗋𝗋}}_{𝒟}\left(\overline{y}\right)$	Order	${\mathrm{𝖾𝗋𝗋}}_{𝒟}\left(\nabla \overline{y}\right)$	Order	${\mathrm{𝖾𝗋𝗋}}_{𝒟}\left(\overline{p}\right)$	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

Table 13. Table 13. Example 3, conforming FEM

$h$	$U_a$	$f_a$	$Y_a$	ni
0.250000	0.0020160	-0.0020160	0.201602$\times {10}^{-6}$	4
0.125000	0.0055595	-0.0055595	0.555952$\times {10}^{-6}$	4
0.062500	-0.0004794	0.0004795	-0.047944$\times {10}^{-6}$	4
0.031250	0.0001470	-0.0001470	0.014705$\times {10}^{-6}$	5
0.015625	-0.0000136	0.0000136	-0.001362$\times {10}^{-6}$	5

Table 14. Table 14. Convergence results, Example 3, conforming FEM

$h$	${\mathrm{𝖾𝗋𝗋}}_{𝒟}\left(\overline{y}\right)$	Order	${\mathrm{𝖾𝗋𝗋}}_{𝒟}\left(\nabla \overline{y}\right)$	Order	${\mathrm{𝖾𝗋𝗋}}_{𝒟}\left(\overline{p}\right)$	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

Table 15. Table 15. Example 3, nc ℙ 1 subscript ℙ 1 \mathbb{P}_{1} FEM

$h$	$U_a$	$f_a$	$Y_a$	ni
0.250000	0.002016	-0.002016	0.0201601$\times {10}^{-5}$	4
0.125000	0.005560	-0.005559	-0.1301803$\times {10}^{-5}$	5
0.062500	-0.000480	0.000479	-0.0017424$\times {10}^{-5}$	5
0.031250	0.000147	-0.000147	0.0011093$\times {10}^{-5}$	5
0.015625	-0.000014	0.000014	-0.0001436$\times {10}^{-5}$	5

Table 16. Table 16. Convergence results, Example 3, nc ℙ 1 subscript ℙ 1 \mathbb{P}_{1} FEM

$h$	${\mathrm{𝖾𝗋𝗋}}_{𝒟}\left(\overline{y}\right)$	Order	${\mathrm{𝖾𝗋𝗋}}_{𝒟}\left(\nabla \overline{y}\right)$	Order	${\mathrm{𝖾𝗋𝗋}}_{𝒟}\left(\overline{p}\right)$	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

Table 17. Table 17. Example 3, HMM

$h$	$U_a$	$f_a$	$Y_a$	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. Table 18. Convergence results, Example 3, HMM

$h$	${\mathrm{𝖾𝗋𝗋}}_{𝒟}\left(\overline{y}\right)$	Order	${\mathrm{𝖾𝗋𝗋}}_{𝒟}\left(\nabla \overline{y}\right)$	Order	${\mathrm{𝖾𝗋𝗋}}_{𝒟}\left(\overline{p}\right)$	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

Equations

$$\underset{u\in {\mathcal{U}}_{\mathrm{ad}}}{\min }J(y,u)\text{ subject to }$$

$$-\mathrm{div}(A\nabla y)=u+f\text{mbox}in\Omega ,$$

$$\displaystyle\quad{A\nabla y\cdot\mathbf{n}_{\Omega}=0\,\,\mbox{ on $\partial\Omega$},\qquad\mathchoice{{\vbox{\hbox{$\textstyle-$ }}\kern-7.83337pt}}{{\vbox{\hbox{$\scriptstyle-$ }}\kern-5.90005pt}}{{\vbox{\hbox{$\scriptscriptstyle-$ }}\kern-4.75003pt}}{{\vbox{\hbox{$\scriptscriptstyle-$ }}\kern-4.25003pt}}\!\int_{\Omega}y({\bm{x}}){\>\rm d}{\bm{x}}=0.}$$

$$J(y,u):=\frac{1}{2}\Vert y-{\overline{y}}_d{\Vert }_{L^2(\Omega )}^2+\frac{\alpha }{2}\Vert u{\Vert }_{L^2(\Omega )}^2$$

$${\mathcal{U}}_{\mathrm{ad}}=\left\{u\in L^2(\Omega ):a\le u\le b\text{mbox}and\text{vbox}-{\int }_{\Omega }u(\mathbf{x})\mathrm{d}\mathbf{x}=0\right\},$$

$$a(y,w)={\int }_{\Omega }uw\mathrm{d}\mathbf{x},$$

$$a(\overline{y},w)=(\overline{u}+f,w)$$

$$a(z,\overline{p})=(\overline{y}-{\overline{y}}_d,z)$$

$$(\overline{p}+\alpha \overline{u},v-\overline{u})\ge 0$$

$${\mathcal{U}}_{\mathrm{ad}}=\left\{u\in L^2(\Omega ):a\le u\le b\text{mbox}and\text{vbox}-{\int }_{\Omega }(u+f)\mathrm{d}\mathbf{x}=0\right\}.$$

$${\mathcal{U}}_{\mathrm{ad}}^{\star }=\left\{u^{\star }\in L^2(\Omega ):a^{\star }\le u^{\star }\le b^{\star }\text{mbox}and\text{vbox}-{\int }_{\Omega }{u}^{\star }\mathrm{d}\mathbf{x}=0\right\}$$

$${\int }_{\Omega }f\mathrm{d}\mathbf{x}+{\int }_{\partial \Omega }g\mathrm{d}s(\mathbf{x})=0.$$

$$a(\overline{y},w)=(\overline{u}+f,w)+(g,\gamma (w){)}_{\partial }$$

$$a(z,\overline{p})=(\overline{y}-{\overline{y}}_d,z)$$

$$(\overline{p}+\alpha \overline{u},v-\overline{u})\ge 0$$

$$\Vert w{\Vert }_{\mathcal{D}}^2:=\Vert {\nabla }_{\mathcal{D}}w{\Vert }^2+{\left|\text{vbox}-{\int }_{\Omega }{\Pi }_{\mathcal{D}}w(\mathbf{x})\mathrm{d}\mathbf{x}\right|}^2$$

$$\left\{\begin{array}[]{llll}-{\rm div}(A\nabla\psi)=F\mbox{ in $\Omega$},\\ A\nabla\psi\cdot\mathbf{n}_{\Omega}=0\mbox{ on $\partial\Omega$}\end{array}\right.$$

$$\mbox{Find $\psi\in H^{1}_{\star}(\Omega)$ such that, for all $w\in H^{1}_{\star}(\Omega)$, }a(\psi,w)=\int_{\Omega}Fw{\>\rm d}{\bm{x}}.$$

$${\psi }_{\mathcal{D}}\in X_{\mathcal{D},\star }$$

$$a_{\mathcal{D}}({\psi }_{\mathcal{D}},w_{\mathcal{D}})={\int }_{\Omega }F{\Pi }_{\mathcal{D}}{w}_{\mathcal{D}}\mathrm{d}\mathbf{x},$$

$${\psi }_{\mathcal{D}}\in X_{\mathcal{D}}$$

$$a_{\mathcal{D}}({\psi }_{\mathcal{D}},w_{\mathcal{D}})+\rho \left(\text{vbox}-{\int }_{\Omega }{\Pi }_{\mathcal{D}}{\psi }_{\mathcal{D}}\mathrm{d}\mathbf{x}\right)\left(\text{vbox}-{\int }_{\Omega }{\Pi }_{\mathcal{D}}{w}_{\mathcal{D}}\mathrm{d}\mathbf{x}\right)={\int }_{\Omega }F{\Pi }_{\mathcal{D}}{w}_{\mathcal{D}}\mathrm{d}\mathbf{x}.$$

$$\nabla_{\mathcal{D}}v=\overline{\nabla}_{K}v+\frac{\sqrt{d}}{d_{K,{\sigma}}}\left(v_{\sigma}-v_{K}-\overline{\nabla}_{K}v\cdot(\overline{{\bm{x}}}_{\sigma}-{\bm{x}}_{K})\right)\mathbf{n}_{K,{\sigma}}\quad\mbox{ on $K$},$$

$${\overline{\nabla }}_{K}v=\frac{1}{|K|}\sum _{\sigma \in {\mathcal{E}}_K}|\sigma |v_{\sigma }{\mathbf{n}}_{K,\sigma }.$$

$$C_{\mathcal{D}}:=\underset{w\in X_{\mathcal{D}}\setminus \{0\}}{\max }\frac{\Vert {\Pi }_{\mathcal{D}}w\Vert }{\Vert w{\Vert }_{\mathcal{D}}}.$$

$$\forall \phi \in H^1(\Omega ),S_{\mathcal{D}}(\phi )=\underset{w\in X_{\mathcal{D}}}{\min }\left(\Vert {\Pi }_{\mathcal{D}}w-\phi \Vert +\Vert {\nabla }_{\mathcal{D}}w-\nabla \phi \Vert \right).$$

$$\forall \mathbf{\phi }\in H_0^{\mathrm{div}}(\Omega ),$$

$$W_{\mathcal{D}}(\mathbf{\phi })=\underset{w\in X_{\mathcal{D}}\setminus \{0\}}{\max }\frac{1}{\Vert w{\Vert }_{\mathcal{D}}}\left|{\int }_{\Omega }({\Pi }_{\mathcal{D}}w\mathrm{div}(\mathbf{\phi })+{\nabla }_{\mathcal{D}}w\cdot \mathbf{\phi })\mathrm{d}\mathbf{x}\right|,$$

$$\displaystyle X\lesssim Y\mbox{ means that }X\leq CY\mbox{ for some $C$ depending}$$

$$\displaystyle\mbox{only on $\Omega$, $A$ and an upper bound of $C_{\mathcal{D}}$}.$$

$$\Vert {\Pi }_{\mathcal{D}}{\psi }_{\mathcal{D}}-\psi \Vert +\Vert {\nabla }_{\mathcal{D}}{\psi }_{\mathcal{D}}-\nabla \psi \Vert \lesssim {\mathrm{WS}}_{\mathcal{D}}(\psi ),$$

$${\mathrm{WS}}_{\mathcal{D}}(\psi )=W_{\mathcal{D}}(A\nabla \psi )+S_{\mathcal{D}}(\psi ).$$

$$\Vert {\Pi }_{\mathcal{D}}{\tilde{\psi }}_{\mathcal{D}}-\tilde{\psi }\Vert +\Vert {\nabla }_{\mathcal{D}}{\tilde{\psi }}_{\mathcal{D}}-\nabla \tilde{\psi }\Vert \lesssim {\mathrm{WS}}_{\mathcal{D}}(\tilde{\psi })$$

$$\Vert {\Pi }_{\mathcal{D}}{\psi }_{\mathcal{D}}-\psi \Vert +\Vert {\nabla }_{\mathcal{D}}{\psi }_{\mathcal{D}}-\nabla \psi \Vert \lesssim {\mathrm{WS}}_{\mathcal{D}}(\tilde{\psi })={\mathrm{WS}}_{\mathcal{D}}(\psi ).$$

$$\Vert {\nabla }_{\mathcal{D}}{\psi }_{\mathcal{D}}\Vert \le \frac{C_{\mathcal{D}}}{\underline{a}}\Vert F\Vert \text{mbox}and\Vert {\Pi }_{\mathcal{D}}{\psi }_{\mathcal{D}}\Vert \le \frac{C_{\mathcal{D}}^2}{\underline{a}}\Vert F\Vert .$$

Full text

1. Introduction

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

[TABLE]

Here, $\Omega\subset{\mathbb{R}}^{d}\;(d\leq 3)$ is a bounded domain with boundary $\partial\Omega$, $\mathbf{n}_{\Omega}$ is the outer unit normal to $\Omega$ and \mathchoice{{\vbox{\hbox{\textstyle- }}\kern-7.83337pt}}{{\vbox{\hbox{\scriptstyle- }}\kern-5.90005pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.75003pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.25003pt}}\!\int_{\Omega}y({\bm{x}}){\>\rm d}{\bm{x}}:=\frac{1}{|\Omega|}\int_{\Omega}y({\bm{x}}){\>\rm d}{\bm{x}} denotes the average value of the function $y$ over $\Omega$ (here and in the rest of the paper, $|E|$ is the Lebesgue measure of a set $E\subset\mathbb{R}^{d}$). The cost functional, dependent on the state variable $y$ and the control variable $u$, is given by

[TABLE]

with $\alpha>0$. The desired state variable $\overline{y}_{d}\in L^{2}(\Omega)$ is chosen to satisfy \mathchoice{{\vbox{\hbox{\textstyle- }}\kern-7.83337pt}}{{\vbox{\hbox{\scriptstyle- }}\kern-5.90005pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.75003pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.25003pt}}\!\int_{\Omega}\overline{y}_{d}({\bm{x}}){\>\rm d}{\bm{x}}=0. The source term $f\in L^{2}(\Omega)$ also satisfies the zero average condition \mathchoice{{\vbox{\hbox{\textstyle- }}\kern-7.83337pt}}{{\vbox{\hbox{\scriptstyle- }}\kern-5.90005pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.75003pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.25003pt}}\!\int_{\Omega}f({\bm{x}}){\>\rm d}{\bm{x}}=0. The diffusion matrix $A:\Omega\rightarrow\mathcal{M}_{d}(\mathbb{R})$ is a measurable, bounded and uniformly elliptic matrix-valued function such that $A({\bm{x}})$ is symmetric for a.e. ${\bm{x}}\in\Omega$. Finally, the admissible set of controls $\mathcal{U}_{\rm ad}$ is the non-empty convex set defined by

[TABLE]

where $a$ and $b$ are constants in $[-\infty,+\infty]$ with $a<0<b$ (this condition is necessary to ensure that $\mathcal{U}_{\rm 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 $\mathcal{O}(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 $\overline{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, nonconforming 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 (BC) (we refer to [20, 30] and the references therein). For Dirichlet BC, 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 $\mathbb{P}_{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 nonlinear 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 BC and the properties needed to prove its convergence. Some classical examples of GDM are presented in Subsection 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 Subsections 3.2 and 3.3. The first super-convergence result provides a nearly quadratic convergence rate for a post-processed control. Under an $L^{\infty}$ 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 Subsection 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]. Subsection 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 Subsection 4.2, we present the results of some numerical experiments. The paper ends with an appendix, Section 5, where we prove the results stated in Subsections 3.2 and 3.3.

Before concluding this introduction, we discuss the optimality conditions for (1.1). For a given $u\in\mathcal{U}_{\rm ad}$, there exists a unique weak solution y(u)\in H^{1}_{\star}(\Omega):=\{w\in H^{1}(\Omega)\,:\mathchoice{{\vbox{\hbox{\textstyle- }}\kern-7.83337pt}}{{\vbox{\hbox{\scriptstyle- }}\kern-5.90005pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.75003pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.25003pt}}\!\int_{\Omega}w({\bm{x}}){\>\rm d}{\bm{x}}=0\} of (1.1b)–(1.1c). That is, for $u\in\mathcal{U}_{\rm ad}$, there exists a unique $y(u)\in H^{1}_{\star}(\Omega)$ such that for all ${w}\in H^{1}_{\star}(\Omega)$,

[TABLE]

where $a(\phi,\psi)=\int_{\Omega}A\nabla\phi\cdot\nabla\psi{\>\rm d}{\bm{x}}$ for all $\phi,\psi\in H^{1}(\Omega)$. The term $y(u)$ is the state associated with the control $u$.

In the following, the norm and scalar product in $L^{2}(\Omega)$ (or $L^{2}(\Omega)^{d}$ for vector-valued functions) are denoted by $\|\cdot\|$ and $(\cdot,\cdot)$. The convex control problem (1.1) has a unique solution $(\overline{y},\overline{u})\in H^{1}_{\star}(\Omega)\times\mathcal{U}_{\rm ad}$ and there exists a co-state $\overline{p}\in H^{1}(\Omega)$ such that the triplet $(\overline{y},\overline{p},\overline{u})\in H^{1}_{\star}(\Omega)\times H^{1}(\Omega)\times\mathcal{U}_{\rm ad}$ satisfies the KKT optimality conditions [28, Chapter 2]:

[TABLE]

Remark 1.1** (Choice of the adjoint state).**

Several co-states satisfy the optimality conditions (1.5), as $\overline{p}$ is only determined up to an additive constant by (1.5). 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:

1. (1)

Design an algorithm (the modified active set algorithm of Subsection 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, 2. (2)

Fix the average of the continuous co-state $\overline{p}$ to be the same as the average of the discrete co-state obtained above.

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

Remark 1.2** (Zero average constraint on the source term and desired state).**

• (i)

If we consider (1.1) without the constraint \mathchoice{{\vbox{\hbox{\textstyle- }}\kern-7.83337pt}}{{\vbox{\hbox{\scriptstyle- }}\kern-5.90005pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.75003pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.25003pt}}\!\int_{\Omega}f({\bm{x}}){\>\rm d}{\bm{x}}=0 on the source term, the set of admissible controls needs to be modified into

[TABLE]

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 $-{\rm div}(A\nabla y)=u^{\star}+f^{\star}$ with u^{\star}=u+\mathchoice{{\vbox{\hbox{\textstyle- }}\kern-7.83337pt}}{{\vbox{\hbox{\scriptstyle- }}\kern-5.90005pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.75003pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.25003pt}}\!\int_{\Omega}f{\>\rm d}{\bm{x}} and f^{\star}=f-\mathchoice{{\vbox{\hbox{\textstyle- }}\kern-7.83337pt}}{{\vbox{\hbox{\scriptstyle- }}\kern-5.90005pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.75003pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.25003pt}}\!\int_{\Omega}f{\>\rm d}{\bm{x}}. Then, \mathchoice{{\vbox{\hbox{\textstyle- }}\kern-7.83337pt}}{{\vbox{\hbox{\scriptstyle- }}\kern-5.90005pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.75003pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.25003pt}}\!\int_{\Omega}f^{\star}{\>\rm d}{\bm{x}}=0 and $u^{\star}\in\mathcal{U}_{\rm ad}^{\star}$ where

[TABLE]

with a^{\star}=a+\mathchoice{{\vbox{\hbox{\textstyle- }}\kern-7.83337pt}}{{\vbox{\hbox{\scriptstyle- }}\kern-5.90005pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.75003pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.25003pt}}\!\int_{\Omega}f{\>\rm d}{\bm{x}} and b^{\star}=b+\mathchoice{{\vbox{\hbox{\textstyle- }}\kern-7.83337pt}}{{\vbox{\hbox{\scriptstyle- }}\kern-5.90005pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.75003pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.25003pt}}\!\int_{\Omega}f{\>\rm d}{\bm{x}}.

• (ii)

If the desired state $\overline{y}_{d}\in L^{2}(\Omega)$ is such that \mathchoice{{\vbox{\hbox{\textstyle- }}\kern-7.83337pt}}{{\vbox{\hbox{\scriptstyle- }}\kern-5.90005pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.75003pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.25003pt}}\!\int_{\Omega}\overline{y}_{d}({\bm{x}}){\>\rm d}{\bm{x}}=:m\neq 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 $\overline{y}_{d}$. In that case, working with $y-m$ and $\overline{y}_{d}-m$ instead of $y$ and $\overline{y}_{d}$ brings back to the original formulation (1.1) with a desired state $\overline{y}_{d}-m$ having a zero average.

Remark 1.3** (Non-homogeneous BCs).**

The study of second order distributed control problem (1.1) with non-homogeneous boundary conditions $A\nabla y\cdot\mathbf{n}_{\Omega}=g$ on $\partial\Omega$ (with $g\in L^{2}(\partial\Omega)$) follows in a similar way. In this case, the source terms and boundary condition are supposed to satisfy the compatibility condition

[TABLE]

The controls are still taken in $\mathcal{U}_{\rm ad}$ defined by (1.3) and the KKT optimality condition is [28]: Seek $(\overline{y},\overline{p},\overline{u})\in H^{1}_{\star}(\Omega)\times H^{1}(\Omega)\times\mathcal{U}_{\rm ad}$ such that

[TABLE]

where $\gamma:H^{1}(\Omega)\rightarrow L^{2}(\partial\Omega)$ is the trace operator and $(\cdot,\cdot)_{\partial}$ is the inner product in $L^{2}(\partial\Omega)$.

2. GDM for elliptic PDE with Neumann BC

The gradient discretisation method (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 BC is given in [15, Definition 3.1]. The following extends this definition by demanding the existence of the element $1_{\mathcal{D}}$ 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 boundary conditions).**

A gradient discretisation (GD) for homogeneous Neumann boundary conditions is given by ${\mathcal{D}}=(X_{{\mathcal{D}}},\Pi_{\mathcal{D}},\nabla_{\mathcal{D}})$ such that

• •

$X_{{\mathcal{D}}}$* is a finite dimensional vector space on $\mathbb{R}$.*

• •

$\Pi_{\mathcal{D}}:X_{{\mathcal{D}}}\rightarrow L^{2}(\Omega)$* and $\nabla_{\mathcal{D}}:X_{{\mathcal{D}}}\rightarrow L^{2}(\Omega)^{d}$ are linear mappings.*

• •

The quantity

[TABLE]

is a norm on $X_{{\mathcal{D}}}$.

• •

There exists $1_{\mathcal{D}}\in X_{\mathcal{D}}$ such that $\Pi_{\mathcal{D}}1_{\mathcal{D}}=1$ on $\Omega$ and $\nabla_{\mathcal{D}}1_{\mathcal{D}}=0$ on $\Omega$.

The flexibility of the GDM analysis framework comes from the wide possible range of choices for $(X_{\mathcal{D}},\Pi_{\mathcal{D}},\nabla_{\mathcal{D}})$. Each of these choices correspond to a particular numerical scheme (see Subsection 2.2 for a few examples). The space $X_{\mathcal{D}}$ represents the degrees of freedom (unknowns) of the method; a vector in $X_{\mathcal{D}}$ gathers values for such unknowns. The operators $\Pi_{\mathcal{D}}$ and $\nabla_{\mathcal{D}}$ reconstruct, from a set of values of these unknowns, a scalar (resp. vector) function on the entire set $\Omega$. 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 ${\mathcal{D}}$.

If $F\in L^{2}(\Omega)$ is such that \mathchoice{{\vbox{\hbox{\textstyle- }}\kern-7.83337pt}}{{\vbox{\hbox{\scriptstyle- }}\kern-5.90005pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.75003pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.25003pt}}\!\int_{\Omega}F({\bm{x}}){\>\rm d}{\bm{x}}=0, the weak formulation of the Neumann boundary value problem

[TABLE]

is given by

[TABLE]

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

[TABLE]

where $\displaystyle a_{{\mathcal{D}}}(\phi_{{\mathcal{D}}},z_{{\mathcal{D}}})=\int_{\Omega}A\nabla_{{\mathcal{D}}}\phi_{{\mathcal{D}}}\cdot\nabla_{{\mathcal{D}}}{z}_{{\mathcal{D}}}{\>\rm d}{\bm{x}}$, for all $\phi_{{\mathcal{D}}},z_{{\mathcal{D}}}\in X_{{\mathcal{D}}}$, and X_{{\mathcal{D}},\star}=\{w_{\mathcal{D}}\in X_{\mathcal{D}}\,:\,\mathchoice{{\vbox{\hbox{\textstyle- }}\kern-7.83337pt}}{{\vbox{\hbox{\scriptstyle- }}\kern-5.90005pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.75003pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.25003pt}}\!\int_{\Omega}\Pi_{\mathcal{D}}w_{\mathcal{D}}{\>\rm d}{\bm{x}}=0\}.

Remark 2.2**.**

As for the continuous formulation (2.3), using the element $1_{\mathcal{D}}\in X_{\mathcal{D}}$ actually enables us to consider in (2.4) test functions $w_{\mathcal{D}}$ in $X_{{\mathcal{D}}}$, rather than just $X_{{\mathcal{D}},\star}$. The simplest technique to achieve this is to use a quadratic penalty method [24, Chapter 11]. For any $\rho>0$, (2.4) can be shown equivalent to

[TABLE]

Indeed, considering $w_{\mathcal{D}}=1_{\mathcal{D}}$ in (2.5) shows that the solution to this problem belongs to $X_{{\mathcal{D}},\star}$ 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 applicables on generic meshes (not just triangular/quadrangular meshes) as encountered for example in reservoir engineering applications.

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

Conforming $\mathbb{P}_{1}$ finite elements. The simplest gradient discretisation is perhaps obtained by considering conforming $\mathbb{P}_{1}$ finite elements. The mesh is made of triangles (in 2D) or tetrahedra (in 3D), with no hanging nodes. Each $v_{\mathcal{D}}\in X_{{\mathcal{D}}}$ is a vector of values at the vertices of the mesh (the standard unknowns of conforming $\mathbb{P}_{1}$ finite elements). $\Pi_{\mathcal{D}}v_{\mathcal{D}}$ is the continuous piecewise linear function on the mesh which takes these values at the vertices, and $\nabla_{\mathcal{D}}v_{\mathcal{D}}=\nabla(\Pi_{\mathcal{D}}v_{\mathcal{D}})$. Then (2.4) is the standard $\mathbb{P}_{1}$ finite element scheme for (2.2).

Non-conforming $\mathbb{P}_{1}$ finite elements. As above, the mesh is made of conforming triangles or tetrahedra. Each $v_{\mathcal{D}}\in X_{{\mathcal{D}}}$ is a vector of values at the centers of mass of the edges/faces, $\Pi_{\mathcal{D}}v_{\mathcal{D}}$ is the piecewise linear function on the mesh which takes these values at these centers of mass, and $\nabla_{\mathcal{D}}v_{\mathcal{D}}=\nabla_{\mathcal{T}}(\Pi_{\mathcal{D}}v_{\mathcal{D}})$ is the broken gradient of $\Pi_{\mathcal{D}}v_{\mathcal{D}}$. In that case, (2.4) gives the non-conforming $\mathbb{P}_{1}$ finite element approximation of (2.2).

Mass-lumped non-conforming $\mathbb{P}_{1}$ finite elements. Still considering a conforming triangular/tetrahedral mesh, the space $X_{\mathcal{D}}$ and gradient reconstruction $\nabla_{\mathcal{D}}$ are identical to those of the non-conforming $\mathbb{P}_{1}$ finite elements described above, but the function reconstruction is modified to be piecewise constant. For each edge/face ${\sigma}$ of the mesh, we consider the diamond $D_{\sigma}$ around ${\sigma}$ constructed from the edge/face and the one or two cell centers on each side (see Figure 1). Then, for $v=(v_{\sigma})_{{\sigma}}\in X_{\mathcal{D}}$, the reconstructed function $\Pi_{\mathcal{D}}v$ is the piecewise constant function on the diamonds, equal to $v_{\sigma}$ on $D_{\sigma}$ for all edge/face ${\sigma}$.

Hybrid mixed mimetic method (HMM). We now consider a generic mesh ${\mathcal{T}}$ (not necessarily triangular/tetrahedral) with one point ${\bm{x}}_{K}$ chosen in each cell $K\in{\mathcal{T}}$ such that $K$ is strictly star-shaped with respect to ${\bm{x}}_{K}$; see Figure 2 for some notations. A vector $v\in X_{{\mathcal{D}}}$ is made of cell $(v_{K})_{K}$ and face $(v_{\sigma})_{\sigma}$ values, and the operator $\Pi_{\mathcal{D}}$ reconstructs a piecewise constant function from the cell values: for any cell $K$, $\Pi_{\mathcal{D}}v=v_{K}$ on $K$. The gradient reconstruction $\nabla_{\mathcal{D}}$ is built in two pieces: a consistent gradient $\overline{\nabla}_{K}$ constant over the cell and stabilisation terms constant over the half-diamonds $D_{K,{\sigma}}$ (and akin to the remainders of first-order Taylor expansions between the cell and face values). For any cell $K$ and any face ${\sigma}$ of $K$, we set

[TABLE]

where $d_{K,{\sigma}}$ is the orthogonal distance between ${\bm{x}}_{K}$ and ${\sigma}$, $\overline{{\bm{x}}}_{\sigma}$ is the center of mass of ${\sigma}$, $\mathbf{n}_{K,{\sigma}}$ is the outer normal to $K$ on ${\sigma}$ and, denoting by ${{{\mathcal{E}}}_{K}}$ the set of faces of $K$,

[TABLE]

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].

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_{\mathcal{D}}$, which ensures the coercivity of the method.

[TABLE]

The second quantity is the interpolation error $S_{\mathcal{D}}$, which measures what is called, in the GDM framework, the GD-consistency of ${\mathcal{D}}$.

[TABLE]

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

[TABLE]

where $H^{{\rm div}}_{0}(\Omega)=\{{\bm{\varphi}}\in L^{2}(\Omega)^{d}\,:\,{\rm div}({\bm{\varphi}})\in L^{2}(\Omega)\,,\;\gamma_{\mathbf{n}}({\bm{\varphi}})=0\}$ with $\gamma_{\mathbf{n}}$ being the normal trace of ${\bm{\varphi}}$ on $\partial\Omega$.

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,

[TABLE]

Theorem 2.3** (Error estimate for the GDM).**

Let ${\mathcal{D}}$ be a GD in the sense of Definition 2.1, let $\psi$ be the solution to (2.3), and let $\psi_{\mathcal{D}}$ be the solution to (2.4). Then

[TABLE]

where

[TABLE]

Remark 2.4** (Rates of convergence).**

For all classical low order methods based on meshes (such as $\mathbb{P}_{1}$ conforming and non-conforming finite element methods, finite volume methods, etc.), if $A$ is Lipschitz continuous and $\psi\in H^{2}(\Omega)$ then $\mathcal{O}(h)$ estimates can be obtained for $W_{\mathcal{D}}(A\nabla\psi)$ and $S_{\mathcal{D}}(\psi)$ [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 $\psi$ and $\Pi_{\mathcal{D}}\psi_{\mathcal{D}}$ with \mathchoice{{\vbox{\hbox{\textstyle- }}\kern-7.83337pt}}{{\vbox{\hbox{\scriptstyle- }}\kern-5.90005pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.75003pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.25003pt}}\!\int_{\Omega}\Pi_{\mathcal{D}}\psi_{\mathcal{D}}{\>\rm d}{\bm{x}}=\mathchoice{{\vbox{\hbox{\textstyle- }}\kern-7.83337pt}}{{\vbox{\hbox{\scriptstyle- }}\kern-5.90005pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.75003pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.25003pt}}\!\int_{\Omega}\psi{\>\rm d}{\bm{x}}. In this case, the estimate (2.10) can be obtained by considering the translation of $\psi_{\mathcal{D}}$ and $\psi$. Set $\widetilde{\psi}_{\mathcal{D}}=\psi_{\mathcal{D}}-c1_{\mathcal{D}}$ and $\widetilde{\psi}=\psi-c1,$ where c=\mathchoice{{\vbox{\hbox{\textstyle- }}\kern-7.83337pt}}{{\vbox{\hbox{\scriptstyle- }}\kern-5.90005pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.75003pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.25003pt}}\!\int_{\Omega}\Pi_{\mathcal{D}}\psi_{\mathcal{D}}{\>\rm d}{\bm{x}}=\mathchoice{{\vbox{\hbox{\textstyle- }}\kern-7.83337pt}}{{\vbox{\hbox{\scriptstyle- }}\kern-5.90005pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.75003pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.25003pt}}\!\int_{\Omega}\psi{\>\rm d}{\bm{x}} and $1$ is the constant function. Using Definition 2.1, we find $\Pi_{\mathcal{D}}\widetilde{\psi}_{\mathcal{D}}=\Pi_{\mathcal{D}}{\psi}_{\mathcal{D}}-c$, $\nabla_{\mathcal{D}}\widetilde{\psi}_{\mathcal{D}}=\nabla_{\mathcal{D}}{\psi}_{\mathcal{D}}$ and $\nabla\widetilde{\psi}=\nabla{\psi}$. This gives \mathchoice{{\vbox{\hbox{\textstyle- }}\kern-7.83337pt}}{{\vbox{\hbox{\scriptstyle- }}\kern-5.90005pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.75003pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.25003pt}}\!\int_{\Omega}\Pi_{\mathcal{D}}\widetilde{\psi}_{\mathcal{D}}{\>\rm d}{\bm{x}}=\mathchoice{{\vbox{\hbox{\textstyle- }}\kern-7.83337pt}}{{\vbox{\hbox{\scriptstyle- }}\kern-5.90005pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.75003pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.25003pt}}\!\int_{\Omega}\widetilde{\psi}{\>\rm d}{\bm{x}}=0. Applying Theorem 2.3,

[TABLE]

which implies

[TABLE]

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

Proposition 2.6** (Stability of the GDM).**

Let $\underline{a}$ be a coercivity constant of $A$. If $\psi_{\mathcal{D}}$ is the solution to the gradient scheme (2.4), then

[TABLE]

Proof.

Take $w_{\mathcal{D}}=\psi_{\mathcal{D}}$ in (2.4) and use the definition of $C_{\mathcal{D}}$ to write

[TABLE]

Since \mathchoice{{\vbox{\hbox{\textstyle- }}\kern-7.83337pt}}{{\vbox{\hbox{\scriptstyle- }}\kern-5.90005pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.75003pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.25003pt}}\!\int_{\Omega}\Pi_{\mathcal{D}}\psi_{\mathcal{D}}{\>\rm d}{\bm{x}}=0, recalling the Definition (2.1) of $\|\cdot\|_{{\mathcal{D}}}$ shows that $\|\psi_{\mathcal{D}}\|_{{\mathcal{D}}}=\|\nabla_{\mathcal{D}}\psi_{\mathcal{D}}\|$ and the proof of first estimate is complete. The second estimate follows from the definition of $C_{\mathcal{D}}$. ∎

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 Subsections 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 Subsection 3.4.

3.1. GDM for the optimal control problem

Let ${\mathcal{D}}$ be a GD as in Definition 2.1. The space $\mathcal{U}_{h}$ is defined as the space of piecewise constant functions on a mesh ${\mathcal{T}}$ of $\Omega$. The space $\mathcal{U}_{{\rm ad},h}=\mathcal{U}_{\rm ad}\cap\mathcal{U}_{h}$ is a finite dimensional subset of $\mathcal{U}_{\rm ad}$. A gradient scheme for (1.5) consists in seeking $(\overline{y}_{{\mathcal{D}}},\overline{p}_{{\mathcal{D}}},\overline{u}_{h})\in X_{{\mathcal{D}},\star}\times X_{{\mathcal{D}}}\times\mathcal{U}_{{\rm ad},h}$, such that

[TABLE]

As in the continuous KKT conditions (1.5), these equations do not define $\overline{p}_{\mathcal{D}}$ uniquely. One possible constraint that fixes a unique $\overline{p}_{\mathcal{D}}$ is described in Lemma 3.10. This particular choice ensures a simple projection relation between $\overline{p}_{{\mathcal{D}}}$ and $\overline{u}_{h}$.

Two projection operators play a major role throughout the paper: the orthogonal projection on piecewise constant functions on ${\mathcal{T}}$, namely $\mathbfcal{P}_{\!{\mathcal{T}}}:L^{1}(\Omega)\to\mathcal{U}_{h}$ and the cut-off function $P_{[a,b]}:\mathbb{R}\to[a,b]$. They are defined as

[TABLE]

[TABLE]

3.2. Basic error estimate for the GDM for the control problem

To state the error estimates, we define the projection error $E_{h}$ by

[TABLE]

The proofs of the following basic error estimates are provided in Section 5. 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 ${\mathcal{D}}$ be a GD, $(\overline{y},\overline{p},\overline{u})$ be a solution to (1.5) and $(\overline{y}_{\mathcal{D}},\overline{p}_{\mathcal{D}},\overline{u}_{h})$ be a solution to (3.1) such that \mathchoice{{\vbox{\hbox{\textstyle- }}\kern-7.83337pt}}{{\vbox{\hbox{\scriptstyle- }}\kern-5.90005pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.75003pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.25003pt}}\!\int_{\Omega}\Pi_{\mathcal{D}}\overline{p}_{\mathcal{D}}{\>\rm d}{\bm{x}}=\mathchoice{{\vbox{\hbox{\textstyle- }}\kern-7.83337pt}}{{\vbox{\hbox{\scriptstyle- }}\kern-5.90005pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.75003pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.25003pt}}\!\int_{\Omega}\overline{p}{\>\rm d}{\bm{x}}. Then, recalling (2.9), (2.11) and (3.4), there exists a constant $C$ depends only on $\alpha$ such that

[TABLE]

Proposition 3.2** (State and adjoint error estimates).**

Let ${\mathcal{D}}$ be a GD, $(\overline{y},\overline{p},\overline{u})$ be a solution to (1.5) and $(\overline{y}_{\mathcal{D}},\overline{p}_{\mathcal{D}},\overline{u}_{h})$ be a solution to (3.1). Assume that \mathchoice{{\vbox{\hbox{\textstyle- }}\kern-7.83337pt}}{{\vbox{\hbox{\scriptstyle- }}\kern-5.90005pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.75003pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.25003pt}}\!\int_{\Omega}\Pi_{\mathcal{D}}\overline{p}_{\mathcal{D}}{\>\rm d}{\bm{x}}=\mathchoice{{\vbox{\hbox{\textstyle- }}\kern-7.83337pt}}{{\vbox{\hbox{\scriptstyle- }}\kern-5.90005pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.75003pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.25003pt}}\!\int_{\Omega}\overline{p}{\>\rm d}{\bm{x}}. Then the following error estimates hold:

[TABLE]

Remark 3.3** (Rates of convergence for the control problem).**

As in Remark 2.4, if $A$ is Lipschitz continuous and $(\overline{y},\overline{p},\overline{u})\in H^{2}(\Omega)^{2}\times H^{1}(\Omega)$ 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 5.

We make here the following assumptions.

• (A1)

[Interpolation operator] For each $w\in H^{2}(\Omega)$, there exists $w_{\mathcal{T}}\in L^{2}(\Omega)$ such that:

i) If $w\in H^{2}(\Omega)$ solves $-{\rm div}(A\nabla w)=g\in H^{1}(\Omega)$, and $w_{\mathcal{D}}$ is the solution to the corresponding GS with \mathchoice{{\vbox{\hbox{\textstyle- }}\kern-7.83337pt}}{{\vbox{\hbox{\scriptstyle- }}\kern-5.90005pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.75003pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.25003pt}}\!\int_{\Omega}\Pi_{\mathcal{D}}w_{\mathcal{D}}{\>\rm d}{\bm{x}}=\mathchoice{{\vbox{\hbox{\textstyle- }}\kern-7.83337pt}}{{\vbox{\hbox{\scriptstyle- }}\kern-5.90005pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.75003pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.25003pt}}\!\int_{\Omega}w{\>\rm d}{\bm{x}}, then

[TABLE]

ii) For any $w\in H^{2}(\Omega)$, it holds

[TABLE]

[TABLE]

• (A2)

The estimate $\|\Pi_{\mathcal{D}}v_{\mathcal{D}}-\mathbfcal{P}_{\!{\mathcal{T}}}(\Pi_{\mathcal{D}}v_{\mathcal{D}})\|\lesssim h\|\nabla_{\mathcal{D}}v_{\mathcal{D}}\|$ holds for any $v_{\mathcal{D}}\in X_{{\mathcal{D}}}$.

• (A3)

[Discrete Sobolev imbedding] For all $v_{\mathcal{D}}\in X_{{\mathcal{D}}}$, it holds

[TABLE]

where $2^{*}$ is a Sobolev exponent of 2, that is, $2^{*}\in[2,\infty)$ if $d=2$, and $2^{*}=6$ if $d=3$.

Let

[TABLE]

and ${\mathcal{T}}_{1}={\mathcal{T}}\setminus{\mathcal{T}}_{2}$. That is, ${\mathcal{T}}_{1}$ is the set of cells where $\overline{u}$ crosses at least one constraint $a$ or $b$. For $i=1,2$, we let $\Omega_{i,{\mathcal{T}}}={\rm int}(\cup_{K\in{\mathcal{T}}_{i}}\overline{K})$. The space $W^{1,\infty}({\mathcal{T}}_{1})$ is the usual broken Sobolev space, endowed with its broken norm. Our last assumption is:

• (A4)

$|\Omega_{1,{\mathcal{T}}}|\lesssim h$ and $\overline{u}_{|\Omega_{1,{\mathcal{T}}}}\in W^{1,\infty}({\mathcal{T}}_{1})$.

Possible choices of the mapping $w\mapsto w_{\mathcal{T}}$, depending on the considered gradient discretisation (that is, the considered numerical method), is discussed in Remark 3.4 below. Note that the assumptions (A1)–(A4) are similar to that in [20] with $X_{{\mathcal{D}},0}$ substituted by $X_{{\mathcal{D}}}$, and an additional average condition in (A1). We also refer to [20] for a detailed discussion on (A1)–(A4).

Assuming $\overline{p}\in H^{2}(\Omega)$ (see Theorem 3.5) and letting $\overline{p}_{\mathcal{T}}$ be defined as in (A1), the post-processed continuous and discrete controls are given by

[TABLE]

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

Remark 3.4** (Choice of $w_{\mathcal{T}}$ in (A1)).**

When $\Pi_{\mathcal{D}}$ is a piecewise linear reconstruction, the super-convergence result (3.8) usually holds with $w_{\mathcal{T}}=w$. This is for example well-known for conforming and non-conforming $\mathbb{P}_{1}$ FE method. When $\Pi_{\mathcal{D}}v_{\mathcal{D}}$ is piecewise constant on ${\mathcal{T}}$ for all $v_{\mathcal{D}}\in X_{{\mathcal{D}}}$, 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_{\mathcal{T}}({\bm{x}})=\mathchoice{{\vbox{\hbox{\textstyle- }}\kern-7.83337pt}}{{\vbox{\hbox{\scriptstyle- }}\kern-5.90005pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.75003pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.25003pt}}\!\int_{K}w({\bm{x}}){\>\rm d}{\bm{x}} for all ${\bm{x}}\in K$ and all $K\in{\mathcal{T}}$ (or, equivalently at order $\mathcal{O}(h^{2})$, $w_{\mathcal{T}}({\bm{x}})=w(\overline{{\bm{x}}}_{K})$ with $\overline{{\bm{x}}}_{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, $\widetilde{u}=P_{[a,b]}(-\alpha^{-1}\overline{p})$.

• •

For mimetic finite difference methods, $\widetilde{u}_{|K}=P_{[a,b]}(-\alpha^{-1}\overline{p}(\overline{{\bm{x}}}_{K}))$ for all $K\in{\mathcal{T}}$.

For $K\in{\mathcal{T}}$ and $\overline{{\bm{x}}}_{K}$ the center of mass of $K$, let $\rho_{K}=\max\{r>0\,:\,B(\overline{{\bm{x}}}_{K},r)\subset K\}$ be the maximal radius of balls centred at $\overline{{\bm{x}}}_{K}$ and included in $K$. Fixing $\eta>0$ such that

[TABLE]

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

[TABLE]

Theorem 3.5** (Super-convergence for post-processed controls I).**

Let ${\mathcal{D}}$ be a GD and ${\mathcal{T}}$ be a mesh. Assume that

• •

(A1)–(A4)* hold,*

• •

$\overline{y}$* and $\overline{p}$ belong to $H^{2}(\Omega)$,*

• •

$\overline{y}_{d}$* and $f$ belong to $H^{1}(\Omega)$,*

and let $\widetilde{u}$, $\widetilde{u}_{h}$ be the post-processed controls defined by (3.11) where $\overline{p}$ and $\overline{p}_{\mathcal{D}}$ are chosen such that \mathchoice{{\vbox{\hbox{\textstyle- }}\kern-7.83337pt}}{{\vbox{\hbox{\scriptstyle- }}\kern-5.90005pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.75003pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.25003pt}}\!\int_{\Omega}\Pi_{\mathcal{D}}\overline{p}_{\mathcal{D}}{\>\rm d}{\bm{x}}=\mathchoice{{\vbox{\hbox{\textstyle- }}\kern-7.83337pt}}{{\vbox{\hbox{\scriptstyle- }}\kern-5.90005pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.75003pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.25003pt}}\!\int_{\Omega}\overline{p}{\>\rm d}{\bm{x}}. Then there exists $C$ depending only on $\alpha$ in (1.2) such that

[TABLE]

[TABLE]

Theorem 3.6** (Super-convergence for post-processed controls II).**

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

[TABLE]

Then there exists $C$ depending only on $\alpha$ and $\delta$ such that

[TABLE]

Remark 3.7**.**

For most methods, assumption (3.14) is satisfied if the mesh is quasi-uniform (see [23] for conforming and non-conforming $\mathbb{P}_{1}$ finite element method, and [20, Theorem 7.1] for HMM methods).

Corollary 3.8** (Super-convergence for the state and adjoint variables).**

Under the assumptions of Theorem 3.5, the following error estimates hold, with $C$ depending only on $\alpha$:

[TABLE]

where $\overline{y}_{\mathcal{T}}$ and $\overline{p}_{\mathcal{T}}$ are defined as in (A1), and $r=2-\frac{1}{2^{*}}$.

Under the assumptions of Theorem 3.6, (3.16) and (3.17) hold with $r=2$ and $C$ depending only $\alpha$ and $\delta$.

3.4. Discussion on post-processed controls

In this section, 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 ${\rm WS}_{\mathcal{D}}(\varphi)\lesssim h$ for all $\varphi\in H^{2}(\Omega)$ (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.23) and (3.27). Hence, the super-convergence result makes sense: since $\overline{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(h^{2})$ convergence result.

Lemma 3.9**.**

Let $-\infty\leq a<0<b\leq\infty$ and $\phi\in L^{1}(\Omega)$. Define $\Gamma:\mathbb{R}\rightarrow\mathbb{R}$ by

[TABLE]

where $P_{[a,b]}$ is given by (3.3). Set $m=a-\mbox{\rm ess sup}(\phi)\in[-\infty,+\infty)$ and $M=b-\mbox{\rm ess inf}(\phi)\in(-\infty,+\infty]$. Then we have the following.

1. (1)

$\Gamma$* is Lipschitz continuous.* 2. (2)

$\lim\limits_{c\to m}\Gamma(c)=a|\Omega|$, $\lim\limits_{c\to M}\Gamma(c)=b|\Omega|$, and there is $c^{\star}\in(m,M)$ such that $\Gamma(c^{\star})=0$. 3. (3)

If $\phi\in H^{1}(\Omega)$, then for any compact interval $Q$ in $(m,M)$, there exists $\rho_{Q}>0$ such that if $c,c^{\prime}\in Q$ with $c<c^{\prime}$, then

[TABLE]

As a consequence, the real number $c^{\star}$ 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 $(c_{n})$ be a sequence in $\mathbb{R}$ such that $c_{n}\rightarrow M$ as $n\rightarrow\infty$. By definition of $M$, this implies $P_{[a,b]}(\phi+c_{n})\rightarrow b$ a.e on $\Omega$. Let $(c_{n})$ be bounded below by $R$ and note that $\phi+R\in L^{1}(\Omega)$. Moreover, for $s\in\mathbb{R}$, $a<0<b$ implies $P_{[a,b]}(s)\geq\min(s,0)$ so $P_{[a,b]}(\phi+c_{n})\geq\min(\phi+c_{n},0)\geq\min(\phi+R,0)\in L^{1}(\Omega)$. By Fatou’s Lemma,

[TABLE]

which gives $b|\Omega|\leq\liminf_{n\rightarrow\infty}\Gamma(c_{n})$. Since $\Gamma(c_{n})\leq b|\Omega|$ (because $P_{[a,b]}(s)\leq b$), we infer that $\lim_{n\to\infty}\Gamma(c_{n})=b|\Omega|$, and thus that $\lim_{c\to M}\Gamma(c)=b|\Omega|$. In a similar way, we deduce that $\lim_{c\to m}\Gamma(c)=a|\Omega|$.

The existence of $c^{\star}$ such that $\Gamma(c^{\star})=0$ then follows from the intermediate value theorem and $\lim_{c\to m}\Gamma(c)=a|\Omega|<0<b|\Omega|=\lim_{c\to M}\Gamma(c)$.

We now assume that $\phi\in H^{1}(\Omega)$ and we turn to Item 3. For a.e $c\in\mathbb{R}$, $\Gamma^{{}^{\prime}}(c)=\int_{\Omega}\mathbbm{1}_{(a,b)}(\phi({\bm{x}})+c){\>\rm d}{\bm{x}},$ where $\mathbbm{1}_{(a,b)}$ is the characteristic function of $(a,b)$. Define $\Theta(c)=\int_{\Omega}\mathbbm{1}_{(a,b)}(\phi({\bm{x}})+c){\>\rm d}{\bm{x}},$ for all $c\in\mathbb{R}$. We claim that

• •

$\Theta$ is lower semi-continuous,

• •

$\forall c\in(m,M)$, $\Theta(c)>0.$

To prove that $\Theta$ is lower semi-continuous, let $c_{n}\rightarrow c$ as $n\rightarrow\infty$. Since $\mathbbm{1}_{(a,b)}$ is lower semi-continuous on $\mathbb{R}$, we have, for all ${\bm{x}}\in\Omega$,

[TABLE]

Applying Fatou’s Lemma,

[TABLE]

Hence, $\Theta$ is lower semi-continuous. We now show that $\Theta>0$ on $(m,M).$ Let $c\in(m,M)$. Then $I=(a-c,b-c)\cap(\mbox{ess inf }\phi,\mbox{ess sup }\phi)$ is an interval of positive length, since $a-c<\mbox{ess sup }\phi$ and $b-c>\mbox{ess inf }\phi$. The set $W_{I,c}=\{{\bm{x}}:\phi({\bm{x}})\in I\}$ has a non-zero measure because $\phi\in H^{1}(\Omega)$ and $\Omega$ is connected. To see this, let $\alpha<\beta$ be the endpoints of $I$ and assume that $\phi\in H^{1}(\Omega)$ takes some values less than $\alpha$ on a non-null set, some values greater than $\beta$ on a non-null set, but that $W_{I,c}$ is a null set. Then $P_{[\alpha,\beta]}(\phi)\in H^{1}(\Omega)$ exactly takes the values $\alpha$ and $\beta$ (outside a set of zero measure). Hence $\nabla P_{[\alpha,\beta]}(\phi)=\mathbbm{1}_{[\alpha,\beta]}(\phi)\nabla\phi=0$ and $P_{[\alpha,\beta]}(\phi)$ should be constant, since $\Omega$ is connected, which is a contradiction. Thus, $W_{I,c}$ has a non-zero measure. Since $\{{\bm{x}}:\phi({\bm{x}})+c\in(a,b)\}\supseteq W_{I,c}$, this gives $\Theta(c)\geq|W_{I,c}|>0$.

Coming back to Item 3, let $Q$ be a compact interval in $(m,M)$. We know that $\Theta>0$ on $Q$ and $\Theta$ is lower semi-continuous. Hence $\Theta$ reaches its minimum on $Q$ and $\inf_{Q}\Theta=\Theta(c_{0})>0$ for some $c_{0}\in Q$. Since $\Gamma^{{}^{\prime}}=\Theta$ a.e, $\Gamma^{{}^{\prime}}\geq\inf_{Q}\Theta$ a.e on $Q$ and, $\Gamma$ being Lipschitz and $[c,c^{\prime}]\subset Q$, we infer

[TABLE]

which establishes (3.18). The uniqueness of $c^{\star}$ such that $\Gamma(c^{\star})=0$ follows from this inequality, which shows that $\Gamma$ is strictly increasing on $(m,M)$. ∎

Lemma 3.10** (Projection formulas for the controls).**

If $\overline{p}\in H^{1}(\Omega)$ is a co-state and $\overline{c}\in\mathbb{R}$ is such that \mathchoice{{\vbox{\hbox{\textstyle- }}\kern-7.83337pt}}{{\vbox{\hbox{\scriptstyle- }}\kern-5.90005pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.75003pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.25003pt}}\!\int_{\Omega}P_{[a,b]}(-\frac{1}{\alpha}\overline{p}({\bm{x}})+\overline{c}){\>\rm d}{\bm{x}}=0, then the continuous optimal control $\overline{u}$ in (1.5) can be expressed in terms of the projection formula

[TABLE]

If ${\mathcal{D}}$ is a GD and $\overline{p}_{\mathcal{D}}$ is chosen such that

[TABLE]

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

[TABLE]

Proof.

Set $\widetilde{p}=\overline{p}-\alpha\overline{c}$. Clearly, $\widetilde{p}\in H^{1}(\Omega)$. From the optimality condition for the control problem (1.5c), we deduce that

[TABLE]

since $\int_{\Omega}\overline{u}{\>\rm d}{\bm{x}}=\int_{\Omega}v{\>\rm d}{\bm{x}}=0$. Set $U=P_{[a,b]}(-\alpha^{-1}\widetilde{p})$ i.e,

[TABLE]

It is then straightforward to see that $U\in\mathcal{U}_{\rm ad}$, i.e, $U\in[a,b]$ and \mathchoice{{\vbox{\hbox{\textstyle- }}\kern-7.83337pt}}{{\vbox{\hbox{\scriptstyle- }}\kern-5.90005pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.75003pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.25003pt}}\!\int_{\Omega}U({\bm{x}})d{\bm{x}}=0 (by choice of $\overline{c}$). Then, using the definitions of $\Omega_{+}$, $\Omega_{0}$ and $\Omega_{-}$, since $v\geq a=U$ on $\Omega_{+}$ and $v\leq b=U$ on $\Omega_{-}$,

[TABLE]

Recall that the optimality condition is nothing but a characterisation of the $L^{2}(\Omega)$ orthogonal projection of $-\alpha^{-1}\widetilde{p}$ on $\mathcal{U}_{\rm ad}$ and, as such, defines a unique element $\overline{u}$ of $\mathcal{U}_{\rm ad}$. We just proved that $U=P_{[a,b]}(-\alpha^{-1}\widetilde{p})$ satisfies this optimality condition, which shows that it is equal to $\overline{u}$. The proof of (3.19) is complete.

The second relation follows in a similar way by noticing that, since controls are piecewise-constants on ${\mathcal{T}}$, (3.1c) is equivalent to $(\mathbfcal{P}_{\!{\mathcal{T}}}(\Pi_{\mathcal{D}}\overline{p}_{\mathcal{D}})+\alpha\overline{u}_{h},v_{h}-\overline{u}_{h})\geq 0$, for all $v_{h}\in\mathcal{U}_{{\rm ad},h}$. We also notice that, by definition of $\mathbfcal{P}_{\!{\mathcal{T}}}$ and the assumption (3.20), $P_{[a,b]}\left(\mathbfcal{P}_{\!{\mathcal{T}}}\left(-\frac{1}{\alpha}\Pi_{\mathcal{D}}\overline{p}_{\mathcal{D}}\right)\right)\in\mathcal{U}_{{\rm ad},h}.$ ∎

Remark 3.11**.**

There is at least one adjoint $\overline{p}_{\mathcal{D}}$ such that (3.20) is satisfied: start from any adjoint $\overline{p}^{0}_{\mathcal{D}}$ and, by applying Lemma 3.9 (Item 2) to $\phi=\mathbfcal{P}_{\!{\mathcal{T}}}(-\alpha^{-1}\Pi_{\mathcal{D}}\overline{p}_{\mathcal{D}}^{0})$ and by noticing that $\phi+c^{\star}=\mathbfcal{P}_{\!{\mathcal{T}}}(-\alpha^{-1}\Pi_{\mathcal{D}}\overline{p}_{\mathcal{D}}^{0})+c^{\star}$, find $c^{\star}$ such that $\overline{p}_{\mathcal{D}}=\overline{p}_{\mathcal{D}}^{0}-\alpha c^{\star}1_{\mathcal{D}}$ satisfies (3.20).

Since the discrete co-state $\overline{p}_{\mathcal{D}}$ is a computable quantity, its average is easier to fix than the average of the non-computable $\overline{p}$. Hence, the projection relation (3.21) is the most natural choice to express the discrete control $\overline{u}_{h}$ in terms of the discrete adjoint variable. This is the choice made in the modified active set strategy presented in Subsection 4.1. Once this choice is made, since $\overline{p}$ must have the same average as $\Pi_{\mathcal{D}}\overline{p}_{\mathcal{D}}$ for $\widetilde{u}$ defined in (3.11) to satisfy super-convergence estimates, it is clear that $P_{[a,b]}(-\frac{1}{\alpha}\overline{p})$ will not have a zero average in general. Hence, if we want to express the continuous control in terms of $\overline{p}$, we need to offset this $\overline{p}$ by the correct $\overline{c}$, as stated in the lemma.

Lemma 3.12** (Stability of the discrete states).**

Let ${\mathcal{D}}$ be a GD, $(\overline{y},\overline{p},\overline{u})$ be a solution to (1.5) and $(\overline{y}_{\mathcal{D}},\overline{p}_{\mathcal{D}},\overline{u}_{h})$ be a solution to (3.1). Assume that \mathchoice{{\vbox{\hbox{\textstyle- }}\kern-7.83337pt}}{{\vbox{\hbox{\scriptstyle- }}\kern-5.90005pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.75003pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.25003pt}}\!\int_{\Omega}\Pi_{\mathcal{D}}\overline{p}_{\mathcal{D}}{\>\rm d}{\bm{x}}=\mathchoice{{\vbox{\hbox{\textstyle- }}\kern-7.83337pt}}{{\vbox{\hbox{\scriptstyle- }}\kern-5.90005pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.75003pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.25003pt}}\!\int_{\Omega}\overline{p}{\>\rm d}{\bm{x}}. Then

[TABLE]

Proof.

Let $\phi\in H^{{\rm div}}_{0}(\Omega)$. Taking $w=0$ in (2.7), we get $S_{\mathcal{D}}(\phi)\leq\|\phi\|+\|\nabla\phi\|.$ By Cauchy-Schwarz inequality, using (2.6) and (2.1), for $w\in X_{\mathcal{D}}$,

[TABLE]

With (2.8), this implies ${\rm WS}_{\mathcal{D}}(\phi)\leq C_{\mathcal{D}}\left(\|{\rm div}(\phi)\|+\|\phi\|\right)$.

Therefore, for $A\nabla\psi\in H^{{\rm div}}_{0}(\Omega)$, recalling the definition (2.11) of ${\rm WS}_{\mathcal{D}}$, we have

[TABLE]

Using Proposition 3.2, Theorem 3.1, this shows that

[TABLE]

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

In the rest of this section, we establish $\mathcal{O}(h)$ estimates between the controls $\overline{u}$, $\overline{u}_{h}$ and their post-processed versions $\widetilde{u}$, $\widetilde{u}_{h}$. These estimates justify that the post-processed controls are indeed meaningful quantities.

Using (3.11), (3.21), the Lipschitz continuity of $P_{[a,b]}$, (A2) and Lemma 3.12, we have the following estimate between $\overline{u}_{h}$ and $\widetilde{u}_{h}$:

[TABLE]

Let us now turn to estimating $\overline{u}-\widetilde{u}$. The co-state $\overline{p}\in H^{1}(\Omega)$ in (1.5) is still taken such that \mathchoice{{\vbox{\hbox{\textstyle- }}\kern-7.83337pt}}{{\vbox{\hbox{\scriptstyle- }}\kern-5.90005pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.75003pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.25003pt}}\!\int_{\Omega}\overline{p}{\>\rm d}{\bm{x}}=\mathchoice{{\vbox{\hbox{\textstyle- }}\kern-7.83337pt}}{{\vbox{\hbox{\scriptstyle- }}\kern-5.90005pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.75003pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.25003pt}}\!\int_{\Omega}\Pi_{\mathcal{D}}\overline{p}_{\mathcal{D}}{\>\rm d}{\bm{x}}. From Lemma 3.9, it follows that there exists a unique constant $\overline{c}\in(m,M)$ such that $\int_{\Omega}P_{[a,b]}(-\frac{1}{\alpha}\overline{p}+\overline{c}){\>\rm d}{\bm{x}}=0$, where $m$ and $M$ are defined as in Lemma 3.9. Using Lemma 3.10 and recalling (3.11),

[TABLE]

Starting from (3.24) and using the Lipschitz continuity of $P_{[a,b]}$, the assumption ${\rm WS}_{\mathcal{D}}(\varphi)\lesssim h$, Corollary 3.8 and Proposition 3.2, we get a constant $C$ depending only on $\alpha$, $f$, $\overline{y}_{d}$, $\overline{p},$ $\overline{y}$ and $\overline{u}$ such that

[TABLE]

To estimate the last term in (3.25), recall the definition of $\Gamma(c)$ from Lemma 3.9 for $\phi=\alpha^{-1}\overline{p}$.

[TABLE]

By choice of $\overline{c}$, $\Gamma(\overline{c})=0$. From Lemma 3.10, the choice of $\overline{p}_{\mathcal{D}}$ shows that

[TABLE]

Let $q_{\mathcal{D}}$ be the solution to (3.1b) with source term $\overline{y}-\overline{y}_{d}$ (that is, the solution to the GS for the equation (1.5b) satisfied by $\overline{p}$) such that \mathchoice{{\vbox{\hbox{\textstyle- }}\kern-7.83337pt}}{{\vbox{\hbox{\scriptstyle- }}\kern-5.90005pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.75003pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.25003pt}}\!\int_{\Omega}\Pi_{\mathcal{D}}q_{\mathcal{D}}{\>\rm d}{\bm{x}}=\mathchoice{{\vbox{\hbox{\textstyle- }}\kern-7.83337pt}}{{\vbox{\hbox{\scriptstyle- }}\kern-5.90005pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.75003pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.25003pt}}\!\int_{\Omega}\Pi_{\mathcal{D}}\overline{p}_{\mathcal{D}}{\>\rm d}{\bm{x}}=\mathchoice{{\vbox{\hbox{\textstyle- }}\kern-7.83337pt}}{{\vbox{\hbox{\scriptstyle- }}\kern-5.90005pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.75003pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.25003pt}}\!\int_{\Omega}\overline{p}{\>\rm d}{\bm{x}}. Using the Lipschitz continuity of $P_{[a,b]}$, the Cauchy–Schwarz inequality, the triangle inequality, Remark 2.5, Proposition 2.6, (A2), Theorem 2.3 and Lemma 3.12, we obtain

[TABLE]

Let $m,M$ be as in Lemma 3.9 for $\phi=\alpha^{-1}\overline{p}$. Relation (3.26) shows that $a|\Omega|<\Gamma(0)<b|\Omega|$ if $h$ is small enough; hence, in this case, $0\in(m,M)$. There is therefore a compact interval $Q$ in $(m,M)$ depending only on $\overline{p}$ such that 0 and $\overline{c}$ belong to $Q$. Without loss of generality, we can assume that $\overline{c}\geq 0$. A use of Lemma 3.9 leads to

[TABLE]

where $\rho_{Q}>0$. This implies $0\leq\overline{c}\lesssim\alpha^{-1}h/\rho_{Q}$, using (3.26) and the fact that $\Gamma(\overline{c})=0$. Combining this with (3.25), we infer that

[TABLE]

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.20) is highlighted in Lemma 3.10: we have the projection relation (3.21) 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 $\overline{p}_{\mathcal{D}}$ such that (3.20) holds, the KKT optimality conditions (3.1) can be rewritten as: Seek $(\overline{y}_{{\mathcal{D}}},\overline{p}_{{\mathcal{D}}},\overline{u}_{h})\in X_{{\mathcal{D}}}\times X_{{\mathcal{D}}}\times\mathcal{U}_{{\rm ad},h}$, such that

[TABLE]

where $\rho>0$ is constant.

Set $\bar{\mu}_{h}=-(\alpha^{-1}\Pi_{\mathcal{D}}\overline{p}_{\mathcal{D}}+\overline{u}_{h})$. As the original active set strategy [32], the modified active set strategy is an iterative algorithm. As initial guesses, two arbitrary functions, $u_{h}^{0},\mu_{h}^{0}$ are chosen. In the $n$th step of the algorithm, we define the set of active and inactive restrictions by

[TABLE]

[TABLE]

If $\max\left(\frac{\|u_{h}^{n}-u_{h}^{n-1}\|_{L^{\infty}(\Omega)}}{\|u_{h}^{n-1}\|_{L^{\infty}(\Omega)}},\frac{\|\Pi_{\mathcal{D}}p_{\mathcal{D}}^{n}-\Pi_{\mathcal{D}}p_{\mathcal{D}}^{n-1}\|_{L^{\infty}(\Omega)}}{\|\Pi_{\mathcal{D}}p_{\mathcal{D}}^{n-1}\|_{L^{\infty}(\Omega)}}\right)\leq 10^{-10}$, then we terminate the algorithm. In this case, we notice that the relative $L^{\infty}$ difference between $\Pi_{\mathcal{D}}y_{\mathcal{D}}^{n-1}$ and $\Pi_{\mathcal{D}}y_{\mathcal{D}}^{n}$ is less than $10^{-6}$ for all examples in Section 4.2. Else we find $y_{\mathcal{D}}^{n},{p}_{\mathcal{D}}^{n}$ and $u_{h}^{n}$ solution to the system

[TABLE]

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

[TABLE]

where

[TABLE]

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

[TABLE]

where $\mathbbm{1}_{a,h}^{n}$ and $\mathbbm{1}_{b,h}^{n}$ denote the characteristic functions of the sets $A^{n}_{a,h}$ and $A^{n}_{b,h}$ respecively.

Remark 4.1**.**

The convergence analysis of the proposed algorithm is a plan for future study. However, if $(\Pi_{\mathcal{D}}y_{\mathcal{D}}^{n},\Pi_{\mathcal{D}}p_{\mathcal{D}}^{n})$ converges weakly to $(\Pi_{\mathcal{D}}\overline{y}_{\mathcal{D}},\Pi_{\mathcal{D}}\overline{p}_{\mathcal{D}})$ in $H^{1}(\Omega)$ and $u_{h}^{n}$ converges to $\overline{u}_{h}$ in $L^{2}(\Omega)$, then the solution to (4.3) converges to the solution of (4.1) as $n\rightarrow\infty$.

4.2. Examples

In this section, 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$\mathbb{P}_{1}$FE) method and hybrid mimetic mixed (HMM) method. All three are GDMs with gradient discretisations with bounds on $C_{\mathcal{D}}$, order $h$ estimate on ${\rm WS}_{\mathcal{D}}$, and satisfying assumptions (A1)–(A4), and (3.14) 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 Subsection 4.1 with zero as an initial guess for both $u$ and $\mu$. Here, $U_{a}$ and $Y_{a}$ denote the average values of the computed control $\overline{u}_{h}$ and the reconstructed state solution $\Pi_{\mathcal{D}}\overline{y}_{\mathcal{D}}$ respectively. Let ni denote the number of iterations required for the convergence of the modified active set algorithm, and $f_{a}$ denote the numerical average of the source term $f$ calculated using the same quadrature rule as in the implementation of the schemes, i.e,

[TABLE]

where $\overline{{\bm{x}}}_{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 $f_{a}$ can be considered to be equal to zero, up to quadrature error. The relative errors are denoted by

[TABLE]

[TABLE]

[TABLE]

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

[TABLE]

where $\overline{c}$ is chosen to ensure that \mathchoice{{\vbox{\hbox{\textstyle- }}\kern-7.83337pt}}{{\vbox{\hbox{\scriptstyle- }}\kern-5.90005pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.75003pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.25003pt}}\!\int_{\Omega}\overline{u}{\>\rm d}{\bm{x}}=0. The matrix-valued function is given by $A={\rm Id}$ unless otherwise specified. The source term $f$ and the desired state $\overline{y}_{d}$ are then computed using

[TABLE]

4.2.1. Example 1 :

$\Omega$ = $\left(0,1\right)^{2}$, $\rho=10^{-4}$, $a=-1$, $b=1$.

We here consider the computational domain $\Omega$ = $\left(0,1\right)^{2}$. We have $\overline{p}(x,y)=-\overline{p}(1-x,y)$ and, since $P_{[-1,1]}$ is odd, $P_{[-1,1]}(-\overline{p})(1-x,y)=-P_{[-1,1]}(-\overline{p})(x,y)$. Integrating this relation over $\Omega$ shows that $P_{[-1,1]}(-\overline{p})$ has a zero average and thus, by Lemma 3.10, that $\overline{c}=0$. We thus see that $\overline{u}=P_{[-1,1]}\left(-\overline{p}\right)$.

Conforming FEM: The discrete solution is computed on a family of uniform grids with mesh sizes $h=\frac{1}{2^{i}},i=2,\ldots,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.

Non-Conforming FEM: For comparison, we compute the solutions of the nc$\mathbb{P}_{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.

HMM scheme: This scheme was tested on a series of regular triangular meshes from [25] where the points $\mathcal{P}$ (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 $L^{2}$ 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 $L^{2}$ norm, together with their orders of convergence, are presented in Table 6.

For all three methods (conforming $\mathbb{P}_{1}$ FE, nc$\mathbb{P}_{1}$ FE and HMM), the theoretical rates of convergence are confirmed by the numerical outputs. Without post-processing, an $\mathcal{O}(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 $\overline{y}_{\mathcal{T}}$ and $\overline{p}_{\mathcal{T}}$ in $\mathsf{err}_{\mathcal{D}}(\overline{y})$ and $\mathsf{err}_{\mathcal{D}}(\overline{p})$).

Remark 4.2**.**

As seen in Table 5, the modified active set algorithm converges in very few iterations if $\rho=10^{-4}$. We however found that, if $\rho=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 $\rho$.

4.2.2. Example 2 :

$A=100Id$, $\Omega$ = $\left(0,1\right)^{2}$, $\rho=10^{-2}$, $a=-1$, $b=1$.

In this subsection, we present some numerical results for the control problem defined on the unit square domain $\Omega=\left(0,1\right)^{2}$ and $A=100Id$. As explained in Example 1, $a=-1$ and $b=1$ imply $\overline{c}=0$.

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.

**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 $\rho=1$.

4.2.3. Example 3 :

$\Omega=(0,1)^{2}$, $\rho=10^{-4},$ $a=-0.5,\,b=1$. In this case, since $P_{[a,b]}$ is no longer odd, $P_{[a,b]}(-\overline{p})$ no longer has a zero average and, to compute $\mathsf{err}_{\mathcal{D}}(\overline{u})$, we need to find $\overline{c}$ such that \mathchoice{{\vbox{\hbox{\textstyle- }}\kern-7.83337pt}}{{\vbox{\hbox{\scriptstyle- }}\kern-5.90005pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.75003pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.25003pt}}\!\int_{\Omega}P_{[a,b]}(-\overline{p}+\overline{c}){\>\rm d}{\bm{x}}=0. This $\overline{c}$ can be found by a bisection method, by computing the averages on a very thin mesh and bisecting until we find a proper $\overline{c}$. Using a mesh of size $h=0.00195$, we find $c\approx-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 $U_{a}$ and $f_{a}$ in Table 13). Here, we observed that the modified active set algorithm converges only when $\rho\leq 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.

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

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 $\rho\leq 10^{-1}$.

5. 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:

[TABLE]

where the co-state $p_{{\mathcal{D}}}(\overline{u})$ is chosen such that \mathchoice{{\vbox{\hbox{\textstyle- }}\kern-7.83337pt}}{{\vbox{\hbox{\scriptstyle- }}\kern-5.90005pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.75003pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.25003pt}}\!\int_{\Omega}\Pi_{\mathcal{D}}p_{\mathcal{D}}(\overline{u}){\>\rm d}{\bm{x}}=\mathchoice{{\vbox{\hbox{\textstyle- }}\kern-7.83337pt}}{{\vbox{\hbox{\scriptstyle- }}\kern-5.90005pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.75003pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.25003pt}}\!\int_{\Omega}\Pi_{\mathcal{D}}\overline{p}_{\mathcal{D}}{\>\rm d}{\bm{x}}=\mathchoice{{\vbox{\hbox{\textstyle- }}\kern-7.83337pt}}{{\vbox{\hbox{\scriptstyle- }}\kern-5.90005pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.75003pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.25003pt}}\!\int_{\Omega}\overline{p}{\>\rm d}{\bm{x}}. For Neumann boundary conditions, this particular choice is essential as it ensures that $p_{\mathcal{D}}(\overline{u})-\overline{p}_{\mathcal{D}}\in X_{{\mathcal{D}},\star}$ can be used as a test function $w_{\mathcal{D}}$ in (3.1a) and (5.1a). Recalling that $\mathbfcal{P}_{\!{\mathcal{T}}}$ is the orthogonal projection on piecewise constant functions on ${\mathcal{T}}$, we obtain $\mathbfcal{P}_{\!{\mathcal{T}}}(\mathcal{U}_{\rm ad})\subset\mathcal{U}_{h}$. Also, for $u\in\mathcal{U}_{\rm ad}$ and $K\in{\mathcal{T}}$, \mathbfcal{P}_{\!{\mathcal{T}}}u_{|K}=\mathchoice{{\vbox{\hbox{\textstyle- }}\kern-7.83337pt}}{{\vbox{\hbox{\scriptstyle- }}\kern-5.90005pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.75003pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.25003pt}}\!\int_{K}u{\>\rm d}{\bm{x}}\in[a,b] and, using (1.3),

[TABLE]

Hence, $\mathbfcal{P}_{\!{\mathcal{T}}}(\mathcal{U}_{\rm ad})\subset\mathcal{U}_{{\rm ad},h}$.

Set $P_{{\mathcal{D}},\alpha}(\overline{u})=\alpha^{-1}\Pi_{\mathcal{D}}p_{{\mathcal{D}}}(\overline{u})$,  $\overline{P}_{{\mathcal{D}},\alpha}=\alpha^{-1}\Pi_{\mathcal{D}}\overline{p}_{{\mathcal{D}}}$ and $\overline{P}_{\alpha}=\alpha^{-1}\overline{p}$. Since $\overline{u}_{h}\in\mathcal{U}_{{\rm ad},h}\subset\mathcal{U}_{\rm ad}$ and $\mathbfcal{P}_{\!{\mathcal{T}}}\overline{u}\in\mathcal{U}_{{\rm ad},h}$, from the optimality conditions ((1.5c) and (3.1c)),

[TABLE]

Adding these two inequalities yields

[TABLE]

By orthogonality property of $\mathbfcal{P}_{\!{\mathcal{T}}}$ we have $(\overline{u}_{h},\overline{u}-\mathbfcal{P}_{\!{\mathcal{T}}}\overline{u})=0$ and $(\mathbfcal{P}_{\!{\mathcal{T}}}\overline{P}_{\alpha},\overline{u}-\mathbfcal{P}_{\!{\mathcal{T}}}\overline{u})=0$. Therefore, the first term in the right-hand side of (5.2) can be re-cast as

[TABLE]

By Cauchy-Schwarz inequality, the first term on the right hand side of (5.3) is estimated as

[TABLE]

Equation (5.1b) shows that $p_{\mathcal{D}}(\overline{u})$ is the solution of the GS corresponding to the adjoint problem (1.5b), whose solution is $\overline{p}$. Therefore, using the fact that \mathchoice{{\vbox{\hbox{\textstyle- }}\kern-7.83337pt}}{{\vbox{\hbox{\scriptstyle- }}\kern-5.90005pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.75003pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.25003pt}}\!\int_{\Omega}\Pi_{\mathcal{D}}p_{\mathcal{D}}(\overline{u}){\>\rm d}{\bm{x}}=\mathchoice{{\vbox{\hbox{\textstyle- }}\kern-7.83337pt}}{{\vbox{\hbox{\scriptstyle- }}\kern-5.90005pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.75003pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.25003pt}}\!\int_{\Omega}\overline{p}{\>\rm d}{\bm{x}} (note that the specific relation between the continuous and discrete co-states is essential here), by Theorem 2.3,

[TABLE]

Hence, using the Cauchy–Schwarz inequality,

[TABLE]

Using the definitions of $C_{\mathcal{D}}$, $\|\cdot\|_{{\mathcal{D}}}$ and the fact that $p_{\mathcal{D}}(\overline{u})-\overline{p}_{\mathcal{D}}\in X_{{\mathcal{D}},\star}$, we find that

[TABLE]

By writing the difference of (5.1b) and (3.1b) we see that $p_{\mathcal{D}}(\overline{u})-\overline{p}_{\mathcal{D}}$ is the solution to the GS (2.4) with source term $F=\overline{y}-\Pi_{\mathcal{D}}\overline{y}_{\mathcal{D}}$. i.e, for all $z_{\mathcal{D}}\in X_{{\mathcal{D}}}$

[TABLE]

Choose $z_{\mathcal{D}}=p_{\mathcal{D}}(\overline{u})-\overline{p}_{\mathcal{D}}$ in the above equality and use it in (5.7) to obtain

[TABLE]

As a consequence,

[TABLE]

Use Theorem 2.3 with $\psi=\overline{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 (5.3) as

[TABLE]

where $C_{1}$ depends only on $\Omega$, $A$ and an upper bound of $C_{\mathcal{D}}$. Plugging (5.4), (5.6) and (5.8) in (5.3) yields

[TABLE]

where $C_{2}$ is the hidden constant in (5.6). Let us turn to the second term in the right-hand side of (5.2). From (3.1b) and (5.1b), for all $z_{\mathcal{D}}\in X_{{\mathcal{D}}}$,

[TABLE]

Also, from (3.1a) and (5.1a), for all $w_{\mathcal{D}}\in X_{{\mathcal{D}},\star}$,

[TABLE]

Choose $z_{\mathcal{D}}=\overline{y}_{\mathcal{D}}-y_{\mathcal{D}}(\overline{u})\in X_{{\mathcal{D}}}$ in (5.10), $w_{\mathcal{D}}=\overline{p}_{\mathcal{D}}-p_{\mathcal{D}}(\overline{u})\in X_{{\mathcal{D}},\star}$ in (5.11), use the symmetry of $a_{\mathcal{D}}(\cdot,\cdot)$, Theorem 2.3 with $\psi=\overline{y}$ and Young’s inequality to obtain

[TABLE]

where $C_{3}$ only depends on $\Omega$, $A$ and an upper bound of $C_{\mathcal{D}}$. The last term in the right hand side of (5.2) can be estimated using (5.5) and Young’s inequality:

[TABLE]

where $C_{4}$ only depends on $\Omega$, $A$, $\alpha$ and an upper bound of $C_{\mathcal{D}}$. Substitute (5.9), (5.12) and (5.13) into (5.2), apply the Young’s inequality and $\sqrt{\sum_{i}a_{i}^{2}}\leq\sum_{i}a_{i}$ to complete the proof. ∎

Proof of Proposition 3.2.

An application of triangle inequality yields

[TABLE]

Subtracting (3.1a) and (5.1a), and using the stability property of GS (Proposition 2.6), the first two terms in the right hand side of the above inequality can be estimated as

[TABLE]

The last two terms on the right hand side of the (5.14) are estimated using Theorem 2.3 as

[TABLE]

A combination of the above two results yields the error estimates (3.6) for the state variable. A use of \mathchoice{{\vbox{\hbox{\textstyle- }}\kern-7.83337pt}}{{\vbox{\hbox{\scriptstyle- }}\kern-5.90005pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.75003pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.25003pt}}\!\int_{\Omega}\Pi_{\mathcal{D}}p_{\mathcal{D}}(\overline{u}){\>\rm d}{\bm{x}}=\mathchoice{{\vbox{\hbox{\textstyle- }}\kern-7.83337pt}}{{\vbox{\hbox{\scriptstyle- }}\kern-5.90005pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.75003pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.25003pt}}\!\int_{\Omega}\Pi_{\mathcal{D}}\overline{p}_{\mathcal{D}}{\>\rm d}{\bm{x}} 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\in L^{2}(\Omega)$, let $p_{\mathcal{D}}^{*}(g)\in X_{{\mathcal{D}}}$ solve

[TABLE]

where $y_{\mathcal{D}}(g)$ is given by (5.1a) with $\overline{u}$ replaced by $g$. We fix $p_{\mathcal{D}}^{*}(g)$ by imposing \mathchoice{{\vbox{\hbox{\textstyle- }}\kern-7.83337pt}}{{\vbox{\hbox{\scriptstyle- }}\kern-5.90005pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.75003pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.25003pt}}\!\int_{\Omega}\Pi_{\mathcal{D}}p_{\mathcal{D}}^{*}(g){\>\rm d}{\bm{x}}=\mathchoice{{\vbox{\hbox{\textstyle- }}\kern-7.83337pt}}{{\vbox{\hbox{\scriptstyle- }}\kern-5.90005pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.75003pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.25003pt}}\!\int_{\Omega}\overline{p}{\>\rm d}{\bm{x}}. This choice is dictated by the pure Neumann boundary condition and will be essential.

For $K\in{\mathcal{T}}$, let $\overline{{\bm{x}}}_{K}$ be the centroid (centre of mass) of $K$. A standard approximation property (see e.g. [19, Lemma A.7] with $w_{K}\equiv 1$) yields

[TABLE]

Define $\widehat{u}$ and $\widehat{p}$ a.e. on $\Omega$ by: For all $K\in{\mathcal{T}}$ and all ${\bm{x}}\in K$, $\widehat{u}({\bm{x}})=\overline{u}(\overline{{\bm{x}}}_{K})$ and $\widehat{p}({\bm{x}})=\overline{p}(\overline{{\bm{x}}}_{K})$. From (3.11) and the Lipschitz continuity of $P_{[a,b]}$, we obtain

[TABLE]

Step 1: estimate of $T_{1}$.

A use of triangle inequality, (1.5b), (5.1b) and (A1)-i) leads to

[TABLE]

We now estimate the last term in this inequality. Use the definitions of $C_{\mathcal{D}}$, $\|\cdot\|_{{\mathcal{D}}}$ and the fact that \mathchoice{{\vbox{\hbox{\textstyle- }}\kern-7.83337pt}}{{\vbox{\hbox{\scriptstyle- }}\kern-5.90005pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.75003pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.25003pt}}\!\int_{\Omega}\Pi_{\mathcal{D}}p_{\mathcal{D}}(\overline{u}){\>\rm d}{\bm{x}}=\mathchoice{{\vbox{\hbox{\textstyle- }}\kern-7.83337pt}}{{\vbox{\hbox{\scriptstyle- }}\kern-5.90005pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.75003pt}}{{\vbox{\hbox{\scriptscriptstyle- }}\kern-4.25003pt}}\!\int_{\Omega}\Pi_{\mathcal{D}}p_{\mathcal{D}}^{*}(\overline{u}){\>\rm d}{\bm{x}} to obtain

[TABLE]

Subtract (5.15) with $g=\overline{u}$ from (5.1b), substitute $z_{\mathcal{D}}=p_{\mathcal{D}}(\overline{u})-p_{\mathcal{D}}^{*}(\overline{u})$, use property (3.9) in (A1)-ii) and Cauchy-Schwarz inequality to obtain

[TABLE]

A use of (5.19) and (A1)-i) leads to $\|\Pi_{\mathcal{D}}p_{\mathcal{D}}(\overline{u})-\Pi_{\mathcal{D}}p_{\mathcal{D}}^{*}(\overline{u})\|\lesssim h^{2}\|\overline{y}\|_{H^{2}(\Omega)}+h^{2}\|f+\overline{u}\|_{H^{1}(\Omega)}$. Plugged into (5.18), this estimate yields

[TABLE]

Step 2: estimate of $T_{2}$.

Subtract the equations (5.15) satisfied by $p_{\mathcal{D}}^{*}(\overline{u})$ and $p_{\mathcal{D}}^{*}(\widehat{u})$ to obtain, for all $z_{\mathcal{D}}\in X_{{\mathcal{D}}}$,

[TABLE]

Since $p_{\mathcal{D}}^{*}(\widehat{u})-p_{\mathcal{D}}^{*}(\overline{u})\in X_{{\mathcal{D}},\star}$, a use of Proposition 2.6 in (5.21) yields

[TABLE]

Choose $z_{\mathcal{D}}=y_{\mathcal{D}}(\overline{u})-y_{\mathcal{D}}(\widehat{u})$ in (5.21), subtract the equations (5.1a) satisfied by $y_{\mathcal{D}}(\overline{u})$ and $y_{\mathcal{D}}(\widehat{u})$, since $p_{\mathcal{D}}^{*}(\overline{u})-p_{\mathcal{D}}^{*}(\widehat{u})\in X_{{\mathcal{D}},\star}$, to obtain

[TABLE]

Set $w_{\mathcal{D}}=p_{\mathcal{D}}^{*}(\overline{u})-p_{\mathcal{D}}^{*}(\widehat{u})$, use orthogonality of $\mathbfcal{P}_{\!{\mathcal{T}}}$, Cauchy-Schwarz inequality and (A2) to infer

[TABLE]

Equation (5.21) and the stability of the GDM (Proposition 2.6) show that

[TABLE]

A use of Holder’s inequality, (A4), (A3), the fact that $w_{\mathcal{D}}\in X_{{\mathcal{D}},\star}$ and (5.24) yields an estimate for the second term on the right hand side of (5.23) as follows:

[TABLE]

Consider now the last term on the right hand side of (5.23). For any $K\in{\mathcal{T}}_{2}$, we have $\overline{u}=a$ on $K$, $\overline{u}=b$ on $K$, or, by (3.19), $\overline{u}=-\alpha^{-1}\overline{p}+\overline{c}$ on $K$. Hence, on $K$, $\mathbfcal{P}_{\!{\mathcal{T}}}\overline{u}-\widehat{u}=0$ or $\mathbfcal{P}_{\!{\mathcal{T}}}\overline{u}-\widehat{u}=\alpha^{-1}\left(\widehat{p}-\mathbfcal{P}_{\!{\mathcal{T}}}\overline{p}\right)$. This leads to $|\mathbfcal{P}_{\!{\mathcal{T}}}\overline{u}-\widehat{u}|\leq\alpha^{-1}|\widehat{p}-\mathbfcal{P}_{\!{\mathcal{T}}}\overline{p}|$ on $\Omega_{2,{\mathcal{T}}}$. Use (5.16), the definition of $C_{\mathcal{D}}$, the fact that $w_{\mathcal{D}}\in X_{{\mathcal{D}},\star}$ and (5.24) to obtain

[TABLE]

Plug (5.24), (5.25) and (5.26) into (5.23) and then in (5.22) to get

[TABLE]

Step 3: estimate of $T_{3}$.

Subtract (3.1b) from (5.15) with $g=\widehat{u}$ and (3.1a) from (5.1a) with $\widehat{u}$ instead of $\overline{u}$, we obtain for all $z_{\mathcal{D}}\in X_{{\mathcal{D}}}$ and $w_{\mathcal{D}}\in X_{{\mathcal{D}},\star}$,

[TABLE]

[TABLE]

Substitute $z_{\mathcal{D}}=p_{\mathcal{D}}^{*}(\widehat{u})-\overline{p}_{\mathcal{D}}\in X_{{\mathcal{D}},\star}$ in (5.28), $w_{\mathcal{D}}=y_{\mathcal{D}}(\widehat{u})-\overline{y}_{\mathcal{D}}\in X_{{\mathcal{D}},\star}$ in (5.29) and use Proposition 2.6 to obtain

[TABLE]

The optimality condition (1.5c) [30, Lemma 3.5] yields for a.e. ${\bm{x}}\in\Omega$,

[TABLE]

Since $\overline{u}$, $\overline{p}$ and $\overline{u}_{h}$ are continuous at the centroid $\overline{{\bm{x}}}_{K}$, we can choose ${\bm{x}}=\overline{{\bm{x}}}_{K}$ and $v(\overline{{\bm{x}}}_{K})=\overline{u}_{h}(\overline{{\bm{x}}}_{K})(=\overline{u}_{h}$ on $K$). All the involved functions being constants over $K$, this gives

[TABLE]

Integrate over $K$ and sum over $K\in{\mathcal{T}}$ to deduce

[TABLE]

Choose $v_{h}=\widehat{u}$ in the discrete optimality condition (3.1c) to establish

[TABLE]

Add the above two inequalities to obtain

[TABLE]

and thus

[TABLE]

Since $\overline{u}_{h}-\widehat{u}$ is piecewise constant on ${\mathcal{T}}$, the orthogonality property of $\mathbfcal{P}_{\!{\mathcal{T}}}$, (5.16) and (3.10) in (A1)-ii) lead to

[TABLE]

By Cauchy–Schwarz inequality, triangle inequality and the notations in (5.17), we obtain

[TABLE]

Subtract the equations (3.1a) and (5.1a) (with $\widehat{u}$ instead of $\overline{u}$) satisfied by $\overline{y}_{\mathcal{D}}$ and $y_{\mathcal{D}}(\widehat{u})$, choose $w_{\mathcal{D}}=p_{\mathcal{D}}^{*}(\widehat{u})-\overline{p}_{\mathcal{D}}$, and use the equations (3.1b) and (5.15) on $\overline{p}_{\mathcal{D}}$ and $p_{\mathcal{D}}^{*}(\widehat{u})$ to arrive at

[TABLE]

A substitution of (5.32)–(5.34) (together with the estimates (5.20) and (5.27) of $T_{1}$ and $T_{2}$) into (5.31) yields an estimate on $\|\overline{u}_{h}-\widehat{u}\|$ which, when plugged into (5.30), gives

[TABLE]

Step 4: conclusion.

A use of (1.2) and the fact that $\overline{u}$ is optimal leads to

[TABLE]

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

[TABLE]

From (3.19) and (3.3), we have

[TABLE]

Note that $|\nabla(-\alpha^{-1}\overline{p}+\overline{c})|=\alpha^{-1}|\nabla\overline{p}|$. Therefore,

[TABLE]

Combine (5.36) and (5.37) to obtain

[TABLE]

Use (5.38) in (5.20), (5.27) and (5.35) and plug the resulting estimates in (5.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 $T_{2}$. This estimate is the only source of the $2-\frac{1}{2^{*}}$ power (instead of 2), and the only place where we used Assumption (A3), here replaced by (3.14). Recall (A4) and use (3.14) in (5.21) satisfied by $p_{\mathcal{D}}^{*}(\overline{u})-p_{\mathcal{D}}^{*}(\widehat{u})$ to write

[TABLE]

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]. ∎

Acknowledgement: 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.