Finite element error analysis for measure-valued optimal control problems governed by a 1D wave equation with variable coefficients

TL;DR

This paper analyzes finite element error estimates for measure-valued optimal control problems governed by a 1D wave equation with variable coefficients, including numerical validation.

Contribution

It introduces three-level bilinear finite element discretizations for measure-valued controls in 1D wave equations and derives associated error estimates.

Findings

Error estimates for the optimal state variable are established.

Numerical results confirm the theoretical error bounds.

The approach handles measure-valued controls effectively.

Abstract

This work is concerned with the optimal control problems governed by a 1D wave equation with variable coefficients and the control spaces $\mathcal M_T$ of either measure-valued functions $L_{w^*}^2(I,\mathcal M(\Omega))$ or vector measures $\mathcal M(\Omega,L^2(I))$. The cost functional involves the standard quadratic tracking terms and the regularization term $\alpha\|u\|_{\mathcal M_T}$ with $\alpha>0$. We construct and study three-level in time bilinear finite element discretizations for this class of problems. The main focus lies on the derivation of error estimates for the optimal state variable and the error measured in the cost functional. The analysis is mainly based on some previous results of the authors. The numerical results are included.

Equations

$$J(u)$$

$$\rho {\partial }_{tt}y-{\partial }_x(\kappa {\partial }_{x}y)=u\text{for}\text{leavevmode}(t,x)\in I\times \Omega =(0,T)\times (0,L)$$

$$u=\sum _{i=1}^n{u}_i(t){\delta }_{x_i},u_i\in L^2(I),\text{leavevmode}{x}_i\in \Omega ,$$

$$u=\sum _{i=1}^n{u}_i(t){\delta }_{x_i(t)},u_i\in L^2(I),x_i:I\to \Omega \text{is measurable},$$

$$\|\bar{y}-\bar{y}_{\tau,h}\|_{L^{2}(I\times\Omega)}=\mathcal{O}\big{(}({\tau+h})^{\alpha}\big{)},\quad|J(\bar{u})-J(\bar{u}_{\tau,h})|=\mathcal{O}\big{(}{(\tau+h)^{2/3}}\big{)}$$

$$J(y,u)=F(y)+\alpha \Vert u{\Vert }_{{\mathcal{M}}_T}\to \underset{u,y}{\min }$$

$$F(y):={\textstyle\frac{1}{2}}\big{(}\left\|y-z_{1}\right\|_{L^{2}(I,{H_{\rho}})}^{2}+\left\|y(T)-z_{2}\right\|_{{H_{\rho}}}^{2}+\left\|\rho\partial_{t}y(T)-z_{3}\right\|_{{\mathcal{V}_{\kappa}^{*}}}^{2}\big{)}$$

$$\{\begin{array}{rlrl}\rho {\partial }_{tt}y-{\partial }_x(\kappa {\partial }_{x}y)& =u& & \text{in}\text{leavevmode}I\times \Omega :=(0,T)\times (0,L)\\ y& =0& & \text{on}\text{leavevmode}I\times \partial \Omega \\ y=y^0,\text{leavevmode}{\partial }_{t}y& =y^1& & \text{in}\text{leavevmode}\{0\}\times \Omega .\end{array}$$

$$\Vert \cdot {\Vert }_V=\Vert {\partial }_x\cdot {\Vert }_H,\Vert \cdot {\Vert }_{V^2}=\Vert {\partial }_{xx}\cdot {\Vert }_H,\Vert \cdot {\Vert }_{V^3}=\Vert {\partial }_x(\kappa {\partial }_x\cdot ){\Vert }_V.$$

$$\Vert w{\Vert }_{H_{\rho }}=\Vert \sqrt{\rho }w{\Vert }_H,\Vert w{\Vert }_{{\mathcal{V}}_{\kappa }}=\Vert \sqrt{\kappa }{\partial }_{x}w{\Vert }_H,\Vert w{\Vert }_{{\mathcal{V}}_{\kappa }^{\ast }}=\underset{\Vert v{\Vert }_{{\mathcal{V}}_{\kappa }}\le 1}{\sup }⟨w,v{⟩}_{\Omega },$$

$$\displaystyle\|\mathbf{z}\|_{\mathcal{Y}}=\big{(}\|z_{1}\|_{L^{2}(I,{H_{\rho}})}^{2}+\|z_{2}\|_{{H_{\rho}}}^{2}+\|z_{3}\|_{{\mathcal{V}_{\kappa}^{*}}}^{2}\big{)}^{1/2},$$

$${\mathcal{C}}_0(\Omega ,L^2(I){)}^{\ast }\cong \mathcal{M}(\Omega ,L^2(I)),L^2(I,{\mathcal{C}}_0(\Omega ){)}^{\ast }\cong L_{w^{\ast }}^2(I,\mathcal{M}(\Omega )),$$

$$⟨u,v{⟩}_{{\mathcal{M}}_T,{\mathcal{C}}_T}:={\int }_{\Omega }{\int }_0^{T}v(x,t)\mathrm{d}u(x)\mathrm{d}t,⟨u,v{⟩}_{{\mathcal{M}}_T,{\mathcal{C}}_T}:={\int }_0^T{\int }_{\Omega }v(t,x)\mathrm{d}u(t)\mathrm{d}t$$

$$\mathcal{M}(\Omega ,L^2(I))\hookrightarrow L_{w^{\ast }}^2(I,\mathcal{M}(\Omega ))\hookrightarrow L^2(I,V^{\ast }).$$

$${B(y,v)}+\big{(}\rho\partial_{t}y(T),v(T)\big{)}_{H}=\int_{I}\langle u,v\rangle_{\Omega}\leavevmode\nobreak\ \mathrm{d}t+\big{(}\rho y^{1},v(0)\big{)}_{H}\leavevmode\nobreak\ \leavevmode\nobreak\ \forall v\in L^{2}(I,V)\cap H^{1}(I,H)$$

$$B(y,v):=-(\rho {\partial }_{t}y,{\partial }_{t}v{)}_{L^2(I\times \Omega )}+(\kappa {\partial }_{x}y,{\partial }_{x}v{)}_{L^2(I\times \Omega )},$$

$${\int }_I⟨\rho {\partial }_{tt}y,v{⟩}_{\Omega }+(\kappa {\partial }_{x}y,{\partial }_{x}v{)}_H\text{leavevmode}\mathrm{d}t={\int }_I⟨u,v{⟩}_{\Omega }\text{leavevmode}\mathrm{d}t\text{leavevmode}\forall v\in L^2(I,V)$$

$$\displaystyle\|y\|_{\mathcal{C}(\bar{I},V)}+\|\partial_{t}y\|_{\mathcal{C}(\bar{I},H)}+\|\partial_{tt}y\|_{L^{2}(I,V^{\ast})}\leq c\,\big{(}\|u\|_{X}+\|\mathbf{y}\|_{V\times H}\big{)}.$$

$$\|\partial_{tt}y\|_{\mathcal{C}(\bar{I},V^{\ast})}\leq c\,\big{(}\|u\|_{H^{1}(I,V^{\ast})}+\|\mathbf{y}\|_{V\times H}\big{)}.$$

$$\|y\|_{\mathcal{C}(\bar{I},V^{2})}+\|\partial_{t}y\|_{\mathcal{C}(\bar{I},V)}+\|\partial_{tt}y\|_{L^{2}(I,H)}\leq c\,\big{(}\|u\|_{X}+\|\mathbf{y}\|_{V^{2}\times V}\big{)}.$$

$$\|\partial_{tt}y\|_{\mathcal{C}(\bar{I},H)}\leq c\,\big{(}\|u\|_{H^{1}(I,H)}+\|\mathbf{y}\|_{V^{2}\times V}\big{)}.$$

$$\int_{I}-(\rho y,\partial_{t}v)_{H}+(\kappa\partial_{x}\mathcal{I}_{t}y,\partial_{x}v)_{H}\leavevmode\nobreak\ \mathrm{d}t+\big{(}\rho y(T),v(T)\big{)}_{H}\\ =\int_{I}\big{\langle}u,\mathcal{I}_{t}^{\ast}v\big{\rangle}_{\Omega}\leavevmode\nobreak\ \mathrm{d}t+\big{(}\rho y^{0},v(0)\big{)}_{H}+\langle\rho y^{1},(\mathcal{I}_{t}^{\ast}v)(0)\rangle_{\Omega}\quad\forall v\in L^{2}(I,V)\cap H^{1}(I,H).$$

$$\|y\|_{\mathcal{C}(\bar{I},H)}+\|\mathcal{I}_{t}y\|_{\mathcal{C}(\bar{I},V)}+\|\partial_{t}y\|_{\mathcal{C}(\bar{I},V^{\ast})}\leq c\,\big{(}\|u\|_{L^{2}(I,V^{\ast})}+\|\mathbf{y}\|_{H\times V^{*}}\big{)}.$$

$$\int_{I}\big{(}y,\rho\partial_{tt}v-\partial_{x}(\kappa\partial_{x}v)\big{)}_{H}\leavevmode\nobreak\ \mathrm{d}t-\big{(}\rho y(T),\partial_{t}v(T)\big{)}_{H}+\langle\rho\partial_{t}y(T),v(T)\rangle_{\Omega}\\ =\int_{I}\langle u,v\rangle_{\Omega}\leavevmode\nobreak\ \mathrm{d}t-\big{(}\rho y^{0},\partial_{t}v(0)\big{)}_{H}+\langle\rho y^{1},v(0)\rangle_{\Omega}$$

$$\{\begin{array}{rlrl}\rho {\partial }_{tt}\tilde{y}-{\partial }_x(\kappa {\partial }_x\tilde{y})& ={\mathcal{I}}_{t}u+\rho y^1& & \text{in}\text{leavevmode}I\times \Omega \\ \tilde{y}& =0& & \text{on}\text{leavevmode}I\times \partial \Omega \\ \tilde{y}=0,\text{leavevmode}{\partial }_t\tilde{y}& =y^0& & \text{in}\text{leavevmode}\{0\}\times \Omega \end{array}$$

$$\partial_{t}y=\partial_{tt}\tilde{y}=(1/\rho)\big{(}\mathcal{I}_{t}u+\partial_{x}(\kappa\partial_{x}\tilde{y})+\rho y^{1}\big{)}\in\mathcal{C}(\bar{I},V^{\ast}).$$

$$\int_{I}(\rho y,\partial_{tt}v)_{H}-\big{(}\kappa\partial_{x}\mathcal{I}_{t}y,\partial_{x}\partial_{t}v\big{)}_{H}\leavevmode\nobreak\ \mathrm{d}t-\big{(}\rho y(T),\partial_{t}v(T)\big{)}_{H}\\ =\int_{I}\big{\langle}u,-\mathcal{I}^{\ast}_{t}(\partial_{t}v)\big{\rangle}_{\Omega}\leavevmode\nobreak\ \mathrm{d}t-\big{(}\rho y^{0},\partial_{t}v(0)\big{)}_{H}+\langle\rho y^{1},-(\mathcal{I}^{\ast}_{t}\partial_{t}v)(0)\rangle_{\Omega}.$$

$$-\int_{I}\big{(}\kappa\partial_{x}\mathcal{I}_{t}y,\partial_{x}\partial_{t}v\big{)}_{H}\leavevmode\nobreak\ \mathrm{d}t=\int_{I}\big{(}\mathcal{I}_{t}y,L\partial_{t}v\big{)}_{H}\leavevmode\nobreak\ \mathrm{d}t\\ =-\int_{I}\big{(}y,Lv\big{)}_{H}\leavevmode\nobreak\ \mathrm{d}t+\big{(}(\mathcal{I}_{t}y)(T),Lv(T)\big{)}_{H}$$

$$\int_{I}\big{(}y,\rho\partial_{tt}v-\partial_{x}(\kappa\partial_{x}v)\big{)}_{H}\leavevmode\nobreak\ \mathrm{d}t-\big{(}\rho y(T),\partial_{t}v(T)\big{)}_{H}+\langle\partial_{x}(\kappa\partial_{x}\mathcal{I}_{t}y)(T),v(T)\rangle_{\Omega}\\ =\int_{I}\big{\langle}u,v\big{\rangle}_{\Omega}\leavevmode\nobreak\ \mathrm{d}t-\langle(\mathcal{I}_{t}u)(T),v(T)\rangle_{\Omega}-\big{(}\rho y^{0},\partial_{t}v(0)\big{)}_{H}+\langle\rho y^{1},v(0)\rangle_{\Omega}-\langle\rho y^{1},v(T)\rangle_{\Omega}.$$

$$⟨\rho {\partial }_{t}y(T),\phi {⟩}_{\Omega }=(({\mathcal{I}}_{t}u)(T),\phi {)}_H+⟨{\partial }_x(\kappa {\partial }_x{\mathcal{I}}_{t}y),\phi {⟩}_{\Omega }+⟨\rho y^1,\phi {⟩}_{\Omega }\forall \phi \in V,$$

$$\int_{I}\big{(}y,-\rho\partial_{t}v)_{H}-(y,\partial_{x}(\kappa\partial_{x}\mathcal{I}_{t}^{\ast}v)\big{)}_{H}\leavevmode\nobreak\ \mathrm{d}t+\big{(}\rho y(T),v(T)\big{)}_{H}\\ =\int_{I}\big{\langle}u,\mathcal{I}_{t}^{\ast}v\big{\rangle}_{\Omega}\leavevmode\nobreak\ \mathrm{d}t+\big{(}\rho y^{0},v(0)\big{)}_{H}+\langle\rho y^{1},(\mathcal{I}_{t}^{\ast}v)(0)\rangle_{\Omega}$$

Full text

The first and the second author were supported by FWF and DFG through the International Research Training Group IGDK 1754 ‘Optimization and Numerical Analysis for Partial Differential Equations with Nonsmooth Structures’. The third has been funded within the framework of the Academic Fund Program at the National Research University Higher School of Economics in 2016-2017 (grant no. 16-01-0054) and by the Russian Academic Excellence Project ‘5-100’. He also thanks the Technical University of Munich for its hospitality in 2014-2015 years.

Philip Trautmann

Department of Mathematics and Scientific Computing

University of Graz

Heinrichstraße 36

8010 Graz, Austria

Boris Vexler

Zentrum Mathematik

Technische Universität München

Boltzmannstraße 3

85748 Garching bei München, Germany

Alexander Zlotnik

Department of Mathematics at Faculty of Economic Sciences

National Research University Higher School of Economics

Myasnitskaya 20

101000 Moscow, Russia

(Communicated by the associate editor name)

1. Introduction

This work is concerned with the discretization and numerical analysis of optimal control problems involving a 1D linear wave equation with variable coefficients and controls taking values in certain measure spaces. The combination of variable coefficients and irregular data leads to significant technical problems.

Motivated by industrial applications as well as applications in the natural sciences, in which one is interested to place actuators in form of point sources in an optimal way, see, e.g., [4, 9] or in the reconstruction of point sources from given measurements, see, e.g., [34, 44], measure valued optimal control problems involving PDEs gained attention in the last years. These problems can be translated into optimization problems in terms of the coordinates and coefficients of the point sources. However, these optimization problem are non-convex since the solution of the state equation (PDE) depends in a non-linear way on the coordinates of the point sources. Thus one has to deal with multiple local minima. Several authors suggested to cast the control problem resp. inverse problem in form of an optimization problem over a suitable measure space $\mathcal{M}_{T}$ involving a convex regularization functional $R$ which favors point sources as solutions. In our case we introduce the following problem formulation involving the 1D wave equation

[TABLE]

with additional initial and boundary conditions. The functional $F$ is given by a quadratic tracking functional involving $y|_{I\times\Omega}$, $y(T,\cdot)|_{\Omega}$ and $\partial_{t}y(T,\cdot)|_{\Omega}$. The regularization functional $R$ and the control space $\mathcal{M}_{T}$ are chosen in a way such that $\mathcal{M}_{T}$ contains point sources of the desired form and $R$ promotes controls of such a form, i.e. linear combinations of point sources with time-dependent intensities or more general controls with a small spatial support. Since problem (1.1) is convex, one does not need to deal with several local minima. However, it is not longer guaranteed that the solution consists of a sum of point sources. We enforce such controls via the regularization functional $R$. Problems of the form (1.1) (also involving other PDEs) have been analysed from theoretical, numerical and algorithmic points of view, see [12, 11, 18, 19, 45, 33, 34, 13, 14, 7, 44, 16, 15]. Optimal control problems governed by the linear wave equation were discussed in several different aspects, see [35, 31, 32, 36, 52, 25, 24, 40, 41, 29, 30]. In our particular case we consider the control spaces $\mathcal{M}_{T}$ of measure-valued functions ${L_{w^{*}}^{2}(I,\mathcal{M}(\Omega))}$ and vector measures $\mathcal{M}(\Omega,L^{2}(I))$ with $R(u)=\alpha\|u\|_{\mathcal{M}_{T}}$. These two different choices imply different structural properties of the optimal controls. A typical non-regular element from the space $\mathcal{M}(\Omega,L^{2}(I))$ is given by

[TABLE]

where $\delta_{x_{i}}$ are the Dirac delta functions. Point sources of such type with fixed positions and time-dependent intensities are of interest in acoustics or geology, see [34, 44]. If one is interested in controls involving moving point sources of the form

[TABLE]

then the control space ${L_{w^{*}}^{2}(I,\mathcal{M}(\Omega))}$ rather than $\mathcal{M}(\Omega,L^{2}(I))$ is more appropriate. The space $\mathcal{M}(\Omega,L^{2}(I))$ and the functional $\|\cdot\|_{\mathcal{M}(\Omega,L^{2}(I))}$ are also related to the term directional sparsity resp. joint sparsity, see [26, 22].

For the discretization of optimal control problem (1.1), we discretize the state equation by space-time finite element method as introduced in [51]. Related methods are also discussed and analyzed in [1, 23], see also [2]. The measure-valued control is not directly discretized, cf. the variational control discretization from [27]. However, there exists optimal controls consisting of Dirac measures in the spatial grid points which can be computed, see also [12, 33]. The numerical analysis of the control problem is based on FEM error estimates for the second order hyperbolic equations from [51] and techniques developed in [12, 33]. It requires to overcome significant technical difficulties caused by non-smoothness of controls and states. To the best of our knowledge, this is the first paper providing such numerical analysis for the studied control problems.

The problem like (1.1) for a parabolic/heat state equation is analyzed for the case $\mathcal{M}_{T}=\mathcal{M}(\Omega,L^{2}(I))$ in [33] and for the case $\mathcal{M}_{T}={L_{w^{*}}^{2}(I,\mathcal{M}(\Omega))}$ in [12]. In particular, in both papers the authors prove existence of optimal controls and derive optimality conditions and FEM error estimates. Our analysis is partly based on these results of [33]. In [34] a problem similar to (1.1) involving the linear wave equation with constant coefficients as state equation is analyzed. In particular, existing regularity results for a Dirac right-hand side are extended to sources from $\mathcal{M}(\Omega,L^{2}(I))$. Based on these regularity results existence of optimal controls is proved as well as optimal conditions are derived in the 3D case.

Now we briefly sum up the contents of this work. First of all we collect and partially prove required existence and regularity results for the linear wave equation in the 1D setting. In particular, we check that the notions of a weaker solution defined in [51] and more commonly used very weak solution, e.g. [38], are equivalent. Most importantly we prove that the solution of the linear wave equation with variable coefficients from $H^{1}(\Omega)$ for any source term $u\in\mathcal{M}(\Omega,L^{2}(I))$ is an element of $\mathcal{C}(\bar{I},H^{1}_{0}(\Omega))\cap\mathcal{C}^{1}(\bar{I},L^{2}(\Omega))$ provided that the initial data have relevant regularity. The proof is based on a non-standard energy type bound in space, not only in time, cf. [37, 21]. In [34] the same result is proved for the wave equation with constant coefficients using duality techniques. This proof in [34] provides also corresponding results for multidimensional case but can not be directly extended for treating variable coefficients. This is due to the fact that it uses estimates of the solution of the wave equation in the whole space with a Dirac measure on the right hand-side which are proven using the Fourier-and Laplace-transformation or explicit solution formulas.

The existence of optimal controls and the derivation of optimality conditions are discussed on the basis of results from [34, 33]. In the case $\mathcal{M}_{T}=\mathcal{M}(\Omega,L^{2}(I))$ we prove that the optimal control $\bar{u}$ belongs to $\mathcal{C}^{1}(\bar{I},\mathcal{M}(\Omega))$.

Further, the FEM discretization of the state equation is introduced. The state variable $y_{h,\tau}$ belongs to the space of bilinear finite elements and is defined by the regularized Galerkin method. The resulting numerical scheme is a three-level method in time (i.e., its main equation relates the approximate solution values at three consecutive grid time levels). Moreover, we pose and prove the FEM error estimates in $\mathcal{C}(\bar{I},L^{2}(\Omega))$ for the discrete state equation which we need for the numerical analysis of the control problem. We base this study mainly on the results from [51] concerning error analysis of FEMs for the second order hyperbolic equations in the classes of the data having integer Sobolev or fractional Nikolskii order of smoothness. Note that their sharpness in a strong sense was stated in [50].

Then we consider a semi-discrete optimal control problem in which the continuous state equation is replaced by its discretized version whereas the controls are not discretized. We prove convergence of the discrete optimal controls to the continuous one and derive optimality conditions based on the Lagrange techniques. Most importantly we derive the discrete adjoint state equation. We can conclude that the first-discretize-then-optimize and first-optimize-then-discretize approaches commute. Therefore an analysis of the discrete adjoint state equation including the error estimates in $\mathcal{C}(\bar{I}\times\bar{\Omega})$ and $L^{2}(I,\mathcal{C}_{0}(\Omega))$ can also be based on techniques from [51]. Then we use results from [33] to represent the numerical error of state variable and of the cost functional in terms of FEM errors of the state equation and the adjoint state equation. Let $\bar{u}$ and $\bar{y}$ be the optimal control and the corresponding optimal state, and the variables $\bar{u}_{\tau,h}$ and $\bar{y}_{\tau,h}$ be their discrete counterparts. As the main result of this paper we prove the error estimates

[TABLE]

where $\tau$ is the step in time, $h$ is the maximal step in space and $\alpha=1/3$ for $\mathcal{M}_{T}={L_{w^{*}}^{2}(I,\mathcal{M}(\Omega))}$ or $\alpha=2/3$ for $\mathcal{M}_{T}=\mathcal{M}(\Omega,L^{2}(I))$. The latter higher order is due to the above mentioned improved regularity results for the state and optimal control. Such estimates are proved for the measure-valued controls in the hyperbolic case for the first time. Similar estimates are impossible in multidimensional settings due to much less fractional Sobolev regularity of optimal states and controls.

Finally we discuss the numerical computation of the discrete control $\bar{u}_{h,\tau}$. Based on a control discretization $u_{h,\tau}$ that given by the sum like (1.2) with $x_{i}$ at the spatial grid points and $u_{i}$ in the space of linear finite elements, a solution of the semi-discrete control problem can be calculated similarly to [33]. For the actual numerical computation of the optimal control we add the term $(\gamma/2)\|u\|_{L^{2}(I\times\Omega)}^{2}$, $\gamma>0$, to (1.1). This regularized problem is solved by a semi-smooth Newton method, see [43]. In a continuation strategy the regularization parameter $\gamma$ is made sufficiently small. We complete this work with a numerical example for $\mathcal{M}_{T}=\mathcal{M}(\Omega,L^{2}(I))$.

The paper is organized in the following way. In Section 2 we introduce the problem setting and the control spaces resp. the regularization functionals. Section 3 is concerned with regularity properties of the linear wave equation with variable coefficients in the 1D setting. In Section 4 the control problem is analyzed from a theoretical point of view. Section 5 deals with discretization of the state equation. Then we obtain stability bounds and error estimates for the discrete state equation in Section 6. Section 7 is concerned with the analysis of the semi-discrete optimal control problem. The next section discusses stability bounds and error estimates for the discrete adjoint sate equation. In Sections 9 resp. 10 error estimates for the optimal state and cost functional are derived being the main theoretical results of the study. Section 11 deals with the time stepping formulation of the discrete state equation. In Section 12 we discuss the control discretization with Dirac measures at the grid points. Then we introduce the $L^{2}(I\times\Omega)$ regularized problem and describe its solutions by a semi-smooth Newton method. Finally Section 13 provides a numerical example.

2. Problem setting

We consider optimal control problems of the following form

[TABLE]

with the parameter $\alpha>0$ and the tracking functional

[TABLE]

using $\mathbf{z}:=(z_{1},z_{2},z_{3})\in\mathcal{Y}:=L^{2}(I\times\Omega)\times L^{2}(\Omega)\times H^{-1}(\Omega)$, subject to the state equation which is an initial-boundary value problem for a 1D linear wave equation with variable coefficients

[TABLE]

Here, in particular, the initial data $\mathbf{y}:=(y^{0},y^{1})\in H^{1}_{0}(\Omega)\times L^{2}(\Omega)$, and $L>0$ and $T>0$. The coefficients $\rho,\kappa\in{H^{1}(\Omega)}$ satisfy $\rho(x)\geq\nu>0$ and $\kappa(x)\geq\nu$ on $\Omega$.

For brevity we denote $H=L^{2}(\Omega)$, $V=H^{1}_{0}(\Omega)$, $V^{2}=H^{2}(\Omega)\cap V$ and $V^{3}=\{v\in V|\partial_{x}(\kappa\partial_{x}v)\in V\}$ equipped with the norms

[TABLE]

Moreover, we utilize the equivalent coefficient-dependent Hilbert norms on $H$, $V$, $V^{\ast}$ and $\mathcal{Y}$

[TABLE]

where $\langle\cdot,\cdot\rangle_{\Omega}$ is the duality relation on $V^{*}\times V$.

For the control space $\mathcal{M}_{T}$ we consider two choices, either the space of vector measures $\mathcal{M}(\Omega,L^{2}(I))$ or the space of weak-star measurable, $\mathcal{M}(\Omega)$-valued functions $L_{w^{*}}^{2}(I,\mathcal{M}(\Omega))$. Recall that $\|u\|_{L_{w^{*}}^{2}(I,\mathcal{M}(\Omega))}=\|\|u(\cdot)\|_{\mathcal{M}(\Omega)}\|_{L^{2}(I)}$ where $\|u(\cdot)\|_{\mathcal{M}(\Omega)}\in L^{2}(I)$ for any $u\in L_{w^{*}}^{2}(I,\mathcal{M}(\Omega))$. Let correspondingly $\mathcal{C}_{T}$ be chosen as $\mathcal{C}_{0}(\Omega,L^{2}(I))$ or $L^{2}(I,\mathcal{C}_{0}(\Omega))$ where $\mathcal{C}_{0}(\Omega)=\{v\in\mathcal{C}(\bar{\Omega})|\,v|_{x=0,L}=0\}$. The following identifications of dual spaces hold

[TABLE]

with the duality pairings respectively

[TABLE]

for any $u\in\mathcal{M}_{T}$ and $v\in\mathcal{C}_{T}$. See [12, 17, 20, 33], where more details on the properties of these spaces can be found. In particular, the following embeddings hold

[TABLE]

3. Existence and regularity of the state

3.1. Weak formulations and preliminary existence, uniqueness and regularity results

In this section we introduce our solution concepts for the state equation (2.1). We begin with defining a weak formulation of (2.1).

Definition 3.1**.**

Let $(u,y^{0},y^{1})\in X\times V\times H$ with $X=L^{2}(I\times\Omega)$ or $H^{1}(I,V^{\ast})$ or $\mathcal{M}(\Omega,L^{2}(I))$. Then $y\in\mathcal{C}(\bar{I},V)\cap\mathcal{C}^{1}(\bar{I},H)$ is called a weak solution of (2.1) if it satisfies the integral identity

[TABLE]

with the indefinite symmetric bilinear form

[TABLE]

and the initial condition $y(0)=y^{0}$.

The right-hand side in (3.1) is well defined for $X=\mathcal{M}(\Omega,L^{2}(I))$ too due to embeddings (2.2).

Remark 1**.**

It is possible (and more common) to suppose that $v(T)=0$ in (3.1) when the last term on the left disappears (for example, see [51]). This leads to an equivalent formulation. To check this, it is enough to replace there $v$ by $v\beta_{\delta}$, where \beta_{\delta}(t)=\min\big{(}1,(T-t)/\delta\big{)}, $0<\delta<T$. Then $\partial_{t}(v\beta_{\delta})=(\partial_{t}v)\beta_{\delta}-(1/\delta)v\chi_{(T-\delta,T)}$, where $\chi_{(T-\delta,T)}$ is the characteristic function of $(T-\delta,T)$. Passing to the limit as $\delta\to 0$ with the help of the dominated convergence theorem and the properties of $y$ and $v$ leads to the result.

Another definition of the weak solution is possible.

Definition 3.2**.**

Let $(u,y^{0},y^{1})\in X\times V\times H$ with $X=L^{2}(I\times\Omega)$ or $H^{1}(I,V^{\ast})$ or $\mathcal{M}(\Omega,L^{2}(I))$. A function $y\in\mathcal{C}(\bar{I},V)\cap H^{2}(I,V^{\ast})\hookrightarrow H^{1}(I,H)$ is called a weak solution of (2.1) if it satisfies

[TABLE]

and $y(0)=y^{0}$ as well as $\partial_{t}y(0)=y^{1}$.

Proposition 1**.**

Definitions 3.1 and 3.2 (up to the property $y\in\mathcal{C}^{1}(\bar{I},H)$) are equivalent.

Proof.

The weak solution from Definition 3.1 has $\partial_{tt}y\in L^{2}(I,V^{*})$ according to the integral identity (3.1). Then the equivalence of (3.3) and (3.1) can be proved using integration by parts in time and the density of $\mathcal{C}^{\infty}(\bar{I},V)$ in $L^{2}(I,V)\cap H^{1}(I,H)$, cf. [38, Chapter 1, Theorem 2.1]. ∎

Proposition 2**.**

1. (1)

Let $(u,y^{0},y^{1})\in X\times V\times H$ with $X=L^{2}(I\times\Omega)$ or $H^{1}(I,V^{\ast})$. Then (2.1) has a unique weak solution satisfying $y\in\mathcal{C}(\bar{I},V)\cap\mathcal{C}^{1}(\bar{I},H)\cap H^{2}(I,V^{\ast})$ and

[TABLE]

Hereafter $c>0$, $c_{1}>0$, etc., are independent of $y$ and the data.

In the case $X=H^{1}(I,V^{\ast})$ there even holds $y\in\mathcal{C}^{2}(\bar{I},V^{\ast})$ as well as

[TABLE] 2. (2)

Let $(u,y^{0},y^{1})\in X\times V^{2}\times V$ with $X=L^{2}(I,V)$ or $H^{1}(I,H)$. Then the weak solution y satisfies $y\in\mathcal{C}(\bar{I},V^{2})\cap\mathcal{C}^{1}(\bar{I},V)\cap H^{2}(I,H)$ and

[TABLE]

In the case $X=H^{1}(I,H)$ there even holds $y\in\mathcal{C}^{2}(\bar{I},H)$ as well as

[TABLE]

Moreover, $y$ satisfies the equation $\rho\partial_{tt}y-\partial_{x}(\kappa\partial_{x}y)=u$ in $L^{2}(I\times\Omega)$, i.e. it is the strong solution.

Proof.

For example, see [51, Propositions 1.1 and 1.3]. ∎

Item 2 ensures the regularity of weak solution for more regular data.

For less regular data $(u,y^{0},y^{1})\in L^{2}(I,V^{\ast})\times H\times V^{\ast}$ one can use other weak formulations. To state the first of them, we define the integration operator $(\mathcal{I}_{t}v)(t):=\int_{0}^{t}v(s)\leavevmode\nobreak\ \mathrm{d}s$ and its adjoint $(\mathcal{I}_{t}^{\ast}v)(t):=\int_{t}^{T}v(s)\leavevmode\nobreak\ \mathrm{d}s$ on $\bar{I}$.

Definition 3.3**.**

Let $(u,y^{0},y^{1})\in L^{2}(I,V^{\ast})\times H\times V^{\ast}$. A function $y\in\mathcal{C}(\bar{I},H)$ with $\mathcal{I}_{t}y\in\mathcal{C}(\bar{I},V)$ is called a weaker solution of (2.1) if it satisfies

[TABLE]

As in the case of Definition 3.1, it is sufficient to take $v(T)=0$ in (3.6), cf. Remark 1.

Proposition 3**.**

Let $(u,y^{0},y^{1})\in L^{2}(I,V^{\ast})\times H\times V^{\ast}$. Then there exists a unique weaker solution $y\in\mathcal{C}(\bar{I},H)\cap\mathcal{C}^{1}(\bar{I},V^{\ast})$ and it satisfies the bound

[TABLE]

Proof.

See [51, Proposition 1.2]. ∎

We infer that there are other weak formulations of (2.1) for solutions $y\in\mathcal{C}(\bar{I},H)\cap\mathcal{C}^{1}(\bar{I},V^{\ast})$. One can use the concept of very weak solutions.

Definition 3.4**.**

Let $(u,y^{0},y^{1})\in L^{2}(I,V^{\ast})\times H\times V^{\ast}$. A function $y\in\mathcal{C}(\bar{I},H)\cap\mathcal{C}^{1}(\bar{I},V^{\ast})$ satisfying

[TABLE]

for any $v\in L^{2}(I,V^{2})\cap H^{2}(I,H)\hookrightarrow H^{1}(I,V)$ is called a very weak solution of (2.1).

Actually, these two last solution concepts are equivalent for the considered data spaces.

Theorem 3.5**.**

Definitions 3.3 and 3.4 are equivalent.

Proof.

First of all, we consider the auxiliary integrated in $t$ problem (2.1):

[TABLE]

for $(u,y^{0},y^{1})\in L^{2}(I,V^{\ast})\times H\times V^{\ast}$. Thus, we have $\mathcal{I}_{t}u\in H^{1}(I,V^{\ast})$. According to Proposition 2 problem (3.8) has a unique weak solution $\tilde{y}\in\mathcal{C}(\bar{I},V)\cap\mathcal{C}^{1}(\bar{I},H)$. Moreover, we set $y=\partial_{t}\tilde{y}$. Thus the weak formulation of (3.8) involving $\tilde{y}$ coincides with the weaker formulation of (2.1) involving $y$. Furthermore there holds $y=\partial_{t}\tilde{y}\in\mathcal{C}(\bar{I},H)$ and

[TABLE]

Now we take any $v\in\mathcal{C}^{\infty}(\bar{I},V^{2})$ and test (3.6) with $-\partial_{t}v$ in the role of $v$:

[TABLE]

Next we rearrange a term on the left integrating by parts in $x$ and $t$:

[TABLE]

with $Lv:=\partial_{x}(\kappa\partial_{x}v)$. Since $\mathcal{I}_{t}y\in\mathcal{C}(\bar{I},V)$, we get

[TABLE]

Since formula (3.9) implies that

[TABLE]

by the density of $\mathcal{C}^{\infty}(\bar{I},V^{2})$ in $L^{2}(I,V^{2})\cap H^{2}(I,H)$ we find that $y$ is a very weak solution of (2.1).

Now let $y\in\mathcal{C}(\bar{I},H)\cap\mathcal{C}^{1}(\bar{I},V^{\ast})$ be a very weak solution of (2.1). Then we take any $v\in\mathcal{C}^{\infty}(\bar{I},V^{2})$ and test (3.7) with $\mathcal{I}_{t}^{\ast}v$. Thus, we get

[TABLE]

and then

[TABLE]

The last equation yields that $L\mathcal{I}_{t}y=-\rho\partial_{t}y+\mathcal{I}_{t}u+\rho y^{1}\in\mathcal{C}(\bar{I},V^{\ast})$. Thus $\mathcal{I}_{t}y\in\mathcal{C}(\bar{I},V)$ and we can transform a term on the left in (3.11) by replacing $v$ by $\mathcal{I}_{t}^{\ast}v$ in (3.10):

[TABLE]

Then the density of $\mathcal{C}^{\infty}(\bar{I},V^{2})$ in $L^{2}(I,V)\cap H^{1}(I,H)$ shows that $y$ is a weaker solution of (2.1). ∎

Moreover, there is the concept of solutions by transposition.

Definition 3.6**.**

Let $(u,y^{0},y^{1})\in L^{2}(I,V^{\ast})\times H\times V^{\ast}$. A solution by transposition $y\in\mathcal{C}(\bar{I},H)\cap\mathcal{C}^{1}(\bar{I},V^{\ast})$ of (2.1) is defined by

[TABLE]

for all $(\phi,p^{0},p^{1})\in L^{2}(I\times\Omega)\times V\times H$ where $p\in\mathcal{C}(\bar{I},V)\cap\mathcal{C}^{1}(\bar{I},H)$ is the weak solution of the adjoint problem

[TABLE]

Proposition 4**.**

Definitions 3.4 and 3.6 are equivalent too.

Proof.

For $\phi\in H^{1}(I,H)$ or $L^{2}(I,V)$, $p^{0}\in V^{2}$ and $p^{1}\in V$ there holds $p\in\mathcal{C}(\bar{I},V^{2})\cap\mathcal{C}^{1}(\bar{I},V)\cap H^{2}(I,H)$, see Proposition 2. Due to the density of $H^{1}(I,H)$ resp. $L^{2}(I,V)$ in $L^{2}(I\times\Omega)$ as well as $V^{2}$ in $V$ and $V$ in $H$ a very weak solution is a solution by transposition. Now let $p\in\mathcal{C}^{\infty}(\bar{I},V^{2})$ and set $\phi=\partial_{tt}p-(1/\rho)\partial_{x}(\kappa\partial_{x}p)\in\mathcal{C}^{\infty}(\bar{I},H)$, $p^{0}=p(T)\in V^{2}$ and $p^{1}=\partial_{t}p(T)\in V^{2}$. Thus $p$ is the solution of (3.13). Then the density of $\mathcal{C}^{\infty}(\bar{I},V^{2})$ in $L^{2}(I,V^{2})\cap H^{2}(I,H)$ implies that a solution by transposition is a very weak solution. ∎

Remark 2**.**

For $(u,y^{0},y^{1})\in L^{2}(I,H)\times V\times H$, the weaker solution coincides with the weak one.

3.2. Existence and regularity of the state

In this section we study the existence, uniqueness and regularity of solution of the state equation for measure valued source terms. We will carry out the analysis for both control spaces. We use the distinct properties of each space in order to show improved regularity of the state.

3.2.1. The control space $\mathcal{M}(\Omega,L^{2}(I))$

The space $\mathcal{M}(\Omega,L^{2}(I))$ is not so broad as ${L_{w^{*}}^{2}(I,\mathcal{M}(\Omega))}$ and contains no moving point sources but contains the standing $\delta$-sources (1.2). Therefore, we expect that the state has better regularity properties in this case and prove that $y\in\mathcal{C}(\bar{I},V)\cap\mathcal{C}^{1}(\bar{I},H)$. The proof will be based on a priori bound and a density argument. First we state the following density result.

Lemma 3.7**.**

Let $u\in\mathcal{M}(\Omega,L^{2}(I))$. Then there exists a sequence $\{u_{n}\}\subset{C_{c}^{\infty}(\Omega,L^{2}(I))}$ such that

[TABLE]

Proof.

We denote by $X$ the locally convex space $\mathcal{M}(\Omega,L^{2}(I))$ endowed with its weak-star topology and define the absolutely convex set

[TABLE]

Assume that (3.14) is wrong. Then there exists $u_{0}\in\mathcal{M}(\Omega,L^{2}(I))$ with $\|u_{0}\|_{\mathcal{M}(\Omega,L^{2}(I))}=1$ such that $u_{0}\not\in\bar{E}$ where $\bar{E}$ is the closure of $E$ in $X$. Owing to the corollary of a theorem on the separation of convex sets [28, Ch. III, Theorem 6] there exists $v\in\mathcal{C}_{0}(\Omega,L^{2}(I))$ such that

[TABLE]

On the other hand, $\mathcal{C}_{c}^{\infty}(\Omega,L^{2}(I))$ is dense in $L^{1}(\Omega,L^{2}(I))$ thus

[TABLE]

that contradicts (3.15). Thus $u_{0}\in\bar{E}$. Since $\mathcal{C}_{0}(\Omega,L^{2}(I))$ is separable, the weak-star topology on $E$ is metrizable. Therefore the closure of $E$ is equal to its sequential closure, see [5, Theorem 3.28, Corollary 3.30]. ∎

Note that clearly $C_{c}^{\infty}(\Omega,L^{2}(I))\subset L^{2}(I,V)$.

Preliminarily we prove the following crucial a priori bound.

Lemma 3.8**.**

Let $(u,y^{0},y^{1})\in L^{1}(\Omega,L^{2}(I))\times V\times H$ and $y$ be the corresponding strong solution of problem (2.1). Then $y$ satisfies the following a priori bound

[TABLE]

Proof.

We first remind the energy equality for problem (2.1)

[TABLE]

After setting

[TABLE]

and c_{0}:=\max\big{(}\|\rho\|_{L^{\infty}(\Omega)},\,\|\kappa\|_{L^{\infty}(\Omega)}\big{)}, the energy equality implies

[TABLE]

We also multiply the equation in (2.1) by $-2\kappa\partial_{x}y$ and integrate over $I$. Integration by parts in $t$ yields the equality

[TABLE]

We define a function $P:=\rho\kappa\|\partial_{t}y\|_{L^{2}(I)}^{2}+\|\kappa\partial_{x}y\|_{L^{2}(I)}^{2}$ on $\Omega$. Since the left-hand side of (3.18) equals \partial_{x}P-\big{(}\partial_{x}(\rho\kappa)\big{)}\|\partial_{t}y\|_{L^{2}(I)}^{2}, taking the modulus and integrating over any $(a,b)\subset\Omega$ we derive

[TABLE]

Let $x_{0}\in\bar{\Omega}$ be such that $\|P\|_{\mathcal{C}(\bar{\Omega})}=P(x_{0})$ hold and let now $[a,b]\ni x_{0}$. Then the mean value theorem for integrals implies

[TABLE]

By the above definitions we clearly have

[TABLE]

Inserting (3.19) into (3.20) and using (3.21), we obtain

[TABLE]

Owing to (3.17) we can write

[TABLE]

Using this in (3.22) and choosing a small enough $(a,b)$ such that

[TABLE]

we derive

[TABLE]

Inserting the last bound in (3.23), we also get

[TABLE]

Finally, this yields bound (3.16). ∎

Remark 3**.**

Lemma 3.8 remains valid for $\rho,\kappa\in W^{1,1}(\Omega)$. Owing to the absolute continuity of the Lebesgue integral we have $\|\partial_{x}(\rho\kappa)\|_{L^{1}(a,b)}\leq\mu(b-a)$, where $\lim_{\theta\to+0}\mu(\theta)=0$, thus one can replace $(b-a)^{1/2}\|\rho\kappa\|_{H^{1}(\Omega)}$ by $\mu(b-a)$ in (3.22) and below in the proof.

Theorem 3.9**.**

Let $(u,y^{0},y^{1})\in\mathcal{M}(\Omega,L^{2}(I))\times V\times H$. Then there exists a unique weak solution $y$ and it satisfies the bound

[TABLE]

Proof.

1. Let first $u=0$. According to Proposition 2 were exists a unique weak solution $y$ of (2.1) for any $\mathbf{y}\in V\times H$ and it satisfies

[TABLE]

2. Now it suffices to consider the case $y^{0}=y^{1}=0$. Let first $u\in\mathcal{M}(\Omega,H^{1}(I))\hookrightarrow H^{1}(I,V^{\ast})$ since $\partial_{t}u\in\mathcal{M}(\Omega,L^{2}(I))\hookrightarrow L^{2}(I,V^{\ast})$. Then according to Proposition 2 there exists a unique weak solution $y\in\mathcal{C}(\bar{I},V)\cap\mathcal{C}^{1}(\bar{I},H)$ of (2.1) and it satisfies bound (3.4). Moreover, it is also a weaker solution.

So it remains to prove the bound

[TABLE]

To this end, according to Lemma 3.7 we approximate $u$ by functions $\{u_{n}\}\subset L^{2}(I,V)$ satisfying (3.14). The strong solution $y_{n}$ of (2.1) corresponding to $u=u_{n}$ satisfies the bound like (3.16) and in particular

[TABLE]

Therefore there exists a subsequence of $\{y_{n}\}$ (not relabeled) and $\tilde{y}\in L^{\infty}(I,V)\cap W^{1,\infty}(I,H)$ such that $y_{n}$ converges to $\tilde{y}$ in the weak-star sense of $L^{\infty}(I,V)\cap W^{1,\infty}(I,H)$. This is sufficient to pass to the limit in the last bound and in (3.6) for $y=y_{n}$, $u=u_{n}$ and $v(T)=0$, see Remark 1. Thus $\tilde{y}$ both satisfies the bound

[TABLE]

and is a weaker solution of (2.1). Due to its uniqueness there holds $\tilde{y}=y$, and bound (3.26) is proved.

3. Let now $u\in\mathcal{M}(\Omega,L^{2}(I))$ and $y$ be the corresponding weaker solution of (2.1), see Proposition 3. The space $\mathcal{M}(\Omega,H^{1}(I))$ is dense in $\mathcal{M}(\Omega,L^{2}(I))$, cf. [34, Proposition 2.1]; this also holds since $\mathcal{M}(\Omega,L^{2}(I))$ is the projective closure of the tensor product between $\mathcal{M}(\Omega)$ and $L^{2}(I)$, see [46], and $H^{1}(I)$ is dense in $L^{2}(I)$. Thus there exists a sequence $\{u_{n}\}\subset\mathcal{M}(\Omega,H^{1}(I))$ such that $u_{n}\to u$ in $\mathcal{M}(\Omega,L^{2}(I))$ as $n\to\infty$. Let $y_{n}\in\mathcal{C}(\bar{I},V)\cap\mathcal{C}^{1}(\bar{I},H)$ be the above weak solution of (2.1) corresponding to $u=u_{n}$. Since $\{u_{n}\}$ is a Cauchy sequence in $\mathcal{M}(\Omega,L^{2}(I))$, $\{y_{n}\}$ is a Cauchy sequence in $\mathcal{C}(\bar{I},V)\cap\mathcal{C}^{1}(\bar{I},H)$ too due to bound (3.26) for $u=u_{n}$. Thus $y_{n}\to\hat{y}$ in $\mathcal{C}(\bar{I},V)\cap\mathcal{C}^{1}(\bar{I},H)$ and

[TABLE]

Then we pass to the limit in (3.1) for $y=y_{n}$, $u=u_{n}$ and $v(T)=0$ and see that $\hat{y}$ is a weak solution of (2.1). Due to uniqueness of the weaker solution we get $\hat{y}=y$, and the proof is complete. ∎

3.2.2. Some function spaces and embeddings

We set

[TABLE]

and introduce the interpolation spaces

[TABLE]

for non-integer $\lambda\in(-1,3)$ using the real $K_{\lambda,q}$-interpolation method of Banach spaces for $q=\infty$, see [3]. Recall that the value $q=\infty$ leads to the broadest intermediate spaces. Their explicit description in terms of the Nikolskii spaces or their subspaces is known, see [42, 48, 49]. In particular,

[TABLE]

where $\|w\|_{\tilde{H}^{1/2,2}(\Omega)}=\|\textsl{o}w\|_{H^{1/2,2}(\tilde{\Omega})}$ and $\textsl{o}w$ is the odd extension of $w$ with respect to $x=0$ and $L$ from $\Omega$ to $\tilde{\Omega}:=(-L,2L)$. Hereafter equalities of Banach spaces are understood up to the equivalence of their norms. It is well known that the space $\tilde{H}^{1/2,2}(\Omega)$ contains discontinuous but piecewise $\mathcal{C}^{1}$-functions.

In addition, let $D_{x}$ be the distributional derivative and, for a Banach space $B(\Omega)\subset L^{1}(\Omega)$, $B_{\perp}(\Omega)$ denote the subspace of $W\in B(\Omega)$ with the mean value

[TABLE]

Define the space $H^{-1/2,2}(\Omega)$ of distributions $w=D_{x}W$ with $W\in H_{\perp}^{1/2,2}(\Omega)$ equipped with the norm $\|w\|_{H^{-1/2,2}(\Omega)}=\|W\|_{H^{1/2,2}(\Omega)}$. Then $H^{(-1/2)}=H^{-1/2,2}(\Omega)$, see Item 3 of the proof of Lemma 3.10 below. (Actually a quite similar result is valid for $H^{(\lambda)}$ for any $-1<\lambda<0$.) Note that, in particular, the Dirac delta-function \delta_{a}(x)=D_{x}\big{(}H(x-a)-(1-a/L)\big{)}\in H^{(-1/2)} for any $a\in\Omega$, where $H(\xi)=0$ for $\xi<0$ and $H(\xi)=1$ for $\xi>0$ is the Heaviside function.

Let $Q=\Omega\times I$ and $\Delta_{h}W(x)=W(x+h)-W(x)$ be the forward difference in $x$. Define the spaces $H^{1/2,0;2}(Q)$ and $SHW^{1/2,1;2}(Q)$ of functions $W\in L^{2}(Q)$ such that respectively $|W|_{H^{1/2,0;2}(Q)}:=\sup_{0<h<L}h^{-1/2}\|\Delta_{h}W\|_{L^{2}((0,L-h)\times I)}<\infty$ and $\partial_{t}W\in H^{1/2,0;2}(Q)$ equipped with the norms

[TABLE]

Here $H^{1/2,0;2}(Q)$ is a particular anisotropic Nikolskii space (of the order $1/2$ in $x$ only) and $SHW^{1/2,1;2}(Q)$ is a particular space of functions having the dominating mixed smoothness (of the order $1/2$ in $x$ in the Nikolskii sense and $1$ in $t$ in the Sobolev sense). Note that $SHW^{1/2,1;2}(Q)\hookrightarrow H^{1/2,0;2}(Q)$.

For a Banach space $B(Q)\subset L^{1}(Q)$, let $B_{\perp}(Q)$ be the subspace of $W\in B(Q)$ such that $\langle W(\cdot,t)\rangle_{\Omega}=0$ on $I$. Define the spaces $H^{-1/2,0;2}(Q)$ and $SHW^{-1/2,1;2}(Q)$ of distributions $w=D_{x}W$ with respectively $W\in H_{\perp}^{1/2,0;2}(Q)$ and $W\in SHW_{\perp}^{1/2,1;2}(Q)$ equipped with the norms

[TABLE]

Note that all the spaces defined above and below in this subsection are Banach ones.

The next technical lemma plays an essential role below.

Lemma 3.10**.**

The following equalities and embeddings hold

[TABLE]

Proof.

1. Define the anisotropic Sobolev spaces $W^{1,0;2}(Q)=\{W\in L^{2}(Q)|\,\partial_{x}W\in L^{2}(Q)\}$ and $W^{0,1;2}(Q)=\{W\in L^{2}(Q)|\,\partial_{t}W\in L^{2}(Q)\}$ equipped with the norms

[TABLE]

The following equalities hold

[TABLE]

for example, see [49, Ch. 1.2]. Recall that the corresponding $\hookleftarrow$-embeddings are proved by the classical techniques of approximation by the Steklov averages and the opposite $\hookrightarrow$-embeddings are rather simple. Moreover, equality (3.31) involving three spaces implies (3.32) since it concerns the closed subspaces of one and the same type for all of these three spaces. The same is valid for pairs of embeddings (3.33)-(3.34) and (3.35) below.

The elements in $L^{2}(I,V^{*})$ and $L^{2}(I,H)=L^{2}(Q)$ can be uniquely identified as distributions $w=D_{x}W$ such that respectively $W\in L_{\perp}^{2}(Q)$ with $\|w\|_{L^{2}(I,V^{*})}=\|W\|_{L^{2}(Q)}$ and $W\in W_{\perp}^{1,0;2}(Q)$ with $\|w\|_{L^{2}(I,H)}=\|\partial_{x}W\|_{L^{2}(Q)}$, where $\|\partial_{x}W\|_{L^{2}(Q)}$ is an equivalent norm in $W_{\perp}^{1,0;2}(Q)$ (in the latter case, of course, $D_{x}W=\partial_{x}W$). In particular, for $w\in L^{2}(Q)$, clearly W(x,t)=\int_{0}^{x}w(\xi,t)\,d\xi-\big{\langle}\int_{0}^{x}w(\xi,t)\,d\xi\big{\rangle}_{\Omega}. Taking into account that one and the same operator establishes the one-to-one correspondence between respectively three spaces involved in equalities (3.32) and (3.27), the latter one is valid too.

2. Define the space $SW^{1,1;2}(Q)=\{W\in W^{1,2}(Q)|\,\partial_{x}\partial_{t}W\in L^{2}(Q)\}$ equipped with the norm $\|W\|_{SW^{1,1;2}(Q)}=\|W\|_{W^{1,2}(Q)}+\|\partial_{x}\partial_{t}W\|_{L^{2}(Q)}$. The following equalities

[TABLE]

can be proved similarly to (3.31)-(3.32).

The elements in $H^{1}(I,V^{*})$ and $H^{1}(I,H)=W^{0,1;2}(Q)$ can be uniquely identified as the distributions $w=D_{x}W$ such that respectively $W\in W_{\perp}^{0,1;2}(Q)$, with the equivalent norms $\|w\|_{H^{1}(I,V^{*})}$ and $\|W\|_{W^{0,1;2}(Q)}$, and $W\in SW_{\perp}^{1,1;2}(Q)$, with the equivalent norms $\|w\|_{H^{1}(I,H)}$ and $\|W\|_{SW^{1,1;2}(Q)}$. Thus equality (3.34) implies (3.28).

3. The following equalities hold

[TABLE]

which are simpler 1D versions of (3.31)-(3.32), for example, see [3, 48] and [49, Ch. 1.2]. The second equality implies the above mentioned one $H^{(-1/2)}=H^{-1/2,2}(\Omega)$.

Let $NBV(\bar{\Omega})$ be the space of normalized functions of bounded variation on $\bar{\Omega}$ that are continuous from the right at $x=0$ and continuous from the left at any $x\in(0,L]$; we equip it with the norm $\|W\|_{NBV(\bar{\Omega})}=\sup_{\bar{\Omega}}W+\operatorname{var}_{\bar{\Omega}}W$. Any $w\in\mathcal{M}(\Omega)$ can be represented as $w=D_{x}W$ with $W\in NBV(\bar{\Omega})$ and $\|w\|_{\mathcal{M}(\Omega)}=\operatorname{var}_{\bar{\Omega}}W$, for example, see [10, Ch. 2]. The representation is clearly unique for $W\in NBV_{\perp}(\Omega)$; in this subspace $\operatorname{var}_{\bar{\Omega}}W$ serves as an equivalent norm.

Notice that the following inequalities hold

[TABLE]

for any $W\in NBV(\bar{\Omega})$; the definition of the Riemann integral implies the latter one. Then for any $W\in NBV_{\perp}(\bar{\Omega})$ we get the inequalities

[TABLE]

(they remain valid for $W\in NBV(\bar{\Omega})$ with $\operatorname{var}_{\bar{\Omega}}W$ replaced by $\|W\|_{NBV(\bar{\Omega})}$). Thus

[TABLE]

Let $w\in L_{w^{*}}^{2}(I,\mathcal{M}(\Omega))$. Then $w\in L_{w}^{2}(I,V^{*})=L^{2}(I,V^{*})$, where the equality is valid due to the classical Pettis theorem [20, Theorem 8.15.2] (since $V^{*}$ is separable), and $w=D_{x}W$ with $W\in L_{\perp}^{2}(Q)$ and $\|W\|_{L^{2}(Q)}\leq c\|w\|_{L_{w^{*}}^{2}(I,\mathcal{M}(\Omega))}$.

Moreover, we have $W(t)\in NBV(\bar{\Omega})$ for a.e. $t\in I$. By applying (3.37) to $W(\cdot,t)$, omitting $\sup_{0<h<L}$ on the left, integrating the squared result over $I$ and taking back $\sup_{0<h<L}$ on the left, we obtain

[TABLE]

that completes the proof of embedding (3.29).

4. Let $w\in\mathcal{C}^{1}(\bar{I},\mathcal{M}(\Omega))$. Then $w=D_{x}W$ and $\partial_{t}w=D_{x}Z$ with $W,Z\in H_{\perp}^{1/2,0;2}(Q)$ and

[TABLE]

according to embedding (3.29).

Moreover, define the forward difference quotients in time $\Delta_{\tau}^{(1)}w(t)=(w(t+\tau)-w(t))/\tau$ for $0\leq t<t+\tau\leq T$. Then for the same $t$ and $\tau$ owing to the first inequality (3.36) we get

[TABLE]

Therefore

[TABLE]

Consequently there exists the derivative $\partial_{t}W=Z\in L^{2}(Q)$, and inequality (3.38) implies embedding (3.30). ∎

We also set $V^{0}(Q)=L^{2}(Q)$ and define the anisotropic Sobolev subspaces

[TABLE]

for $\ell=1,2$ and the anisotropic Nikolskii subspaces

[TABLE]

equipped with the norm $\|w\|_{\tilde{H}^{\ell+1/2,0;2}(Q)}=\|\partial_{x}^{\ell}\textsl{o}w\|_{H^{1/2,0;2}(\tilde{Q})}$ for $\ell=0,1$, where $\tilde{Q}=\tilde{\Omega}\times I$. Then the following equality holds

[TABLE]

which is similar to equality (3.31).

4. Analysis of the control problem

According to Theorem 3.9 and Proposition 3 the state equation (2.1) is uniquely solvable for any $u$ in either $\mathcal{M}(\Omega,L^{2}(I))$ or ${L_{w^{*}}^{2}(I,\mathcal{M}(\Omega))}$ and the solution $y$ depends continuously on the data. Therefore, we can introduce the linear and bounded operator $\hat{S}\colon(u,y^{0},y^{1})\mapsto(y,y(T),{\rho}\partial_{t}y(T))$. The control-to-state mapping

[TABLE]

is given by $Su=\hat{S}(u,0,0)+\hat{S}(0,y^{0},y^{1})$ for fixed $y^{0}$ and $y^{1}$ and it is an affine and bounded operator. So we can rewrite the original control problem ($\mathcal{P}$) in its reduced form

[TABLE]

Proposition 5**.**

Problem ($\mathcal{P}$) has a unique solution $\bar{u}\in\mathcal{M}_{T}$.

Proof.

The control-to-state operator $S$ is weak-star-to-strong sequential continuous, i.e., if $\{u_{n}\}\subset\mathcal{M}_{T}$ and $u_{n}\rightharpoonup^{\ast}u$ in $\mathcal{M}_{T}$, then $Su_{n}\rightarrow Su$ in $\mathcal{Y}$. The proof of this continuity property is similar to [34, Lemma 6.1] in the case of solutions by transposition resp. very weak solutions. The strong continuity follows from the compact embeddings and well known Aubin-Lions-Lemma. Then the direct method of calculus of variations combined with the sequential Banach-Alaoglu theorem ($\mathcal{C}_{T}$ is separable) can be applied to show existence of an optimal control. Additionally the control is unique since the control-to-state operator $S$ is injective and the data tracking functional is strictly convex. ∎

Owing to Proposition 3 the optimal control $\bar{u}\in\mathcal{M}_{T}$ satisfies the inequalities

[TABLE]

and thus

[TABLE]

Hereafter $C>0$ depends on the norms of data.

Next we discuss first order optimality conditions. We introduce the adjoint control-to-solution operator $S^{\star}\colon\mathcal{Y}\rightarrow C(\bar{I},V)\hookrightarrow\mathcal{C}_{T}$, $(\phi,p^{1},p^{0})\mapsto p$ where $p$ is a weak solution of (3.13). This operator is well defined and bounded according to Proposition 2.

We also need the operator $A^{-1}\colon V^{\ast}\rightarrow V$, $f\mapsto w$ where $w\in V$ is the unique solution of

[TABLE]

The next result provides the necessary and sufficient optimality condition for the optimal pair $(\bar{p},\bar{u})$.

Proposition 6**.**

An element $\bar{u}\in\mathcal{M}_{T}$ is an optimal control of ($\mathcal{P}$) if and only if

[TABLE]

or equivalently

[TABLE]

where \bar{p}=S^{*}\big{(}\bar{y}-z_{1},-(\bar{y}(T)-z_{2}),A^{-1}(\rho\partial_{t}\bar{y}-z_{3})\big{)} with $(\bar{y},\bar{y}(T),{\rho}\partial_{t}\bar{y}(T))=\hat{S}(\bar{u},y^{0},y^{1})$.

Proof.

For $\mathcal{M}_{T}=\mathcal{M}(\Omega,L^{2}(I))$ the proof of [34, Theorem 7.1] remains valid; for $\mathcal{M}_{T}={L_{w^{*}}^{2}(I,\mathcal{M}(\Omega))}$ it is similar to [12, Theorem 3.2]. ∎

To discuss further the properties of the optimal control $\bar{u}$, we introduce the Jordan decomposition of a signed measure $\mu\in\mathcal{M}(\Omega)$, see [6]. There exists unique elements $\mu^{\pm}\in\mathcal{M}(\Omega)^{+}$ such that $\mu=\mu^{+}-\mu^{-}$. Moreover, we recall the polar decomposition of a vector measure $\mu\in\mathcal{M}(\Omega,L^{2}(I))$: $\mathrm{d}\mu=\mu^{\prime}\mathrm{d}|\mu|$, where $\mu^{\prime}$ is the Radon-Nikodym-derivative of $\mu$ with respect to $|\mu|$.

The subgradient condition in Proposition 6 implies the following conditions.

Proposition 7**.**

Let $\bar{u}\in\mathcal{M}_{T}$ be the optimal control of ($\mathcal{P}$) and $\bar{p}\in\mathcal{C}_{T}$ be the corresponding adjoint state. Then there holds $\|\bar{p}\|_{\mathcal{C}_{T}}\leq\alpha$.

In the cases $\mathcal{M}_{T}={L_{w^{*}}^{2}(I,\mathcal{M}(\Omega))}$ and $\mathcal{M}_{T}=\mathcal{M}(\Omega,L^{2}(I))$ there respectively hold

[TABLE]

and

[TABLE]

Proof.

A detailed discussion of the proof of these results can be found in [12, 33]. ∎

The regularity of the adjoint state $\bar{p}$ is now applied to show improved regularity of the optimal control $\bar{u}$.

Theorem 4.1**.**

Let $\mathcal{M}_{T}=\mathcal{M}(\Omega,L^{2}(I))$, $\mathbf{z}\in\mathcal{Y}^{1}:=L^{2}(I,V)\times V\times H$, $\mathbf{y}\in V\times H$ and $\bar{u}$ be the optimal control of ($\mathcal{P}$). Then $\bar{u}\in\mathcal{C}^{1}(\bar{I},\mathcal{M}(\Omega))$ and the following bound holds

[TABLE]

Proof.

There holds $\bar{y}\in{\mathcal{C}(\bar{I},V)\cap\mathcal{C}^{1}(\bar{I},H)}$ according to Theorem 3.9. Thus, the optimal adjoint state has the following regularity $\bar{p}\in\mathcal{C}(\bar{I},V^{2})\cap\mathcal{C}^{1}(\bar{I},V)$ by Proposition 2. We have $\bar{u}=-{\alpha^{-1}}\bar{p}\,|\bar{u}|$ according to (4.6). Moreover, we define the function

[TABLE]

and show that it serves the time derivative of $\bar{u}$. For any $t_{0},t\in\bar{I}$ and $t_{0}\neq t$, we define the difference quotient \bar{u}{(t_{0},t)}=\big{(}\bar{u}(t)-\bar{u}(t_{0})\big{)}/(t-t_{0}). Then we consider

[TABLE]

as $t\to t_{0}$ since $\bar{p}\in\mathcal{C}^{1}(\bar{I},V)$. Next, quite similarly we get

[TABLE]

as $t\to t_{0}$. Consequently $\partial_{t}\bar{u}=w\in\mathcal{C}(\bar{I},\mathcal{M}(\Omega))$. Finally, we bound $\partial_{t}\bar{u}$ as follows

[TABLE]

owing to Proposition 2 and Theorem 3.9. Utilizing bound (4.2) for $\bar{u}$, we complete the proof. ∎

5. Discretization of the state equation

We introduce the uniform grid $t_{m}=m\tau$ in time with the step $\tau=T/M$ and a non-uniform grid $0=x_{0}<x_{1}<\ldots<x_{N}=L$ in space with the steps $h_{j}=x_{j}-x_{j-1}$, where $M\geq 2$ and $N\geq 2$. Let also $h=\max_{j=1,\ldots,N}h_{j}$, $h_{\rm\min}=\min_{j=1,\ldots,N}h_{j}$ and $\vartheta=(\tau,h)$. We assume that the space grid is quasi-uniform, i.e., $h\leq c_{1}h_{\rm\min}$. Hereafter $c,c_{1},C$, etc., are grid-independent.

Let $V_{\tau}\subset H^{1}(I)$ and $V_{h}\subset V$ be the spaces of piecewise linear finite elements with respect to the introduced grids on $\bar{I}$ and $\bar{\Omega}$.

We approximate the state variable $y$ by $y_{\vartheta}\in V_{\vartheta}:=V_{\tau}\otimes V_{h}\subset H^{1}(I,V)$ and additionally $\partial_{t}y(T)$ by $y_{Th}^{1}\in V_{h}$. For $(u,y^{0},y^{1})\in\mathcal{M}_{T}\times H\times V^{*}$ the discrete state equation has the following form

[TABLE]

involving the indefinite symmetric bilinear form

[TABLE]

with the grid independent parameter $\sigma$, cf. (3.1). This definition follows [51] but notice carefully that normally $y_{\vartheta}$ is uniquely defined by (5.1) with $v(T)=0$ and (5.2). To treat general $v$, we need $y_{Th}^{1}$.

Remark 4**.**

The second term on the right hand-side of (5.3) regularizes the Galerkin (i.e. projection) method with respect to bilinear form (3.2). It is included to ensure unconditional stability for suitable values of $\sigma$. Moreover, the term

[TABLE]

is the error term of the compound trapezoidal rule applied for the calculation of the temporal integral in $(\kappa\partial_{x}y,\partial_{x}v)_{L^{2}(I\times\Omega)}$. So that, in particular, for $\sigma=0$ in (5.3) this temporal integral is calculated using this rule whereas for $\sigma=1/6$ it is not approximated.

Next we recall the inverse inequality

[TABLE]

where the least constant satisfies $c_{1}h^{-1}\leq\alpha_{h}\leq c_{2}h^{-1}$ for the quasi-uniform grid. For $\sigma\leq 1/4$ we need to state conditions linking the temporal and spatial grids to ensure stability of the numerical method.

Assumption 1**.**

In what follows, let

[TABLE]

Remark 5**.**

The parameters $\varepsilon_{0}$ and $\varepsilon_{1}$ can be chosen arbitrarily small but then constants in the stability and error estimates for our FEM can tend to infinity.

Remark 6**.**

As we see below in Section 11, the method is related to well known time-stepping methods, in particular, to the explicit Leap-Frog-method for $\sigma=0$. Then conditions (5.5) and (5.6) reduce to a CFL-type one $\tau\alpha_{h}\leq 2\sqrt{1-\varepsilon_{0}^{2}}$. For $\sigma=1/4$ the method is related to the Crank-Nicolson scheme and is unconditionally stable but in a weaker norm than we need to derive our error estimates so that we impose a very weak CFL-type condition $\tau\alpha_{h}\leq 2/\varepsilon_{1}$.

Below in proofs we utilize the auxiliary squared norms

[TABLE]

for $\varphi\in V_{h}$ and $y\in V_{\tau}\otimes V_{h}$. We need to bound them by standard norms.

Lemma 5.1**.**

Under conditions (5.5) and (5.6) the following inequalities hold

[TABLE]

with $\varepsilon_{0}:=1$ for $\sigma\geq 1/4$ and $\varepsilon_{1}:=\sqrt{4\sigma-1}$ for $\sigma>1/4$.

Proof.

For $\sigma\geq 1/4$, the first inequality is obvious; for $\sigma<1/4$ it can be checked by a direct calculation using (5.4). The proof of the second inequality is covered in [51, Corollary 2.1]. ∎

Now we discuss some properties of $y_{Th}^{1}$ and $\partial_{t}y(T)$ that are essential below.

Proposition 8**.**

Let $(y_{\vartheta},{y_{Th}^{1}})\in V_{\vartheta}\times V_{h}$ be the solution of (5.1)-(5.2). Then there holds

[TABLE]

Proof.

This is proved by testing (5.1) with time constant functions $v=\varphi\in V_{h}$. ∎

The non-local in time identity (5.8) is convenient for our error analysis but not for the implementation; for the latter issue see Section 11. Identities similar to (5.8) also hold on the continuous level.

Proposition 9**.**

1. (1)

Let $y\in\mathcal{C}(\bar{I},V)\cap\mathcal{C}^{1}(\bar{I},H)$ be the weak solution of (2.1) for $\mathcal{M}_{T}=\mathcal{M}(\Omega,L^{2}(I))$. Then there holds

[TABLE] 2. (2)

Let $y\in\mathcal{C}(\bar{I},H)\cap\mathcal{C}^{1}(\bar{I},V^{\ast})$ be the weaker (very weak) solution of (2.1) for $\mathcal{M}_{T}={L_{w^{*}}^{2}(I,\mathcal{M}(\Omega))}$. Then there holds

[TABLE]

Proof.

For $\mathcal{M}_{T}=\mathcal{M}(\Omega,L^{2}(I))$ identity (5.9) is proved by testing (3.1) with time constant function $v=\varphi\in V$. For $\mathcal{M}_{T}={L_{w^{*}}^{2}(I,\mathcal{M}(\Omega))}$ we test (3.7) with any $\varphi\in V^{2}$ and get

[TABLE]

According to Proposition 3 we have $\mathcal{I}_{t}y\in\mathcal{C}(\bar{I},V)$. Thus there holds

[TABLE]

The density of $V^{2}$ in $V$ implies (5.10). ∎

For our analysis, we need some projection and interpolation operators. We introduce the standard projectors $\pi_{h}^{0}$: ${H_{\rho}}\to V_{h}$ and $\pi_{h}^{1}$: ${\mathcal{V}_{\kappa}}\to V_{h}$ defined by

[TABLE]

Clearly $\|\pi_{h}^{0}w\|_{H_{\rho}}\leq\|w\|_{H_{\rho}}$ and $\|\pi_{h}^{1}w\|_{\mathcal{V}_{\kappa}}\leq\|w\|_{\mathcal{V}_{\kappa}}$. Identity (5.2) means that $y_{\vartheta}(0)=\pi_{h}^{0}y^{0}$.

Moreover the following property holds

[TABLE]

Following [51], we also introduce the regularized ${H_{\rho}}$ projector $\pi_{h,\sigma_{0}}^{0}$: $V\to V_{h}$ defined by

[TABLE]

with the grid independent parameter $\sigma_{0}\geq\sigma-1/4$. Clearly $\pi_{h,\sigma_{0}}=\pi_{h}^{0}$ for $\sigma_{0}=0$.

Let $i_{\tau}$: $\mathcal{C}(\bar{I})\to V_{\tau}$ be the interpolation operator such that $i_{\tau}w(t_{m})=w(t_{m})$ for all $m=0,\ldots,M$.

Next we define the operator $A_{h}^{-1}\colon V^{\ast}\rightarrow V_{h}$, $f\mapsto w_{h}$ where $w_{h}\in V_{h}$ is the unique solution of

[TABLE]

Clearly $A_{h}^{-1}=\pi_{h}^{1}A^{-1}$, see (4.3) with $w=A^{-1}f$, and the norm in $\mathcal{V}_{\kappa}^{\ast}$ and its discrete counterpart can be written as

[TABLE]

Moreover, we set $r_{h}A^{-1}:=A^{-1}-A_{h}^{-1}=A^{-1}-\pi_{h}^{1}A^{-1}$. First we note that

[TABLE]

Then by the standard FEM error analysis [8] and operator interpolation theory we have

[TABLE]

[TABLE]

6. Stability and error estimates for the discrete state equation

In this section we present error estimates for the state equation. We begin with an auxiliary result.

Lemma 6.1**.**

For $\sigma_{0}\geq\sigma-1/4\geq 0$, the following estimate holds

[TABLE]

Proof.

We recall the well known estimates

[TABLE]

which are valid using the inverse inequality (5.4). We also remind inequality (5.7) and notice also that for $\sigma_{0}\geq 0$ the following additional inequality holds

[TABLE]

Let $w\in V$ and $\varphi\in V_{h}$. We apply identities (5.11) and (5.14) and get

[TABLE]

Now we set $\varphi=\pi_{h,\sigma_{0}}^{0}{w}-\pi_{h}^{0}{w}$, and from the former equality and estimate (6.2) for $\lambda=1$ as well as the latter equality and estimate (6.3) for $\lambda=2$ we obtain the estimate

[TABLE]

By using the $K_{\lambda,\infty}$-method, we complete the proof. ∎

Now we get a stability bound and error estimates in $\mathcal{C}(\bar{I},H)\times\mathcal{V}_{\kappa,h}^{*}$ for the discrete state equation.

Proposition 10**.**

Let $y$ and $(y_{\vartheta},y_{Th}^{1})$ be the solutions to the state equation (2.1) and the discrete state equation (5.1)-(5.2).

1. (1)

For $(u,y^{0},y^{1})\in L^{2}(I,V^{\ast})\times V\times V^{*}$, the following stability bound holds:

[TABLE] 2. (2)

For $(u,y^{0},y^{1})\in{H^{-1/2,0;2}(Q)}\times V\times H^{(-1/2)}$, the following error estimate holds:

[TABLE] 3. (3)

For $(u,y^{0},y^{1})\in{SHW^{-1/2,1;2}(Q)}\times V\times H$, the higher order error estimate holds:

[TABLE]

Proof.

1. According to [51, Theorem 2.1 (1)], the bound

[TABLE]

is valid for any $y_{\vartheta}(0)\in V_{h}$. We have $y_{\vartheta}(0)=\pi_{h}^{0}y^{0}$. In the case $\sigma\leq 1/4$, there clearly holds

[TABLE]

For $\sigma>1/4$, we alternatively get using (6.1) for $\lambda=1$

[TABLE]

for any ${\sigma_{0}\geq\sigma-1/4}$.

We proceed with the bound for $y_{Th}^{1}$. Identity (5.8) and bound (6.8) together with the generalized Minkowski inequality

imply

[TABLE]

Finally we derive bound (6.5).

2. Let $\tilde{y}_{\vartheta}$ be the solution of equation (5.1) for $\tilde{y}_{\vartheta}(0)=\pi_{h,\sigma_{0}}^{0}y^{0}$. Bound (6.8) together with (6.1) for $\lambda=1$, the bound in Proposition 3 and the stability of $\pi_{h}^{1}$ in $V$ imply

[TABLE]

Owing to [51, Theorem 4.1] the following error estimate holds

[TABLE]

Using the $K_{1/2,\infty}$-method and equality (3.27) we get the intermediate error estimate

[TABLE]

In the case $\sigma\leq 1/4$ we can choose $\sigma_{0}=0$, then $y_{\vartheta}(0)=\pi_{h,\sigma_{0}}y^{0}=\pi_{h}^{0}y^{0}$ and $\tilde{y}_{\vartheta}=y_{\vartheta}$. In the case $\sigma\geq 1/4$ we can use the stability bound (6.8) and estimate (6.1) to get

[TABLE]

Then by subtracting (5.8) from (5.10) and applying identity (5.12) we find

[TABLE]

consequently

[TABLE]

Thus we obtain (6.6).

3. Once again we apply [51, Theorem 4.1] and first get the estimate

[TABLE]

Combining it together with (6.11), we derive

[TABLE]

In this proof, we apply this estimate in the case $u=0$ only (but in general case below).

In the remaining case $\mathbf{y}=0$, from [51, Theorem 4.1] we also get the higher order error estimate

[TABLE]

Moreover owing to Proposition 3 and bound (6.5) (both for $\mathbf{y}=0$) we have

[TABLE]

The last bound and estimate (6.14) imply by the $K_{1/2,\infty}$-method and equality (3.28):

[TABLE]

for any $u\in SHW^{-1/2,1;2}(Q)$. Applying inequality (6.12) we complete the proof. ∎

Remark 7**.**

A priori stability bound (6.5) implies the unique solvability of the discrete state equation (5.1)-(5.2).

Remark 8**.**

According to the given proof, for $\tilde{y}_{\vartheta}$ in place of $y_{\vartheta}$ the norms of $\mathbf{y}$ in (6.5) and (6.6) can be weakened down to respectively $\|\mathbf{y}\|_{H\times V^{*}}$ and $\|\mathbf{y}\|_{H^{(1/2)}\times H^{(-1/2)}}$. For $\sigma\leq 1/4$, we have $\tilde{y}_{\vartheta}=y_{\vartheta}$. The same can be shown for $y_{\vartheta}$ also for $\sigma>1/4$ provided that $\tau\alpha_{h}\leq c_{0}$ with any $c_{0}>0$.

7. Discrete control problem

First we introduce the discrete mapping

[TABLE]

and the discrete affine linear control-to-state mapping

[TABLE]

defined by $S_{\vartheta}{u}=\hat{S}_{\vartheta}(u,0,0)+\hat{S}_{\vartheta}(0,y_{0},y_{1})$, with $\rho\times V_{h}=\{\rho\varphi;\varphi\in V_{h}\}$. The mapping $S_{\vartheta}$ is a composition of

[TABLE]

where $\{e_{m,n}^{\vartheta}\}$ is a basis in $V_{\vartheta}$, and $\vec{u}\mapsto(y_{\vartheta},y_{\vartheta}(T),{\rho y_{Th}^{1}})$. The former mapping is bounded due to $e_{m,n}^{\vartheta}\in\mathcal{C}_{T}$ and the latter one is finite dimensional. Thus $S_{\vartheta}$ is a bounded operator. Then we consider the following semi-discrete optimal control problem

[TABLE]

with the squared semi-norm corresponding to the inner product

[TABLE]

Using the similar argument as in the continuous case it can be shown that ($\mathcal{P}_{\vartheta}$) has a solution $\bar{u}_{\vartheta}$ which is not unique in general, and due to the optimality, the stability bound (6.5) and property (5.16) (for $\lambda=-1$) one gets

[TABLE]

cf. (4.1), and consequently

[TABLE]

Theorem 7.1**.**

Let $\mathbf{z}\in\mathcal{Y}$, $\mathbf{y}\in V\times H^{(-1/2)}$ and $\bar{u},\bar{u}_{\vartheta}\in\mathcal{M}_{T}$ be the optimal controls of respectively problems ($\mathcal{P}$) and ($\mathcal{P}_{\vartheta}$). Then there holds

[TABLE]

Proof.

Owing to (7.1) there exists a sequence $\{\vartheta_{n}\}$, $\vartheta_{n}\rightarrow 0$, and $u\in\mathcal{M}_{T}$ such that $\bar{u}_{\vartheta_{n}}\rightharpoonup^{\ast}u$ in $\mathcal{M}_{T}$ as $n\rightarrow\infty$. This implies the limit relation

[TABLE]

To prove it, we write the chain of inequalities

[TABLE]

The first term on the right in the last inequality converges to zero according to the error estimate (6.6). The convergence of the second term follows from the weak-star-to-strong continuity of $S\colon\mathcal{M}_{T}\rightarrow{\mathcal{Y}}$ and the stability of $\pi_{h}^{1}$ in $V$. Finally, property (5.13) for $\tilde{w}=w$ implies the convergence of the last term.

Then (7.2) and the weak-star lower semicontinuity of $\|\cdot\|_{\mathcal{M}_{T}}$ in $\mathcal{M}_{T}$ implies

[TABLE]

Thus, the uniqueness of $\bar{u}$ means that $u=\bar{u}$ and in addition implies the convergence of the whole sequence $\bar{u}_{\vartheta}\rightharpoonup^{\ast}\bar{u}$ in $\mathcal{M}_{T}$ as $\vartheta\rightarrow 0$. Moreover, we have $j_{\vartheta}(\bar{u}_{\vartheta})\rightarrow j(\bar{u})$. This and (7.2) lead to $\|\bar{u}_{\vartheta}\|_{\mathcal{M}_{T}}\rightarrow\|\bar{u}\|_{\mathcal{M}_{T}}$. ∎

For convenience we set $F_{h}({\mathbf{z}})={(1/2)}\|{\mathbf{z}}\|_{{\mathcal{Y}_{h}}}^{2}$. In the following the directional derivative of a functional $g\colon\mathcal{M}_{T}\rightarrow\mathbb{R}$ at $u\in\mathcal{M}_{T}$ in direction $\delta u\in\mathcal{M}_{T}$ is denoted by $Dg(u)\delta u$. In the case $Dg(u)\in\mathcal{M}_{T}^{\ast}$, $g$ is the Gateaux differentiable in $u$. Moreover, we make use of the convex subdifferential of $\|\cdot\|_{\mathcal{M}_{T}}$. Let $\hat{u}\in\mathcal{M}_{T}$ and $p\in\mathcal{C}_{T}$. Then there holds $p\in\partial\|\hat{u}\|_{\mathcal{M}_{T}}$ if and only if

[TABLE]

An element $\bar{u}_{\vartheta}\in\mathcal{M}_{T}$ is an optimal solution of ($\mathcal{P}_{\vartheta}$) if and only if $-D((F_{h}\circ S_{\vartheta})(\bar{u}_{\vartheta}))\in\alpha\partial\|\bar{u}_{\vartheta}\|_{\mathcal{M}_{T}}$. To calculate $D((F_{h}\circ S_{\vartheta})(u))$ for $u\in\mathcal{M}_{T}$, we apply the Lagrange technique and define the Lagrange functional by

[TABLE]

with $(p_{\vartheta},{p_{0h}^{1}})\in V_{\vartheta}\times V_{h}$ (where we base on identities (5.1)-(5.2)). We obviously have

[TABLE]

Thus there holds

[TABLE]

provided that $(p_{\vartheta},{p_{0h}^{1})}\in{V}_{\vartheta}\times V_{h}$ is the solution of the discrete problem

[TABLE]

and

[TABLE]

Therefore the discrete optimality system consists of the discrete state equation

[TABLE]

the discrete adjoint state equation

[TABLE]

and the discrete variational inequality

[TABLE]

8. Stability and error estimates for the discrete adjoint state equation

We define the general discrete adjoint state equation

[TABLE]

Here $y$ is the solution to the state equation (2.1). Clearly identity (8.2) means simply that $p_{\vartheta}(T)=A_{h}^{-1}q_{T}=\pi_{1}^{h}A^{-1}q_{T}$ with $q_{T}:=\rho\partial_{t}y(T)-z_{3}$.

Now we get a stability bound and error estimates in $\mathcal{C}(\bar{I},H)\times\mathcal{V}_{\kappa,h}^{*}$ and $\mathcal{C}_{T}$ for the discrete adjoint state equation.

Proposition 11**.**

Let p=S^{*}\big{(}y-z_{1},-(y(T)-z_{2}),A^{-1}(\rho\partial_{t}y(T)-z_{3})\big{)} and $(p_{\vartheta},p_{0h}^{1})$ be the solution of the corresponding general discrete adjoint state equation (8.1)-(8.2).

1. (1)

If $y\in\mathcal{C}(\bar{I},H)\cap\mathcal{C}^{1}(\bar{I},V^{*})$ and $\mathbf{z}\in\mathcal{Y}$, then the following stability bound holds

[TABLE] 2. (2)

If $u\in L^{2}(I,V^{*})$, $\mathbf{z}\in\mathcal{Y}$ and $\mathbf{y}\in H\times V^{*}$, then the following error estimate holds

[TABLE] 3. (3)

If $u\in{H^{-1/2,0;2}(Q)}$, $\mathbf{z}\in\mathcal{Y}^{1/2}:={\tilde{H}^{1/2,0;2}(Q)}\times H^{(1/2)}\times H^{(-1/2)}$ and $\mathbf{y}\in H^{(1/2)}\times H^{(-1/2)}$, then the following error estimate holds

[TABLE] 4. (4)

If $u\in{SHW^{-1/2,1;2}(Q)}$, $\mathbf{z}\in\mathcal{Y}^{3/2}:={\tilde{H}^{3/2,0;2}(Q)}\times H^{(3/2)}\times H^{(1/2)}$ and $\mathbf{y}\in H^{(3/2)}\times H^{(1/2)}$, then the following higher order error estimate holds

[TABLE]

Proof.

1. According to [51, Theorem 2.1 (2)] the following energy bound hold

[TABLE]

for any $p_{\vartheta}(T)\in V_{h}$. Using (6.2), $A_{h}^{-1}=\pi_{h}^{1}A^{-1}$ and (5.16) we get

[TABLE]

By applying also the counterpart of inequalities (6.9) we derive bound (8.3).

2. The counterpart of the error estimate (6.13) for the adjoint state equation case and bound (8.7) give

[TABLE]

Owing to inequality (6.12) and Proposition 3 we obtain estimate (8.4).

3. Below we need the multiplicative inequalities

[TABLE]

Let $\check{p}_{\vartheta}$ be the auxiliary solution to (8.1) for $\check{p}_{\vartheta}(T)=\pi_{h}^{0}A^{-1}q_{T}$. Owing to inequality (8.8) and the stability bounds [51, Theorem 2.1] we get

[TABLE]

Consequently, for $q_{T}\in H^{(\alpha-2)}$, by (6.2), (5.17) and (5.18) the following chain of inequalities hold

[TABLE]

for $1\leq\alpha\leq 2$. Thus it is enough to prove error estimates (8.5) and (8.6) for $\check{p}_{\vartheta}$ instead of $p_{\vartheta}$.

According to [51, Theorem 5.3 and estimate (5.18)] we have the error estimate

[TABLE]

for $\alpha=1,2$. We emphasize that due to [51, Theorem 4.3 (2) (e)] and (6.1) this estimate holds for $\check{p}_{\vartheta}(T)=\pi_{h}^{0}A^{-1}q_{T}$.

Inequality (8.8), Proposition 2 (applied to the adjoint state problem) and property (5.16) imply the following error estimate for the time interpolation

[TABLE]

for $\alpha=1,2$. Owing to estimates (8.10) and (8.11) and Propositions 3 and 2 we get

[TABLE]

for $\leavevmode\nobreak\ \alpha=1,2$, where $\mathcal{Y}^{(0)}:=\mathcal{Y}$ and $\mathcal{Y}^{(1)}:=L^{2}(I,H)\times V\times H$.

Applying the $K_{1/2,\infty}$-method together with equalities (3.27) and (3.39) for $\ell=0$, we get (8.5) for $\check{p}_{\vartheta}$ in the role of $p_{\vartheta}$.

4. First notice that the multiplicative inequality (8.9), Proposition 2 (2) (applied for the adjoint state problem) and property (5.16) imply another error estimate for the time interpolation

[TABLE]

Then Proposition 2 (1) leads to

[TABLE]

Next we derive the error estimate

[TABLE]

According to [51, Theorem 5.3 and estimate (5.18)] and equality (3.39) for $\ell=1$ together with Propositions 3 and 2 the following three estimates hold

[TABLE]

[TABLE]

[TABLE]

and for $\check{p}_{\vartheta}(T)=\pi_{h}^{0}A^{-1}q_{T}$ (for the same reason as above). Then applying the $K_{1/2,\infty}$-method to the two last estimates and using equality (3.28) we get

[TABLE]

By combining this estimate and (8.15) we obtain (8.14).

Estimates (8.13) and (8.14) imply

[TABLE]

that completes the proof of (8.6) for $\check{p}_{\vartheta}$ in the role of $p_{\vartheta}$. ∎

Remark 9**.**

A priori stability bound (6.5) (taken for $y=0$) implies the unique solvability of the general discrete adjoint state equation (8.1)-(8.2).

9. Error estimates for the state variable

We introduce the discrete adjoint control-to-state operator $S^{\star}_{\vartheta}\colon L^{2}(I\times\Omega)\times V\times H\rightarrow V_{\vartheta}$, $(\phi,p^{1},p^{0})\mapsto p_{\vartheta}$ defined by

[TABLE]

with $p_{\vartheta}(T)=\pi_{h}^{0}p^{0}$. Similarly to bound (8.3) and Remark 9 it is well defined and satisfies

[TABLE]

Let for brevity $W,W_{h}\colon\mathcal{Y}\rightarrow\mathcal{Y}^{\ast}$ be the duality mappings defined by

[TABLE]

for any $(y_{1},y_{2},y_{3})\in\mathcal{Y}$. With this notation, the function

[TABLE]

solves the general discrete adjoint state equation (8.1)-(8.2).

Proposition 12**.**

Let $\mathbf{z}\in\mathcal{Y}$ and $\mathbf{y}\in V\times V^{*}$. Then the following estimate holds

[TABLE]

Proof.

We recall that $\bar{p}=S^{\star}W(S\bar{u}-\mathbf{z})$ and $\bar{p}_{\vartheta}=S^{\star}_{\vartheta}W_{h}(S_{\vartheta}\bar{u}_{\vartheta}-\mathbf{z})$ and test the continuous subgradient condition (4.5) with the discrete optimal control $\bar{u}_{\vartheta}$ and the discrete subgradient condition (7.5) with the continuous optimal control $\bar{u}$. Then we subtract the first inequality from the second one and get

[TABLE]

We define $\hat{p}_{\vartheta}:=S^{\star}_{\vartheta}W_{h}(S\bar{u}-\mathbf{z})$, insert it between $\bar{p}$ and $\bar{p}_{\vartheta}$ and obtain

[TABLE]

For convenience we introduce the variables $(\hat{y}_{\vartheta},\hat{y}_{\vartheta}(T),{\rho\hat{y}_{Th}^{1}})=S_{\vartheta}\bar{u}$ and remark that the state equations for $(\bar{y}_{\vartheta},{\bar{y}_{Th}^{1}})$ and $(\hat{y}_{\vartheta},{\hat{y}_{Th}^{1}})$ have the same initial data. With the help of them we rewrite the second term on the right in (9.2) taking first the difference of the discrete state equations (7.3) and (5.1) (taken for $(\hat{y}_{\vartheta},\hat{y}_{Th}^{1})$) for $v=\hat{p}_{\vartheta}-\bar{p}_{\vartheta}$, next the difference of the discrete adjoint state equations (7.4) and (8.1)-(8.2) (taken for $\hat{p}_{\vartheta}$) for $v=\bar{y}_{\vartheta}-\hat{y}_{\vartheta}$ and $\varphi=\bar{y}_{Th}^{1}-\hat{y}_{Th}^{1}$ and finally using (5.15)

[TABLE]

Further we easily get

[TABLE]

Thus (9.2) implies

[TABLE]

Finally by applying bounds (4.2) and (7.1) we derive (9.1). ∎

This proposition is important since it allows one to derive estimates for $\bar{y}-\bar{y}_{\vartheta}$ with the help of the above error estimates for the discrete state and adjoint state equations.

Theorem 9.1**.**

1. (1)

Let $\mathcal{M}_{T}={L_{w^{*}}^{2}(I,\mathcal{M}(\Omega))}$, $\mathbf{z}\in\mathcal{Y}^{1/2}$ and $\mathbf{y}\in V\times H$. Then the following error estimate holds

[TABLE] 2. (2)

Let $\mathcal{M}_{T}=\mathcal{M}(\Omega,L^{2}(I))$, $\mathbf{z}\in\mathcal{Y}^{3/2}$ and $\mathbf{y}\in H^{(3/2)}\times H^{(1/2)}$. Then the following higher order error estimate holds

[TABLE]

Proof.

1. Let us base on Proposition 12. First, Proposition 10 (4) implies

[TABLE]

Second, Proposition 11 (3) leads to

[TABLE]

Now owing to Proposition 12, embedding (3.29) and bound (4.2) for $\bar{u}$ error estimate (9.3) is proved.

2. First, Proposition 10 (3) implies

[TABLE]

Second, Proposition 11 (4) leads to

[TABLE]

Now owing to Proposition 12, embedding (3.30) and Theorem 4.1 for $\bar{u}$ error estimate (9.4) is proved too. ∎

Remark 10**.**

Note that our error bounds could be better provided that one would improve the last term on the right in (9.1) by increasing the power $1/2$. But this seems a complicated problem.

10. Error estimate for the cost functional

In this section we derive error estimate for the cost functional. We first observe the inequalities

[TABLE]

which can be equivalently rewritten in the form

[TABLE]

Therefore, to bound $|j(\bar{u})-j_{\vartheta}(\bar{u}_{\vartheta})|$ below we apply the following result.

Proposition 13**.**

Let $\mathbf{y}\in V\times H$. Then for any $u\in\mathcal{M}_{T}$

[TABLE]

with $(y,y(T),{\rho}\partial_{t}y(T))=Su$ and the same $p$ and $(p_{\vartheta},{p}_{0h}^{1})$ as in Proposition 11.

Proof.

Let $u\in\mathcal{M}_{T}$. According to the definitions of the continuous and discrete cost functionals and property (5.13) for $\tilde{w}=w$ and $\tilde{w}_{h}=w_{h}$ we get

[TABLE]

We set $p_{Th}:=A_{h}^{-1}(\rho\partial_{t}y(T)-z_{3})$.

Owing to the adjoint problem (3.12) with

[TABLE]

we have

[TABLE]

Similarly owing to the general discrete adjoint state equation (8.1)-(8.2) for $v=y_{\vartheta}$ and the discrete state equation (5.1)-(5.2) for $v=p_{\vartheta}$ and $\varphi=p_{0h}^{1}$ we get

[TABLE]

In addition owing to the definitions (8.2) of $p_{\vartheta}(T)$ and (5.15) of $A_{h}^{-1}$, we can write

[TABLE]

Consequently we obtain

[TABLE]

In addition using property (5.13) we derive

[TABLE]

Next, for the term $(\rho y^{0},\partial_{t}p(0)-{p_{0h}^{1}})_{H}$ in (10.4) we have

[TABLE]

due to the bounds $\|y^{0}-\pi_{h}^{0}y^{0}\|_{H_{\rho}}\leq\|y^{0}-\pi_{h}^{1}y^{0}\|_{H_{\rho}}$, (5.18) and (6.2). Clearly also $|(\rho y^{1},p(0)-p_{\vartheta}(0))_{H}|\leq\|y^{1}\|_{H}\|p(0)-p_{\vartheta}(0)\|_{H}$. Finally from (10.3)-(10.6) we derive (10.2). ∎

Now we prove for the cost functional a higher order error estimate than (9.3) for the state variable in the case $\mathcal{M}_{T}={L_{w^{*}}^{2}(I,\mathcal{M}(\Omega))}$.

Theorem 10.1**.**

Let $\mathcal{M}_{T}={L_{w^{*}}^{2}(I,\mathcal{M}(\Omega))}$, $\mathbf{z}\in\mathcal{Y}^{1/2}$ and $\mathbf{y}\in V\times{H}$. Then the following error estimate for the cost functional holds

[TABLE]

Proof.

Let us base on Proposition 13 and take any $u\in{L_{w^{*}}^{2}(I,\mathcal{M}(\Omega))}$. Owing to Proposition 10 (2) we have

[TABLE]

Proposition 11 (3) leads to

[TABLE]

Owing to Propositions 2(1) (applied to the adjoint state problem) and 3 we have

[TABLE]

(like in estimates (8.11)-(8.12) for $\alpha=1$). By using estimate (5.17) for $\lambda=-1/2$ we obtain

[TABLE]

By collecting all these estimates together with embedding (3.29), Proposition 11 (2) to bound $\|\rho(\partial_{t}{p(0)-{p}_{0h}^{1}})\|_{\mathcal{V}_{\kappa,h}^{*}}$ and applying Proposition 13, we derive

[TABLE]

Owing to inequalities (10.1) together with bounds (4.2) for $\bar{u}$ and (7.1) for $\bar{u}_{\vartheta}$ the proof is complete. ∎

Remark 11**.**

In the case $\mathcal{M}_{T}=\mathcal{M}(\Omega,L^{2}(I))$ we know that $\bar{u}\in{\mathcal{C}^{1}(\bar{I},\mathcal{M}(\Omega))}$ (cf. Theorem 4.1). The lack of the corresponding bound at least $\|\bar{u}_{\vartheta}\|_{SHW^{-1/2,1;2}(Q)}\leq C$ at the discrete level does not allow us to prove the error estimate $|j(\bar{u})-j_{\vartheta}(\bar{u}_{\vartheta})|\leq C(\tau+h)^{4/3}$. The estimate $|j(\bar{u})-j_{\vartheta}(\bar{u}_{\vartheta})|\leq C(\tau+h)^{2/3}$ follows directly from (9.4).

11. Time-stepping formulation

In this section we discuss the time-stepping formulation of the discrete state equation (5.1)-(5.2) and the discrete adjoint state equation (7.4).

We introduce the piecewise-linear “hat” functions such that $e^{\tau}_{m}(t_{k})=\delta_{m,k}$ for any $k,m=0,\ldots,M$, where $\delta_{m,k}$ is the Kroneker delta. We recall that $e^{\tau}_{m}$ are “half” hat functions for $m=0,M$. There holds $V_{\tau}=\operatorname{span}\{e_{0}^{\tau},\ldots,e_{M}^{\tau}\}$. Similarly, we introduce the spatial hat functions such that $e^{h}_{j}(x_{k})=\delta_{j,k}$ for any $j=1,\ldots,N-1$ and $k=0,\ldots,N$; then $V_{h}=\operatorname{span}\{e^{h}_{1},\ldots,e^{h}_{N-1}\}$.

Then the approximate state variable $y_{\vartheta}\in V_{\vartheta}$ can be represented in the following forms

[TABLE]

for $(t,x)\in\bar{I}\times\bar{\Omega}$ with $y_{m,j}\in\mathbb{R}$, $y_{m}^{h}\in V_{h}$ and $y_{j}^{\tau}\in V_{\tau}$.

We also define the forward and backward difference quotients and the average in time operator

[TABLE]

We define the self-adjoint positive-definite operators $B_{h}$ and $L_{h}$ acting in $V_{h}$ (in other words, the mass and stiffness matrices) such that

[TABLE]

For $w\in V^{\ast}$ and ${u}\in L^{2}(I,V^{\ast})$ we define the vectors $w^{h}=\{\langle w,e_{j}^{h}\rangle_{\Omega}\}_{j=1}^{N-1}$ and

[TABLE]

We recall the form of the discrete state (11.1).

The forward time-stepping is implemented as follows. The integral identities (5.1)-(5.2) are equivalent to the operator equations

[TABLE]

followed by the counterpart of (11.3) at time $T$ for $y_{Th}^{1}$:

[TABLE]

Next the adjoint (backward) time-stepping is implemented in a similar manner. Namely, the integral identities (7.4) are equivalent to the operator equations

[TABLE]

followed by the counterpart of (11.5) for $p_{0h}^{1}$:

[TABLE]

Remark 12**.**

For $\sigma=1/4$ the three-level time stepping scheme (11.2)-(11.5) is closely related to the well-known two-level Crank-Nicolson method applied to the first order in time system

[TABLE]

see [51, Section 8] for details, as well as to the Petrov-Galerkin method described in [32]. After the mass lumping, for $\sigma=0$ our method becomes explicit and is related to the Leap-Frog method; moreover, for any $\sigma$ it becomes close to three-level finite-difference schemes with such weight in time, eg. see [47].

12. Control discretization. Solution process and $L^{2}(I\times\Omega)$-regularization

Now we discuss in more detail solving of the semi-discrete optimization problem ($\mathcal{P}_{\vartheta}$) in the case $\mathcal{M}_{T}=\mathcal{M}(\Omega,L^{2}(I))$.

An important point is that we can seek its solution in the form

[TABLE]

To show that, let $\pi_{\tau}^{0}$ be the projector in $L^{2}(I)$ on $V_{\tau}$. Note that, for $\eta\in L^{2}(I)$, it satisfies

[TABLE]

Then we define $\Pi_{h}$: $\mathcal{M}(\Omega)\to\mathcal{M}_{h}$ by $\Pi_{h}w:=\sum_{j=1}^{N-1}\langle w,e_{j}^{h}\rangle_{\Omega}\delta_{x_{j}}$ and $\Pi_{\vartheta}=\pi_{\tau}^{0}\Pi_{h}$. The following identity holds

[TABLE]

with the interpolation operator $i_{h}$: $\mathcal{C}_{0}(\Omega)\to V_{h}$ such that $i_{h}w(x_{j})=w(x_{j})$ for all $j=0,\ldots,N$. In particular, if $v\in V_{\vartheta}$, then

[TABLE]

and consequently (like in [33, Lemma 3.11]) we have $S_{\vartheta}=S_{\vartheta}\circ\Pi_{\vartheta}$ as well as $\|\Pi_{\vartheta}u\|_{\mathcal{M}_{T}}\leq\|u\|_{\mathcal{M}_{T}}$. Thus for each solution $\tilde{u}_{\vartheta}$ of problem ($\mathcal{P}_{\vartheta}$), the discrete control $\Pi_{\vartheta}\tilde{u}_{\vartheta}$ satisfies

[TABLE]

Therefore $\Pi_{\vartheta}\tilde{u}_{\vartheta}$ is also a solution of ($\mathcal{P}_{\vartheta}$). This is a justification for solving the fully discrete problem

[TABLE]

in order to get a solution of ($\mathcal{P}_{\vartheta}$).

The direct solution of (12.1) by means of a generalized Newton type method is a challenging problem since a proper globalization strategy is needed, see [39]. Thus we propose a solution strategy based on an additional $L^{2}(I\times\Omega)$-regularization of (12.1) with a parameter $\gamma>0$ and a continuation method. For high values of $\gamma$ the corresponding Newton type method converges independently of the initial guess in numerical practice. Thus the continuation strategy can be seen as simple globalization strategy.

On the continuous level we consider the following regularized problem

[TABLE]

It is possible to formulate a semi-smooth Newton method for this problem on the continuous level which is based on the following necessary and sufficient optimality condition

[TABLE]

with $\bar{p}=S^{\star}{W_{h}(S\bar{u}_{\gamma}-\mathbf{z}})$. Moreover, this semi-smooth Newton method is superlinear convergent. Let $\bar{u}_{\gamma}$ and $\bar{u}$ be the unique solutions of (12.2) and ($\mathcal{P}$). Then we have $\bar{u}_{\gamma}\rightharpoonup^{\ast}\bar{u}$ in $\mathcal{M}(\Omega,L^{2}(I))$, see [33, 43, 26]. This justifies the use of a continuation strategy in $\gamma$. The control discretization described above can not be used for (12.2). Instead we propose to use discrete controls from $V_{\vartheta}$, i.e.,

[TABLE]

cf. (11.1). In particular, we solve the following fully discrete regularized problem

[TABLE]

with

[TABLE]

where $D=\operatorname{diag}(d_{1},\ldots,d_{N-1})$ is the lumped mass matrix. Moreover, the operator $l_{\vartheta}$ is defined by

[TABLE]

The use of $D$ allows us to derive the following optimality conditions for (12.4)

[TABLE]

for all $m$ and $j$, with $\bar{p}_{\vartheta}=S^{\star}_{\vartheta}{W_{h}}(S_{\vartheta}\bar{u}_{\vartheta}-{\mathbf{z}})$, cf. (12.3). Based on (12.5) we can set up a semi-smooth Newton method. Since problem (12.4) is a discretization of (12.2), we can expect that this method behaves mesh independently. Let $\bar{u}_{\vartheta}^{\gamma}=\sum_{j=1}^{N-1}u_{j}(t)e^{h}_{j}$ be the solution of (12.4) and we define

[TABLE]

As $\gamma\rightarrow 0$ the control $\tilde{u}_{\vartheta}^{\gamma}$ tends to a solution of (12.1) justifying the use of this control discretization and the continuation strategy. For more details see [43].

13. Numerical results

In this section, we present results of numerical experiments and consider two examples both involving zero initial data $y^{0}=y^{1}=0$, the control space $\mathcal{M}_{T}=\mathcal{M}(\Omega,L^{2}(I))$ and the tracking functional

[TABLE]

with the time independent desired state $z$ which is a Gaussian centered at $x=\lambda$. We choose $\rho=0.1$ and $\lambda$ as an irrational parameter.

For sufficiently large $\alpha$ ($\alpha=0.1$), we expect that the optimal control $\bar{u}$ consists of one point source with a position close to $\lambda$. If the Gaussian would move through the domain, a point source shaped $\bar{u}$ is not able to follow the center of the Gaussian since $\mathcal{M}(\Omega,L^{2}(I))$ contains no moving point sources. The optimal control would rather consist of some additional fixed point sources. This would not lower the regularity of the state whereas a moving point source can cause it.

The domain $\Omega$ and the time segment $\bar{I}$ are discretized by the uniform grids for $N=2^{r_{h}}$ and $M=2^{r_{\tau}}$ where $r_{\tau},r_{h}=2,3,\ldots,r^{\max}$ with $r^{\max}=10$. The stability parameter is fixed to its lowest value $\sigma=1/4$ ensuring unconditional stability of the time-stepping method. The discrete control problem is solved for $r_{h}=2,3,\ldots,{r^{\max}}$ and the fixed $r_{\tau}=r^{\max}$ and then vice versa. The solution process has been described above in Section 12. Numerically the desired state $z$ is replaced by $i_{h}z$ for simplicity, moreover the corresponding error $\mathcal{O}(h^{2})$ is negligible. Since the optimal pairs $(\bar{u},\bar{y})$ are not known in our examples, we replace them by reference solutions $(\hat{u},\hat{y})$ which are taken as the approximate solutions on the finest grid level.

Example 1. We first take the constant coefficients $\rho\equiv 1$ and $\kappa\equiv 1$ and set $\lambda=\pi/20$. We depict the reference solution $(\hat{u},\hat{y})$ in Figure 1.

As expected, the optimal control $\hat{u}$ consists only of one point source positioned in the vicinity of $\lambda$. Thus, the state $\hat{y}$ has a kink at this position. Due to reflections at the boundary, $\hat{y}$ has also kinks at other positions.

Next, we discuss the convergence results. In Figure 2, we see the convergence rate of $\|\bar{y}_{\sigma}-\hat{y}\|_{L^{2}(I\times\Omega)}$ (left) and the objective functional (right) as $h$ refines. The state error behaves mostly in a linear way and the rate for the functional is close to two; as usual the latter is approximately the doubled rate of the former, and fortunately both are better than the above proved theoretical rates.

In Figure 3, we see the similar results as $\tau$ refines. The error of the functional stagnates at the last $\tau$ refinement that is caused by a too coarse space grid. Nevertheless, we observe reduced rates for $\hat{y}$ much less than two caused by its reduced regularity (kinks).

Example 2. Now we take the variable coefficient

[TABLE]

and set $\lambda=\pi/6$. Our analysis does not cover discontinuous coefficients, but they are of great importance in applications, for example, in seismic tomography. A jump discontinuity in $\kappa$ translates to a jump in the wave speed which can be related to two different material characteristics changing at the point of discontinuity. Note that the point of discontinuity is a grid point for all grid levels.

The reference solution $(\hat{u},\hat{y})$ is displayed in Figure 4. Once again $\hat{u}$ consists of one Dirac measure with a time-dependent intensity located in the vicinity of $\lambda=\pi/6$ and thus $\hat{y}$ has a kink at this position. Moreover, we can clearly see that at $x=0.25$ the wave speed changes and the wave propagation becomes slower.

In Figure 5 we observe that the error of the state variable converges in a linear way whereas the error measured in the objective functional behaves quadratically.

Finally in Figure 6 we study the error behaviour for $\tau$-refinement and find the similar rates of convergence. So somewhat surprisingly we find that the convergence behavior of the error is comparable to the previous Example 1 with $\kappa\equiv 1$ that stimulates further possible studies.