A finite element method for Dirichlet boundary control of elliptic partial differential equations

TL;DR

This paper develops a new variational finite element approach for Dirichlet boundary control of elliptic PDEs, proving well-posedness and convergence, and validating with numerical examples.

Contribution

It introduces a novel variational formulation that relates the state and adjoint state through boundary control, with proven well-posedness and convergence analysis for finite element approximations.

Findings

Finite element approximations converge at order k-1/2.

The variational formulation is well-posed in H^1 spaces.

Numerical results confirm theoretical convergence rates.

Abstract

This paper introduces a new variational formulation for Dirichlet boundary control problem of elliptic partial differential equations, based on observations that the state and adjoint state are related through the control on the boundary of the domain, and that such a relation may be imposed in the variational formulation of the adjoint state. Well-posedness (unique solvability and stability) of the variational problem is established in the $H^{1}(\Omega)\times H_{0}^{1}(\Omega)$ space for the respective state and adjoint state. A finite element method based on this formulation is analyzed. It is shown that the conforming $k-$th order finite element approximations to the state and the adjoint state, in the respective $L^{2}$ and $H^{1}$ norms converge at the rate of order $k-1/2$ on quasi-uniform mesh for conforming element of order $k$. Numerical examples are presented to validate the theory.

Tables

Table 1. Table 1: Numerical data of γ = 1 𝛾 1 \gamma=1 for Example 1: h ℎ h – maximum size of quasi-uniform mesh; ‖ ∇ ( y − y h ) ‖ norm ∇ 𝑦 subscript 𝑦 ℎ ||\nabla(y-y_{h})|| – numerical error for the state variable y 𝑦 y ; order y – the speed of convergence for y 𝑦 y ; ‖ ∇ ( z − z h ) ‖ norm ∇ 𝑧 subscript 𝑧 ℎ ||\nabla(z-z_{h})|| – numerical error for the adjoint state variable z 𝑧 z ; order z – the speed of convergence for z 𝑧 z ; ‖ u − u h ‖ 0 , Γ subscript norm 𝑢 subscript 𝑢 ℎ 0 Γ ||u-u_{h}||_{0,\Gamma} – numerical error for the control variable u 𝑢 u ; order u – the speed of convergence for u 𝑢 u .

$h$	$\Vert \nabla \left(y-y_h\right)\Vert$	${\mathrm{order}}_y$	$\Vert \nabla \left(z-z_h\right)\Vert$	${\mathrm{order}}_z$	${\Vert u-u_h\Vert }_{0,\Gamma }$	${\mathrm{order}}_u$
0.7071	0.7187	–	0.1069	–	0.1901	–
0.3536	0.3603	0.9964	0.0539	0.9881	0.0663	1.5200
0.1768	0.1928	0.9021	0.0278	0.9552	0.0345	0.9424
0.0884	0.0898	1.1023	0.0140	0.9897	0.0154	1.1637
0.0442	0.0446	1.0097	0.0070	1.000	0.0066	1.2224

Table 2. Table 2: Numerical data of γ = 1 𝛾 1 \gamma=1 for Example 1: h ℎ h – maximum size of quasi-uniform mesh; ‖ y − y h ‖ norm 𝑦 subscript 𝑦 ℎ ||y-y_{h}|| – numerical error for the state variable y 𝑦 y ; order y subscript order 𝑦 {\rm order}_{y} – the speed of convergence for y 𝑦 y in L 2 superscript 𝐿 2 L^{2} norm; ‖ z − z h ‖ norm 𝑧 subscript 𝑧 ℎ ||z-z_{h}|| – numerical error for the adjoint state variable z 𝑧 z ; order z subscript order 𝑧 {\rm order}_{z} – the speed of convergence for z 𝑧 z in L 2 superscript 𝐿 2 L^{2} norm.

$h$	0.7071	0.3536	0.1768	0.0884	0.0442	0.0221
$\Vert y-y_h\Vert$	0.0897	0.0250	0.0078	0.0025	7.86e-004	2.55e-004
${\mathrm{order}}_y$	–	1.8436	1.6804	1.6415	1.6686	1.6259
$\Vert z-z_h\Vert$	0.0181	0.0054	0.0014	3.55e-004	8.74e-005	2.10e-005
${\mathrm{order}}_z$	–	1.7453	1.9475	1.9775	2.0234	2.0604

Table 3. Table 3: Numerical data of γ = 0.01 𝛾 0.01 \gamma=0.01 for Example 1: h ℎ h – maximum size of quasi-uniform mesh; ‖ y − y h ‖ norm 𝑦 subscript 𝑦 ℎ ||y-y_{h}|| – numerical error for the state variable y 𝑦 y ; order y – the speed of convergence for y 𝑦 y in L 2 superscript 𝐿 2 L^{2} norm; ‖ ∇ ( z − z h ) ‖ norm ∇ 𝑧 subscript 𝑧 ℎ ||\nabla(z-z_{h})|| – numerical error for the adjoint state variable z 𝑧 z ; order z – the speed of convergence for z 𝑧 z ; ‖ u − u h ‖ 0 , Γ subscript norm 𝑢 subscript 𝑢 ℎ 0 Γ ||u-u_{h}||_{0,\Gamma} – numerical error for the control variable u 𝑢 u ; order u – the speed of convergence for u 𝑢 u .

$h$	$\Vert y-y_h\Vert$	${\mathrm{order}}_y$	$\Vert \nabla \left(z-z_h\right)\Vert$	${\mathrm{order}}_z$	${\Vert u-u_h\Vert }_{0,\Gamma }$	${\mathrm{order}}_u$
0.7071	2.9637	–	0.1134	–	9.0693	–
0.3536	0.7594	1.9649	0.0561	1.0156	2.5525	1.8295
0.1768	0.2101	1.8538	0.0279	1.0077	0.8519	1.5832
0.0884	0.0663	1.6640	0.0140	0.9948	0.3384	1.3320
0.0442	0.0191	1.7954	0.0070	1.000	0.1187	1.5114

Table 4. Table 4: Numerical data of γ = 1 𝛾 1 \gamma=1 for Example 2: h ℎ h – maximum size of quasi-uniform mesh; ‖ y − y h ‖ norm 𝑦 subscript 𝑦 ℎ ||y-y_{h}|| – numerical error for the state variable y 𝑦 y in L 2 superscript 𝐿 2 L^{2} norm; order y – the speed of convergence for y 𝑦 y ; ‖ z − z h ‖ norm 𝑧 subscript 𝑧 ℎ ||z-z_{h}|| – numerical error for the adjoint state variable z 𝑧 z in L 2 superscript 𝐿 2 L^{2} norm; order z – the speed of convergence for z 𝑧 z ; ‖ u − u h ‖ 0 , Γ subscript norm 𝑢 subscript 𝑢 ℎ 0 Γ ||u-u_{h}||_{0,\Gamma} – numerical error for the control variable u 𝑢 u ; ‖ ∇ ( y − y h ) ‖ norm ∇ 𝑦 subscript 𝑦 ℎ ||\nabla(y-y_{h})|| – numerical error for the state variable y 𝑦 y in H 1 superscript 𝐻 1 H^{1} seminorm.

$h$	$\Vert y-y_h\Vert$	${\mathrm{order}}_y$	$\Vert z-z_h\Vert$	${\mathrm{order}}_z$	${\Vert u-u_h\Vert }_{0,\Gamma }$	$\left|\right|\nabla (y-y_h|\left|\right)$
0.1768	0.0117	–	7.28e-005	–	0.3212	0.4462
0.0884	0.0034	1.7833	2.51e-005	1.5398	0.2999	0.3283
0.0442	0.0011	1.6280	7.95e-006	1.6554	0.2896	0.3167
0.0221	3.61e-004	1.6078	2.10e-006	1.9218	0.2794	0.3048
0.0111	1.18e-004	1.6185	5.29e-007	1.9875	0.2693	0.2977

Equations

$$\min J(u),J(u)=\frac{1}{2}||y-y_d|{|}_{L^2(\Omega )}^2+\frac{\gamma }{2}||u|{|}_{L^2(\Gamma )}^2,$$

$$-△y=f\mathrm{in}\Omega ,$$

$$y=u\mathrm{on}\Gamma .$$

$$-△y=f\mathrm{in}\Omega ,$$

$$y=u\mathrm{on}\Gamma ,$$

$$-△z=y-y_d\mathrm{in}\Omega ,$$

$$z=0\mathrm{on}\Gamma ,$$

$$u=\frac{1}{\gamma }\frac{\partial z}{\partial \mathbf{n}}\mathrm{on}\Gamma ,$$

$${\int }_{\Omega }\nabla y\cdot \nabla \psi d\mathbf{x}={\int }_{\Omega }f\psi d\mathbf{x}.$$

$${\int }_{\Omega }\nabla z\cdot \nabla \varphi d\mathbf{x}-{\int }_{\Gamma }\frac{\partial z}{\partial \mathbf{n}}\varphi ds={\int }_{\Omega }(y-y_d)\varphi d\mathbf{x}.$$

$$\displaystyle\frac{\partial z}{\partial{\bf n}}\big{|}_{\Gamma}=\gamma y|_{\Gamma}.$$

$${\int }_{\Omega }\nabla z\cdot \nabla \varphi d\mathbf{x}-{\int }_{\Gamma }\gamma y\varphi ds={\int }_{\Omega }(y-y_d)\varphi d\mathbf{x}.$$

$$(\nabla y,\nabla \psi )=(f,\psi )\forall \psi \in H_0^1(\Omega ),$$

$$(\nabla z,\nabla \varphi )-(\gamma y,\varphi {)}_{\Gamma }-(y,\varphi )=-(y_d,\varphi )\forall \varphi \in H^1(\Omega ).$$

$$p_{\ast }^{\Omega }=2{\omega }_{\max }/(2{\omega }_{\max }-\pi ),$$

$$(y,u,z)\in H^{3/2-1/p}(\Omega )\times H^{1-1/p}(\Gamma )\times (H_0^1(\Omega )\cap W_p^2(\Omega )),2\le p<p_{\ast }.$$

$${\gamma }^{1/2}||u|{|}_{0,\Gamma }+||y||+||\nabla z||\le C_{\gamma }\left(||f|{|}_{-1}+||y_d||\right).$$

$$\begin{array}[]{lll}||\nabla z||^{2}&=&(\gamma y,z)_{\Gamma}+(y,z)-(y_{d},z)\vspace{2mm}\\ &\leq&\gamma||y||_{-1/2,\Gamma}||z||_{1/2,\Gamma}+||y||_{-1}||z||_{1}+||y_{d}||_{-1}||z||_{1}\vspace{2mm}\\ &\leq&C(\gamma||y||_{0,\Gamma}||z||_{1}+\|y\|\,\|z\|_{1}+\|y_{d}\|\,\|z\|_{1}),\end{array}$$

$$||z|{|}_1\le C\left(\gamma ||y|{|}_{0,\Gamma }+||y||+||y_d||\right).$$

$$\begin{array}[]{lll}\gamma||y||_{0,\Gamma}^{2}+||y||^{2}&=&(\nabla z,\nabla y)+(y_{d},y)\vspace{2mm}\\ &=&(f,z)+(y_{d},y)\leq||f||_{-1}||z||_{1}+\|y_{d}\|\,\|y\|\vspace{2mm}\\ &\leq&C||f||_{-1}\left(\gamma||y||_{0,\Gamma}+||y||+||y_{d}||\right)+\|y_{d}\|\,\|y\|\vspace{2mm}\\ &\leq&C\left(||f||_{-1}+||y_{d}||\right)\left(\gamma||y||_{0,\Gamma}+||y||\right)+C\left(||f||_{-1}^{2}+||y_{d}||^{2}\right),\end{array}$$

$$\gamma ||y|{|}_{0,\Gamma }^2+||y|{|}^2\le C_{\gamma }\left(||f|{|}_{-1}^2+||y_d|{|}^2\right).$$

$$||\nabla y||\le C_{\gamma }\left(||f|{|}_{-1}+||y_d||\right).$$

$$\begin{array}[]{lll}||\nabla y||&\leq&C\left(||f|_{-1}+||u||_{1/2,\Gamma}\right)\vspace{2mm}\\ &\leq&\displaystyle C\left(||f|_{-1}+\frac{1}{\gamma}\Big{\|}\frac{\partial z}{\partial{\bf n}}\Big{\|}_{1/2,\Gamma}\right).\end{array}$$

$$\begin{array}[]{lll}\displaystyle||z||_{2}+\Big{\|}\frac{\partial z}{\partial{\bf n}}\Big{\|}_{1/2,\Gamma}&\leq&C||y-y_{d}||\vspace{2mm}\\ &\leq&C\left(||y||+||y_{d}||\right).\end{array}$$

$${\gamma }^{1/2}||u|{|}_{0,\Gamma }+||y|{|}_1+||z|{|}_1\le C_{\gamma }\left(||f|{|}_{-1}+||y_d|\right).$$

$$V_h:=\{v_h\in C(\overline{\Omega }):v_h{|}_K\in P_k(K),\forall K\in {\mathcal{T}}_h\}$$

$$(\nabla y_h,\nabla {\psi }_h)=(f,{\psi }_h)\forall {\psi }_h\in V_h^0,$$

$$(\nabla z_h,\nabla {\varphi }_h)-(\gamma y_h,{\varphi }_h{)}_{\Gamma }-(y_h,{\varphi }_h)=-(y_d,{\varphi }_h)\forall {\varphi }_h\in V_h.$$

$$(\nabla y_h,\nabla {\psi }_h)=0\forall {\psi }_h\in V_h^0,$$

$$(\nabla z_h,\nabla {\varphi }_h)-(\gamma y_h,{\varphi }_h{)}_{\Gamma }-(y_h,{\varphi }_h)=0\forall {\varphi }_h\in V_h.$$

$$(\nabla z_h,\nabla y_h)-(\gamma y_h,y_h{)}_{\Gamma }-||y_h|{|}^2=0.$$

$${\int }_{\Gamma }\gamma y_h^{2}ds+||y_h|{|}^2=0,$$

$$(\nabla z_h,\nabla {\varphi }_h)=0\forall {\varphi }_h\in V_h,$$

$$||\nabla {\theta }_1^b||\lesssim h^{-1/2}||{\theta }_1^b|{|}_{L^2(\Gamma )}.$$

Full text

1 Introduction

Let $\Omega\subset\mathbb{R}^{d},d\geq 2$, be a bounded polygonal or polyhedral domain with Lipshitz boundary $\Gamma=\partial\Omega$. Consider the following Dirichlet boundary control problem of elliptic partial differential equations (PDEs):

[TABLE]

where the regularization parameter $\gamma>0$ and $y$ is the solution of the Poisson equation with nonhomogeneous Dirichlet boundary conditions

[TABLE]

After the pioneering works of Falk Falk and Geveci Geveci , there were some efforts on the error estimates for finite element approximation to control problems governed by PDEs. Arada et al. in Arada ; Casas2002 derived error estimates for the control in the $L^{\infty}$ and $L^{2}$ norms for semilinear elliptic control problem. The articles Gunzburger1996 ; Fursikov studied the error estimates of finite element approximation for some important flow control problems. Casas et al. Casas2005 carried out the study of the Neumann boundary control problem.

However, the works mentioned above are mainly contributions to the distributed control. Since the Dirichlet boundary control plays an important role in many applications such as flow control problems and has been a hot topic for decades. It is well known that Dirichlet control problems are difficult theoretically and numerically, because the Dirichlet boundary data does not directly enter a standard variational setting for the PDEs. On the one hand, the traditional finite element method (see, e.g., Casa2006 ; Vexler2007 ; Deckelnick2009 ; May2013 ; Apel2017 ) deals with the state variable ($y$) using its weak formulation, e.g., allowing for solutions $y\in L^{2}(\Omega)$; on the other hand, the attempt of the first order optimality condition involves the normal derivative of the adjoint state ($z$) on the boundary of the domain. Therefore, it is crucial to obtain this normal derivative numerically by using additional equation. But in doing so the problem becomes complicated in both theoretical analysis and numerical practice. Note that the regularity of the solution and error estimates for finite element approximates have been studied in Apel2015 ; Mateos2016 ; Apel2017 .

To avoid the difficulty described above, there are two ways to deal with the control variable. One is in Of2015 replaced the $L^{2}$ norm in the cost functional with the $H^{1/2}$ norm and attained a priori estimate of the numerical error of the control by using piecewise linear elements, and the other approximate the nonhomogeneous Dirichlet boundary condition with a Robin boundary condition or weak boundary penalization. However, the former changed the problem and the latter had to deal with the penalization which is computationally expensive. Recently, techniques similar to Of2015 were used in Gunzburger1991 ; Gunzburger1992 ; Chowdhury2017 .

Recently, Gong et al. considered the mixed finite element method in Gong2011 , where the optimal control and the adjoint state were involved in a variational form in a natural sense. This makes its theoretical analysis easier, but the corresponding fluxes of the two states $y$ and $z$ are required to be introduced. It points out that the mixed finite method obtained the same rate of convergence as the regularity of the control on boundary. Apel et al. Apel2017 have considered a standard finite element method on a special class of meshes and guaranteed a superlinear convergence rate for the control. Very recently, Hu et al. Hu2017 considered a hybridizable discontinuous Galerkin method to obtain optimal a priori error estimates for the control by solving an algebraic system of seven unknown functions.

Based on both the facts that the control $u$ is equal to the restriction of the state $y$ on the boundary (see the original equation (3)), and that the restriction of an approximation of the state $y$ on the boundary is also an approximation to the control $u$, we realize that the restriction of the numerical error for the state $y$ on the boundary can be used to measure the numerical error of the control $u$ in the $L^{2}(\Gamma)$-norm. This idea is done in the way that the state $y$ and its adjoint state $z$ will be coupled by the original equation (3) and an extra equation (8) as well as by the right-hand side term $y$ of the equation (6), i.e., the control $u$ and the normal derivative of the adjoint state $z$ along the boundary are cancelled. This is different from the idea in literatures e.g., Casa2006 ; Vexler2007 ; Deckelnick2009 ; May2013 ; Apel2017 ; Of2015 ; Gunzburger1991 ; Gunzburger1992 ; Chowdhury2017 , where both the original equation (3) and an extra equation (8) were taken into account in variational formulation.

This paper introduces a new variational formulation for Dirichlet boundary control problem of elliptic partial differential equations, based on observations that the state and adjoint state are related through the control on the boundary of the domain, and that such a relation may be imposed in the variational formulation of the adjoint state, i.e., one can substitute the control by the control law (the control is the normal derivative of the adjoint state (up to a factor)). Well-posedness (unique solvability and stability) of the variational problem is established in the $H^{1}(\Omega)\times H_{0}^{1}(\Omega)$ space for the respective state and adjoint state. A finite element method based on this formulation is analyzed. It is shown that the conforming $k-$th order finite element approximations to the state and the adjoint state, in the respective $L^{2}$ and $H^{1}$ norms converge at the rate of order $k-1/2$ on quasi-uniform mesh for conforming element of order $k$. Numerical examples are presented to validate the theory.

This paper is organized as follows. In Section 2, we introduce a new variational setting based on an observation. Section 3 is devoted to the unique solvability and stability of the variational problem. In Section 4, we introduce finite element approximation to the variational setting and prove a preliminary result, which will prepare us for the a priori error estimation on an approximation of the conforming element of order $k$ over quasi-uniform mesh in Section 5. In Section 6, we analyze the stability of the discrete control in $L^{2}(\Gamma)$ norm and $H^{1/2}(\Gamma)$ norm in the sense that the restriction of the discrete state on the boundary is considered as an approximation of the control. Finally numerical tests are provided in Section 7 to support our theory.

2 A variational formulation

For any bounded open subset $\omega$ of $\Omega$ with Lipschitz boundary $\gamma$, let $L^{2}(\gamma)$ and $H^{m}(\omega)$ be the standard Lebesgue and Sobolev spaces equipped with standard norms $\|\cdot\|_{\gamma}=\|\cdot\|_{L^{2}(\gamma)}$ and $\|\cdot\|_{m,\omega}=\|\cdot\|_{H^{m}(\omega)}$, $m\in\mathbb{N}$. Note that $H^{0}(\omega)=L^{2}(\omega)$. Denote by $|\cdot|_{m,\omega}$ the semi-norm in $H^{m}(\omega)$. Similarly, denote by $(\cdot,\cdot)_{\gamma}$ and $(\cdot,\cdot)_{\omega}$ the $L^{2}$ inner products on $\gamma$ and $\omega$, respectively. We shall omit the symbol $\Omega$ in the notations above if $\omega=\Omega$.

It is well known that the Dirichlet boundary control problem in (1)-(3) is equivalent to the optimality system

[TABLE]

where ${\bf n}$ is the unit outer normal to $\Gamma$. Note that these equations must be understood in a weak sense.

To see the idea of variational setting, we consider the following several cases under an assumption that the domain and known data are respectively satisfied with these cases:

Case one: the state $y\in H^{1/2}(\Omega)$, so $y|_{\Gamma}$ belongs to $L^{2}(\Gamma)$. We know $u\in L^{2}(\Gamma)$ from (5). The equation (8) implies that $\partial z/\partial{\bf n}|_{\Gamma}\in L^{2}(\Gamma)$, this needs the adjoint state to satisfy $z\in H^{3/2}(\Omega)$.

Case two: $y\in H^{1}(\Omega),y|_{\Gamma}=u\in H^{1/2}(\Gamma)$, the equation (8) means that $\partial z/\partial{\bf n}|_{\Gamma}\in H^{1/2}(\Gamma)$, this requires $z\in H^{2}(\Omega)$.

Case three: $y\in L^{2}(\Omega),y|_{\Gamma}=u\in H^{-1/2}(\Gamma)$ (the dual space of $H^{1/2}(\Gamma)$), the equation (8) requires $\partial z/\partial{\bf n}|_{\Gamma}\in H^{-1/2}(\Gamma)$, this indicates $z\in H^{1}(\Omega)$.

Case four: $y\in H^{3/2}(\Omega),y|_{\Gamma}=u\in H^{1}(\Gamma)$, the equation (8) needs $\partial z/\partial{\bf n}|_{\Gamma}\in H^{1}(\Gamma)$, this requires $z\in H^{5/2}(\Omega)$.

Since natural functional analytical setting of this problem uses $L^{2}(\Gamma)$ as a “control space”, Case one is an ideal choice for the control $u$, state $y$, and adjoint state $z$. However, it is difficult to bring this characteristics of $y$ and $z$ into their respective variational formulation if (5) and (8) are regarded as two independent equations. For Cases two and three, it is convenient to incorporate the spaces of $y$ and $z$ into their respective variational formulation, but doing so expands or narrows down the space of the control $u$, and can not provide the variational formulation of $u$ if (3) and (8) are still regarded as two independent equations. Case four further does not only enlarges the space of $u$, but also requires a higher regularity on $y$ and $z$, and bring an unexpected difficulty to variation and computation.

These cases show that it is difficult to keep the compatibility of the spaces of $u,y$ and $z$ and incorporate their respective space into their respective variational formulation. We realize that the state $y$ and adjoint state $z$ are connected by the control $u$ on the boundary (see (5) and (8)) as well as by the state $y$ being the right-hand side term of the equation of the adjoint state (see (6)), and that it can overcome these difficulties referred above to cancel the control $u$ and to absorb \displaystyle\frac{1}{\gamma}\frac{\partial z}{\partial{\bf n}}\big{|}_{\Gamma}=y|_{\Gamma} into the variational formulation of $z$ as a boundary condition.

Based on the above observation, the Dirichlet boundary condition in (5) indicates that the control $u$ is equal to the restriction of the state $y$ on the boundary $\Gamma$. Therefore, we simultaneously obtain the control $u$ if the state $y$ is got. The equation (8) is an additional equation with respect to the adjoint state $z$. Here, we do not regard (8) as an additional equation, instead we understand it as a boundary condition, through which the state $y$ and its adjoint state $z$ will be coupled. So the control $u$ can be cancelled in form, but it can be reflected by the state $y$ in essence. It is pointed out that the right hand term of (6) includes the state variable $y$, through which the adjoint state $z$ is coupled over the whole domain.

Based on this idea, multiplying both sides of (4) by $\psi\in H_{0}^{1}(\Omega)$, and applying integration by parts, we attain

[TABLE]

Similarly, multiplying both sides of (6) by $\phi\in H^{1}(\Omega)$, an integration by parts leads to

[TABLE]

Cancelling $u$ from a combination of (5) and (8) yields to

[TABLE]

Substituting (11) into (10), we get

[TABLE]

Collecting (9) and (12) gives the following variational formulation: Find $(y,z)\in H^{1}(\Omega)\times H_{0}^{1}(\Omega)$ such that

[TABLE]

In what follows, we clarify the unique solvability of the variational problem in (13)-(14). For a 2D convex polygonal domain, we recall a regularity result of May et al. in May2013 below, which gives conditions on the domain and data to guarantee the regularity of the solution. To this end, let $\omega_{\rm max}$ be the maximum interior angle of the polygonal domain $\Omega$, and denote $p_{*}^{\Omega}$ by

[TABLE]

including the special case $p_{*}^{\Omega}=\infty$ for $\omega_{\rm max}=\pi/2$. For a higher dimensional convex polygonal domain, we do not attempt to provide condition on the regularity of the solution, because we put an emphasis on a variational setting and the corresponding finite element approximation. Of course, the regularity theory is more complicated in three-dimensional case.

Lemma 1

(May2013 Lemma 2.9).* Suppose that $f\in L^{2}(\Omega)$ and $y_{d}\in L^{p_{*}^{d}}(\Omega),p_{*}^{d}>2$, and that $\Omega\subset\mathbb{R}^{2}$ is a bounded convex domain with polygonal boundary $\Gamma$. Let $p_{*}^{\Omega}\geq 2$ be defined by (15) and $p_{*}:=\min(p_{*}^{d},p_{*}^{\Omega})$. Then, the solution $(y,u)$ of the optimization problem (1)-(3) and the associated adjoint state determined by (6) have the regularity properties*

[TABLE]

Owing to the optimal system (4)-(8) equivalent to the problem (1)-(3), the regularity of the solution for the system (13)-(14) is guaranteed in terms of Lemma 1 in case of $p=2$.

3 Unique solvability and stability

This section establishes unique solvability for the variational problem in (13)-(14), and stability estimate of the control and the state and the adjoint state variables.

Theorem 3.1

For $f\in H^{-1}(\Omega)$ and $y_{d}\in L^{2}(\Omega)$, the system (13)-(14) is uniquely solvable, and is stable in the sense that there exists a positive constant $C_{\gamma}$, depending on $\gamma$, such that

[TABLE]

Proof

We first prove that the variational problem in (13)-(14) is solvable. Since the existence of the solution for the optimization problem in (1)-(3) has been proven by introducing the so-called “solution operator” and using convex analysis (see Lemma 2.4 in May2013 ), and the first order optimal condition shows that the solution of the optimization problem in (1)-(3) satisfies (4)-(8). Hence, the system (4)-(8) is solvable. Due to the solution of (4)-(8) satisfying the variational problem in (13)-(14), then (13)-(14) has a solution.

In what follows, we prove the stability of the system in (13)-(14). By (14) with $\psi=z$, we obtain

[TABLE]

which, together with the Poincáre inequality, implies

[TABLE]

It follows from (14) with $\psi=y$, (13), the Cauchy-Schwarz and Young inequalities, and (17) that

[TABLE]

which implies

[TABLE]

Now, (16) is a direct consequence of (17) and the fact that $u=y$ on $\Gamma$, and the uniqueness of the solution follows from (16) immediately, since the corresponding homogeneous system has vanishing solution. This completes the proof of the theorem.

Theorem 3.2

Assume that the domain $\Omega$ is convex with Lipshitz boundary. For $f\in H^{-1}(\Omega)$ and $y_{d}\in L^{2}(\Omega)$, there exists a positive constant $C_{\gamma}$ dependent on $\gamma$ such that

[TABLE]

Proof

By the standard $H^{1}(\Omega)$ a priori estimate of the problem in (4)-(5), and equation (8), we have

[TABLE]

The standard $H^{2}(\Omega)$ a priori estimate (see, e.g., May2013 ; Casas2009 ) of the problem in (5)-(6) gives

[TABLE]

Now, (18) is a direct consequence of (19), (20), and (16).

Remark 1

Due to $z\in H_{0}^{1}(\Omega)$, the Pioncaré inequality implies $||z||_{1}\leq C||\nabla z||$. Hence, under the assumption of Theorem 3.2, it holds the following stable estimate:

[TABLE]

4 Finite element approximation and preliminary result

We introduce the discrete formulation of (13)-(14). To this end, let $\mathcal{T}_{h}$ be a partition of $\Omega$ into triangles (tetrahedra for $d=3$) or parallelograms (parallelepiped for $d=3$). With each element $K\in\mathcal{T}_{h}$, we associate two parameters $\rho(K)$ and $\sigma(K)$, where $\rho(K)$ denotes the diameter of the set $K$, and $\sigma(K)$ is the diameter of the largest ball contained in $K$. Let us define the size of the mesh by $h=\max_{K\in\mathcal{T}_{h}}\rho(K)$. About the partition, we also assume that there exists a constant $\rho>0$ such that $h/\rho(K)\leq\rho$ for all $K\in\mathcal{T}_{h}$ and $h>0$, i.e., the mesh $\mathcal{T}_{h}$ is quasi-uniform.

Denote $P_{k}(K)$ be the space of polynomials of total degree at most $k$ if $K$ is a simplex, or the space of polynomials with degree at most $k$ for each variable if $K$ is a parallelogram/parallelepiped. Define the finite element space by

[TABLE]

Furthermore, denote $V_{h}^{0}=V_{h}\cap H_{0}^{1}(\Omega)$.

In the rest of this paper, we denote by $C$ a constant independent of mesh size with different context in different occurrence, and also use the notation $A\lesssim F$ to represent $A\leq CF$ with a generic constant $C>0$ independent of mesh size. In addition, $A\approx F$ abbreviates $A\lesssim F\lesssim A$.

The discrete form reads: Find $(y_{h},z_{h})\in V_{h}\times V_{h}^{0}$ such that

[TABLE]

Theorem 4.1

The discrete variational problem in (21)-(22) exists a unique solution $(y_{h},z_{h})\in V_{h}\times V_{h}^{0}$.

Proof

Since the existence of the solution is equivalent to its uniqueness for a finite-dimensional system, it is sufficient to prove that the corresponding homogeneous system has trivial solution. To this end, let $f=0$ and $y_{d}=0$ in (21) and (22), respectively, we get

[TABLE]

Taking $\phi_{h}=y_{h}$ in (24) leads to

[TABLE]

Noticing $z_{h}\in V_{h}^{0}$ gives $(\nabla z_{h},\nabla y_{h})=0$. Combining this with (25) yields to

[TABLE]

which, altogether with the assumption $\gamma>0$, results in $y_{h}=0$.

(24) with $y_{h}=0$ gives

[TABLE]

which, in turn, yields to $(\nabla z_{h},\nabla z_{h})=0$, by choosing $\phi_{h}=z_{h}$. Noticing $z_{h}\in V_{h}^{0}$, we get $z_{h}=0$. Thus, the corresponding homogeneous system has vanishing solution.

Lemma 2

Assume that $\theta_{1}^{b}\in V_{h}$ and vanishes at all interior nodes of the mesh, and let $h$ be the size of the quasi-uniform mesh $\mathcal{T}_{h}$. It holds the following estimate

[TABLE]

Proof

Denote $\Omega_{h}^{b}$ the set of element with at least one vertex on the boundary. Since $\theta_{1}^{b}$ vanishes at any node of an element whose vertices completely contained in the interior of the domain $\Omega$, it’s restriction on the element is zero. This implies that

[TABLE]

For the sake of simplicity, we consider only triangular element in two dimensions as an example, since the similar proof is easily extended to the other types of element and three-dimensional case. There are three cases as following:

(1). Two vertices of an element $K$ lie on the boundary, i.e., $K$ has an edge $E$ contained in $\Gamma$ ($E=\partial K\cap\Gamma$) (see (Case 1) in Figure 1). Assuming $||\nabla\theta_{1}^{b}||_{K}=0$ indicates that $\theta_{1}^{b}$ is a constant over the element $K$. And since $\theta_{1}^{b}$ vanishes at internal node of $K$ (there exists at least an internal node such as internal vertex), this shows that $\theta_{1}^{b}$ is zero over $K$, and that $||\nabla\theta_{1}^{b}||_{K}$ is a norm of $\theta_{1}^{b}$ over $K$. Further assuming $||\theta_{1}^{b}||_{E}=0$, this leads that $\theta_{1}^{b}$ vanishes over $E$, and that $\theta_{1}^{b}$ vanishes at nodes of $E$, and that $\theta_{1}^{b}$ is zero over $K$. Therefore, $||\theta_{1}^{b}||_{E}$ is another norm of $\theta_{1}^{b}$ over $K$. Since any two norms are equivalent to each other over a finite-dimension space, we attain

[TABLE]

where the positive constant $C_{K}$ depends on the size $h_{K}$ of $K$ (and number of dimensions of $V_{h}|_{K}$). To see the dependence on the size of $K$, we apply scaling argument.

To this end, for any element $K\in\mathcal{T}_{h}$ there exists a bijection $F_{K}:\hat{K}\rightarrow K$, where $\hat{K}$ is the reference element. Denote by $DF_{K}$ the Jacobian matrix and let $J_{K}=|det(DF_{K})|$. It is easy to see that for all element types, the mapping definition and shape-regularity and quasi-uniformity of the grids imply that

[TABLE]

which, results in

[TABLE]

Let $E$ be an edge (side) of $K$, and $\hat{E}$ be an edge (side) of $\hat{K}$ with respect to $E$. Similarly, we have

[TABLE]

We have from (29)

[TABLE]

A combination of (30), (31), and (32) yields to

[TABLE]

(2). Three vertices of an element $K$ lie on the boundary $\Gamma$, i.e., $K$ has two edges contained in $\Gamma$ (see (Case 2) in Figure 1). Suppose that one can always choose an element $K^{\prime}$ that has an internal vertex and a common edge with $K$. Now consider $\theta_{1}^{b}$ over $K\cup K^{\prime}$. Repeating the proof of Case (1), we have

[TABLE]

where $C_{K\cup K^{\prime}}$ relies on the size $h_{K\cup K^{\prime}}$, of $K\cup K^{\prime}$. Using the scaling argument again, we easily obtain

[TABLE]

(34) indicates that $||\hat{\nabla}\hat{\theta}_{1}^{b}||_{\widehat{K\cup K^{\prime}}}\approx||\hat{\theta}_{1}^{b}||_{\widehat{\partial K\cap\Gamma}}$. Hence, we obtain from a combination of (35) and (36)

[TABLE]

(3). Only one vertex ${\bf x}_{j}$, of $K$ lies on the boundary $\Gamma$ (see (Case 3) in Figure 1). Suppose that one can always choose an element $K^{\prime}$ such that $\partial(K\cup K^{\prime})$ contains an boundary edge $E$ and $K^{\prime}$ has a common edge with $K$, i.e., $E\subset\partial(K\cup K^{\prime})\cap\Gamma$. Similarly to Case (2) or (1), we easily obtain

[TABLE]

In fact, in this case, we can also consider $\theta_{1}^{b}$ over the patch $\omega_{{\bf x}_{j}}$ (the set of element shared ${\bf x}_{j}$ with $K$), of ${\bf x}_{j}$. By using the scaling argument, we can obtain

[TABLE]

Collecting (33) and (37)-(39), we obtain from (28)

[TABLE]

which results in the desired estimate (27).

5 Analysis of error

Since the control $u$ is equal to the restriction of the state $y$ on the boundary, i.e., $u=y|_{\Gamma}$, it is natural that the restriction of an approximation $y_{h}$ of $y$ on the boundary is also an approximation of $u$. This shows that $||y-y_{h}||_{0,\Gamma}$ can be used to measure the numerical error of the control, in this sense we write $||u-u_{h}||_{0,\Gamma}=||y-y_{h}||_{0,\Gamma}$.

Theorem 5.1

Assume that $(y,z)\in H^{1}(\Omega)\times H_{0}^{1}(\Omega)$ and $(y_{h},z_{h})\in V_{h}\times V_{h}^{0}$ be the solutions to (13)-(14) and (21)-(22), respectively. For $y\in H^{k+1}(\Omega),z\in H^{k+1}(\Omega)\cap H_{0}^{1}(\Omega)$, and for the numerical error of the state variable $y$, there exists a positive constant $C_{\gamma}$ depending on $\gamma$ such that

[TABLE]

Proof

Denote $R_{h}:H^{1}(\Omega)\rightarrow V_{h}$ the Ritz projection operator by

[TABLE]

Recalling the properties of the Ritz projection Brezz1991 ; Brenner1994 as following

[TABLE]

Setting $\eta_{1}=y-R_{h}y$ and $\theta_{1}=R_{h}y-y_{h}$ gives $y-y_{h}=\eta_{1}+\theta_{1}$. We have from triangle inequality and (42)

[TABLE]

The trace inequality and the properties, (42), of the Ritz projection imply that

[TABLE]

(43) and (44) indicates that we only need to estimate $||\theta_{1}||$ and $||\gamma^{1/2}\theta_{1}||_{0,\Gamma}$ in order to estimate $||y-y_{h}||+||\gamma^{1/2}(y-y_{h})||_{0,\Gamma}$. To this end, let $R_{h}^{0}:H_{0}^{1}(\Omega)\rightarrow V_{h}^{0}$ be the Ritz projection operator by

[TABLE]

Again recalling the properties of the Ritz projection Brezz1991 ; Brenner1994 as following

[TABLE]

Setting $\eta_{2}=z-R_{h}^{0}z,\theta_{2}=R_{h}^{0}z-z_{h}$ gives $z-z_{h}=\eta_{2}+\theta_{2}$. From (14) and (22), we obtain the following orthogonality

[TABLE]

Especially taking $\phi_{h}=\theta_{1}\in V_{h}$ in (46) yields to

[TABLE]

which results in

[TABLE]

From (13) and (21), we get the following orthogonal property

[TABLE]

Taking $\psi_{h}=\theta_{2}\in V_{h}^{0}$ in (48) yields to

[TABLE]

In the second step above, we apply the orthogonal property of the Ritz projection, because of $\theta_{2}\in V_{h}^{0}\subset V_{h}$. Combining (47) with (49), we attain

[TABLE]

In what follows, we estimate each term on the right-hand side of (50). In terms of the proof of (43) and (44), we immediately obtain the estimates of the last two terms on the right-hand side of (50)

[TABLE]

and

[TABLE]

To estimate the first term of on the right-hand side of (50), we decompose $\theta_{1}$ into $\theta_{1}^{i}$ and $\theta_{1}^{b}$, where the value of $\theta_{1}^{i}$ at the internal node equals to the one of $\theta_{1}$ at the corresponding node, and the value of $\theta_{1}^{i}$ at the boundary node is zero; the value of $\theta_{1}^{b}$ at the internal node is zero, and the value of $\theta_{1}^{b}$ at boundary node equals to the one of $\theta_{1}$ at the corresponding node. Obviously, $\theta_{1}=\theta_{1}^{i}+\theta_{1}^{b}$.

Noticing $\theta_{1}^{i}\in V_{h}^{0},\theta_{1}^{b}\in V_{h}$, we have from the definition of the Ritz projection

[TABLE]

We further derive from Lemma 2, together with $\theta_{1}^{b}=\theta_{1}$ on the boundary $\Gamma$

[TABLE]

By combining (50)-(52) with (54), and applying the properties, (45), of the Ritz projection, and Young inequality, we obtain

[TABLE]

which, implies

[TABLE]

Collecting (43), (44), and (55), we get

[TABLE]

which results in the desired estimate (40).

Theorem 5.2

Assume that $(y,z)\in H^{1}(\Omega)\times H_{0}^{1}(\Omega)$ and $(y_{h},z_{h})\in V_{h}\times V_{h}^{0}$ be the solutions to (13)-(14) and (21)-(22), respectively. For $y\in H^{k}(\Omega),z\in H^{k+1}(\Omega)\cap H_{0}^{1}(\Omega)$, and for the numerical error of the state variable $y$, there exists a positive constant $C_{\gamma}$ depending on $\gamma$ such that

[TABLE]

Proof

Following the idea of the proof in Theorem 5.1, using $y\in H^{k}(\Omega)$ instead of $y\in H^{k+1}(\Omega)$ while concerning the Ritz projection of $y$, we obtain the desired estimate (56).

Remark 2

As pointed at the beginning of this section, we understand $||u-u_{h}||_{L^{2}(\Gamma)}$ as $||y-y_{h}||_{0,\Gamma}$. For $y\in H^{k+1}(\Omega),z\in H_{0}^{1}(\Omega)\cap H^{k+1}(\Omega)$, Theorem 5.1 gives the control an estimate

[TABLE]

For $y\in H^{k}(\Omega),z\in H_{0}^{1}(\Omega)\cap H^{k+1}(\Omega)$, Theorem 5.2 gives the control an estimate

[TABLE]

Theorem 5.3

Assume that $(y,z)\in H^{1}(\Omega)\times H_{0}^{1}(\Omega)$ and $(y_{h},z_{h})\in V_{h}\times V_{h}^{0}$ be the solutions to (13)-(14) and (21)-(22), respectively. For $y\in H^{k+1}(\Omega),z\in H^{k+1}(\Omega)\cap H_{0}^{1}(\Omega)$, the numerical errors of the adjoint state $z$ is bounded by

[TABLE]

For $y\in H^{k}(\Omega),z\in H^{k+1}(\Omega)\cap H_{0}^{1}(\Omega)$, the numerical errors of the adjoint state $z$ is bounded by

[TABLE]

Proof

Recalling the decomposition of the error $z-z_{h}$ in the proof of Theorem 5.1, we obtain from the orthogonal property of the Ritz projection

[TABLE]

which, together with the property (45) of the Ritz projection, results in,

[TABLE]

The inequality (60) means that it is sufficient to only estimate $||\nabla\theta_{2}||$ in order to estimate $||\nabla(z-z_{h})||$.

Taking $\phi_{h}=\theta_{2}\in V_{h}^{0}$ in (46) yields to

[TABLE]

which, together with the orthogonal relation $(\nabla\eta_{2},\nabla\theta_{2})=0$, results in,

[TABLE]

Applying the Poincaré inequality, we obtain from (61)

[TABLE]

Combing (60) with (62), we obtain

[TABLE]

which, together with (40) and (56), respectively, results in the desired estimates (57) and (58).

Remark 3

Lemma 1 suggests that the $H^{2}$ regularity for the state cannot be reached on polygonal$/$polyhedral domain. This makes these estimates restricted to the case of $k=1$. However, the $H^{k}$ regularity for the state can be reached for domains with sufficiently smooth boundary. Since the Dirichlet boundary control problem is completely different from the Dirichlet boundary value problem, it is non-trivial to generalize analytical technique for high order element (including isoparametric-equivalent element) for the Dirichlet boundary value problem to the Dirichlet boundary control problem. Here are two remedies in two dimensional case for the sake of simplicity.

In first case, let $\Omega\subset\mathbb{R}^{2}$ be a bounded domain with smooth boundary and $\mathcal{T}_{h}$ be a “triangulation” of $\Omega$, where each triangle at the boundary has at most one curved side. The finite element spaces $V_{h}$ and $V_{h}^{0}$ are defined by

[TABLE]

By using standard interpolation error estimates, we can easily verify that the properties of the Ritz projection on $V_{h}$ (and $V_{h}^{0}$) are still true. Assume that the “triangulation” $\mathcal{T}_{h}$ guarantees Lemma 2. Indeed, this is easily realised by assuming that there exists $\rho>0$ such that for each triangle $T\in\mathcal{T}_{h}$ one can find two concentric circular discs $D_{1}$ and $D_{2}$ such that

[TABLE]

Since $\partial\Omega$ is smooth, for h small enough, we have $h_{e}<2{\rm diam}T<2{\rm diam}D_{2}$ (curved side $e\subset\partial T$, $h_{e}$ denote the arc length of $e$). This indicates Lemma 2 is still valid. Therefore, the results of Theorems 5.1-5.3 are applicable to high order curved-triangle Lagrange element.

In the second case, recall that we have a polyhedral approximation, $\Omega_{h}$ to $\Omega$, and an isoparametric mapping $F^{h}$ such that $F^{h}(\Omega_{h})$ closely approximates to $\Omega$, and denote $\tilde{V}_{h}$ a base finite element space defined on $\Omega_{h}$, the resulting space,

[TABLE]

is an isoparametric-equivalent finite element space (we refer to Ciarlet on details). Let $V_{h}^{0}=V_{h}\cap H_{0}^{1}(F^{h}(\Omega_{h}))$. If we impose the control rule on $\partial(F^{h}(\Omega_{h}))$, i.e., $\displaystyle\frac{1}{\gamma}\frac{\partial z}{\partial{\bf n}}=u$ on $\partial(F^{h}(\Omega_{h}))$, this shows that we are considering the problem (1)-(3) on the domain $F^{h}(\Omega_{h})$. The only difference is that we substitute the domain $\Omega$ in the precious context with $F^{h}(\Omega_{h})$. Since the corresponding Ritz projection $R_{h}:H^{1}(F^{h}(\Omega_{h}))\rightarrow V_{h}$ ($R_{h}^{0}:H_{0}^{1}(F^{h}(\Omega_{h}))\rightarrow V_{h}^{0}$) still possesses the same approximation properties as (42)((45)), and since the result of Lemma 2 can be achieved by the similar proof. Therefore, by repeating the proof of Theorem 5.1, we can obtain the following estimate

[TABLE]

under the assumption that $y\in H^{k}(F^{h}(\Omega_{h})),z\in H^{k+1}(F^{h}(\Omega_{h}))$.

Furthermore, we will assume there is auxiliary mapping $F:\Omega_{h}\rightarrow\Omega$ and that $F_{i}^{h}=I^{h}F_{i}$ for each component of the mapping. Here $I^{h}v$ denotes the isoparametric interpolation by $I^{h}v(F^{h}({\bf x}))=\tilde{I}^{h}\tilde{v}({\bf x})$ for all ${\bf x}\in\Omega_{h}$ where $\tilde{v}({\bf x})=v(F^{h}({\bf x}))$ for all ${\bf x}\in\Omega_{h}$ and $\tilde{I}_{h}$ is the global interpolation for the base finite element space, $\tilde{V}_{h}$ (we refer to Ciarlet on details). Thus, the mapping $\Phi^{h}:\Omega\rightarrow F^{h}(\Omega_{h})$ defined by $\Phi^{h}({\bf x})=F^{h}(F^{-1}({\bf x}))$, suggests that $y$ regarded as a function in $F^{h}(\Omega_{h})$ possesses the same regularity as $J_{\Phi^{h}}^{-1}$ (inverse matrix of the Jacobian $J_{\Phi^{h}}$) when $y$ is smooth enough in the domain $\Omega$, this can easily be observed by the chain rule. Therefore, the key is the regularity of the inverse mapping $F^{-1}$ , because the regularity of $F^{h}$ may be reached by using isoparametric interpolation operator of high order. Unfortunately, in mapping a polyhedral domain to a smooth domain, a $C^{1}$ mapping is inappropriate. However, since $F^{h}(\Omega_{h})$ closely approximates to $\Omega$, and $\partial(F^{h}(\Omega_{h}))$ consists of curved sides, it is certain that $y$ regarded as a function in $F^{h}(\Omega_{h})$ has higher regularity than $y$ regarded as a function in $\Omega_{h}$ . This shows that the regularity of $y$ in the domain $F^{h}(\Omega_{h})$ can be reached asymptotically. Of course, the construction of such a mapping $F$ is non-trivial , but it is done in Lenoir .

It is well known that the $L^{2}$ norm of numerical error is controlled by the $H^{1}$ norm for conforming finite element approximation to the standard Laplacian equation, and that the $L^{2}$ norm of numerical error is of order one higher than the $H^{1}$ norm. The following Theorem 5.4 shows that $||y-y_{h}||$ is still controlled by $||y-y_{h}||_{1}$, but isn’t of order one higher than $||y-y_{h}||_{1}$. This will be testified by numerical experiments in Section7.

Theorem 5.4

Assume that $(y,z)\in H^{1}(\Omega)\times H_{0}^{1}(\Omega)$ and $(y_{h},z_{h})\in V_{h}\times V_{h}^{0}$ be the solutions to (13)-(14) and (21)-(22), respectively. It holds

[TABLE]

Proof

Consider the following Neumann boundary-value problem

[TABLE]

The continuous weak formulation for the problem (65) reads: Find $w\in H^{1}(\Omega)$ such that

[TABLE]

We get the following orthogonality from a combination of (14) and (22)

[TABLE]

Owing to $v_{h}=1\in V_{h}$, the above identity implies that

[TABLE]

This shows that the problem (65) satisfies the consistent condition. Therefore, the weak formulation (66) has a unique solution in the sense that the solutions differ by a constant, and satisfies the following estimate

[TABLE]

Taking $\psi=y$ and $\psi=y_{h}$, respectively, in (66) yields to

[TABLE]

and

[TABLE]

A combination of (68) and (69) leads to

[TABLE]

We obtain from (67)

[TABLE]

A combination (70) and (71) yields to

[TABLE]

which, results in

[TABLE]

we complete the proof of (64).

Remark 4

In terms of the proof of Theorem 5.4, for a function $v\in H^{1}(\Omega)$ satisfying

[TABLE]

it holds an analogue of the Poincaré inequality

[TABLE]

6 Stability for discrete solution

Since the control is firstly concerned in practice for the optimal control problem, this section specially devotes to an analysis of the stability for the control in the sense that the restriction of the discrete state $y_{h}$ on the boundary is an approximation of the control $u$. To this end, let $V_{h}^{\partial}$ be the trace space corresponding to $V_{h}$, i.e., $V_{h}^{\partial}=V_{h}|_{\Gamma}$. Recall the following “inverse estimate” for finite element functions $\chi_{h}\in V_{h}^{\partial}$:

[TABLE]

Indeed, this can be found in May2013 or be proven by combining estimates in Ciarlet ; Brenner1994 with standard results from interpolation theory. We define the $L^{2}$ projection $P_{h}^{\partial}:L^{2}(\Gamma)\rightarrow V_{h}^{\partial}$ by

[TABLE]

By standard results for finite element elements we have the error estimate (see Ciarlet ; Brenner1994 ; Casa2006 )

[TABLE]

Theorem 6.1

Assume that $f\in H^{-1}(\Omega),y_{d}\in L^{2}(\Omega)$, the domain $\Omega$ is convex, and its boundary $\Gamma$ is Lipschitz continuous. There exists a positive constant $C_{\gamma}$ depending on $\gamma$ such that

[TABLE]

Proof

Taking $\phi_{h}=y_{h}$ and $\psi_{h}=z_{h}$ in (22) and (21), respectively, gives

[TABLE]

Noticing $z-z_{h}\in H_{0}^{1}(\Omega)$, we obtain from the Poincaré inequality, (58) with $k=1$, and (20)

[TABLE]

Combining (75) with (76), together with Young inequality , gives

[TABLE]

Applying the stable estimates (16) and (18) of the state $y$ and adjoint state $z$, respectively, we get

[TABLE]

which, results in the desired estimate (74).

Theorem 6.2

Under the assumption of Theorem 6.1, the discrete solutions admit the uniform bound

[TABLE]

Proof

From triangle inequality, “inverse estimate” (72), and the property, (73), of the $L^{2}$ projection operator $P_{h}^{\partial}$, we get

[TABLE]

From (56) with $k=1$ and (20), we have

[TABLE]

A combination (79) and (80) yields to

[TABLE]

The desired estimate (78) follows from a combination of (81), (16) and (18).

7 Numerical experiments

In this section, we test the performance of finite element approximation to the variational formulation developed in this paper with two model problems. The actual solution of the first model problem is known, and the true solution of the second example is unknown, and two settings of the regularization parameter $\gamma$ will be considered here, We are thus able to study the convergence rate of the state $y$ and adjoint state $z$, as well as the control variable $u$ over quasi-uniform mesh, and to study the relation between the singularity of the actual solution and the regularization parameter in Example two. Note that we shall employ piecewise linear element in both examples. Let $\{\psi_{i}\}$ and $\{\phi_{j}\}$ be respectively the basis of $V_{h}^{0}$ and $V_{h}$, then the algebraic system with respect to (21)-(22) has the following form

[TABLE]

7.1 Example one

We consider the problem (1)-(2) over a unit square $\Omega=(0,1)\times(0,1)$ with

[TABLE]

The exact solutions are given by

[TABLE]

It is easy to verify that the control $u$, state $y$, and adjoint state $z$ satisfy

[TABLE]

Here, we consider two settings, $\gamma=1$ and $\gamma=0.01$, of regularization parameter.

We start with an initial mesh consisting of 8 congruent right triangles. Figure 2 reports an approximation solution of the state variable $y$ over the mesh with 8192 elements, which generated by uniform refinement of iterations 5 for regularization parameter $\gamma=1$ (left), and over the mesh with 32768 elements, which generated by uniform refinement of iterations 6 for regularization parameter $\gamma=0.01$ (right). In Figure 3, we depict the pictures of an approximation solution of the adjoint state $z$ over the mesh with 8192 elements for $\gamma=1$ (left) and over the mesh with 32768 elements for $\gamma=0.01$ (right). Figure 4 shows an restriction (which is regarded as an approximation solution of the control variable $u$) of $y_{h}$ on the boundary over the mesh with 32768 elements for $\gamma=1$ (left) and for $\gamma=0.01$ (right).

Table 1 shows respectively the exact errors $\|\nabla(y-y_{h})\|,\|\nabla(z-z_{h})\|$ and $\|u-u_{h}\|_{0,\Gamma}$ for the regularization parameter $\gamma=1$. It is observed that they have the rate of convergence of order one for linear element, which is order half higher than theoretical results. Table 2 reports the true errors of the state and adjoint state in $L^{2}$ norm for $\gamma=1$. It can be seen that $||y-y_{h}||$ has the rate of convergence of order 1.5 at least, and that the speed of convergence of $||z-z_{h}||$ is close to 2. Table 3 provides the exact errors of $||y-y_{h}||,||\nabla(z-z_{h})||$ and $||u-u_{h}||_{0,\Gamma}$ for the regularization parameter $\gamma=0.01$, and the similar rate of convergence to $\gamma=1$ can be observed.

In addition, comparing Table 1 with Table 2, we can see that the speed of convergence of $||y-y_{h}||$ is order half higher than $||\nabla(y-y_{h})||$, and that the rate of convergence of $||z-z_{h}||$ is order one higher than $||\nabla(z-z_{h})||$

7.2 Example two

We consider a 2D example over a square domain $\Omega=(0,1/4)\times(0,1/4)\subset\mathbb{R}^{2}$. The data is chosen as

[TABLE]

where $s=10^{-5}$. Since we do not have an explicit expression for the exact solution, the “reference solution” has been calculated over a fine mesh with 131072 elements. Here, we also consider two settings, $\gamma=1$ and $\gamma=0.01$, of regularization parameter.

We still start with an initial mesh consisting of 8 congruent right triangles. Figures 5 and 6 show an approximation solution to the state $y$ and adjoint state $z$ over the mesh generated by uniform refinement of iteration 6 (with 32768 element) for different regularization parameter $\gamma=1$ (left) and $\gamma=0.01$ (right). Figure 7 reports the restriction of an approximation of the state on the boundary, i.e., an approximation solution of the control $u$, over the mesh generated by uniform refinement of iteration 7 (with 131072 element) for different regularization parameter $\gamma=1$ (left) and $\gamma=0.01$ (right).

From Figures 5 and 7, we observe that the control changes quickly at the four corners of the boundary $\Gamma$. Furthermore, we remark that the control for the regularization parameter $\gamma=0.01$ changes more sharply at the four corners of the boundary than for $\gamma=1$, and that the singularity of the exact solution for $\gamma=0.01$ is stronger than for $\gamma=1$.

From Table 4, we observe that the numerical error of the state $y$ in $L^{2}$ norm has the speed of convergence of order 1.6, and that the rate of convergence of the numerical error for the adjoint state $z$ is still close to order 2. However, the numerical errors $||u-u||_{0,\Gamma}$ and $||\nabla(y-y_{h})||$ have a very slow speed of convergence, this is due to the very low regularity of the exact solutions. In fact, the exact control $u$ has strong singularity at four corners of the boundary. This indicates that adaptive mesh based on a posteriori error estimator is efficient to this type of problems, we refer to the articles cai1 ; Boundary ; Ani ; Babuska1978 ; NUMER ; Verfurth1996 ; Liu2008 ; Li2002 ; Kohls2012 ; Kohls2014 ; Schneider2016 ; Becker2000 ; Hint2008 about adaptive finite element methods on the base of a posteriori error estimates.