Error estimates for optimal control problems involving the Stokes system and Dirac measures

TL;DR

This paper derives a priori error estimates for finite element methods applied to optimal control problems involving the Stokes system with Dirac measures, addressing reduced regularity and providing numerical validation.

Contribution

It introduces new finite element error estimates for control problems with Dirac measures and reduced regularity in Stokes systems, supported by numerical experiments.

Findings

Error estimates for velocity and control variables

Numerical validation in 2D and 3D cases

Handling reduced regularity due to Dirac measures

Abstract

The aim of this work is to derive a priori error estimates for finite element discretizations of control--constrained optimal control problems that involve the Stokes system and Dirac measures. The first problem entails the minimization of a cost functional that involves point evaluations of the velocity field that solves the state equations. This leads to an adjoint problem with a linear combination of Dirac measures as a forcing term and whose solution exhibits reduced regularity properties. The second problem involves a control variable that corresponds to the amplitude of forces modeled as point sources. This leads to a solution of the state equations with reduced regularity properties. For each problem, we propose a finite element solution technique and derive a priori error estimates. Finally, we present numerical experiments, in two and three dimensions, that illustrate our theoretical developments.

Equations

$$J(\mathbf{y},\mathbf{u}):=\frac{1}{2}\sum _{t\in \mathcal{Z}}|\mathbf{y}(t)-{\mathbf{y}}_t{|}^2+\frac{\lambda }{2}\Vert \mathbf{u}{\Vert }_{{\mathbf{L}}^2(\Omega )}^2,$$

$$\left\{\begin{array}[]{rcll}-\Delta\mathbf{y}+\nabla p&=&\mathbf{u}&\text{ in }\quad\Omega,\\ \text{div }\mathbf{y}&=&0&\text{ in }\quad\Omega,\\ \mathbf{y}&=&\mathbf{0}&\text{ on }\quad\partial\Omega,\end{array}\right.$$

$$\mathbf{u}\in {\mathbb{U}}_{ad},{\mathbb{U}}_{ad}:=\{\mathbf{v}\in {\mathbf{L}}^2(\Omega ):\mathbf{a}\le \mathbf{v}\le \mathbf{b}\text{ a.e. in }\Omega \},$$

$$\left\{\begin{array}[]{rcll}-\Delta\mathbf{z}-\nabla r&=&\sum_{t\in\mathcal{Z}}(\mathbf{y}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(t)}-\mathbf{y}_{t})\delta_{t}&\text{ in }\Omega,\\ \text{div }\mathbf{z}&=&0&\text{ in }\Omega,\\ \mathbf{z}&=&\mathbf{0}&\text{ on }\partial\Omega.\end{array}\right.$$

$$\mathfrak{J}(\mathbf{y},\mathcal{U}):=\frac{1}{2}\Vert \mathbf{y}-{\mathbf{y}}_{\Omega }{\Vert }_{{\mathbf{L}}^2(\Omega )}^2+\frac{\lambda }{2}\sum _{t\in \mathcal{D}}|{\mathbf{u}}_t{|}^2,\mathcal{U}=({\mathbf{u}}_1,\dots ,{\mathbf{u}}_l),$$

$$-\Delta \mathbf{y}+\nabla p=\sum _{t\in \mathcal{D}}{\mathbf{u}}_t{\delta }_t\text{ in }\Omega ,\text{div }\mathbf{y}=0\text{ in }\Omega ,\mathbf{y}=0\text{ on }\partial \Omega ,$$

$${\mathfrak{U}}_{ad}:=\{\mathcal{V}=({\mathbf{v}}_1,\dots ,{\mathbf{v}}_l)\in [{\mathbb{R}}^d{]}^l:{\mathbf{a}}_t\le {\mathbf{v}}_t\le {\mathbf{b}}_t\text{ for all }t\in \mathcal{D}\},$$

$$(\mathbf{w},\mathbf{v}{)}_{{\mathbf{L}}^2(G)}={\int }_G\mathbf{w}\cdot \mathbf{v},\Vert \mathbf{v}{\Vert }_{{\mathbf{L}}^2(G)}=(\mathbf{v},\mathbf{v}{)}_{{\mathbf{L}}^2(G)}^{\frac{1}{2}}\forall \mathbf{w},\mathbf{v}\in {\mathbf{L}}^2(G).$$

$$C_{\omega }=\underset{B}{\sup }\left(\frac{1}{|B|}{\int }_B\omega \right)\left(\frac{1}{|B|}{\int }_B{\omega }^{-1}\right)<\infty ,$$

$$\Vert v{\Vert }_{H^1(\omega ,G)}:={\left(\Vert v{\Vert }_{L^2(\omega ,G)}^2+\Vert \nabla v{\Vert }_{L^2(\omega ,G)}^2\right)}^{\frac{1}{2}}.$$

$$\Vert \nabla \mathbf{v}{\Vert }_{{\mathbf{L}}^2(\omega ,G)}={\left(\sum _{i=1}^d\Vert \nabla v_i{\Vert }_{L^2(\omega ,G)}^2\right)}^{\frac{1}{2}}.$$

$$\left\{\begin{array}[]{rcll}-\Delta\bm{\Phi}+\nabla\zeta&=&\sum_{t\in\mathcal{E}}\mathbf{F}_{t}\delta_{t}&\text{ in }\quad\Omega,\\ \text{div}\>\bm{\Phi}&=&0&\text{ in }\quad\Omega,\\ \bm{\Phi}&=&\mathbf{0}&\text{ on }\quad\partial\Omega,\end{array}\right.$$

$$|\nabla \mathbf{\Phi }(x)|\approx |x-t{|}^{1-d},|\zeta (x)|\approx |x-t{|}^{1-d}.$$

$${\int }_{B(t^{\ast },r)}|x-t^{\ast }{|}^{\alpha }|\nabla \mathbf{\Phi }(x){|}^2\mathrm{d}x<\infty ,{\int }_{B(t^{\ast },r)}|x-t^{\ast }{|}^{\alpha }|\zeta (x){|}^2\mathrm{d}x<\infty ,$$

$$d_{\mathcal{E}}:=\left\{\begin{array}[]{ll}\mathrm{dist}(\mathcal{E},\partial\Omega),&\mbox{if }\#\mathcal{E}=1,\\ \min\left\{\mathrm{dist}(\mathcal{E},\partial\Omega),\min\{|t-t^{\prime}|:t,t^{\prime}\in\mathcal{E},\ t\neq t^{\prime}\}\right\},&\mbox{otherwise}.\end{array}\right.$$

$$\rho (x)={\mathsf{d}}_t^{\alpha }(x),$$

$$\rho (x)=\{\begin{array}{ll}{\mathsf{d}}_t^{\alpha }(x),& \exists t\in \mathcal{E}:{\mathsf{d}}_t(x)<\frac{d_{\mathcal{E}}}{2},\\ 1,& {\mathsf{d}}_t(x)\ge \frac{d_{\mathcal{E}}}{2}\forall t\in \mathcal{E},\end{array}$$

$$\mathbb{X}:$$

$$\mathbb{Y}:$$

$$a:\mathbb{X}\times \mathbb{X}\to \mathbb{R},a(\mathbf{w},\mathbf{v}):={\int }_{\Omega }\nabla \mathbf{w}:\nabla \mathbf{v},$$

$$b:\mathbb{X}\times \mathbb{Y},b(\mathbf{v},q):=-{\int }_{\Omega }q\text{div}\mathbf{v}.$$

$$\left\{\!\!\begin{array}[]{rcll}a(\bm{\Phi},\mathbf{v})+b(\mathbf{v},\zeta)&=&\sum_{t\in\mathcal{E}}\langle\mathbf{F}_{t}\delta_{t},\mathbf{v}\rangle&\quad\forall\mathbf{v}\in\mathbf{H}_{0}^{1}(\rho^{-1},\Omega),\\ b(\bm{\Phi},q)&=&0&\quad\forall q\in L^{2}(\rho^{-1},\Omega)/\mathbb{R},\end{array}\right.$$

$$\Vert \nabla \mathbf{\Phi }{\Vert }_{{\mathbf{L}}^2(\rho ,\Omega )}+\Vert \zeta {\Vert }_{L^2(\rho ,\Omega )/\mathbb{R}}\lesssim \sum _{t\in \mathcal{E}}|{\mathbf{F}}_t|\Vert {\delta }_t{\Vert }_{H_0^1({\rho }^{-1},\Omega {)}^{\prime }}.$$

$$\Vert \mathbf{v}{\Vert }_{{\mathbf{L}}^2(\Omega )}\lesssim \Vert \nabla \mathbf{v}{\Vert }_{{\mathbf{L}}^2(\rho ,\Omega )}\forall \mathbf{v}\in {\mathbf{H}}_0^1(\rho ,\Omega ),$$

$${\mathcal{N}}_T:=\left\{T^{\prime }\in \mathcal{T}:T\cap T^{\prime }\ne \varnothing \right\}.$$

$$\begin{array}{rl}{Q}_h& =\left\{q_h\in C(\bar{\Omega }):q_h{|}_T\in {\mathbb{P}}_1(T)\forall T\in {\mathcal{T}}_h\right\}\cap L^2(\Omega )/\mathbb{R},\\ {\mathbf{V}}_h& =\left\{{\mathbf{v}}_h\in \mathbf{C}(\bar{\Omega }):{\mathbf{v}}_h{|}_T\in [{\mathbb{P}}_1(T)\oplus \mathbb{B}(T){]}^d\forall T\in {\mathcal{T}}_h\right\}\cap {\mathbf{H}}_0^1(\Omega ),\end{array}$$

$$\begin{array}{rl}{Q}_h& =\left\{q_h\in C(\bar{\Omega }):q_h{|}_T\in {\mathbb{P}}_1(T)\forall T\in {\mathcal{T}}_h\right\}\cap L^2(\Omega )/\mathbb{R},\\ {\mathbf{V}}_h& =\left\{{\mathbf{v}}_h\in \mathbf{C}(\bar{\Omega }):{\mathbf{v}}_h{|}_T\in [{\mathbb{P}}_2(T){]}^d\forall T\in {\mathcal{T}}_h\right\}\cap {\mathbf{H}}_0^1(\Omega ).\end{array}$$

$$\Vert \mathbf{\Phi }-{\mathbf{\Phi }}_h{\Vert }_{{\mathbf{L}}^2(\Omega )}\lesssim h^{2-d/2}\sum _{t\in \mathcal{E}}|{\mathbf{F}}_t|\Vert {\delta }_t{\Vert }_{\mathcal{M}(\Omega )},$$

$$\min \{J(\mathbf{y},\mathbf{u}):(\mathbf{y},\mathbf{u})\in {\mathbf{H}}_0^1(\Omega )\times {\mathbb{U}}_{ad}\},$$

$$\left\{\begin{array}[]{rcll}a(\mathbf{y},\mathbf{v})+b(\mathbf{v},p)&=&(\mathbf{u},\mathbf{v})_{{\mathbf{L}}^{2}(\Omega)}&\quad\forall\mathbf{v}\in\mathbf{H}_{0}^{1}(\Omega),\\ b(\mathbf{y},q)&=&0&\quad\forall q\in L^{2}(\Omega)/\mathbb{R}.\end{array}\right.$$

$$\Vert \nabla \mathbf{y}{\Vert }_{{\mathbf{L}}^2(\Omega )}+\Vert p{\Vert }_{L^2(\Omega )}\lesssim \Vert \mathbf{u}{\Vert }_{{\mathbf{L}}^2(\Omega )};$$

Full text

{AMS}

35Q35, 35R06, 49M25, 65N15, 65N30.

1 Introduction

Let $d\in\{2,3\}$ and $\Omega\subset\mathbb{R}^{d}$ be an open and bounded polytopal domain with Lipschitz boundary $\partial\Omega$. The purpose of this work is the study of a priori error estimates for finite element discretizations that approximate the solution to optimal control problems involving the Stokes equations and Dirac measures; control–constraints are also considered. We consider two illustrative examples, which we proceed to describe in what follows.

1.1 Optimization with Point Observations

Let $\mathcal{Z}\neq\emptyset$ be a finite ordered subset of $\Omega$ with cardinality $\#\mathcal{Z}=m$. Given a set of desired states $\{{\bf{y}}_{t}\}_{t\in{\mathcal{Z}}}\subset\mathbb{R}^{d}$, a regularization parameter $\lambda>0$, and the cost functional

[TABLE]

we are interested in finding $\min J(\mathbf{y},\mathbf{u})$ subject to the Stokes system

[TABLE]

and the control constraints

[TABLE]

with $\mathbf{a},\mathbf{b}\in\mathbb{R}^{d}$ satisfying $-\bm{\infty}<\mathbf{a}<\mathbf{b}<\bm{\infty}$. We immediately comment that, throughout this work, vector inequalities must be understood componentwise and that $|\cdot|$ will denote the Euclidean norm in $\mathbb{R}^{d}$.

In contrast to standard PDE–constrained optimization problems, the cost functional (1) involves point evaluations of the state velocity field. This leads to a subtle formulation of the adjoint problem:

[TABLE]

Here, $\delta_{t}$ corresponds to the Dirac delta supported at the interior point $t$.

The optimal control problem (1)–(3) finds relevance in some applications where the observations are carried out at specific locations of the domain. For instance, in the active control of sound [27, 10] and in the active control of vibrations [31, 41]; see also [13, 15] for other applications. We immediately comment that, since the domain $\Omega$ is Lipschitz and $\mathbb{U}_{ad}\subset\mathbf{L}^{\infty}(\Omega)$, the point observations $\mathbf{y}(t)$, with $t\in\mathcal{Z}$, are well–defined; see [14, Theorem 2.9] and [21, Lemma 12]. These point observations tend to enforce the state velocity field $\mathbf{y}$ to have the fixed vector value $\mathbf{y}_{t}$ at the point $t\in\mathcal{Z}$. Consequently, problem (1)–(3) can be understood as a penalty version of a PDE–constrained optimization problem where the velocity field that solves the state equation is constrained at a collection of points.

There are several works that provide a priori error estimates for the optimal control problem (1)–(3) but when the state equations are a Poisson problem. Under the fact that the associated adjoint variable belongs to $W_{0}^{1,r}(\Omega)$, for $r\in(1,d/(d-1))$, the authors of [15] obtain, for $d\in\{2,3\}$, a priori and a posteriori error estimates for the so-called variational discretization scheme when applied for discretizing the underlying optimal control problem; the state and adjoint equations are discretized on the basis of standard piecewise linear finite element functions. The derived a priori error estimates for the control, the state, and adjoint state variables are optimal in terms of regularity; see [15, Theorem 3.2]. Later, the authors of [13] propose and analyze a fully discrete scheme that discretizes the optimal state, adjoint, and control variables with piecewise linear functions. For this scheme, the authors provide error estimates when $d=2$ [13, Theorem 5.1]; the control and the state are discretized using meshes of size $\mathcal{O}(h^{2})$ and $\mathcal{O}(h)$, respectively. The authors of [13] also analyze the variational discretization scheme and derive an optimal error estimate, in terms of regularity, for the control variable [13, Theorem 5.2]. In [7], the authors invoke the theory of Muckenhoupt weights and Muckenhoupt–weighted Sobolev spaces to provide error estimates for a numerical scheme that discretizes the control variable with piecewise constant functions; the state and adjoint equations are discretized with piecewise linear finite elements. To be precise, the authors derive a priori error estimates for the error approximation of the optimal control variable when $d\in\{2,3\}$; the one for $d=2$ being nearly–optimal in terms of approximation [7, Theorem 4.3]. However, the estimate for $d=3$ is suboptimal in terms of approximation; it behaves as $\mathcal{O}(h^{1/2}|\log h|)$. This has been recently improved in [9, Theorem 6.6].

1.2 Optimization with Singular Sources

Let $\mathcal{D}\neq\emptyset$ be a finite ordered subset of $\Omega$ with cardinality $\#\mathcal{D}=l$. Given a desired state ${\bf{y}}_{\Omega}\in\mathbf{L}^{2}(\Omega)$, a regularization parameter $\lambda>0$, and the cost functional

[TABLE]

the problem under consideration reads as follows: Find $\min\mathfrak{J}(\mathbf{y},\mathcal{U})$ subject to

[TABLE]

and the control constraints $\mathcal{U}=(\mathbf{u}_{1},\ldots,\mathbf{u}_{l})\in\mathfrak{U}_{ad}$, where

[TABLE]

with $\mathbf{a}_{t},\mathbf{b}_{t}\in\mathbb{R}^{d}$ satisfying $\mathbf{a}_{t}<\mathbf{b}_{t}$ for all $t\in\mathcal{D}$.

The optimization problem (5)–(7) is of relevance in applications where it can be specified a control at finitely many prespecified points. We observe that, in view of the particular structure of the control variable $\mathcal{U}$, the state equation (6) corresponds to a Stokes system that has a linear combination of Dirac measures on the right–hand side of the momentum equation.

There are a few works that consider the numerical approximation of problem (5)–(7) when the Stokes equations are replaced by a Poisson problem. In [35], the authors use the variational discretization concept to derive a priori error estimates. Their technique is based on the fact that the state belongs to $W_{0}^{1,r}(\Omega)$ for $r\in(1,d/(d-1))$. An approach involving weighted estimates has also been considered in [7], where the authors obtain the following rates of convergence for the error approximation of the control variable: $\mathcal{O}(h^{2-\epsilon})$ in two dimensions and $\mathcal{O}(h^{1-\epsilon})$ in three dimensions, where $\epsilon>0$.

Since the Stokes system with a linear combination of Dirac measures in the momentum equation appears as the adjoint system (4) of problem (1)–(3) and as the state equation of problem (5)–(7), it is of importance to understand the regularity properties of the involved solution and the development of numerical techniques to approximate it. The main difficulty in the study of the aforementioned problem is that it does not admit a solution in the classical Hilbert space $\mathbf{H}^{1}_{0}(\Omega)\times L^{2}(\Omega)/\mathbb{R}$. In spite of this fact, there are a few works that consider the numerical approximation of this system. In [49] the author presents a numerical method in two and three–dimensional bounded domains. However, convergence properties are not investigated. Later, the authors of [11] derive quasi–optimal local convergence results in $\mathbf{H}^{1}\times L^{2}$; the error is analyzed on a subdomain which does not contain the singularity of the solution. The authors operate in two dimensions and consider the mini element and Taylor–Hood schemes. On the other hand, on the basis of the fact that the solution to the Stokes system with singular sources can be seen as an element of a weighted space, a priori and a posteriori error estimates, for classical low–order inf–sup stable finite element approximations, have been developed in [26] and [4], respectively. The a priori error estimates derived in [26, Corollary 5.4] have been recently improved in [57].

In spite of these advances and to the best of our knowledge, this is the first work that provides a priori error estimates for the optimal control problems (1)–(3) and (5)–(7). We discretize the adjoint and state equations with classical low–order inf–sup stable finite element schemes and the control variable with piecewise constant functions. We derive error estimates for problem (1)–(3), in two and three–dimensions, and problem (5)–(7), in two–dimensions. In particular, in two–dimensions, and for both problems, the derived error estimates for the discretization of the optimal control variable are nearly–optimal in terms of approximation. In three–dimensions, and for problem (1)–(3), the obtained error estimate for the control variable is suboptimal in terms of approximation. However, numerical experiments seem to indicate that such an estimate is attained in practice.

The rest of the paper is organized as follows. In Section 2 we introduce the notation and functional framework we shall work with. We also briefly review, in Section 2.3, the well–posedness of the Stokes system with singular sources. Section 3 contains the numerical analysis for problem (1)–(3). In Section 4 we derive error estimates for the optimal control problem (5)–(7). We conclude, in Section 5, with a series of numerical examples that illustrate the developed theory.

2 Notation and Preliminaries

Let us fix the notation and conventions in which we will operate.

2.1 Notation

Throughout this work $d\in\{2,3\}$ and $\Omega\subset\mathbb{R}^{d}$ is an open, bounded, and convex polytopal domain. If $\mathcal{X}$ and $\mathcal{Y}$ are normed vector spaces, we write $\mathcal{X}\hookrightarrow\mathcal{Y}$ to denote that $\mathcal{X}$ is continuously embedded in $\mathcal{Y}$. We denote by $\mathcal{X}^{\prime}$ and $\|\cdot\|_{\mathcal{X}}$ the dual and the norm of $\mathcal{X}$, respectively. Given a Lebesgue measurable subset $A$ of $\mathbb{R}^{d}$, we denote by $|A|$ its Lebesgue measure.

To denote vector-valued functions we shall use lower–case bold letters, whereas to denote function spaces we shall use upper-case bold letters. For a bounded domain $G\subset\mathbb{R}^{d}$, if $X(G)$ corresponds to a function space over $G$, we shall denote $\mathbf{X}(G)=[X(G)]^{d}$. In particular, we denote $\mathbf{L}^{2}(G)=[L^{2}(G)]^{d}$, which is equipped with the following inner product and norm, respectively:

[TABLE]

Finally, the relation $a\lesssim b$ indicates that $a\leq Cb$, with a positive constant that does not depend on $a$, $b$ nor the discretization parameter. The value of $C$ might change at each occurrence.

2.2 Weighted Sobolev Spaces

We start this section with a notion which will be fundamental for further discussions, that of a weight. A weight is a nonnegative locally integrable function on $\mathbb{R}^{d}$ that takes values in $(0,\infty)$ almost everywhere. We will be particularly interested in the weights belonging to the so-called Muckenhoupt class $A_{2}(\mathbb{R}^{d})$ [23, 29, 55, 61].

Definition 2.1** (Muckenhoupt class $A_{2}(\mathbb{R}^{d})$).**

Let $\omega$ be a weight. We say that $\omega\in A_{2}(\mathbb{R}^{d})$, or that $\omega$ is an $A_{2}(\mathbb{R}^{d})$–weight, if there exists a positive constant $C_{\omega}$ such that

[TABLE]

*where the supremum is taken over all balls $B$ in $\mathbb{R}^{d}$. *

Let $\omega\in A_{2}(\mathbb{R}^{d})$ and $G\subset\mathbb{R}^{d}$ be an open and bounded domain. We define the weighted Lebesgue space L^{2}(\omega,{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}G}) as the space of square–integrable functions with respect to the measure $\omega\textrm{d}x$. We also define the weighted Sobolev space $H^{1}(\omega,G):=\{v\in L^{2}(\omega,G):|\nabla v|\in L^{2}(\omega,G)\},$ which we equip with the norm

[TABLE]

It is remarkable that most of the properties of classical Sobolev spaces have a weighted counterpart. This is not because of the specific form of the weight but rather due to the fact that the weight $\omega\in A_{2}(\mathbb{R}^{d})$. In particular, $L^{2}(\omega,G)$ and $H^{1}(\omega,G)$ are Hilbert spaces and C^{\infty}({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}G}) is dense in $H^{1}(\omega,G)$; see, for instance, [61, Proposition 2.1.2, Corollary 2.1.6] and [34, Theorem 1]. Define $H^{1}_{0}(\omega,G)$ as the closure of C_{0}^{\infty}({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}G}) in $H^{1}(\omega,G)$. In view of a weighted Poincaré inequality, that follows from [29, Theorem 1.3] for {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}G} being a ball and [17, 44] for more general domains, we conclude that, in $H^{1}_{0}(\omega,G)$, the seminorm $\|\nabla v\|_{L^{2}(\omega,G)}$ is equivalent to the norm $\|v\|_{H^{1}(\omega,G)}$.

Finally, we define the vector space $\mathbf{H}^{1}_{0}(\omega,G):=[H^{1}_{0}(\omega,G)]^{d}$, which, in view of the aforementioned Poincaré inequality, we equip with the norm

[TABLE]

2.3 The Stokes System with Dirac Sources

In this section we review well–posedness results in weighted spaces for the Stokes system with a linear combination of Dirac measures as a forcing term in the momentum equation. We also comment on the finite element approximation of such a problem. The review is motivated since the adjoint equation (4) for the pointwise tracking optimal control problem of Section 1.1 and the state equation (6) for the optimization with singular sources problem of Section 1.2 are particular instances of the aforementioned singular setting.

Let $\mathcal{E}$ be a finite ordered subset of $\Omega$ with cardinality $\#\mathcal{E}$. Given $\{\mathbf{F}_{t}\}_{t\in\mathcal{E}}\subset\mathbb{R}^{d}$, we consider the following boundary value problem: Find $(\bm{\Phi},\zeta)$ such that

[TABLE]

where $\delta_{t}$ denotes the Dirac delta supported at ${t}\in\Omega$. Let us assume that $\Omega=\mathbb{R}^{d}$. If this is the case, the results of [32, Section IV.2] yield the following asymptotic behavior near the points $t\in\mathcal{E}$:

[TABLE]

Consequently, $(\bm{\Phi},\zeta)\not\in\mathbf{H}^{1}(\Omega)\times L^{2}(\Omega)/\mathbb{R}$. However, if $B(t^{*},r)$ denotes a ball of radius $r$ and center $t^{*}$ with $t^{*}\in\mathcal{E}$ and $r>0$ such that $B(t^{*},r)\cap\mathcal{E}=\{t^{*}\}$, then

[TABLE]

for $\alpha\in(d-2,\infty)$. This heuristic argument suggests to study problem (8) on weighted spaces.

2.3.1 Weak Formulation

Define

[TABLE]

Since $\mathcal{E}\subset\Omega$ and $\mathcal{E}$ is finite, we thus have that $d_{\mathcal{E}}>0$. With this notation, we define the weight $\rho$ as follows: if $\#\mathcal{E}=1$, then

[TABLE]

otherwise

[TABLE]

where $\mathsf{d}_{t}^{\alpha}(x):=|x-t|^{\alpha}$ and $\alpha\in(-d,d)$. Since $\alpha\in(-d,d)$, owing to [1, Theorem 6] and [30, Lemma 2.3 (v)], we have that the function $\rho$ belongs to the Muckenhoupt class $A_{2}(\mathbb{R}^{d})$. Define the spaces

[TABLE]

Finally, we define the bilinear forms

[TABLE]

and

[TABLE]

With all these ingredients at hand, we present the following weak formulation of problem (8) [4, Section 3]: Find $(\bm{\Phi},\zeta)\in\mathbf{H}_{0}^{1}(\rho,\Omega)\times L^{2}(\rho,\Omega)/\mathbb{R}$ such that

[TABLE]

where $\langle\cdot,\cdot\rangle$ denotes the duality pairing between $\mathbf{H}_{0}^{1}(\rho^{-1},\Omega)^{\prime}$ and $\mathbf{H}_{0}^{1}(\rho^{-1},\Omega)$. We immediately mention that in order to guarantee that $\mathbf{F}_{t}\delta_{t}\in\mathbf{H}_{0}^{1}(\rho^{-1},\Omega)^{\prime}$, and thus that $\langle\mathbf{F}_{t}\delta_{t},\mathbf{v}\rangle$ is well defined for $\mathbf{v}\in\mathbf{H}_{0}^{1}(\rho^{-1},\Omega)$, the parameter $\alpha$ should be restricted to $(d-2,d)$; see [47, Lemma 7.1.3], [40, Remark 21.19], and [26, Proposition 5.2] for details.

It can be proved that problem (11) admits a unique solution; see [56, Theorem 14]. Moreover, the following estimate can be obtained [56, Theorem 14]:

[TABLE]

We conclude this section with the following embedding result.

Theorem 2.2** ($\mathbf{H}_{0}^{1}(\rho,\Omega)\hookrightarrow\mathbf{L}^{2}(\Omega)$).**

If $\alpha\in(d-2,2)$, then $\mathbf{H}_{0}^{1}(\rho,\Omega)\hookrightarrow\mathbf{L}^{2}(\Omega)$. Moreover, the following weighted Poincaré inequality holds

[TABLE]

*where the hidden constant depends only on $\Omega$ and $d_{\mathcal{E}}$. *

Proof. The proof follows from [3, Lemmas 1 and 2].

2.3.2 Finite Element Approximation and Error Estimates

We start the discussion by introducing some standard finite element notation [12, 18, 28]. We denote by $\mathscr{T}_{h}=\{T\}$ a conforming partition, or mesh, of $\bar{\Omega}$ into closed simplices $T$ with size $h_{T}=\text{diam}(T)$. Define $h:=\max_{T\in\mathscr{T}_{h}}h_{T}$. We denote by $\mathbb{T}=\{\mathscr{T}_{h}\}_{h>0}$ a collection of conforming and quasi–uniform meshes. For $T\in\mathscr{T}_{h}$, we define the patch $\mathcal{N}_{T}$ associated with the element $T$ as

[TABLE]

In an abuse of notation, in what follows, by $\mathcal{N}_{T}$ we will indistinctively denote either this set or the union of the triangles that comprise it.

Given a mesh $\mathscr{T}_{h}\in\mathbb{T}$, we denote by $\mathbf{V}_{h}$ and $Q_{h}$ the finite element spaces that approximate the velocity field and the pressure, respectively, constructed over $\mathscr{T}_{h}$. In this work we will consider the following popular finite element discretizations:

• (a)

The mini element [28, Section 4.2.4]:

[TABLE]

where $\mathbb{B}(T)$ denotes the space spanned by local bubble functions.

• (b)

The classical Taylor–Hood elements [28, Section 4.2.5]:

[TABLE]

We now present an error estimate.

Lemma 2.3** (Error estimate for Stokes system with Dirac sources).**

Let $\Omega$ be convex. Let $(\bm{\Phi},\zeta)\in\mathbf{H}_{0}^{1}(\rho,\Omega)\times L^{2}(\rho,\Omega)/\mathbb{R}$ be the solution of (11). Let $(\bm{\Phi}_{h},\zeta_{h})\in\mathbf{V}_{h}\times Q_{h}$ be the finite element approximation of $(\bm{\Phi},\zeta)$ on the basis of the discrete spaces (14) or (15). Then, we have the error estimate

[TABLE]

*with a hidden constant that does not depend on $(\bm{\Phi},\zeta)$, $(\bm{\Phi}_{h},\zeta_{h})$, nor $h$. Here, $\mathcal{M}(\Omega)$ denotes the space of bounded Radon measures. *

Proof. See [57, Theorem 7].

3 The Pointwise Tracking Optimal Control Problem

In this section we analyze a weak version of the optimal control problem (1)–(3), which reads as follows: Find

[TABLE]

subject to

[TABLE]

Standard arguments that rely on the coercivity of $a$ on $\mathbf{H}_{0}^{1}(\Omega)$ and an inf–sup condition for $b$ yield the existence of a unique solution $(\mathbf{y},p)\in\mathbf{H}_{0}^{1}(\Omega)\times L^{2}(\Omega)/\mathbb{R}$ to problem (18). In addition, the pair $(\mathbf{y},p)$ satisfies the following stability estimate:

[TABLE]

see, for instance, [28, Theorem 4.3]. We also present the following regularity result for the pair $(\mathbf{y},p)$ that solves (18).

Proposition 3.1** (regularity).**

Let $\Omega\subset\mathbb{R}^{d}$ be a convex polytope and $\mathbf{u}\in\mathbf{L}^{2}(\Omega)$. If $(\mathbf{y},p)$ solves (18), then $\mathbf{y}\in\mathbf{H}^{2}(\Omega)\cap\mathbf{H}_{0}^{1}(\Omega)$, $p\in H^{1}(\Omega)\cap L^{2}(\Omega)/\mathbb{R}$, and

[TABLE]

*with a hidden constant that is independent of $(\mathbf{y},p)$ and $\mathbf{u}$.

Proof. See [45, Theorem 2] and [37] for $d=2$, and [19] and [48] for $d=3$; see also [54, Corollary 1.8].

Note that, since the control variable $\mathbf{u}\in\mathbb{U}_{ad}$ and $\Omega$ is convex, the results of Proposition 3.1 guarantee that $\mathbf{y}\in\mathbf{H}^{2}(\Omega)\hookrightarrow\mathbf{C}(\Omega)$. As a consequence, point evaluations of the velocity field $\mathbf{y}$ on the cost functional $J$ are well defined.

3.1 Optimality Conditions

Since $\lambda>0$ and the underlying control-to-state operator $\mathbf{S}$ is linear and continuous, standard arguments yield the existence of a unique solution $(\bar{\mathbf{y}},\bar{\mathbf{u}})\in\mathbf{H}_{0}^{1}(\Omega)\times\mathbb{U}_{ad}$ to the optimal control problem (17)–(18). To present optimality conditions, we introduce the adjoint pair $(\mathbf{z},r)$ as the unique solution to the following problem: Find $(\mathbf{z},r)\in\mathbf{H}_{0}^{1}(\rho,\Omega)\times L^{2}(\rho,\Omega)/\mathbb{R}$ such that

[TABLE]

where $\mathbf{y}$ solves (18). The weight $\rho$ is defined as in (9)–(10) with obvious modifications. The following first–order sufficient and necessary optimality condition follows from [2, Theorem 8]: $(\bar{\mathbf{y}},\bar{\mathbf{u}})$ is optimal for (17)–(18) if and only if $\bar{\mathbf{u}}$ satisfies

[TABLE]

Here, $(\bar{\mathbf{z}},\bar{r})\in\mathbf{H}^{1}_{0}(\rho,\Omega)\times L^{2}(\rho,\Omega)/\mathbb{R}$ denotes the optimal adjoint state, which solves (21) with $\mathbf{y}$ replaced by $\bar{\mathbf{y}}$. The well–posedness of problem (21) follows from the results elaborated in Section 2.3.

We now recall the so-called projection formula for $\bar{\mathbf{u}}$: The optimal control $\bar{\mathbf{u}}$ satisfies (22) if and only if

[TABLE]

where $\Pi_{[\mathbf{a},\mathbf{b}]}:\mathbf{L}^{1}(\Omega)\rightarrow\mathbb{U}_{ad}$ is such that $\Pi_{[\mathbf{a},\mathbf{b}]}(\mathbf{v}):=\min\{\mathbf{b},\max\{\mathbf{v},\mathbf{a}\}\}.$ This projection formula leads to the following regularity result for $\bar{\mathbf{u}}$.

Proposition 3.2**.**

*(Regularity of $\bar{\mathbf{u}}$)

If $\bar{\mathbf{u}}$ is optimal for problem (17)–(18), then $\bar{\mathbf{u}}\in\mathbf{H}^{1}(\rho,\Omega)$. Moreover, the following estimate holds:*

[TABLE]

Proof. Note that, in view of (23), $\bar{\mathbf{u}}$ can be written as

[TABLE]

The regularity result thus follows directly from [46, Theorem A.1]. This concludes the proof.

To summarize, the pair $(\bar{\mathbf{y}},\bar{\mathbf{u}})$ is optimal for the pointwise tracking optimal control problem (17)–(18) if and only if $(\bar{\mathbf{y}},\bar{p},\bar{\mathbf{z}},\bar{r},\bar{\mathbf{u}})\in\mathbf{H}_{0}^{1}(\Omega)\times L^{2}(\Omega)/\mathbb{R}\times\mathbf{H}_{0}^{1}(\rho,\Omega)\times L^{2}(\rho,\Omega)/\mathbb{R}\times\mathbb{U}_{ad}$ solves the optimality system (18), (21), and (22).

3.2 A Fully Discrete Scheme: Error Estimates

In this section, we design and analyze a fully discrete scheme for the numerical approximation of problem (17)–(18). To begin with, we define the discrete admissible set

[TABLE]

where \mathbf{U}_{h}:=\{\mathbf{u}\in\mathbf{L}^{2}({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Omega}):\mathbf{u}|_{T}\in[\mathbb{P}_{0}(T)]^{d}\ \forall T\in\mathscr{T}_{h}\}.

The discrete counterpart of the continuous problem (17)–(18) thus reads as follows: Find min $J(\mathbf{y}_{h},\mathbf{u}_{h})$ subject to the discrete state equations

[TABLE]

and the discrete control constraints $\mathbf{u}_{h}\in\mathbb{U}_{ad}^{h}$. Standard arguments reveal the existence of a unique optimal pair ({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(\bar{\mathbf{y}}_{h},\bar{p}_{h})},\bar{\mathbf{u}}_{h}). In addition, the pair ({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(\bar{\mathbf{y}}_{h},\bar{p}_{h})},\bar{\mathbf{u}}_{h}) is optimal for the aforementioned discrete optimal control problem if and only if {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(\bar{\mathbf{y}}_{h},\bar{p}_{h})} solves (24), with $\mathbf{u}_{h}$ replaced by $\bar{\mathbf{u}}_{h}$, and $\bar{\mathbf{u}}_{h}$ satisfies the variational inequality

[TABLE]

where $(\bar{\mathbf{z}}_{h},\bar{r}_{h})$ solves

[TABLE]

with ${\mathbf{y}}_{h}$ replaced by $\bar{\mathbf{y}}_{h}$.

3.2.1 Auxiliary Problems

We introduce two auxiliary problems that will be instrumental to derive error estimates for the proposed discrete scheme.

The first problem reads as follows: Find $(\hat{\mathbf{y}}_{h},\hat{p}_{h})\in\mathbf{V}_{h}\times Q_{h}$ such that

[TABLE]

The second auxiliary problem is: Find $(\hat{\mathbf{z}}_{h},\hat{r}_{h})\in\mathbf{V}_{h}\times Q_{h}$ such that

[TABLE]

Before providing error estimates, we present the following auxiliary result.

Lemma 3.3** (Discrete pointwise stability).**

Let $(\bm{\xi},\theta)\in\mathbf{H}_{0}^{1}(\Omega)\times L^{2}(\Omega)/\mathbb{R}$ be the solution to

[TABLE]

and let $(\bm{\xi}_{h},\theta_{h})\in\mathbf{V}_{h}\times Q_{h}$ be its Galerkin approximation on the basis of the discrete spaces (14) or (15). Then,

[TABLE]

*where the hidden constant is independent of $(\bm{\xi},\theta)$, $(\bm{\xi}_{h},\theta_{h})$, $\bar{\mathbf{u}}$, $\bar{\mathbf{u}}_{h}$, and $h$. *

Proof. We begin the proof by noticing that, since $\bar{\mathbf{u}}-\bar{\mathbf{u}}_{h}\in\mathbf{L}^{2}(\Omega)$, then $(\bm{\xi},\theta)\in\mathbf{H}^{2}(\Omega)\cap\mathbf{H}_{0}^{1}(\Omega)\times H^{1}(\Omega)\cap L^{2}(\Omega)/\mathbb{R}$. This results follows from Proposition 3.1. Let us denote by $\mathbf{I}_{h}:\mathbf{C}(\bar{\Omega})\rightarrow\mathbf{V}_{h}$ the Lagrange interpolation operator. An application of the triangle inequality in conjunction with a standard inverse estimate yield

[TABLE]

To control the first term on the right hand side of (30) we invoke the continuous Sobolev embedding $\mathbf{H}^{2}\hookrightarrow\mathbf{C}(\bar{\Omega})$ and the regularity estimate (20) to arrive at

[TABLE]

For the remaining terms in the right hand side of (30) we utilize standard interpolation error estimates for $\mathbf{I}_{h}$ and error estimates for the finite element approximation of problem (29). This concludes the proof.

3.2.2 Error Estimates

To perform an a priori error analysis, it is useful to introduce the $\mathbf{L}^{2}(\Omega)$–orthogonal projection onto $[\mathbb{P}_{0}(\mathscr{T}_{h})]^{d}$, which is defined by

[TABLE]

where $T\in\mathscr{T}_{h}$. Note that, in view of the weighted Poincaré inequality of Theorem 2.2, $\Pi_{\mathbf{L}^{2}}$ is well defined over $\mathbf{H}_{0}^{1}(\rho,\Omega)$.

To obtain error estimates for $\Pi_{\mathbf{L}^{2}}$, we proceed with the help of regularity estimates for a suitable Green’s function. For $x\in\Omega$, we define the Green’s function $\mathcal{G}:\Omega\times\Omega\rightarrow\mathbb{R}$ as the solution (in the sense of distributions) to

[TABLE]

The following properties of the Green’s function are essential in what follows.

Proposition 3.4** (Properties of $\mathcal{G}$).**

Let $\Omega\subset\mathbb{R}^{d}$ be a Lipschitz polytope. The Green’s function $\mathcal{G}$ satisfies the following properties:

$\nabla\mathcal{G}(x,\cdot)\in L^{\frac{d}{d-1},\infty}(\Omega)$. 2. 2.

$\mathcal{G}(x,\cdot)\in W^{1,\frac{d}{d-1}}(\Omega\setminus B_{R})$* and $\mathcal{G}(x,\cdot)\in W^{2,1}(\Omega\setminus B_{R})$, where $B_{R}$ denotes the ball of radius $R$ and center $x\in\Omega$. In addition, we have*

[TABLE]

*Proof. A proof of $\nabla\mathcal{G}\in L^{d/(d-1),\infty}(\Omega)$ can be found in [38, Theorem 1.1] and [43, Theorem 4.1], for $d=3$, and [53, Theorem 1.1] and [22, Section 6], for $d=2$. Here, $L^{s,\infty}(\Omega)$, with $0<s<\infty$, corresponds to the so–called weak Lebesgue space; see definition 1.15 in [36]. A proof of the estimate (31) can be found in [20, Theorem 1]. *

We now provide an error estimate for the $\Pi_{L^{2}}$ projection.

Theorem 3.5**.**

Let $\Omega$ be an open, bounded and Lipschitz polytope. Let $x_{0}\in\Omega$ and let $u$ be the weak solution to the Poisson problem $-\Delta u=\delta_{x_{0}}$ in $\Omega$ and $u=0$ on $\partial\Omega$. Thus, we have the error estimate

[TABLE]

*with a hidden constant that is independent of $h$.

Proof. Let us consider a partition of $\mathscr{T}_{h}$ into the sets $\mathcal{N}_{x_{0}}:=\{T\in\mathscr{T}_{h}:x_{0}\in\mathcal{N}_{T}\}$ and $\mathcal{M}_{x_{0}}:=\mathscr{T}_{h}\setminus\mathcal{N}_{x_{0}}$ and utilize the regularity results of Proposition 3.4 accordingly; we recall that $\mathcal{N}_{T}$ is defined in (13). We proceed in two steps.

Step 1. Let $T\in\mathcal{M}_{x_{0}}$. Standard approximation results yields the estimate

[TABLE]

Since the mesh $\mathscr{T}_{h}$ is quasi–uniform, we thus obtain

[TABLE]

where, to obtain the last estimate, we have used the finite overlapping property of stars and (31) to conclude that $\|\nabla u\|_{L^{d/(d-1)}(\Omega\setminus B_{R})}\lesssim|\log h|$, with $B_{R}$, in this case, being the largest ball such that $B_{R}\subset\cup\{T\in\mathscr{T}_{h}:x_{0}\in T\}$.

Step 2. Let $T\in\mathcal{N}_{x_{0}}$. Observe that

[TABLE]

Let us use now the fact that $\nabla u\in L^{d/(d-1),\infty}(\Omega)$. Since $\nabla u\in L^{s,\infty}(\Omega)$ for $s\leq d/(d-1)$, we have that the distribution function of $|\nabla u|$, i.e., $d_{\nabla u}(\gamma):=|\{x\in\Omega:|\nabla u(x)|>\gamma\}|$, satisfies $d_{\nabla u}(\gamma)\lesssim\gamma^{-s}$ for every $\gamma>0$; see [36, (1.1.8) and (1.1.9)]. Here, $|A|$ denotes the Lebesgue measure of a Lebesgue measurable subset $A$ of $\mathbb{R}^{d}$. We recall that

[TABLE]

Let $C_{T}:=\|\nabla u\|_{L^{s,\infty}(T)}$. Invoke [36, Proposition 1.1.4] to obtain, for $s\in(1,d/(d-1)]$, that

[TABLE]

Set $\sigma=C_{T}[|T|(s-1)]^{-1/s}$ to arrive at $\int_{T}|\nabla u|\mathrm{d}x\lesssim C_{T}h_{T}^{d(1-1/s)}$. Set now $s=d/(d-1)$. This yields $\int_{T}|\nabla u|\mathrm{d}x\lesssim C_{T}h_{T}$. Consequently,

[TABLE]

with a hidden constant independent of $u$ and $h$. The previous estimate implies

[TABLE]

where we have used that, since $\{\mathscr{T}_{h}\}$ is shape–regular, $\#\mathcal{N}_{x_{0}}$ is uniformly bounded. This also yields the bound $\max_{T\in\mathcal{N}_{x_{0}}}\{C_{T}^{2}\}\lesssim 1$, with a hidden constant independent of $h$.

A collection of the bounds (33) and (34) yields the desired estimate (32). This concludes the proof.

Remark 3.6** (Green’s function for Stokes problem and error estimate for $\Pi_{\mathbf{L}^{2}}$).**

The Dirichlet Green’s matrix associated with the Stokes system can be defined as follows: $\mathcal{G}=(\mathcal{Q},\mathcal{R})$, where $\mathcal{Q}$ is a $d\times d$ matrix–valued function and $\mathcal{R}$ is a $1\times d$ vector–valued function that satisfy system (7.8) in [53]. Theorem 7.1 in [53] yields, for $d\in\{2,3\}$,

[TABLE]

see also [58, Theorem 1.1] and [16, Theorem 2.4]. Estimate (31) does not appear directly in the literature to our knowledge. In the case of Poisson’s problem, we stress that (31) can be found in [20, Lemma 2]. We shall thus assume estimate (31) in the case of the Stokes system. With this assumption and the regularity estimate (35) at hand, the arguments elaborated in Theorem 3.5 yield the error estimate

[TABLE]

where the hidden constant depends on the $\mathbf{L}^{d/(d-1),\infty}(\Omega)$-norm of $\bar{\mathbf{z}}$ but, more importantly, is independent of $h$. In view of the projection formula (23), [46, Theorem A.1], and the definition of the weak Lebesgue space $\mathbf{L}^{s,\infty}(\Omega)$, with $0<s<\infty$, we can conclude that the regularity estimates needed to obtain (36) are inherited to the control variable $\bar{\mathbf{u}}$. Consequently,

[TABLE]

*with, again, a hidden constant independent of $h$. *

We are ready to prove the main result of this section.

Theorem 3.7** (Rates of convergence for $\bar{\mathbf{u}}$).**

Let $(\bar{\mathbf{y}},\bar{p},\bar{\mathbf{z}},\bar{r},\bar{\mathbf{u}})\in\mathbf{H}_{0}^{1}(\Omega)\times L^{2}(\Omega)/\mathbb{R}\\ \times\mathbf{H}_{0}^{1}(\rho,\Omega)\times L^{2}(\rho,\Omega)/\mathbb{R}\times\mathbb{U}_{ad}$ be the solution to the optimality system (18), (21), and (22) and $(\bar{\mathbf{y}}_{h},\bar{p}_{h},\bar{\mathbf{z}}_{h},\bar{r}_{h},\bar{\mathbf{u}}_{h})\in\mathbf{V}_{h}\times Q_{h}\times\mathbf{V}_{h}\times Q_{h}\times\mathbb{U}_{ad}^{h}$ its numerical approximation given as the solution to (24)–(26). Then,

[TABLE]

*with $\vartheta>d$. The hidden constant is independent of the continuous and discrete solutions, the size of the elements in the mesh $\mathscr{T}_{h}$, and $\#\mathscr{T}_{h}$. The constant, however, blows up as $\lambda\downarrow 0$. *

Proof. We proceed in four steps.

Step 1. Let us consider $\mathbf{u}=\bar{\mathbf{u}}_{h}$ in (22) and $\mathbf{u}_{h}=\Pi_{\mathbf{L}^{2}}\bar{\mathbf{u}}$ in (25). Adding the obtained inequalities we arrive at

[TABLE]

Step 2. The goal of this step is to bound the term $(\bar{\mathbf{z}}-\bar{\mathbf{z}}_{h},\bar{\mathbf{u}}_{h}-\bar{\mathbf{u}})_{\mathbf{L}^{2}(\Omega)}$. To accomplish this task, we add and subtract the auxiliary term $\hat{\mathbf{z}}_{h}$, where $(\hat{\mathbf{z}}_{h},\hat{r}_{h})$ corresponds to the solution to (28), to obtain

[TABLE]

Let us concentrate on $\mathbf{I}$. Note that $(\bar{\mathbf{y}}_{h}-\hat{\mathbf{y}}_{h},\bar{p}_{h}-\hat{p}_{h})\in\mathbf{V}_{h}\times Q_{h}$ solves

[TABLE]

for all $\mathbf{v}_{h}\in\mathbf{V}_{h}$ and $q_{h}\in Q_{h}$, and that $(\bar{\mathbf{z}}_{h}-\hat{\mathbf{z}}_{h},\bar{r}_{h}-\hat{r}_{h})\in\mathbf{V}_{h}\times Q_{h}$ solves

[TABLE]

for all $\mathbf{w}_{h}\in\mathbf{V}_{h}$ and $s_{h}\in Q_{h}$. Set $\mathbf{w}_{h}=\bar{\mathbf{y}}_{h}-\hat{\mathbf{y}}_{h}\in\mathbf{V}_{h}$ in (41) and $\mathbf{v}_{h}=\bar{\mathbf{z}}_{h}-\hat{\mathbf{z}}_{h}\in\mathbf{V}_{h}$ in (40) to conclude

[TABLE]

To estimate $\mathbf{II}$ we proceed as follows. First, we use Young’s inequality to obtain

[TABLE]

Second, to estimate the term $\|\hat{\mathbf{z}}_{h}-\bar{\mathbf{z}}\|_{\mathbf{L}^{2}(\Omega)}^{2}$, we introduce $(\tilde{\mathbf{z}}_{h},\tilde{r}_{h})\in\mathbf{V}_{h}\times Q_{h}$ as the solution to

[TABLE]

Third, add and subtract $\tilde{\mathbf{z}}_{h}$ and use the triangle inequality to arrive at

[TABLE]

We now analyze $\|\tilde{\mathbf{z}}_{h}-\bar{\mathbf{z}}\|_{\mathbf{L}^{2}(\Omega)}$. Since $(\tilde{\mathbf{z}}_{h},\tilde{r}_{h})$ corresponds to the Galerkin approximation of $(\bar{\mathbf{z}},\bar{r})$, an application of Lemma 2.3 yields

[TABLE]

It thus remains to estimate $\|\hat{\mathbf{z}}_{h}-\tilde{\mathbf{z}}_{h}\|_{\mathbf{L}^{2}(\Omega)}$. To accomplish this, we invoke the weighted Poincaré inequality of Theorem 2.2 to obtain the estimate

[TABLE]

This, in view of the stability of the discrete Stokes system in weighted spaces [26, Theorem 4.1], allows us to obtain

[TABLE]

Now, let us recall that $\hat{\mathbf{y}}_{h}$ is the Galerkin approximation of $\bar{\mathbf{y}}$. In addition, since $d\in\{2,3\}$, $\Omega$ is a convex polytope, and $\bar{\mathbf{u}}\in\mathbf{L}^{\infty}(\Omega)$, we have \bar{\mathbf{y}}\in{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mathbf{W}}^{1,\infty}(\Omega) [52] (see also [39, 33]). Therefore, the pointwise error estimates of [25, Theorem 4.1], combined with the weighted estimates of [24] for $d=3$, yield

[TABLE]

with $\vartheta>d$. To obtain the last estimate, we have used [33, estimate (1.17)]. Recall that $\mathbf{I}\leq 0$. We thus replace (47) into (46) and combine the obtained estimate with (45) to conclude that

[TABLE]

with $\vartheta>d$, and $C_{1}$ and $C_{2}$ being positive constants.

Step 3. The goal of this step is to bound the remaining term in (39). Note that, by adding and subtracting the term $\lambda\bar{\mathbf{u}}$ and the adjoint variables $\bar{\mathbf{z}}$ and $\hat{\mathbf{z}}_{h}$, we obtain the following identity:

[TABLE]

We now bound the terms $\mathbf{III}_{1}$, $\mathbf{III}_{2}$, $\mathbf{III}_{3}$, and $\mathbf{III}_{4}$. To estimate $\mathbf{III}_{1}$ we invoke the error estimates (36) and (37). We thus proceed as follows:

[TABLE]

Now we estimate $\mathbf{III}_{2}$. We begin with a simple application of the Cauchy–Schwarz inequality and the Poincaré inequality of Theorem 2.2 to write

[TABLE]

Invoking the stability of the discrete Stokes system in weighted spaces of [26, Theorem 4.1], we obtain

[TABLE]

To estimate the right hand side of the previous expression we introduce the auxiliary variables $(\hat{\mathbf{y}},\hat{p})\in\mathbf{H}_{0}^{1}(\Omega)\times L^{2}(\Omega)/\mathbb{R}$ as the solution to

[TABLE]

Since the pair $(\bar{\mathbf{y}}-\hat{\mathbf{y}},\bar{p}-\hat{p})$ solves the Stokes system with the term $\bar{\mathbf{u}}-\bar{\mathbf{u}}_{h}$ in the right hand side of the momentum equation and $(\hat{\mathbf{y}}_{h}-\bar{\mathbf{y}}_{h},\hat{p}_{h}-\bar{p}_{h})$ corresponds to the Galerkin approximation of $(\bar{\mathbf{y}}-\hat{\mathbf{y}},\bar{p}-\hat{p})$, it follows from Lemma 3.3 that

[TABLE]

where we have considered, in the notation of Lemma 3.3, $(\bm{\xi},\theta)=(\bar{\mathbf{y}}-\hat{\mathbf{y}},\bar{p}-\hat{p})\in\mathbf{H}_{0}^{1}(\Omega)\times L^{2}(\Omega)/\mathbb{R}$ and $(\bm{\xi}_{h},\theta_{h})=(\hat{\mathbf{y}}_{h}-\bar{\mathbf{y}}_{h},\hat{p}_{h}-\bar{p}_{h})\in\mathbf{V}_{h}\times Q_{h}$. We thus invoke (50), (52), the error estimate (37), and Young’s inequality to conclude that

[TABLE]

where $C$ denotes a positive constant that is independent of $h$ and $\lambda$.

We bound the term $\mathbf{III}_{3}$ by using Young’s inequality and the estimates provided in (44)–(47). These arguments reveal that

[TABLE]

with $\vartheta>d$. To estimate $\mathbf{III}_{4}$ we use Young’s inequality to immediately arrive at

[TABLE]

where $C$ denotes a positive constant that is independent of $h$ and $\lambda$.

Step 4. The proof concludes by gathering (39), (48), (49), (53), (54), and (55).

Remark 3.8** (Rates of convergence for $\bar{\mathbf{u}}$).**

The error estimate of Theorem 3.7 reads as follows:

[TABLE]

*for $d=2$ and $d=3$, respectively. The two–dimensional error estimate is nearly–optimal in terms of approximation (nearly because of the presence of the $\log$-term). This rate of convergence is dictated by the regularity properties of $\bar{\mathbf{u}}$, namely, $\bar{\mathbf{u}}\in\mathbf{L}^{\infty,d/(d-1)}(\Omega)$ and the polynomial degree that is used for its approximation.

The following result establishes rates of convergence for the errors $\bar{\mathbf{y}}-\bar{\mathbf{y}}_{h}$, $\bar{p}-\bar{p}_{h}$, and $\bar{\mathbf{z}}-\bar{\mathbf{z}}_{h}$.

Theorem 3.9** (Rates of convergence).**

Let $(\bar{\mathbf{y}},\bar{p},\bar{\mathbf{z}},\bar{r},\bar{\mathbf{u}})\in\mathbf{H}_{0}^{1}(\Omega)\times L^{2}(\Omega)/\mathbb{R}\times\mathbf{H}_{0}^{1}(\rho,\Omega)\times L^{2}(\rho,\Omega)/\mathbb{R}\times\mathbb{U}_{ad}$ be the solution to the optimality system (18), (21), and (22) and $(\bar{\mathbf{y}}_{h},\bar{p}_{h},\bar{\mathbf{z}}_{h},\bar{r}_{h},\bar{\mathbf{u}}_{h})\in\mathbf{V}_{h}\times Q_{h}\times\mathbf{V}_{h}\times Q_{h}\times\mathbb{U}_{ad}^{h}$ its numerical approximation given as the solution to (24)–(26). Then,

[TABLE]

[TABLE]

and

[TABLE]

*with $\vartheta>d$. The hidden constants are independent of the continuous and discrete solutions, the size of the elements in the mesh $\mathscr{T}_{h}$, and $\#\mathscr{T}_{h}$. The constants, however, blow up as $\lambda\downarrow 0$. *

Proof. We first control the error $\bar{\mathbf{y}}-\bar{\mathbf{y}}_{h}$. To accomplish this task, we invoke the pair $(\hat{\mathbf{y}}_{h},\hat{p}_{h})$, defined as the solution to (27), and use the triangle inequality to write

[TABLE]

The first term on the right hand side of the previous expression is bounded in (47). In view of (52), the second term can be bounded by $\|\bar{\mathbf{u}}-\bar{\mathbf{u}}_{h}\|_{\mathbf{L}^{2}(\Omega)}$. The desired estimate (56) follows by collecting the previous estimates and the one obtained in Theorem 3.7.

We now control $\bar{p}-\bar{p}_{h}$. We invoke, again, the auxiliary variable $(\hat{\mathbf{y}}_{h},\hat{p}_{h})$, defined as the solution to (27), and the triangle inequality to obtain

[TABLE]

Since $(\hat{\mathbf{y}}_{h},\hat{p}_{h})$ corresponds to the Galerkin approximation of $(\bar{\mathbf{y}},\bar{p})$, then [28, Proposition 4.16] yields the error estimate

[TABLE]

where to obtain the last estimate we have used (20). On the other hand, since $\mathbf{V}_{h}\times Q_{h}$ satisfy the so–called inf-sup condition [28, Proposition 4.13], we can thus invoke the discrete problem that $(\hat{\mathbf{y}}_{h}-\bar{\mathbf{y}}_{h},\hat{p}_{h}-\bar{p}_{h})$ solves to obtain

[TABLE]

The result of Theorem 3.7 allows us to control $\|\bar{\mathbf{u}}-\bar{\mathbf{u}}_{h}\|_{\mathbf{L}^{2}(\Omega)}$.

Finally, we bound $\bar{\mathbf{z}}-\bar{\mathbf{z}}_{h}$. We begin with the estimate

[TABLE]

where $(\tilde{\mathbf{z}}_{h},\tilde{r}_{h})$ solves (43). The term $\|\bar{\mathbf{z}}-\tilde{\mathbf{z}}_{h}\|_{\mathbf{L}^{2}(\Omega)}$ can be estimated by using the result of Lemma 2.3. The remaining term is bounded by using the result of Theorem 2.2 in conjunction with the stability, in weighted spaces, of the discrete Stokes system [26, Theorem 4.1]:

[TABLE]

The proof concludes by invoking (56).

3.3 A Semi Discrete Scheme: Error Estimates

In this section, we propose a semidiscrete scheme for the pointwise tracking optimal control problem that is based on the so-called variational discretization approach [42]. This approach discretizes only the state space (the control space $\mathbb{U}_{ad}$ is not discretized) and induces a discretization of the optimal control variable by projecting the optimal discrete adjoint state into the admissible control set.

The semidiscrete scheme reads as follows: Find min $J(\mathbf{y}_{h},\mathbf{q})$ subject to the discrete state equations

[TABLE]

and the control constraints $\mathbf{q}\in\mathbb{U}_{ad}$. Standard arguments yield the existence and uniqueness of an optimal solution $((\bar{\mathbf{y}}_{h},\bar{p}_{h}),\bar{\mathbf{q}})$. In addition, $((\bar{\mathbf{y}}_{h},\bar{p}_{h}),\bar{\mathbf{q}})$ is optimal for the semidiscrete scheme if and only if $(\bar{\mathbf{y}}_{h},\bar{p}_{h})$ solves (59), with $\mathbf{q}$ replaced by $\bar{\mathbf{q}}$, and $\bar{\mathbf{q}}$ satisfies the variational inequality

[TABLE]

where $(\bar{\mathbf{z}}_{h},\bar{r}_{h})$ solves (26) with ${\mathbf{y}}_{h}$ replaced by $\bar{\mathbf{y}}_{h}$.

In the next result we provide error estimates for the semidiscrete scheme.

Theorem 3.10** (Rates of convergence for $\bar{\mathbf{u}}$).**

Let $(\bar{\mathbf{y}},\bar{p},\bar{\mathbf{z}},\bar{r},\bar{\mathbf{u}})\in\mathbf{H}_{0}^{1}(\Omega)\times L^{2}(\Omega)/\mathbb{R}\\ \times\mathbf{H}_{0}^{1}(\rho,\Omega)\times L^{2}(\rho,\Omega)/\mathbb{R}\times\mathbb{U}_{ad}$ be the solution to the optimality system (18), (21), and (22) and $(\bar{\mathbf{y}}_{h},\bar{p}_{h},\bar{\mathbf{z}}_{h},\bar{r}_{h},\bar{\mathbf{q}})\in\mathbf{V}_{h}\times Q_{h}\times\mathbf{V}_{h}\times Q_{h}\times\mathbb{U}_{ad}$ its numerical approximation given as the solution to the semidiscrete scheme. Then, we have

[TABLE]

*The hidden constant is independent of the continuous and discrete solutions, the size of the elements in the mesh $\mathscr{T}_{h}$, and $\#\mathscr{T}_{h}$. The constant, however, blows up as $\lambda\downarrow 0$. *

Proof. Set $\mathbf{u}=\bar{\mathbf{q}}$ and $\mathbf{q}=\bar{\mathbf{u}}$ in the variational inequalities (22) and (60), respectively. Add the obtained inequalities to arrive at

[TABLE]

The desired estimate (61) thus follows from the arguments elaborated in the proof of Theorem 3.7.

Remark 3.11** (variational discretization).**

There is an improvement over the fully discrete scheme of section 3.2. The error estimate of Theorem 3.10 reads

[TABLE]

*for $d=2$ and $d=3$, respectively. Notice that, for $d=2$, this error estimate coincides with the one obtained in Theorem 3.7. In particular, it is nearly–optimal in terms of approximation. When $d=3$, (61) is optimal in terms of regularity and improves upon (38). *

4 The Optimal Control Problem with Singular Sources

In this section, we precisely describe and analyze the optimal control problem with point sources (5)–(7) introduced in Section 1.2. We begin by defining the weight $\rho$ as in (9)–(10) with obvious modifications that basically entails replacing $\mathcal{E}$ by $\mathcal{D}$. We recall that the cost functional $\mathfrak{J}$ and the set of admissible controls $\mathfrak{U}_{ad}$ are defined by (5) and (7), respectively. Before we proceed with our analysis, we comment that when deriving a priori error estimates for suitable finite element approximations of problem (5)–(7), it will be essential to assume the existence of $d_{\mathcal{D}}>0$ such that

[TABLE]

The weak version of the optimal control problem with point sources reads as follows: Find

[TABLE]

subject to the following weak formulation of the state equation (6): Find $(\mathbf{y},p)\in\mathbf{H}_{0}^{1}(\rho,\Omega)\times L^{2}(\rho,\Omega)/\mathbb{R}$ such that

[TABLE]

We recall that $\langle\cdot,\cdot\rangle$ denotes the duality pairing between $\mathbf{H}_{0}^{1}(\rho^{-1},\Omega)^{\prime}$ and $\mathbf{H}_{0}^{1}(\rho^{-1},\Omega)$. If $\alpha\in(d-2,d)$ problem (64) is well–posed; see Section 2.3.1 for details. Finally, we mention that in view of the continuous embedding of Theorem 2.2, $\mathfrak{J}$ is well defined over $\mathbf{H}_{0}^{1}(\rho,\Omega)\times\mathfrak{U}_{ad}$. This further restricts $\alpha$ to belong to the interval $(d-2,2)\subset(d-2,d)$.

To analyze the optimal control problem with point sources we introduce the so-called control-to-state operator

[TABLE]

where $\mathbf{y}=\mathcal{C}\mathcal{U}$ solves problem (64). Since $\alpha\in(d-2,2)$, the map $\mathcal{C}$ is well defined. We can thus define the reduced cost functional

[TABLE]

We immediately conclude that $\mathfrak{j}$ is weakly lower semicontinuous and strictly convex $(\lambda>0)$. This, combined with the fact that $\mathfrak{U}_{ad}$ is compact, allow us to conclude the existence and uniqueness of an optimal control $\bar{\mathcal{U}}\in\mathfrak{U}_{ad}$ and an optimal state $\bar{\mathbf{y}}=\mathcal{C}\bar{\mathcal{U}}\in\mathbf{H}_{0}^{1}(\rho,\Omega)$ that satisfies (64) [60, Theorem 2.14]. In addition, the control variable $\bar{\mathcal{U}}$ is optimal for our optimal control problem if and only if [60, Lemma 2.21]

[TABLE]

To explore this variational inequality, we introduce the adjoint pair $(\mathbf{z},r)\in\mathbf{H}_{0}^{1}(\Omega)\times L^{2}(\Omega)/\mathbb{R}$, which satisfies

[TABLE]

Since $\mathbf{y}-\mathbf{y}_{\Omega}\in\mathbf{L}^{2}(\Omega)$, the well–posedness of problem (66) is immediate. Moreover, since $\Omega$ is convex, the results of Proposition 3.1 yield $(\mathbf{z},r)\in\mathbf{H}^{2}(\Omega)\times H^{1}(\Omega)$. This combined with [26, Proposition 2.3] reveal that the adjoint pair $(\mathbf{z},r)\in\mathbf{H}_{0}^{1}(\rho^{-1},\Omega)\times L^{2}(\rho^{-1},\Omega)$. This result is important because it allows to set $(\mathbf{v},q)=(\mathbf{z},r)$ as a test function in problem (64).

With these ingredients at hand, we proceed to further explore optimality conditions for problem (63)–(64). Note that, since $(\mathbf{y},p)\in\mathbf{H}_{0}^{1}(\rho,\Omega)\times L^{2}(\rho,\Omega)/\mathbb{R}\setminus\mathbf{H}_{0}^{1}(\Omega)\times L^{2}(\Omega)/\mathbb{R}$, we are not allow to set $(\mathbf{y},p)$ as a test function in problem (66). We thus have to proceed on the basis of different arguments.

Theorem 4.1** (Optimality conditions).**

Let $\alpha\in(d-2,2)$. The pair $(\bar{\mathbf{y}},\bar{\mathcal{U}})\in\mathbf{H}^{1}_{0}(\rho,\Omega)\times\mathfrak{U}_{ad}$ is optimal for problem (63)–(64) if and only if $\bar{\mathbf{y}}=\mathcal{C}\bar{\mathcal{U}}$ and the optimal control $\bar{\mathcal{U}}$ satisfies the variational inequality

[TABLE]

*where $(\bar{\mathbf{z}},\bar{r})\in\mathbf{H}_{0}^{1}(\Omega)\times L^{2}(\Omega)/\mathbb{R}$ corresponds to the optimal adjoint state, which solves (66) with $\mathbf{y}$ replaced by $\bar{\mathbf{y}}=\mathcal{C}\bar{\mathcal{U}}$. *

Proof. A simple computation shows that the variational inequality (65) can be rewritten as

[TABLE]

for all $\mathcal{U}=(\mathbf{u}_{1},\ldots,\mathbf{u}_{l})\in\mathfrak{U}_{ad}$. In what follows, to simplify the presentation of the material, we let $\mathbf{\mathbf{y}}=\mathcal{C}\mathcal{U}$. Let us concentrate on the first term of the left hand side of the previous expression. To study such a term, we note that $(\mathbf{y}-\bar{\mathbf{y}},p-\bar{p})$ solves

[TABLE]

for all $\mathbf{v}\in\mathbf{H}_{0}^{1}(\rho^{-1},\Omega)$ and $q\in L^{2}(\rho^{-1},\Omega)/\mathbb{R}$, respectively. Since the variable $\bar{\mathbf{z}}\in\mathbf{H}_{0}^{1}(\rho^{-1},\Omega)$, we are allowed to set $\mathbf{v}=\bar{\mathbf{z}}$ and $q=0$ in (69). This yields

[TABLE]

With this identity at hand, a density argument allows us to conclude

[TABLE]

In fact, let $\{\mathsf{y}_{n}\}_{n\in\mathbb{N}}\subset\mathbf{C}^{\infty}_{0}(\Omega)$ be such that $\mathsf{y}_{n}\rightarrow\mathbf{y}-\bar{\mathbf{y}}$ in $\mathbf{H}_{0}^{1}(\rho,\Omega)$. Since, for $n\in\mathbb{N}$, $\mathsf{y}_{n}$ is smooth, we can set $\mathbf{w}=\mathsf{y}_{n}$ and $s=0$ in (66). This yields

[TABLE]

Now, observe that, on the basis of Theorem 2.2, we have

[TABLE]

On the other hand, since $\bar{r}\in L^{2}(\rho^{-1},\Omega)$, we can set $q=\bar{r}$ in (69) to arrive at $b(\mathbf{y}-\bar{\mathbf{y}},\bar{r})=0$. This and the continuity of the bilinear form $b$ on $\mathbf{H}_{0}^{1}(\rho,\Omega)\times L^{2}(\rho^{-1},\Omega)$ imply that $b(\mathsf{y}_{n},\bar{r})$ converges to [math] as $n\uparrow\infty$. Finally, since $\bar{\mathbf{z}}\in\mathbf{H}_{0}^{1}(\rho^{-1},\Omega)$, the continuity of bilinear form $a$ on $\mathbf{H}_{0}^{1}(\rho,\Omega)\times\mathbf{H}_{0}^{1}(\rho^{-1},\Omega)$ allows us to conclude that $a((\mathbf{y}-\bar{\mathbf{y}})-\mathsf{y}_{n},\bar{\mathbf{z}})$ tends to [math] as $n\uparrow\infty$. The collection of these arguments yield the required identity (71).

The proof concludes upon using (68), (70), and (71).

We now introduce, for each $t\in\mathcal{D}$, the projection operator

[TABLE]

With this operator at hand, similar arguments to the ones elaborated in the proof of [60, Lemma 2.26] reveal that $\bar{\mathcal{U}}=(\bar{\mathbf{u}}_{1},\ldots,\bar{\mathbf{u}}_{l})$ satisfies (67) if and only if

[TABLE]

To summarize, the pair $(\bar{\mathbf{y}},\bar{\mathcal{U}})$ is optimal for problem (5)–(7) if and only if $(\bar{\mathbf{y}},\bar{p},\bar{\mathbf{z}},\bar{r},\bar{\mathcal{U}})\in\mathbf{H}_{0}^{1}(\rho,\Omega)\times L^{2}(\rho,\Omega)/\mathbb{R}\times\mathbf{H}_{0}^{1}(\Omega)\times L^{2}(\Omega)/\mathbb{R}\times\mathfrak{U}_{ad}$ solves (64), (66), and (67).

4.1 Discretization and Error Estimates

We begin by introducing the discrete counterpart of (5)–(7), which reads as follows: Find min $\mathfrak{J}(\mathbf{y}_{h},\mathcal{U}_{h})$ subject to the discrete state equations

[TABLE]

and the control constraints $\mathcal{U}_{h}\in\mathfrak{U}_{ad}$. The spaces $\mathbf{V}_{h}$ and $Q_{h}$ are given by (15). We restrict ourselves to consider such a pair of finite element spaces in order to apply the recent local error estimates of [8]. We comment that no discretization is needed for the optimal control variable, since the admissible set $\mathfrak{U}_{ad}$ is a subset of a finite dimensional space.

Standard arguments reveal the existence of a unique optimal pair $(\bar{\mathbf{y}}_{h},\bar{\mathcal{U}}_{h})$. In addition, the pair $(\bar{\mathbf{y}}_{h},\bar{\mathcal{U}}_{h})$ is optimal for the previously defined discrete optimal control problem if and only if $\bar{\mathbf{y}}_{h}$ solves (72) and $\bar{\mathcal{U}}_{h}$ satisfies the variational inequality

[TABLE]

where $(\bar{\mathbf{z}}_{h},\bar{r}_{h})\in\mathbf{V}_{h}\times Q_{h}$ solves

[TABLE]

To provide an error analysis for the previous scheme we introduce the following auxiliary problem: Find $(\hat{\mathbf{y}}_{h},\hat{p}_{h})\in\mathbf{V}_{h}\times Q_{h}$ such that

[TABLE]

To simplify the presentation of the material, we define, for $\mathcal{U}=(\mathbf{u}_{1},...,\mathbf{u}_{l})\in\mathfrak{U}_{ad}$ and $\mathcal{V}=(\mathbf{v}_{1},...,\mathbf{v}_{l})\in\mathfrak{U}_{ad}$,

[TABLE]

If $\mathbf{w}\in\mathbf{C}(\bar{\Omega})$ and \mathcal{V}=(\mathbf{v}_{1},...,\mathbf{v}_{l})\in{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}[\mathbb{R}^{d}]^{l}}, $\langle\mathbf{w},\mathcal{V}\rangle_{\mathcal{D}}:=\sum_{t\in\mathcal{D}}\mathbf{w}(t)\cdot\mathbf{v}_{t}.$

With the discrete system (72)–(74) at hand, we are in conditions to present the main result of this section. Our arguments are inspired by ideas developed in [35, 50] and yield a nearly–optimal error estimate in terms of approximation.

Theorem 4.2** (Rates of convergence for $\mathcal{\bar{U}}$).**

Let $(\bar{\mathbf{y}},\bar{p},\bar{\mathbf{z}},\bar{r},\bar{\mathcal{U}})\in\mathbf{H}_{0}^{1}(\rho,\Omega)\times L^{2}(\rho,\Omega)/\mathbb{R}\times\mathbf{H}_{0}^{1}(\Omega)\times L^{2}(\Omega)/\mathbb{R}\times\mathfrak{U}_{ad}$ be the solution to the optimality system (64), (66), and (67) and $(\bar{\mathbf{y}}_{h},\bar{p}_{h},\bar{\mathbf{z}}_{h},\bar{r}_{h},\bar{\mathcal{U}}_{h})\in\mathbf{V}_{h}\times Q_{h}\times\mathbf{V}_{h}\times Q_{h}\times\mathfrak{U}_{ad}$ its numerical approximation given as the solution to (72)–(74). Let $\Omega_{0}$, $\Omega_{1}\subset\Omega$ be such that $\mathcal{D}\Subset\Omega_{0}\subset\Omega_{1}$ with $\mathrm{dist}(\bar{\Omega}_{0},\partial\Omega_{1})\geq\mathfrak{d}\geq\mathfrak{c}h$ for $\mathfrak{c}$ sufficiently large. If $d=2$ and $\mathbf{y}_{\Omega}\in\mathbf{L}^{\kappa}(\Omega)$, for every $\kappa\in(2,\infty)$, then

[TABLE]

*The hidden constant is independent of the continuous and discrete solutions, the size of the elements in the mesh $\mathscr{T}_{h}$, and $\#\mathscr{T}_{h}$. The constant, however, blows up as $\lambda\downarrow 0$. *

Proof. We proceed in $4$ steps.

Step 1. Set $\mathcal{U}=\bar{\mathcal{U}}_{h}$ in (67) and $\mathcal{U}=\bar{\mathcal{U}}$ in (73). Adding the obtained inequalities we arrive at the basic estimate

[TABLE]

Step 2. Define the discrete auxiliary variables $(\hat{\mathbf{z}}_{h},\hat{r}_{h})$ and $(\tilde{\mathbf{z}}_{h},\tilde{r}_{h})$ as the solutions to

[TABLE]

and

[TABLE]

respectively, where $(\hat{\mathbf{y}}_{h},\hat{p}_{h})$ solves (75).

We now invoke (76) and add and subtract $\hat{\mathbf{z}}_{h}$ and $\tilde{\mathbf{z}}_{h}$ to obtain

[TABLE]

Similar arguments to those elaborated in Step 2 of Theorem 3.7 allow us to conclude that $\mathbf{III}=-\|\bar{\mathbf{y}}_{h}-\hat{\mathbf{y}}_{h}\|_{\mathbf{L}^{2}(\Omega)}^{2}\leq 0$. As a result, we obtain

[TABLE]

Step 3. We now estimate the term $\mathbf{I}$. To accomplish this task, we first note that, since $d=2$, the velocity field that solves (64) satisfies $\bar{\mathbf{y}}\in\mathbf{W}^{1,\nu}(\Omega)$ for every $\nu<2$. A standard Sobolev embedding result thus implies that $\bar{\mathbf{y}}\in\mathbf{L}^{\kappa}(\Omega)$ for every $\kappa<\infty$. Consequently, $\bar{\mathbf{y}}-\mathbf{y}_{\Omega}\in\mathbf{L}^{\kappa}(\Omega)$ for every $\kappa<\infty$. We can thus apply the interior regularity results of [32, Theorem IV 4.1] to conclude that $(\bar{\mathbf{z}},\bar{r})\in\mathbf{W}^{2,\kappa}(\Omega_{1})\times W^{1,\kappa}(\Omega_{1})$ for every $\kappa<\infty$ together with the bound

[TABLE]

where we have also used the first item in Section 5.5 of [51]. One of the main ingredients in the proof of [32, Theorem IV 4.1] are the global regularity estimates of [32, Theorem IV 2.1], which in turn follow from the Calderón–Zygmund theorem. The constant $C_{\kappa}$ can be traced and behaves as the constant involved in Calderón–Zygmund theorem: $C_{\kappa}\leq C\kappa$, with $C$ being independent of $\kappa$; see [32, Remark II.11.2]. On the other hand, since $\Omega$ is convex we have at hand the global regularity result $(\bar{\mathbf{z}},\bar{r})\in\mathbf{H}^{2}(\Omega)\times H^{1}(\Omega)$; see Proposition 3.1. We are thus in position to apply [8, Theorem 6.3] combined with [8, Remark 2.18] to arrive at

[TABLE]

where the infimum is taken over the whole space $\mathbf{V}_{h}\times Q_{h}$. Utilize the aforementioned regularity results and standard interpolation error estimates to obtain

[TABLE]

for any $\kappa<\infty$. Inspired by [59, page 3], we thus set $\kappa=|\log h|$ to conclude

[TABLE]

where we have also used a stability estimate for the problem that $(\bar{\mathbf{y}},\bar{p})$ solves and the regularity results of Proposition 3.1.

Step 4. We conclude by estimating the term $\mathbf{II}$ in (77). To accomplish this task, we proceed on the basis of a duality argument. Let us define the pair $(\bm{\varphi},\pi)\in\mathbf{H}_{0}^{1}(\Omega)\times L^{2}(\Omega)/\mathbb{R}$ as the solution to

[TABLE]

where the pair $(\hat{\mathbf{y}}_{h},\hat{p}_{h})$ solves (75). Since $\|\text{sgn}(\bar{\mathbf{y}}-\hat{\mathbf{y}}_{h})\|_{\mathbf{L}^{\infty}(\Omega)}\leq 1$, we can apply, again, the results of [32, Theorem IV 4.1] to conclude that $(\bm{\varphi},\pi)\in\mathbf{W}^{2,\infty}(\Omega_{1})\times W^{1,\infty}(\Omega_{1})$. We are thus in position to invoke, again, [8, Theorem 6.3] combined with [8, Remark 2.18] to conclude that

[TABLE]

Here, $(\bm{\varphi}_{h},\pi_{h})$ corresponds to the finite element approximation of $(\bm{\varphi},\pi)$. Notice that, since $\Omega$ is convex, we immediately have that $(\bm{\varphi},\pi)\in\mathbf{H}^{2}(\Omega)\times H^{1}(\Omega)$. To obtain the last inequality, we have used that $\|\text{sgn}(\bar{\mathbf{y}}-\hat{\mathbf{y}}_{h})\|_{\mathbf{L}^{\infty}(\Omega)}\leq 1$.

On the other hand, we can apply [26, Proposition 2.3] to obtain that $(\bm{\varphi},\pi)\in\mathbf{H}^{2}(\Omega)\cap\mathbf{H}_{0}^{1}(\Omega)\times H^{1}(\Omega)\cap L^{2}(\Omega)/\mathbb{R}\hookrightarrow\mathbf{H}_{0}^{1}(\rho^{-1},\Omega)\times L^{2}(\rho^{-1},\Omega)$. Since the pair $(\bm{\varphi},\pi)$ solves (81), similar arguments to the ones that allowed us to derive (71) yield

[TABLE]

Note that $(\hat{\mathbf{y}}_{h},\hat{p}_{h})$ corresponds to the finite element approximation, within the space $\mathbf{V}_{h}\times Q_{h}$, of $(\bar{\mathbf{y}},\bar{p})$, the solution to problem (64). We thus invoke Galerkin orthogonality, twice, and set $\mathbf{v}=\bm{\varphi}-\bm{\varphi}_{h}\in\mathbf{H}_{0}^{1}(\rho^{-1},\Omega)$ in (64) to arrive at

[TABLE]

Finally, apply (82) to obtain the error estimate

[TABLE]

With the previous estimates at hand, we can thus bound $\|\tilde{\mathbf{z}}_{h}-\hat{\mathbf{z}}_{h}\|_{\mathbf{L}^{\infty}(\Omega)}$. To accomplish this task, we invoke a standard inverse estimate [12, Lemma 4.9.2], the problem that $\tilde{\mathbf{z}}_{h}-\hat{\mathbf{z}}_{h}$ solves, and estimate (83). In fact, we have

[TABLE]

The proof concludes by gathering the estimates (78), (80), and (84).

Theorem 4.3** (Rates of convergence).**

Let $(\bar{\mathbf{y}},\bar{p},\bar{\mathbf{z}},\bar{r},\bar{\mathcal{U}})$ be the solution to the optimality system (64), (66), and (67), in $\mathbf{H}_{0}^{1}(\rho,\Omega)\times L^{2}(\rho,\Omega)/\mathbb{R}\times\mathbf{H}_{0}^{1}(\Omega)\times L^{2}(\Omega)/\mathbb{R}\times\mathfrak{U}_{ad}$, and $(\bar{\mathbf{y}}_{h},\bar{p}_{h},\bar{\mathbf{z}}_{h},\bar{r}_{h},\bar{\mathcal{U}}_{h})\in\mathbf{V}_{h}\times Q_{h}\times\mathbf{V}_{h}\times Q_{h}\times\mathfrak{U}_{ad}$ its numerical approximation given as the solution to (72)–(74). Let $\Omega_{0}$, $\Omega_{1}\subset\Omega$ be such that $\mathcal{D}\Subset\Omega_{0}\subset\Omega_{1}$ with $\mathrm{dist}(\bar{\Omega}_{0},\partial\Omega_{1})\geq\mathfrak{d}\geq\mathfrak{c}h$ for $\mathfrak{c}$ sufficiently large. If $d=2$ and $\mathbf{y}_{\Omega}\in\mathbf{L}^{\kappa}(\Omega)$ for every $\kappa\in(2,\infty)$, then

[TABLE]

and

[TABLE]

*with a similar estimate for $\|\bar{r}-\bar{r}_{h}\|_{L^{2}(\Omega)}$. The hidden constants are independent of the continuous and discrete solutions, the size of the elements in the mesh $\mathscr{T}_{h}$, and $\#\mathscr{T}_{h}$. The constants, however, blow up as $\lambda\downarrow 0$. *

Proof. We first estimate $\bar{\mathbf{y}}-\bar{\mathbf{y}}_{h}$. To accomplish this task, we invoke the pair $(\hat{\mathbf{y}}_{h},\hat{p}_{h})$, defined as the solution to (75), and write

[TABLE]

Apply now the Poincaré inequality of Theorem 2.2 and the stability of the discrete Stokes system in weighted spaces [26, Theorem 4.1] to obtain

[TABLE]

On the other hand, since $(\hat{\mathbf{y}}_{h},\hat{p}_{h})$ corresponds to the Galerkin approximation of $(\bar{\mathbf{y}},\bar{p})$, we estimate the second term on the right hand side of (85) in view of the error estimate (16). Collect the derived estimate to arrive at the desired estimate for $\|\bar{\mathbf{y}}-\bar{\mathbf{y}}_{h}\|_{\mathbf{L}^{2}(\Omega)}$.

Similar arguments can be used to estimate the terms $\|\nabla(\bar{\mathbf{z}}-\bar{\mathbf{z}}_{h})\|_{\mathbf{L}^{2}(\Omega)}$ and $\|\bar{r}-\bar{r}_{h}\|_{L^{2}(\Omega)}$. This concludes the proof.

Remark 4.4** (Rates of convergence for $\bar{\mathcal{U}}$ and $\bar{\mathbf{y}}$).**

The error estimate of Theorem 4.2 reads

[TABLE]

*This estimate is nearly–optimal in terms of approximation (nearly because of the presence of the $\log$-term) and improves the one derived in [7, Theorem 5.1], for the Poisson problem, which behaves as $\mathcal{O}(h^{2-\epsilon})$, for every $\epsilon>0$. The error estimate obtained in Theorem 4.3 for the discretization of the state velocity field reads as follows: $\|\bar{\mathbf{y}}-\bar{\mathbf{y}}_{h}\|_{\mathbf{L}^{2}(\Omega)}\lesssim h$. We notice that this estimate is optimal in terms of regularity.

5 Numerical Examples

In this section, we conduct a series of numerical examples that illustrate the performance of the discrete schemes (24)–(26) and (72)–(74) when approximating the solutions to the optimization problems described in Sections 3 and 4, respectively.

5.1 Implementation

All the experiments have been carried out with the help of a code that we implemented using C++. All matrices have been assembled exactly. The right hand sides as well as the approximation errors are computed with the help of a quadrature formula that is exact for polynomials of degree $19$ for two dimensional domains and degree $14$ for three dimensional domains. The global linear systems were solved using the multifrontal massively parallel sparse direct solver (MUMPS) [5, 6].

In all the examples we set $\lambda=1$ and $\Omega=(0,1)^{d}$ with $d\in\{2,3\}$. For a given partition $\mathscr{T}_{h}$ of $\Omega$, for the first problem, we seek $(\bar{\mathbf{y}}_{h},\bar{p}_{h},\bar{\mathbf{z}}_{h},\bar{r}_{h},\bar{\mathbf{u}}_{h})\,\in\mathbf{V}_{h}\times Q_{h}\times\mathbf{V}_{h}\times Q_{h}\times\mathbb{U}_{ad}^{h}$ that solves the discrete optimality system (24)–(26), while for the second problem, we seek $(\bar{\mathbf{y}}_{h},\bar{p}_{h},\bar{\mathbf{z}}_{h},\bar{r}_{h},\bar{\mathcal{U}}_{h})\,\in\mathbf{V}_{h}\times Q_{h}\times\mathbf{V}_{h}\times Q_{h}\times\mathfrak{U}_{ad}$ that solves the discrete optimality system (72)–(74). In all the numerical examples we make use of the Taylor–Hood element defined in (15). To solve the associated minimization problems, we use the Newton–type primal–dual active set strategy as described in [60, Section 2.12.4].

We consider problems with inhomogeneous Dirichlet boundary conditions whose exact solutions are known. Note that this violates the assumption of homogeneous Dirichlet boundary conditions which is needed for the analysis that we have performed, but it retains its essential difficulties and singularities and allows us to evaluate the rates of convergences. In both problems we construct exact solutions in terms of fundamental solutions of the Stokes equations [32, Section IV.2]:

[TABLE]

where $\{\mathbf{e}_{i}\}_{i=1}^{d}$ denotes the canonical basis of $\mathbb{R}^{d}$ and

[TABLE]

with $\mathbf{r}_{t}=x-t$. Here, $\mathbb{I}_{d}$ denotes the identity matrix in $\mathbb{R}^{d\times d}$.

Finally, we define $\mathbf{e}_{\mathbf{y}}:=\bar{\mathbf{y}}-\bar{\mathbf{y}}_{h}$, $e_{p}:=\bar{p}-\bar{p}_{h}$, $\mathbf{e}_{\mathbf{z}}:=\bar{\mathbf{z}}-\bar{\mathbf{z}}_{h}$, $e_{r}:=\bar{r}-\bar{r}_{h}$, $\mathbf{e}_{\mathbf{u}}:=\bar{\mathbf{u}}-\bar{\mathbf{u}}_{h}$, and $\mathbf{e}_{\mathcal{U}}:=\bar{\mathcal{U}}-\bar{\mathcal{U}}_{h}$.

5.2 Optimization with Point Observations

The finite sequence of vectors $\{\mathbf{y}_{t}\}_{t\in\mathcal{D}}$ is computed from the constructed solutions in such a way that the adjoint system (21) holds. A straightforward computation reveals that, for $t\in\mathcal{D}$, $\mathbf{y}_{t}=\bar{\mathbf{y}}(t)-(\mathbf{e}_{1}+\cdots+\mathbf{e}_{d}).$ In order to simplify the construction of exact solutions, we have incorporated, in the momentum equation of (18), an extra forcing term $\mathbf{f}\in\mathbf{L}^{\infty}(\Omega)$. With such a modification, the right hand side of the momentum equation reads as follows: $(\mathbf{f}+\mathbf{u},\mathbf{v})_{\mathbf{L}^{2}(\Omega)}$. Finally, we will denote the total number of degrees of freedom as $\mathsf{Ndof}=2\dim(\mathbf{V}_{h})+2\dim(Q_{h})+\dim(\mathbf{U}_{h}).$

 

Example 1. We let $\Omega=(0,1)^{2}$, $\alpha=1.5$, $\mathbf{a}=(-5,-5)^{T}$, $\mathbf{b}=(5,5)^{T}$, and $\mathcal{Z}=\{(0.25,0.25),(0.25,0.75),(0.75,0.25),(0.75,0.75)\}.$ We define the exact optimal state as

[TABLE]

and $\bar{p}(x_{1},x_{2})=x_{1}x_{2}(1-x_{1})(1-x_{2})-{1}/{36}$, while the exact optimal adjoint state is taken to be as in (86).

 

Example 2. We set $\Omega=(0,1)^{3}$, $\mathbf{a}=(-10,-10,-10)^{T}$, $\mathbf{b}=(2,2,2)^{T}$, $\alpha=1.99$, and $\mathcal{Z}=\left\{(0.5,0.5,0.5)\right\}.$ The exact optimal state is given by

[TABLE]

and $\bar{p}(x_{1},x_{2},x_{3})=x_{1}x_{2}x_{3}-{1}/{8}$. The optimal adjoint state is as in (86).

We observe, in Fig. 2, that when approximating the optimal control variable $\bar{\mathbf{u}}$ and the adjoint velocity field $\bar{\mathbf{z}}$, the obtained experimental rates of convergence are in agreement with the estimates provided in (38) and (58), respectively.

5.3 Optimization with Singular Sources

We now explore the performance of the discrete scheme (72)–(74) with $d=2$. In this case, the number of degrees of freedom is given by $\mathsf{Ndof}=2\dim(\mathbf{V}_{h})+2\dim(Q_{h})+2l$, where $l=\#\mathcal{D}$.

 

Example 3. We let $\Omega=(0,1)^{2}$, $\alpha=1.99$, $\mathbf{a}_{t}=(0,0)^{T}$, and $\mathbf{b}_{t}=(2,2)^{T}$ for all $t\in\mathcal{D}$ and $\mathcal{D}=\{(0.75,0.25)\}.$ We define the exact optimal adjoint state as follows

[TABLE]

and $\bar{r}(x_{1},x_{2})=x_{1}x_{2}(1-x_{1})(1-x_{2})-{1}/{36}$. The exact optimal state is as in (86).

We observe, in Fig. 3, that when approximating the state velocity field $\bar{\mathbf{y}}$, the experimental rate of convergence for this variable is in agreement with the estimate provided in Theorem 4.3. We also observe that the experimental rate of convergence for the error approximation of the optimal control variable $\bar{\mathcal{U}}$ is in agreement with Theorem 4.2. Finally we observe improved experimental rates of convergence for the adjoint variables $\bar{\mathbf{z}}$ and $\bar{r}$.