1 Introduction↩︎

Let \(\Omega \subset \mathbb{R}^2\) be a bounded, convex open domain with a Lipschitz boundary \(\Gamma\). In the present work, we are interested in the numerical analysis of the Dirichlet boundary control problem \[\label{p1} \underset{u \in {U^{\mathrm{ad}}}}{ minimize } \; \frac{1}{2}\|y-y^d\|^{2}_{0,\Omega} + \frac{\omega}{2} \|u\|^{2}_{0,\Gamma} =: J(y,u)\tag{1}\] subject to \[\label{p2} \begin{align} \nabla \cdot ( - \epsilon \nabla y + \beta y ) + \alpha y =& \, f \qquad \qquad in\Omega, \\ y=& \, u \qquad \qquad on\Gamma, \end{align}\tag{2}\] where the admissible control set \({U^{\mathrm{ad}}}\) is specified by \[\begin{align} \label{p3} {U^{\mathrm{ad}}}:= \{ u \in L^{2}(\Gamma): \; u_{a} \leq u(x) \leq u_{b} \;\; a.e. on \;\; \Gamma\}, \end{align}\tag{3}\] with the real numbers \(u_a\) and \(u_b\) satisfying \(u_a \leq u_b\). For the discretization of this control problem, we consider variational discretization [1], as well as a discontinuous Galerkin approximation of the control variable with piecewise linear finite elements. In both cases, we employ the local discontinuous Galerkin (LDG) method for the discretization of the state. For the notation and conventions, we refer to Sections 2 and 3.

Contribution:

For both control schemes, we present a general estimate for the control and state errors. For the variational discrete scheme, in Theorem 2 for \(u\in H^{s}(\Gamma)\) \((s\in [0,1])\) we obtain the error relation \[\| u - u_h \|_{0,\Gamma} + \|y-y_h\|_{0,\Omega} \lesssim \|y-y_h(u)\|_{0,\Omega} + \|({\mathbf{p}}_h(u) - {\mathbf{p}})\cdot{\mathbf{n}}\|_{0,\Gamma} + \|z_h(u) - z\|_{0,\Gamma} =:E(y,{\mathbf{p}},z),\] which only involves Galerkin approximation errors for the primal potential \(y\), as well as for the adjoint flux \({\mathbf{p}}\) and the adjoint potential \(z\). For the fully discrete case, in Theorem 3 we obtain the error relation \[\begin{gather} \| u - u_h \|_{0,\Gamma} + \|y-y_h\|_{0,\Omega} \lesssim E(y,{\mathbf{p}},z) + \\ +\|\pi_hu-u\|_{0,\Gamma} +\big(\|({\mathbf{p}}_h - {\mathbf{p}}(y_h))\cdot{\mathbf{n}}\|_{0,\Gamma}+\|z_h - z(y_h)\|_{0,\Gamma}\big)^{1/2}\|\pi_hu-u\|_{0,\Gamma}^{1/2} + \|\pi_h u - u\|_{-s,\Gamma}^{1/2}, \end{gather}\] which, in addition to the Galerkin approximation errors for \(y,z\), and \({\mathbf{p}}\), also contains approximation errors for the auxiliary adjoint flux \({\mathbf{p}}(y_h)\) and the potential \(z(y_h)\), as well as projection errors associated with the Carstensen quasi-interpolation operator \(\pi_h\). It is noted that here and throughout the paper, the notation \(A \lesssim B\) means that there exists a positive constant \(C\), independent of the relevant discretization parameters, such that \(A \leq C B\). We would like to emphasize that our numerical analysis is inspired by the techniques introduced in [2] and does not rely on the stability and continuity properties of the discrete solution operator associated with the LDG approximation scheme, thus differing from common practices in the literature; compare, e.g., [3], [4]. In the fully discrete case, for the error estimation we only rely on the uniform boundedness of the discrete optimal states \(y_h\) in the \(L^2(\Omega)\) norm, which can be easily deduced from the structure of the cost functional together with the boundedness of the continuous solution operator. For our LDG approach with linear finite elements for the approximation of the discrete fluxes and potentials, we obtain \[\| u - u_h \|_{0,\Gamma} + \|y-y_h\|_{0,\Omega} \lesssim h^{1/2}\] for both the fully discrete and the variational discrete approaches. This is in accordance with the results reported in [5] for the case where piecewise constant fluxes, piecewise linear potentials, and piecewise linear controls for the approximation of the control problem are utilized and convergence with \(h^{1/2}\) is shown, regardless of the smoothness of the involved states and adjoint variables. One conclusion from the analysis here can be that in the context of DG methods, stability can only be achieved at the expense of accuracy.

Finite element approximation of optimal control problems has been an active research area in engineering design for a long time and has been extensively studied in many scientific and engineering applications. We refer to the monographs [2], [6], [7], as well as the references therein for the theory of optimal control problems and for the development of numerical methods. Compared to distributed control problems, in which the control acts throughout the domain \(\Omega\), there has been relatively less research on boundary control problems. Moreover, the majority of existing studies on boundary control problems focus on Neumann boundary control problems (see, e.g., [8]–[11]), where the control is applied through Neumann boundary conditions rather than through conditions of the form 2 . In contrast to distributed and Neumann control problems, Dirichlet boundary control problems pose additional challenges in both theoretical and numerical analysis, as the Dirichlet boundary data cannot be directly incorporated into a standard variational formulation. As a result, various alternative formulations have been proposed in the literature.

The first approach involves approximating the nonhomogeneous Dirichlet boundary condition using either a Robin boundary condition or a weak boundary penalization method; see, e.g., [12], [13]. The second one, which is adopted in the present study, employs the \(L^2(\Gamma)\) norm. In this case, the traditional finite element method treats the state variable in a very weak sense, as the nonhomogeneous Dirichlet boundary condition is prescribed only in \(L^2(\Gamma)\). Extensive numerical studies have been conducted on elliptic Dirichlet boundary control problems using this formulation. We refer to [4], [14]–[18] for a priori error estimates and to [19] for the regularity analysis. In [15], a control-constrained optimal control problem governed by a semilinear elliptic equation on a convex polygonal domain was investigated, where error estimates of order \(h^s\) for both the control and state variables were derived with \(s < \min(1, \pi/(2\theta))\) and \(\theta\) denoting the largest interior angle of the domain. Subsequently, an enhanced error estimate of order \(h^{3/2}\) for the control variable in the case of smooth domains was presented in [16] by exploiting the superconvergence properties of regular triangulations. Then, improved convergence rates for the state variable in the unconstrained case on convex domains were obtained in [17] through a duality argument combined with a control estimate in a norm weaker than \(L^2(\Gamma)\). Later, error estimates of order \(h^s\), with \(s < \min(1, \pi/\theta - 1/2)\) for general meshes and \(s < \min(3/2, \pi/\theta - 1/2)\) for superconvergent meshes, were established for the control variable on general polygonal domains in [14]. Another common approach to formulate Dirichlet boundary control problems is the energy space method, where the \(L^2\)-norm in the cost functional is replaced by the \(H^{1/2}\)-norm. This choice of control space allows for a standard weak formulation of the state equation and yields optimal estimates on a convex polygonal domain; however, it introduces an additional operator, namely the Steklov–Poincaré operator, through the harmonic extension of the given Dirichlet boundary data; see, e.g., [20], [21]. In a similar fashion, the \(H^1\)-norm, which naturally results in harmonic control, is utilized in [22], [23] to regularize the control variable. Further, it is noted in [24] that a fourth-order formulation of the adjoint variable, derived by eliminating the control and state variables, can be employed to address Dirichlet boundary control problems.

Unlike continuous finite element approximations of the state variables, alternative discretization techniques, such as the mixed finite element method [3] and the hybridizable discontinuous Galerkin (HDG) method [25], have been employed to overcome the challenges arising from the variational formulation of continuous finite element methods. In [3], error estimates of order \(h^{1-1/s}\) with \(s \geq 2\) were established for polygonal domains and of order \(h|\ln h|\) for smooth domains, where the classical \(RT_0\) mixed finite element method was used for the approximation. In [25], an HDG method with discrete fluxes in \(P_k\) and potentials in \(P_{k+1}\) for \(k\ge 0\) was applied. Error estimates of order \(h^s\) for \(s < \min(3/2, \pi/\theta -1/2)\) for \(k\ge 1\), and of order \(h^{1/2}\) in the case \(k=0\), were reported. Here and throughout the paper, \(k\) denotes the polynomial degree of the finite element space.

All of the aforementioned studies concentrate on Dirichlet boundary control problems for elliptic equations. However, such control problems are also of significant importance in various applications governed by fluid dynamics models, including the Stokes [26], [27] and Navier–Stokes equations [28], [29]. Prior to analyzing the numerical behavior of complex optimal control problems governed by fluid dynamics models, it is essential to first study optimization problems constrained by convection–diffusion equations. Numerical studies for unconstrained Dirichlet boundary control problems governed by convection-diffusion equations have been carried out in [5], [30], [31] using an HDG method, and in [32] employing a symmetric interior penalty Galerkin approach. In the present work, we focus on the numerical analysis of Dirichlet boundary control problems governed by a convection–diffusion equation with \(L^2\)–boundary controls subject to pointwise bounds on the control. The local discontinuous Galerkin (LDG) method is employed as a discretization technique to address the variational challenges that arise when the control space is taken as \(L^2(\Gamma)\) and to leverage the strengths of discontinuous Galerkin methods for convection–diffusion equations. The LDG method is one of several discontinuous Galerkin approaches that have been extensively studied, particularly for convection–diffusion problems, due to its versatility in handling a wide range of applications and advantageous features such as local conservativity and strong locality. Like the HDG method introduced in [33], the LDG method rewrites higher-order PDEs as systems of first-order equations and employs numerical fluxes to weakly couple discontinuous solutions across element interfaces. However, the two methods differ fundamentally in how these fluxes are constructed and how global coupling is achieved. In LDG, numerical fluxes are defined directly between neighboring elements using information from both sides of each interface, allowing considerable flexibility (e.g., alternating or upwind choices) and making the method particularly robust for a wide range of problems, including convection–diffusion systems. This direct coupling avoids the introduction of additional globally coupled trace variables, as required in HDG, and keeps the formulation conceptually simpler. In contrast, HDG defines its numerical fluxes through a single-valued hybrid variable on element interfaces. This hybrid variable becomes the only globally coupled unknown, and the volumetric element degrees of freedom are recovered locally via static condensation, significantly reducing the size and bandwidth of the global system at the cost of increased formulation complexity. For studies on the use of discontinuous Galerkin methods in convection-dominated PDEs, we refer the reader to [34]–[38]. In the present work, we use the results from [39]–[41] to estimate the approximation errors of the state and adjoint variables. Finally, let us mention that optimal control problems governed by convection–diffusion equations were also investigated in [42]–[46].

The rest of this paper is organized as follows. In the following section, we discuss the regularity of the solutions and present the optimality system. Section 3 is devoted to the formulation of the LDG scheme for approximating the Dirichlet boundary control problem. In Section 4, we establish a priori error estimates for the LDG approximation on a convex polygonal domain. These estimates are derived for two distinct discretization strategies used in the optimal control problem: the variational discretization approach and the piecewise linear discretization approach. Section 5 presents numerical experiments to support and validate the theoretical findings. Finally, concluding remarks and further discussion are provided in Section 6.

2 Regularity and optimality system↩︎

Throughout this paper, we adopt the standard notation \(W^{m,p}(\Omega)\) (see, e.g., [47]) for Sobolev spaces equipped with the norm \(\| \cdot \|_{m,p,\Omega}\) and seminorm \(| \cdot |_{m,p,\Omega}\) for \(m \geq 0\) and \(1 \leq p \leq \infty\) on an open bounded domain \(\Omega\) in \(\mathbb{R}^2\) with boundary \(\Gamma=\partial \Omega\). We denote \(W^{m,2}(\Omega)\) by \(H^m(\Omega)\) with norm \(\| \cdot \|_{m,\Omega}\) and seminorm \(| \cdot |_{m,\Omega}\). The spaces of square-integrable functions over \(\Omega\) and \(\Gamma\) are represented by \(L^2(\Omega)\) and \(L^2(\Gamma)\), respectively, with norms denoted by \(\| \cdot \|_{0,\Omega}\) and \(\| \cdot \|_{0,\Gamma}\). The inner products on \(L^2 (\Omega)\) and \(L^2 (\Gamma)\) are defined by \[(v, w)_{\Omega} = \int_{\Omega} v \, w \, dx \qquad \forall v, w \in L^2(\Omega), \quad \text{and} \quad \langle v, w \rangle_{\Gamma} = \int_{\Gamma} v \, w \, ds \qquad \forall v, w \in L^2(\Gamma),\] respectively. Note that \(H^1_0(\Omega) = \{ v \in H^1(\Omega): v=0 \; on \; \Gamma \}\) and \(H(\mathrm{div},\Omega) := \left\{ \mathbf{v} \in (L^2(\Omega))^2: \nabla \cdot \mathbf{v} \in L^2(\Omega) \right\}\). Moreover, \(C\) denotes a generic positive constant independent of the mesh size \(h\) and may differ in various estimates. For clarity and convenience, we assume the following conditions on the domain \(\Omega\) and the given data, which are valid throughout the paper:

Owing to the low a priori regularity of the control \(u\in U:= L^2(\Gamma)\) in the optimal control problem 1 , the state equation 2 is understood in a very weak sense, i. e., we seek a solution \(y \in L^2(\Omega)\) satisfying \[\begin{align} \label{p4} -\epsilon (y, \Delta v )_\Omega - (y, \beta \cdot \nabla v )_\Omega + (\alpha y, v)_\Omega = (f,v)_\Omega - \epsilon \langle u, \partial_n v \rangle_{\Gamma} \end{align}\tag{4}\] for all \(v \in H_0^1 (\Omega) \cap H^2(\Omega)\) with \(\partial_n v := \nabla v \cdot {\mathbf{n}}\) and \({\mathbf{n}}\) representing the unit outward normal vector to \(\Gamma\). It is well known that problem 4 for \(f\in L^2(\Omega)\) and \(u \in L^2(\Gamma)\) admits a unique solution \(y\in H^{1/2}(\Omega)\), which satisfies \[\label{p5} \|y\|_{H^{1/2}(\Omega)} \leq C \big( \|f\|_{0,\Omega} + \|u\|_{0,\Gamma} \big),\tag{5}\] see, e.g., [4], [48], [49].

Based on the very weak formulation of the state equation 4 , the optimal control problem (1 )-(2 ) can be stated as follows

\[\label{p8} \underset{(y,u) \in L^2(\Omega) \times U^{ad}}{ minimize } \; J(y,u)subject to\eqref{p4}.\tag{6}\] Using the standard arguments in the present setting, it can be shown that problem 6 admits a unique solution \((y,u)\in H^1(\Omega)\times (H^{1/2}(\Gamma)\cap {U^{\mathrm{ad}}})\), which, together with the unique adjoint state \(z\in H^2(\Omega)\cap H^1_0(\Omega)\), weakly solves the following optimality system: \[\tag{7} \begin{equation}\tag{8} \begin{align} \nabla \cdot ( - \epsilon \nabla y + \beta y ) + \alpha y =& \, f \qquad \qquad \quad in\Omega, \\ y=& \,u \qquad \qquad \quad on\Gamma, \end{align} \end{equation} \begin{equation}\tag{9} \begin{align} \nabla \cdot ( - \epsilon \nabla z - \beta z ) + \alpha z = \,& y-y^d \qquad \quad in\Omega, \\ z=& \,0 \qquad \qquad \quad on\Gamma, \end{align} \end{equation} and \begin{equation}\tag{10} \langle \omega u - \epsilon \partial_n z, w -u \rangle_{\Gamma} \geq 0 \qquad \quad \quad \;\;\;\;\; \forall w \in {U^{\mathrm{ad}}}. \end{equation}\] Since \(z\in H^2(\Omega)\cap H^1_0(\Omega)\), it follows for convex polygonal domains from, e.g., [4] that \(\partial_n z \in H^{1/2}(\Gamma)\). The variational inequality 10 is equivalent to \(u = \mathcal{P}_{{U^{\mathrm{ad}}}} \left(\frac{\epsilon}{\omega} \, \partial_n z \right) \in H^{1/2}(\Gamma)\); compare, e.g., [49]. Here, \(\mathcal{P}_{U^{\mathrm{ad}}}\) denotes the orthogonal projection in \(L^2(\Gamma)\) onto the admissible set \({U^{\mathrm{ad}}}\). Since 8 admits a unique very weak solution \(y(u)\) for every \(u\in{U^{\mathrm{ad}}}\), the reduced cost functional \(\widehat J(u):= J(y(u),u)\) is well-defined. Furthermore, \(\widehat J'(u) = \omega u-\epsilon \partial_n z\), where \(z\) denotes the unique adjoint state satisfying 9 with \(y=y(u)\).

To apply the local discontinuous Galerkin (LDG) method for approximating the solution of the Dirichlet boundary control problem 1 –2 , we employ a mixed formulation of the optimality system 8 –10 . This is possible since the optimal control of problem 6 satisfies \(u\in H^{1/2}(\Gamma)\), and in this case it follows from 4 that \(\nabla y \in H(\mathrm{div},\Omega)\), and consequently \(\nabla y\cdot {\mathbf{n}}\in H^{-1/2}(\Gamma)\); compare [3], [25]. One now has that with the unique solution \((y,z,u)\) of the optimality system 8 –10 the functions \(y,{\mathbf{q}},z,{\mathbf{p}},u\), with \({\mathbf{q}}:=-\epsilon^{{\frac{1}{2}}} \nabla y\), \({\mathbf{p}}:= \epsilon^{{\frac{1}{2}}} \nabla z\), \({\mathbf{q}}, {\mathbf{p}}\in H(\mathrm{div},\Omega)\), form the unique weak solution of the mixed system given by \[\tag{11} \begin{align}\tag{12} \begin{aligned} \nabla \cdot (\beta\, y + \epsilon^{{\frac{1}{2}}} \mathbf{q}) + \alpha\, y &= f &&\text{in }\Omega,\\[2mm] \mathbf{q} &= -\,\epsilon^{{\frac{1}{2}}} \nabla y &&\text{in }\Omega,\\[2mm] y &= u &&\text{on }\Gamma . \end{aligned} \end{align} \begin{align}\tag{13} \begin{aligned} \nabla \cdot (-\beta\, z - \epsilon^{{\frac{1}{2}}} \mathbf{p}) + \alpha\, z &= y - y^{d} &&\text{in }\Omega,\\[2mm] \mathbf{p} &= \epsilon^{{\frac{1}{2}}} \nabla z &&\text{in }\Omega,\\[2mm] z &= 0 &&\text{on }\Gamma . \end{aligned} \end{align} \begin{align}\tag{14} \langle\, \omega u - \epsilon^{{\frac{1}{2}}} \mathbf{p}\cdot {\mathbf{n}},\, w - u\, \rangle_{\Gamma} \;\ge\; 0, \qquad \forall w \in {U^{\mathrm{ad}}}. \end{align}\]

To ensure the approximation properties of the LDG method, specific regularity assumptions on the involved states and fluxes are required. For the solution of our optimal Dirichlet boundary control problem on polygonal domains, a detailed study can be found in [19]. In this respect, we obtain the following theorem from e.g., [30], which is also valid in our control constrained case; compare the investigations in [19] related to the convex case.

3 Local discontinuous Galerkin formulation↩︎

Since we assume that the domain \(\Omega\) is polygonal, its boundary is exactly represented by triangle edges. We denote \(\{ {\mathcal{T}}_h\}_h\) as a family of shape-regular simplicial triangulations of \(\Omega\). Each mesh \({\mathcal{T}}_h\) consists of closed triangles such that \(\overline{\Omega} = \bigcup_{K \in {\mathcal{T}}_h} K\) holds. We assume that the mesh is regular in the following sense: for different triangles \(K_i, K_j \in {\mathcal{T}}_h\), \(i \not= j\), the intersection \(K_i \cap K_j\) is either empty, a vertex, or an edge, i.e., hanging nodes are not allowed. The diameter of an element \(K\) and the length of an edge \(E\) are denoted by \(h_{K}\) and \(h_E\), respectively, with \(h = \max \limits_{K \in {\mathcal{T}}_h} h_K\). The set of all edges \({\mathcal{E}}_h\) is split into the set \({\mathcal{E}}^0_h\) of interior edges and the set \({\mathcal{E}}^{\partial}_h\) of boundary edges, so that \({\mathcal{E}}_h={\mathcal{E}}^{0}_h \cup {\mathcal{E}}^{\partial}_h\). For the outward unit normal vector \({\mathbf{n}}\) on \(\Gamma\), the inflow and outflow parts of \(\Gamma\) are denoted by \(\Gamma^-\) and \(\Gamma^+\), respectively, \[\Gamma^- = \left\{{x \in \Gamma}\,:~{\beta\cdot {\mathbf{n}}< 0}\right\}, \quad \Gamma^+ = \left\{{x \in \Gamma}\,:~{\beta \cdot {\mathbf{n}}\geq 0}\right\}.\] Similarly, the boundaries of an element \(K \in {\mathcal{T}}_h\) can be categorized into inflow and outflow boundaries: \(\partial K^-\) and \(\partial K^+\).

Let \(E\) be a common edge for two elements \(K\) and \(K^e\). For a piecewise continuous scalar function \(y\), there are two traces of \(y\) along \(E\), denoted by \(y|_E\) from inside \(K\) and \(y^e|_E\) from inside \(K^e\). Then, the jump and average of \(y\) across the edge \(E\) are defined by: \[\begin{align} \left[\!\left[ y \right]\!\right]&=y|_E{\mathbf{n}}_{K}+y^e|_E{\mathbf{n}}_{K^e}, \quad \left\{\!\!\left\{ y \right\}\!\!\right\}=\frac{1}{2}\big( y|_E+y^e|_E \big), \end{align}\] where \({\mathbf{n}}_K\) (resp. \({\mathbf{n}}_{K^e}\)) denotes the outward unit normal vector to \(\partial K\) (resp. \(\partial K^e\)). Similarly, for a piecewise continuous vector field \(\mathbf{q}\), the jump and average across an edge \(E\) are defined, respectively, by \[\begin{align} \left[\!\left[ \mathbf{q} \right]\!\right]&=\mathbf{q}|_E \cdot {\mathbf{n}}_{K} + \mathbf{q}^e|_E \cdot {\mathbf{n}}_{K^e}, \quad \left\{\!\!\left\{ \mathbf{q} \right\}\!\!\right\}=\frac{1}{2}\big(\mathbf{q}|_E+\mathbf{q}^e|_E \big). \end{align}\] For a boundary edge \(E \in K \cap \Gamma\), we set \(\left\{\!\!\left\{ \mathbf{q} \right\}\!\!\right\}=\mathbf{q}\) and \(\left[\!\left[ y \right]\!\right]=y{\mathbf{n}}\), where \({\mathbf{n}}\) is the outward unit normal vector on \(\Gamma\). Note that the jump in \(y\) is vector-valued, whereas the jump in \(\mathbf{q}\) is scalar-valued which only involves the normal component of \(\mathbf{q}\).

To adapt the continuous mixed formulation to the LDG setting, we note that the weak solution of 12 also satisfies the following elementwise system on each element \(K \in {\mathcal{T}}_h\) \[\tag{15} \begin{eqnarray} -(\epsilon^{\frac{1}{2}}{\mathbf{q}}+ \beta y, \nabla v)_K + (\alpha y, v)_K + \langle ( \epsilon^{\frac{1}{2}}{\mathbf{q}}+ \beta y ) \cdot {\mathbf{n}}, v \rangle_{\partial K} &=& (f,v)_K \quad \;\; \forall v \in V, \tag{16} \\ ( {\mathbf{q}}, \mathbf{r} )_K - ( \epsilon^{\frac{1}{2}}y, \nabla \cdot \mathbf{ r})_K +\langle \epsilon^{\frac{1}{2}}u, \mathbf{r} \cdot {\mathbf{n}}\rangle_{\partial K} &=&0 \qquad \qquad \forall \mathbf{r} \in \mathbf{W}, \tag{17} \end{eqnarray}\] where \[\begin{align} \mathbf{W} :=& \, \left\{{\mathbf{w} \in \left( L^2(\Omega) \right)^2}\,:~{\mathbf{w}_K \in \left( H^1(K) \right)^2, \;\; \forall K \in {\mathcal{T}}_h}\right\}, \\ V :=& \, \left\{{v\in L^2(\Omega)}\,:~{v_K \in H^1(K), \;\; \forall K \in {\mathcal{T}}_h}\right\}. \end{align}\]

Next, we seek to approximate the state solution \((y, {\mathbf{q}})\) by functions \((y_h, {\mathbf{q}}_h)\) in the following finite element spaces \({\mathbf{W}}_h \times V_h \subset {\mathbf{W}}\times V\): \[\label{DG1} \begin{align} \mathbf{W}_h &= \left\{{\mathbf{w} \in \big(L^2(\Omega)\big)^2}\,:~{ \mathbf{w}\mid_{K}\in \big( \mathbb{S}^1(K) \big)^2, \quad \forall K \in {\mathcal{T}}_h}\right\}, \\[4pt] V_h &= \left\{{v \in L^2(\Omega)}\,:~{ v \mid_{K}\in \mathbb{S}^1(K), \quad \forall K \in {\mathcal{T}}_h}\right\}, \\[4pt] U_h &= \left\{{u \in L^2(\Gamma)}\,:~{ u \mid_{E}\in \mathbb{S}^1(E), \quad \forall E \in {\mathcal{E}}_h^{\partial}}\right\}, \end{align}\tag{18}\] where \(\mathbb{S}^1(K)\) (resp. \(\mathbb{S}^1(E)\)) denotes the local finite element space, which consists of linear polynomials on each element \(K\) (resp. on \(E\)). Then, for all \((v, {\mathbf{r}}) \in V_h \times {\mathbf{W}}_h\) the approximate solution \((y_h,{\mathbf{q}}_h)\) of the state solution \((y,{\mathbf{q}})\) satisfies \[\label{DG2} \begin{eqnarray} -(\epsilon^{\frac{1}{2}}{\mathbf{q}}_h + \beta y_h, \nabla v)_{K} + (\alpha y_h, v)_{K} + \langle ( \epsilon^{\frac{1}{2}}\widehat{{\mathbf{q}}}_h + \beta \widetilde{y}_h ) \cdot {\mathbf{n}}, v \rangle_{\partial K} &=& (f,v)_{K}, \\ ( {\mathbf{q}}_h, {\mathbf{r}})_{K} - (\epsilon^{\frac{1}{2}}y_h, \nabla \cdot {\mathbf{r}})_{K} + \langle \epsilon^{\frac{1}{2}}\widehat{y}_h, {\mathbf{r}}\cdot {\mathbf{n}}\rangle_{ \partial K} &=&0. \end{eqnarray}\tag{19}\] The numerical fluxes, denoted by \(\widehat{{\mathbf{q}}}_h, \widetilde{y}_h\), and \(\widehat{y}_h\), must be chosen appropriately to ensure the stability of the method and to improve its accuracy.

We are now ready to introduce the expressions that define the numerical fluxes. The numerical traces of \(y_h\) associated with the diffusion and convection terms are characterized, respectively, by \[\label{DG3} \widehat{y}_h = \left\{ \begin{array}{ll} \left\{\!\!\left\{ y_h \right\}\!\!\right\} + \mathbf{C}_{12} \cdot \left[\!\left[ y_h \right]\!\right], & E \in {\mathcal{E}}^{0}_h, \\ u_h, & E \in {\mathcal{E}}^{\partial}_h, \end{array} \right. \quad and \quad \widetilde{y}_h =\left\{ \begin{array}{ll} u_h, & E \in \Gamma^{-}, \\ \left\{\!\!\left\{ y_h \right\}\!\!\right\} + \mathbf{D}_{11} \cdot \left[\!\left[ y_h \right]\!\right], & E \in {\mathcal{E}}^{0}_h, \\ y_h, & E \in \Gamma^+. \end{array} \right.\tag{20}\] From 20 , it follows that the numerical trace \(\widetilde{y}_h\) aligns with the traditional upwinding trace. Moreover, the numerical flux \(\widehat{{\mathbf{q}}}_h\) is given by \[\begin{align} \label{DG4} \widehat{{\mathbf{q}}}_h &=& \left\{ \begin{array}{ll} \left\{\!\!\left\{ {\mathbf{q}}_h \right\}\!\!\right\} - C_{11}\left[\!\left[ y_h \right]\!\right] - \mathbf{C}_{12} \left[\!\left[ {\mathbf{q}}_h \right]\!\right], & E \in {\mathcal{E}}^{0}_h, \\ {\mathbf{q}}_h - C_{11} (y_h -u_h) {\mathbf{n}}, & E \in {\mathcal{E}}^{\partial}_h. \end{array} \right. \end{align}\tag{21}\] Here, we set \(C_{11} = \epsilon/h_E\) for each \(E \in {\mathcal{E}}_h\), choose \(\mathbf{C}_{12}\) such that \(\mathbf{C}_{12} \cdot {\mathbf{n}}= \frac{1}{2} \mathop{\mathrm{sign}}\big( {\mathbf{n}}\cdot \mathbf{v} \big)\) for a nonzero arbitrary but fixed vector \(\mathbf{v}\), and define the vector-valued function \(\mathbf{D}_{11}\) by \[\mathbf{D}_{11} \cdot {\mathbf{n}}= \frac{1}{2} \mathop{\mathrm{sign}}\big( {\mathbf{n}}\cdot \beta \big).\] Note that the auxiliary variable \(\mathbf{v}\) need not be associated with the convective velocity \(\boldsymbol{\beta}\); it may be chosen independently without affecting the convergence properties of the method. Nevertheless, when \(\boldsymbol{\beta}\) is not identically zero on any element \(K \in \mathcal{T}_h\), it is natural and often convenient to relate \(\mathbf{v}\) to \(\boldsymbol{\beta}\). For a more detailed discussion, we refer the reader to [39], [40], [46], [50].

Inserting the numerical fluxes in 20 –21 into 19 and summing over all elements, we obtain \[\begin{align} &-\sum\limits_{K \in {\mathcal{T}}_h} \int_{K} (\epsilon^{{\frac{1}{2}}} {\mathbf{q}}_h + \beta y_h) \cdot \nabla v \, dx + \sum\limits_{E \in {\mathcal{E}}_h^0} \int_E \epsilon^{{\frac{1}{2}}} \left( \left\{\!\!\left\{ {\mathbf{q}}_h \right\}\!\!\right\} - C_{11} \left[\!\left[ y_h \right]\!\right] - \mathbf{C}_{12} \left[\!\left[ {\mathbf{q}}_h \right]\!\right] \right) \cdot \left[\!\left[ v \right]\!\right] \, ds \\ &\quad + \sum\limits_{K \in {\mathcal{T}}_h} \int_{K} \alpha y_h \, v \, dx + \sum\limits_{E \in {\mathcal{E}}_h^0} \int_E \big( \left\{\!\!\left\{ y_h \right\}\!\!\right\} + \mathbf{D}_{11} \cdot \left[\!\left[ y_h \right]\!\right] \big) \beta \cdot \left[\!\left[ v \right]\!\right] \, ds + \sum\limits_{E \subset \Gamma^+} \int_E ({\mathbf{n}}\cdot \beta) y_h \, v \, ds \\ &\quad + \sum\limits_{E \in {\mathcal{E}}_h^{\partial}} \int_E \epsilon^{{\frac{1}{2}}} \big( {\mathbf{q}}_h \cdot {\mathbf{n}}- C_{11} y_h \big) v \, ds = \int_{\Omega} f \, v \, dx - \sum\limits_{E \in {\mathcal{E}}_h^{\partial}} \int_E \epsilon^{{\frac{1}{2}}} C_{11} u_h \, v \, ds + \sum\limits_{E \subset \Gamma^-} \int_E |\beta \cdot {\mathbf{n}}| u_h \, v \, ds, \end{align}\]

and

\[\begin{align} &\int_{\Omega} {\mathbf{q}}_h \cdot {\mathbf{r}}\, dx - \sum\limits_{K \in {\mathcal{T}}_h} \int_{K} \epsilon^{{\frac{1}{2}}} y_h \, \nabla \cdot {\mathbf{r}}\, dx + \sum\limits_{E \in {\mathcal{E}}_h^0} \int_E \epsilon^{{\frac{1}{2}}} \big( \left\{\!\!\left\{ y_h \right\}\!\!\right\} + \mathbf{C}_{12} \cdot \left[\!\left[ y_h \right]\!\right] \big) \left[\!\left[ {\mathbf{r}} \right]\!\right] \, ds \\ &\quad = -\sum\limits_{E \in {\mathcal{E}}_h^{\partial}} \int_E \epsilon^{{\frac{1}{2}}} u_h \, {\mathbf{r}}\cdot {\mathbf{n}}\, ds. \end{align}\] It is convenient to introduce the following bi(linear) forms: \[\begin{align} a_h({\mathbf{q}}, {\mathbf{r}}) :=& \int_{\Omega} {\mathbf{q}}\cdot {\mathbf{r}}\, dx, \\[0.3em] b_h(y, {\mathbf{r}}) :=& -\sum\limits_{K \in {\mathcal{T}}_h} \int_{K} \epsilon^{{\frac{1}{2}}} y \, \nabla \cdot {\mathbf{r}}\, dx + \sum\limits_{E \in {\mathcal{E}}_h^0} \int_E \epsilon^{{\frac{1}{2}}} \big( \left\{\!\!\left\{ y \right\}\!\!\right\} + \mathbf{C}_{12} \cdot \left[\!\left[ y \right]\!\right] \big) \left[\!\left[ {\mathbf{r}} \right]\!\right] \, ds, \\[0.3em] c_h(y, v) :=& \sum\limits_{K \in {\mathcal{T}}_h} \int_{K} \big( \alpha y \, v - y \, \beta \cdot \nabla v \big) \, dx + \sum\limits_{E \in {\mathcal{E}}_h^0} \int_E \big( \left\{\!\!\left\{ y \right\}\!\!\right\} + \mathbf{D}_{11} \cdot \left[\!\left[ y \right]\!\right] \big) \beta \cdot \left[\!\left[ v \right]\!\right] \, ds \\ & - \sum\limits_{E \in {\mathcal{E}}_h^0} \int_E \epsilon^{{\frac{1}{2}}} C_{11} \left[\!\left[ y \right]\!\right] \cdot \left[\!\left[ v \right]\!\right] \, ds + \sum\limits_{E \subset \Gamma^+} \int_E ({\mathbf{n}}\cdot \beta) y \, v \, ds - \sum\limits_{E \in {\mathcal{E}}_h^{\partial}} \int_E \epsilon^{{\frac{1}{2}}} C_{11} y \, v \, ds, \end{align}\] \[\begin{align} m_{h,1}(u, {\mathbf{r}}) :=& -\sum\limits_{E \in {\mathcal{E}}_h^{\partial}} \int_E \epsilon^{{\frac{1}{2}}} u \, {\mathbf{r}}\cdot {\mathbf{n}}\, ds, \\[0.3em] m_{h,2}(u, v) :=& -\sum\limits_{E \in {\mathcal{E}}_h^{\partial}} \int_E \epsilon^{{\frac{1}{2}}} C_{11} u \, v \, ds + \sum\limits_{E \subset \Gamma^-} \int_E |\beta \cdot {\mathbf{n}}| u \, v \, ds, \\[0.3em] F(v) :=& \int_{\Omega} f \, v \, dx. \end{align}\] By applying integration by parts to the first term in \(b_h(\cdot,\cdot)\), we obtain \[\begin{align} b_h(y, {\mathbf{r}}) =& -\sum\limits_{K \in {\mathcal{T}}_h} \int_{K} \epsilon^{{\frac{1}{2}}} y \, \nabla \cdot {\mathbf{r}}\, dx + \sum\limits_{E \in {\mathcal{E}}_h^0} \int_E \epsilon^{{\frac{1}{2}}} \big( \left\{\!\!\left\{ y \right\}\!\!\right\} + \mathbf{C}_{12} \cdot \left[\!\left[ y \right]\!\right] \big) \left[\!\left[ {\mathbf{r}} \right]\!\right] \, ds \\ =& \sum\limits_{K \in {\mathcal{T}}_h} \int_{K} \epsilon^{{\frac{1}{2}}} \nabla y \cdot {\mathbf{r}}\, dx - \sum\limits_{K \in {\mathcal{T}}_h} \int_{\partial K} \epsilon^{{\frac{1}{2}}} y \, {\mathbf{r}}\cdot {\mathbf{n}}\, ds \\ &\quad + \sum\limits_{E \in {\mathcal{E}}_h^0} \int_E \epsilon^{{\frac{1}{2}}} \big( \left\{\!\!\left\{ y \right\}\!\!\right\} + \mathbf{C}_{12} \cdot \left[\!\left[ y \right]\!\right] \big) \left[\!\left[ {\mathbf{r}} \right]\!\right] \, ds. \end{align}\] Then, by the following straightforward computation \[\sum\limits_{K \in {\mathcal{T}}_h} \int_{\partial K} \epsilon^{{\frac{1}{2}}} y \, {\mathbf{r}}\cdot {\mathbf{n}}\, ds = \sum\limits_{E \in {\mathcal{E}}_h^0 \cup {\mathcal{E}}_h^{\partial}} \int_E \left\{\!\!\left\{ {\mathbf{r}} \right\}\!\!\right\} \cdot \left[\!\left[ \epsilon^{{\frac{1}{2}}} y \right]\!\right] \, ds + \sum\limits_{E \in {\mathcal{E}}_h^0} \int_E \left[\!\left[ {\mathbf{r}} \right]\!\right] \cdot \left\{\!\!\left\{ \epsilon^{{\frac{1}{2}}} y \right\}\!\!\right\} \, ds,\] we get \[\begin{align} b_h(y, {\mathbf{r}}) =& \sum\limits_{K \in {\mathcal{T}}_h} \int_{K} \epsilon^{{\frac{1}{2}}} \nabla y \cdot {\mathbf{r}}\, dx - \sum\limits_{E \in {\mathcal{E}}_h^0} \int_E \epsilon^{{\frac{1}{2}}} \big( \left\{\!\!\left\{ {\mathbf{r}} \right\}\!\!\right\} - \mathbf{C}_{12} \left[\!\left[ {\mathbf{r}} \right]\!\right] \big) \cdot \left[\!\left[ y \right]\!\right] \, ds \\ &\quad - \sum\limits_{E \in {\mathcal{E}}_h^{\partial}} \int_E \epsilon^{{\frac{1}{2}}} y \, {\mathbf{r}}\cdot {\mathbf{n}}\, ds. \end{align}\]

For given right-hand side \(f\) and boundary condition \(u\), the LDG approximation of the unique solution \((y,{\mathbf{q}})\) to 12 is given by \((y_h,{\mathbf{q}}_h) \in V_h\times {\mathbf{W}}_h\) satisfying \[\label{LDG-approx} \begin{align} a_h({\mathbf{q}}_h, {\mathbf{r}}) + b_h(y_h, {\mathbf{r}}) &= m_{h,1}(u, {\mathbf{r}}) && \forall {\mathbf{r}}\in {\mathbf{W}}_h, \\ -b_h(v, {\mathbf{q}}_h) + c_h(y_h, v) &= m_{h,2}(u, v) + F(v) && \forall v \in V_h. \end{align}\tag{22}\] For given \(u\in L^2(\Gamma)\) and \(f\in L^2(\Omega)\), this system admits a unique solution, which can be established using arguments similar to those in [39].

Hence, the LDG approximation scheme of the Dirichlet boundary control problem 1 –2 reads as \[\label{DG6} \underset{u_h \in {U^{\mathrm{ad}}_h}, \, (y_h, {\mathbf{q}}_h) \in V_h \times {\mathbf{W}}_h}{minimize} \; J(y_h,u_h) = \frac{1}{2}\|y_h - y^d\|^{2}_{0,\Omega} + \frac{\omega}{2} \|u_h\|^{2}_{0,\Gamma}\tag{23}\] subject to \[\label{DG5} \begin{align} a_h({\mathbf{q}}_h, {\mathbf{r}}) + b_h(y_h, {\mathbf{r}}) &= m_{h,1}(u_h, {\mathbf{r}}) && \forall {\mathbf{r}}\in {\mathbf{W}}_h, \\ -b_h(v, {\mathbf{q}}_h) + c_h(y_h, v) &= m_{h,2}(u_h, v) + F(v) && \forall v \in V_h. \end{align}\tag{24}\]

Here, the set of admissible controls is chosen either as \({U^{\mathrm{ad}}_h}:= {U^{\mathrm{ad}}}\), which is the case of variational discretization [1], or as \({U^{\mathrm{ad}}_h}:={U^{\mathrm{ad}}}\cap U_h\), which we refer to as the case of full discretization. Since \({U^{\mathrm{ad}}_h}\) is nonempty, closed, and convex and \(J\) is uniformly convex, problem 23 24 admits a unique solution \(((y_h,{\mathbf{q}}_h),u_h)\in (V_h\times {\mathbf{W}}_h)\times {U^{\mathrm{ad}}_h}\). Moreover, it follows from standard arguments that \(((y_h,{\mathbf{q}}_h),u_h)\) is the unique solution of the optimal control problem 23 24 if and only if a pair \(({\mathbf{p}}_h,z_h)\in {\mathbf{W}}_h\times V_h\) of discrete adjoint flux and discrete adjoint potential exists, which together with \(((y_h,{\mathbf{q}}_h),u_h,({\mathbf{p}}_h,z_h))\) solves the following discrete optimality conditions \[\tag{25} \begin{align} a_h({\mathbf{q}}_h, {\mathbf{r}}) + b_h(y_h, {\mathbf{r}}) &= m_{h,1}(u_h, {\mathbf{r}}) && \forall {\mathbf{r}}\in {\mathbf{W}}_h, \tag{26} \\[0.3em] -b_h(v, {\mathbf{q}}_h) + c_h(y_h, v) &= m_{h,2}(u_h, v) + F(v) && \forall v \in V_h, \tag{27} \\[0.3em] a_h({\mathbf{p}}_h, \psi) - b_h(z_h, \psi) &= 0 && \forall \psi \in {\mathbf{W}}_h, \tag{28} \\[0.3em] b_h(\phi, {\mathbf{p}}_h) + c_h(\phi, z_h) &= (y_h - y^d, \phi)_\Omega && \forall \phi \in V_h, \tag{29} \\[0.3em] \langle \omega u_h + m'_{1,h}(u_h, {\mathbf{p}}_h) + m'_{2,h}(u_h, z_h), \, w - u_h \rangle_\Gamma &\ge 0 && \forall w \in {U^{\mathrm{ad}}_h}. \tag{30} \end{align}\] The variational inequality in 30 in terms of \({\mathbf{q}}_h\) and \(z_h\) reads \[\langle \omega u_h - \epsilon^{{\frac{1}{2}}} {\mathbf{p}}_h \cdot {\mathbf{n}} + \kappa_z z_h, \, w - u_h \rangle_{\Gamma} \ge 0 \qquad \forall w \in {U^{\mathrm{ad}}_h},\] where \(\kappa_z = -\epsilon^{{\frac{1}{2}}} C_{11} + \chi_{\Gamma^-} |\beta \cdot {\mathbf{n}}|\) with \(\chi_{\Gamma^-}\) denoting the indicator function of \(\Gamma^-\).

If the optimal solution \(\big( y, {\mathbf{q}}, z, {\mathbf{p}}, u \big) \in V \times {\mathbf{W}}\times V \times {\mathbf{W}}\times {U^{\mathrm{ad}}}\) is smooth enough, it also satisfies the system \[\label{sys:optimality95cont} \begin{align} a_h({\mathbf{q}}, {\mathbf{r}}) + b_h(y, {\mathbf{r}}) &= m_{h,1}(u, {\mathbf{r}}) && \forall {\mathbf{r}}\in {\mathbf{W}}, \\[0.3em] -b_h(v, {\mathbf{q}}) + c_h(y, v) &= m_{h,2}(u, v) + F(v) && \forall v \in V, \\[0.3em] a_h({\mathbf{p}}, \psi) - b_h(z, \psi) &= 0 && \forall \psi \in {\mathbf{W}}, \\[0.3em] b_h(\phi, {\mathbf{p}}) + c_h(\phi, z) &= (y - y^d, \phi)_\Omega && \forall \phi \in V, \\[0.3em] \langle \omega u + m'_{1,h}(u, {\mathbf{p}}) + m'_{2,h}(u, z), \, w - u \rangle_\Gamma &\ge 0 && \forall w \in {U^{\mathrm{ad}}}. \end{align}\tag{31}\]

4 Error estimates↩︎

In this section, we present error estimates for the LDG approximation of the Dirichlet boundary control problem posed on a convex polygonal domain. The analysis is carried out for two different discretization strategies applied to the optimal control problem: the variational discretization approach and the piecewise linear discretization approach.

Let us introduce the (linear) control-to-state operator \(S : L^2(\Gamma) \to H^{{\frac{1}{2}}}(\Omega),\) defined by \(u \mapsto S(u) := y\), where for a given \(u\) the state \(y\) satisfies 4 . Using this, the reduced cost functional can be written as \[\widehat J : U \to \mathbb{R}, \qquad \widehat J(u) = J(S u, u).\] Hence, the optimal control problem (1 )-(2 ) can be reformulated as \[\label{eqn:reduced95problem} \min_{u \in {U^{\mathrm{ad}}}} \widehat J(u).\tag{32}\] Following standard arguments in [51], the existence and uniqueness of an optimal solution \(u \in {U^{\mathrm{ad}}}\) can be established by invoking the continuity and strict convexity of \(\widehat J\). For a control \(u \in U\) and a direction \(\delta u \in U\), the directional derivatives are given by \[\label{continuous95derivative} \widehat J'(u)(\delta u) = (y - y^d, y(\delta u))_\Omega + \omega \langle u, \delta u \rangle_{\Gamma} = \langle\, \omega u - \epsilon^{{\frac{1}{2}}} \mathbf{p}\cdot {\mathbf{n}},\, \delta u\, \rangle_{\Gamma},\tag{33}\] where \(y\) is the state associated with the control \(u\), \(y(\delta u)\) denotes the state assigned to the control \(\delta u\), and \({\mathbf{p}}\) is the adjoint flux associated with the observation \(y-y^d\). The necessary optimality condition of the reduced problem 32 is given as a variational inequality: \[\label{reduced95cont95variational} \widehat J'(u)(\delta u - u) \geq 0 \quad \forall \delta u \in {U^{\mathrm{ad}}}.\tag{34}\]

To discretize the reduced optimal control problem 32 , we define the discrete reduced cost functional as \[\widehat{J}_h : U \to \mathbb{R}, \qquad \widehat{J}_h(u) = J(S_h u,u),\] where \(S_h : u \mapsto y_h(u)\) represents the discrete solution operator, which, according to [39], is well-defined. The corresponding discretized reduced problem is then given by \[\label{eqn:discrete95reduced95problem} \min_{u \in {U^{\mathrm{ad}}_h}} \; \widehat{J}_h(u).\tag{35}\] Analogous to the continuous case, the discrete reduced cost functional \(\widehat{J}_h\) is quadratic. Its directional derivative is given by \[\label{discrete95derivative} \widehat{J}_h'(u)(\delta u) = (y_h(u) - y^d,\, y_h(\delta u))_{\Omega} + \omega \langle u,\, \delta u \rangle_{\Gamma} = \langle \omega u - \epsilon^{{\frac{1}{2}}} {\mathbf{p}}_h \cdot {\mathbf{n}} + \kappa_z z_h, \, \delta u \rangle_{\Gamma},\tag{36}\] where \(y_h(u)\), \(y_h(\delta u)\) denote discrete states associated with \(u\) and \(\delta u\), respectively, and \({\mathbf{p}}_h\) and \(z_h\) denote the adjoint flux and potential associated with \(y_h(u)-y^d\).

In this work, we consider two approaches for the discretization of the control variable: the variational discretization proposed in [1], and the piecewise linear discretization. For the linear discretization, we employ the discrete space defined in 18 , while in the case of variational discretization, the discrete control space is chosen as \(U_h = U = L^2(\Gamma)\). In both settings, the discrete admissible control set is given by \({U^{\mathrm{ad}}_h}= U_h \cap {U^{\mathrm{ad}}}\). Furthermore, for the variational discretization approach, the discrete optimal control \(u_h\) is characterized by \[\label{reduced95disc95var95variational} \widehat{J}_h'(u_h)(\delta u - u_h) \ge 0, \qquad \forall \delta u \in {U^{\mathrm{ad}}},\tag{37}\] whereas for the piecewise linear discretization approach, the discrete optimal control \(u_h\) satisfies \[\label{reduced95disc95variational} \widehat{J}_h'(u_h)(\delta u_h - u_h) \ge 0, \qquad \forall \delta u_h \in {U^{\mathrm{ad}}_h}.\tag{38}\] In the analysis, we require an interpolation/projection operator that preserves admissibility; however, the \(L^2\)-projection does not, in general, preserve the admissibility of functions in \({U^{\mathrm{ad}}}\). Hence, we employ the quasi-interpolation operator \(\pi_h : L^1(\partial \Omega) \to \mathcal{E}_h^\partial\) inspired by [52], and used in the context of Dirichlet boundary control, see, e.g., [4]. The operator \(\pi_h\) is defined as \[\label{def:quasi95interpolation} \pi_h u = \sum_{n \in \mathcal{N}_h^\partial} \frac{(u, \xi_n)}{(1, \xi_n)} \, \xi_n,\tag{39}\] where \(\mathcal{N}_h^\partial\) denotes the set of all mesh nodes of \(\mathcal{T}_h\) lying on the boundary \(\Gamma\), and \(\xi_n \in U_h\) represents the local (nodal) basis function associated with node \(n\) on the boundary. By construction, it follows that \(\pi_h u \in {U^{\mathrm{ad}}_h}\) for every \(u \in {U^{\mathrm{ad}}}\). For \(u \in H^{s}(\Gamma)\) with \(0 \leq s \leq 1\), the following estimates also hold [4] \[\tag{40} \begin{align} \| u - \pi_h u \|_{0,\Gamma} &\leq C h^{s} \| u \|_{s,\Gamma}, \tag{41} \\[4pt] \| u - \pi_h u \|_{-s,\Gamma} &\leq C h^{2s} \| u \|_{s,\Gamma}. \tag{42} \end{align}\]

In the numerical error analysis, it will be convenient to use the adjoint-based representations of the respective derivatives \(\widehat J'\) and \(\widehat J_h'\). In particular, we write \[\label{adj-deriv-z} \widehat J'(u)(\delta u) = \langle\, \omega u - \epsilon^{{\frac{1}{2}}} \mathbf{p}\cdot {\mathbf{n}}+ \kappa_z z,\, \delta u\, \rangle_{\Gamma}.\tag{43}\] This representation is valid since \(z=0\) on \(\Gamma\).

To facilitate the analysis, we also introduce the following auxiliary problem: find \((z_h(y), {\mathbf{p}}_h(y)) \in V_h \times {\mathbf{W}}_h\) such that \[\label{aux2} \begin{align} a_h({\mathbf{p}}_h(y), \psi) - b_h(z_h(y), \psi) &= 0, && \forall \psi \in {\mathbf{W}}_h, \\[4pt] b_h(\phi, {\mathbf{p}}_h(y)) + c_h(\phi, z_h(y)) &= (y - y^d, \phi)_\Omega, && \forall \phi \in V_h. \end{align}\tag{44}\] Here, \((z_h(y), {\mathbf{p}}_h(y))\) denotes the LDG approximation of \((z, {\mathbf{p}})\).

We begin by establishing an error estimate for the variational discretization.

For the fully discrete approach we have

4.1 Discussion of Theorem 2 and Theorem 3:↩︎

In both cases, namely the variational discrete and the fully discrete setting, to obtain error estimates for the optimal control and the associated state we need LDG error estimates for the state, the adjoint flux as well as for the adjoint potential under natural regularity requirements on the continuous optimal control \(u\) and the data \(y^d\) (and \(f\)).

Let us first recall what is known in the literature with respect to the approximation properties of LDG approximations. In the first instance we here use results of [39]. In particular, from Theorem 2.1 and Tables 2.1 and 4.1 therein, one obtains for the case \(k=1\):

If one inspects the proof of this lemma, one recognizes that higher regularity of \((v,{\mathbf{w}})\in H^{s+2}(\Omega)\times H^{s+2}(\Omega)^2\) for \(s\in (0,1)\) in our case \(k=1\) does not affect the error estimate. If we look at the variables of our optimization problem we observe that for \(y^d\in L^2(\Omega)\) (and \(u\in H^{1/2}(\Gamma)\), which is the minimal regularity we obtain for the control \(u\)), the adjoint potential satisfies \(z\in H^2(\Omega)\cap H^1_0(\Omega)\) with flux \({\mathbf{p}}\in H^1(\Omega,\mathbb{R}^2)\). In this case, we conclude that \[\kappa_z\|z(u)-z_h(u)\|_{0,\Gamma} \lesssim h^{1/2}, \quad \|{\mathbf{p}}(u)\cdot {\mathbf{n}}-{\mathbf{p}}_h(u) \cdot {\mathbf{n}}\|_{0,\Gamma} \lesssim h^{1/2},\] since in our LDG method \(C_{11} \sim \epsilon h^{-1}\). Moreover, with Lemma 1 we for higher regularity of the involved variables according to Theorem 1 may not expect an improvement for the case \(k=1\). But we only have \(y\in H^1(\Omega)\) with flux \({\mathbf{q}}\in H(\mathrm{div},\Omega)\). With this, Lemma 1 is not applicable. However, in Lemma 4.3 of the recent paper [41] treating stability and convergence of the HDG method in approximating solutions to elliptic PDEs, it among other things is shown that for \(s\in [0,1]\) \[\label{L2est} \|y-y_h(u)\|_{0,\Omega} \le Ch^{s+1} |y|_{H^{s+1}(\Omega)}\tag{52}\] holds. Assuming that this estimate can also be established within the LDG framework, one may then expect a corresponding error bound \[\| u - u_h \|_{0,\Gamma} + \|y-y_h\|_{0,\Omega} \lesssim h + h^{1/2} + h^{1/2} \lesssim h^{1/2}.\] This compares very well with the results of [5] for the case \(k=0\), that is, piecewise constant flux and piecewise linear potential approximation.

In the fully discrete case, additional terms arise that need to be estimated. For the error associated with the quasi-interpolation operator \(\pi_h\), one has for \(u\in H^s(\Gamma)\) and \(0 \le s \le 1\) from 41 42 \[\|u-\pi_hu\|_{0,\Gamma} \lesssim h^s \quad \text{ and } \quad \|u-\pi_hu\|_{-s,\Gamma} \lesssim h^{2s}.\] For the generic case \(s=\frac{1}{2}\), this gives the convergence behavior \(h^{1/2}\) and \(h^1\), respectively. We also need to estimate the approximation errors for \({\mathbf{p}}(y_h)\) and \(z(y_h)\), where \(y_h\) denotes the discrete optimal potential. Since \(y_h\in L^2(\Omega)\) with \[\|y_h\|_{0,\Omega} \le C\] uniformly in \(h\) (this directly follows from \(\|y_h\|_{0,\Omega} \le \sqrt{4\widehat J_h(0) + 2\|y^d\|_{0,\Omega}^2}\)) we conclude from Lemma 1 \[\|({\mathbf{p}}_h - {\mathbf{p}}(y_h))\cdot{\mathbf{n}}\|_{0,\Gamma}+\kappa_z\|z_h - z(y_h)\|_{0,\Gamma} \lesssim h^{1/2}.\]

Finally, since \({\mathbf{p}}\cdot {\mathbf{n}}\in H^{1/2}(\Gamma)\) (see, e.g., [19]) we may apply 42 with \(s=\frac{1}{2}\). All together, we for \(u\in H^{1/2}(\Gamma)\) also in the fully discrete case obtain the error estimate \[\| u - u_h \|_{0,\Gamma} + \|y-y_h\|_{0,\Omega} \lesssim h^{1/2}.\]

To conclude, for both control discretization strategies in the case \(k=1\), we only obtain convergence estimates with leading factor \(h^{1/2}\), so that the variational discretization concept in fact simplifies the numerical analysis, but does not lead to improved convergence estimates, as it is frequently observed in PDE-constrained optimization [2].

5 Numerical Experiments↩︎

In this section, we provide numerical experiments to demonstrate the performance of the local discontinuous Galerkin discretization and to underpin the theoretical findings established for the Dirichlet boundary control problem. All numerical results are obtained using piecewise linear approximations for the state \(y\), the state flux \({\mathbf{q}}\), the adjoint \(z\), the adjoint flux \({\mathbf{p}}\), and the control \(u\). Unless otherwise stated, the initial mesh is constructed by partitioning the domain \(\Omega\) into a \(2 \times 2\) array of uniform squares, each of which is subsequently subdivided into two triangles. The resulting triangulation is then uniformly refined by subdividing each triangle into four congruent sub-triangles; see Fig. 1. To solve the discretized control constrained problem, the primal-dual active set algorithm is applied as a semismooth Newton method; see, e.g., [53]. The optimization procedure is terminated when two consecutive active sets coincide. Moreover, the experimental order of convergence is computed using \[rate = \frac{1}{\ln 2} \ln \left( \frac{\|e(h)\|_{0,\Omega}}{\|e(h/2)\|_{0,\Omega}}\right),\] where \(e(h)\) denotes the error associated with a triangulation of mesh size \(h\). Although not discussed in Section 4.1, for the convenience of the reader, in all tables in addition to \(\|y-y_h\|_{0,\Omega}\) and \(\|u-u_h\|_{0,\Gamma}\) we also report the numerically observed convergence orders for \(\|z-z_h\|_{0,\Gamma}\) and \(\|({\mathbf{p}}-{\mathbf{p}}_h)\cdot \mathbf{n}\|_{0,\Gamma}\).

5.1 Example 1↩︎

Our first example is a modified form of the elliptic problem in [17] and represents an unconstrained problem with analytical solutions. The problem data are chosen as \[\Omega=(0,1) \times (0,1), \quad \beta = (1,1)^T, \quad \alpha=1, \quad \omega=1.\] The source function \(f(x_1, x_2)\) and the desired state \(y^d(x_1, x_2)\) are generated so that the analytical solutions for the state \(y\), adjoint \(z\), and control \(u\) are given by \[\begin{align} y(x_1,x_2) &=& - \frac{\epsilon^{1/2}}{\omega} \left( x_1 (1-x_1) + x_2 (1-x_2) \right), \\ z(x_1,x_2) &=& \epsilon^{-1/2}x_1 x_2(1-x_1)(1-x_2), \\ u(x_1,x_2) &=& - \frac{\epsilon^{1/2}}{\omega} \left( x_1 (1-x_1) + x_2 (1-x_2) \right), \end{align}\] respectively.

In the present example we have a polynomial exact solution of the optimal control problem. For the LDG approximation errors of these functions, we may expect the best possible order in terms of the mesh size of the finite element mesh. Inspection of the errors on the right-hand side of estimate ?? shows that then, the limiting factors in the case \(\epsilon=1\) are the terms \(\|({\mathbf{p}}-{\mathbf{p}}_h(u))\cdot \mathbf{n}\|_{0,\Gamma}\) and \(\kappa_z\|z-z_h(u)\|_{0,\Gamma}\), where \(\kappa_z = \frac{-\epsilon^{3/2}}{h}+\chi_{\Gamma^-}|\beta\cdot{\mathbf{n}}|\). Numerical tests (not explicitly displayed here) performed for \(\epsilon \in [10^{-8},1]\) show that \(\|({\mathbf{p}}-{\mathbf{p}}_h(u))\cdot \mathbf{n}\|_{0,\Gamma}\) converges with order one, and \(\|z-z_h(u)\|_{0,\Gamma}\) with order two. This behavior is also observed for the approximation of the state \(y\) and the flux \({\mathbf{q}}\). It is also expected that the projection errors of \(u\) converge with best possible order. This altogether for \(\epsilon =1\) could explain the convergence order of one for the \(L^2(\Gamma)-\)error in the optimal control. The \(L^2-\)error in the optimal state however seems to converge with the higher order of \(3/2\). If one compares to the expected rates in the classical finite element approach taken in e.g., [17], we observe that our LDG errors in the optimal control with order one and optimal state with order 3/2 behave as the convergence orders proven there. We observe this behaviour of the optimal LDG approximation also in the subsequent examples. With decreasing \(\epsilon\) the convergence behaviour of control and state firstly start to improve and then for further decrease of \(\epsilon\) become irregular. This may be explained by the fact that for very small \(\epsilon\), one approximates zero functions. In contrast, the convergence behaviour of the adjoint state remains stable.

Fig. 2 shows the computed solutions of the state \(y\), adjoint \(z\), and control \(u\) on a fine mesh with 524288 elements. Tab. 1 and Tab. 2 present the \(L^2(\Omega)\) errors of the state variable \(y\), as well as the \(L^2(\Gamma)\) errors of the control \(u\), the adjoint state \(z\), and the adjoint flux \({\mathbf{p}}\), corresponding to \(\epsilon = 1\) and \(\epsilon = 10^{-6}\), respectively.

5.2 Example 2↩︎

The following example, posed on the unit square \(\Omega=(0,1) \times (0,1)\), is adopted from [3], [15] and does not admit an explicit analytical solution. The remaining problem data are specified as \[f=0, \quad y^d = \frac{1}{(x_1^2 + x_2^2)^{1/3}}, \quad \beta=(1,1)^T, \quad \alpha=1, \quad \omega=1.\] The admissible control set is defined by \[{U^{\mathrm{ad}}}= \bigl\{ u \in L^2(\Gamma) : 0 \leq u(x) \leq 0.2 \;\; \text{a.e. } x \in \Gamma \bigr\}.\] Here, the largest interior angle is \(\theta=\frac{\pi}{2}\) and \(y^d \in H^{1/3 - \eta}(\Omega)\) for any \(\eta >0\) due to the singularity on the boundary. Consequently, from Theorem 1, it follows that \(s \in [\frac{1}{2}, \frac{5}{6})\).

The problem is solved numerically on a fine mesh consisting of 524288 elements (that is, \(h=2^{-8}/\sqrt{2}\) and 1572864 degrees of freedom), and the resulting solution is used as a reference for comparison with solutions computed on coarser meshes. Fig. 3 displays the numerical solutions computed on the reference mesh for \(\epsilon=1\). The numerical results reported in Tab. 3 for \(\epsilon=1\) again show higher convergence rates for \(\|u-u_h\|_{0,\Gamma}\) and \(\|y-y_h\|_{0,\Omega}\). Let us note that similar observations are also reported for a related example in [5], which corresponds to our numerical setting. As for the previous example, we observe that the convergence rates are in good agreement with those proved for the classical finite element approach; compare e.g., [17]. However, from the numerical results reported in Tab. 4 for \(\epsilon=10^{-4}\) it is difficult to deduce an error behavior from the estimate ?? , since we expect that the exact state and adjoint develop boundary layers, so that besides the explicit appearance of \(\epsilon\) in this estimate one also has to take into account the growth of the respective \(H^2\)-seminorms of \(y\) and \(z\) in the estimation of \(\|y-y_h(u)\|_{0,\Omega}\) and \(\|z-z_h(u)\|_{0,\Gamma}\), which in the case of boundary layers grow critically with decreasing \(\epsilon\).

5.3 Example 3↩︎

Our last example, adapted from [17], is formulated on a polygonal domain with maximum interior angle \(\theta = \frac{5}{6} \pi\), as depicted in Fig. 4. The remaining problem data are \[y^d = \begin{cases} -1, & 0 \leq x_2 < 0.5, \\ 1, & 0.5 \leq x_2 < 1, \end{cases} \quad f = 1, \quad \beta = (1,0)^T, \quad \alpha = 2, \quad \omega = 1.\] The admissible control set is defined by \[{U^{\mathrm{ad}}}= \bigl\{ u \in L^2(\Gamma) : 0 \leq u(x) \;\; \text{a.e. } x \in \Gamma \bigr\}.\] In this example, \(y^d \in H^{1/2-\eta}(\Omega)\) for any \(\eta>0\), and the maximum interior angle is \(\theta = \frac{5\pi}{6}\). Consequently, it follows that \(s \in [\frac{1}{2},\frac{7}{10})\) from Theorem 1.

The reference solutions have been computed on a fine mesh with 458752 elements (that is, 1376256 degrees of freedom); see Fig. 5 for the case \(\epsilon=1\). In this case, we again observe better convergence rates than predicted by our theory. For the state we observe a convergence rate close to 1.1, which is the rate proven in [17] for the classical finite element approach. According to this reference, the control error for the classical approach should converge with order 0.6. Here, however, we observe an order between 1 and 1.2. Concerning the behaviour for small \(\epsilon\), a discussion similar to that of Example 5.2 applies. The numerical results are reported in Tab. 5 for \(\epsilon=1\) and Tab. 6 for \(\epsilon=10^{-4}\).

6 Conclusions↩︎

In this work, we have investigated Dirichlet boundary control of a convection–diffusion equation with \(L^2\)-boundary controls subject to pointwise constraints, discretized using the LDG method. The LDG framework naturally incorporates the Dirichlet boundary conditions into the variational formulation, even when the control space is chosen as \(L^2(\Gamma)\). We have derived general a priori error estimates for both the fully discrete and the variational discrete approximation of the optimal control problem in convex polygonal domains and presented numerical results that for moderate sizes of \(\epsilon\) underpin the theoretical predictions for the numerical approximation of Dirichlet boundary control problems. As future work, the development of a posteriori error estimates and adaptive discontinuous Galerkin methods could help to further robustify the numerical treatment in the convection-dominated case.

References↩︎