Carleman contraction mapping for a 1D inverse scattering problem with experimental time-dependent data

We consider in this paper a 1D Coefficient Inverse Problem (CIP) for a hyperbolic PDE and construct a globally convergent numerical method for it. Unlike all previous works of this group on the convexification method for CIPs, which are cited below in this section, we construct in this paper a sequence of linearized boundary value problems with overdetermined boundary data. Each of these problems is solved using a weighted version of the Quasi-Reversibility Method (QRM). The weight is the Carleman Weight Function (CWF), i.e. this is the function, which is involved as the weight in the Carleman estimate for the corresponding PDE operator, see, e.g. [7, 4, 5, 22, 18, 31, 38, 53] for Carleman estimates. Thus, we call this method “Carleman Quasi-Reversibility Method” (CQRM). The key convergence estimate for this sequence is similar with the one of the conventional contraction mapping principle, see the first item of Remarks 8.1 in section 8. This explains the title of our paper. That estimate implies the global convergence of that sequence to the true solution of our CIP. Furthermore, that estimate implies a rapid convergence. As a result, computations are done here in real time, also, see Remark 11.1 in section 11 as well as section 12. The real time computation is obviously important for our target application, which is described below. The numerical method of this paper can be considered as the second generation of the convexification method for CIPs.

The convexification method ends up with the minimization of a weighted globally strictly convex Tikhonov functional, in which the weight is the CWF. The convexification has been developed by our research group since the first two publications in 1995 [20] and 1997 [21]. Only the theory  was presented in those two originating papers. The start for numerical studies of the convexification was given in the paper [1]. Thus, in follow up publications [24, 25, 30, 16, 14, 15, 24, 25, 27, 28, 29, 32, 33, 30, 40, 50, 51] the theoritical results for the convexification are combined with the numerical ones. In particular, some of these papers work with experimentally collected backscattering data [14, 15, 26, 32, 33]. We also refer to the recently published book of Klibanov and Li [31], in which many of these results are described.

The CIP of our paper has a direct application in the problem of the standoff detection and identification of antipersonnel land mines and improvised explosive devices (IEDs). Thus, in the numerical part of this publication, we present results of the numerical performance of our method for both computationally simulated and experimentally collected data for targets mimicking antipersonnel land mines and IEDs. The experimental data were collected in the field, which is more challenging than the case of the data collected in a laboratory [14, 15, 26, 33]. The data used in this paper were collected by the forward looking radar of the US Army Research Laboratory [46]. Since these data were described in the previous publications of our research group [13, 23, 24, 25, 30, 35, 36, 51], then we do not describe them here. Those previous publications were treating these experimental data by several different numerical methods for 1D CIPs. We also note that the CIP of our paper was applied in [32, 33] to the nonlinear synthetic aperture (SAR) imaging, including the cases of experimentally collected data.

Our goal is to compute approximate values of dielectric constants of targets. We point out that our experimental data are severely underdetermined. Indeed, any target is a 3D object. On the other hand, we have only one experimentally measured time resolved curve for each target. Therefore, we can compute only a sort of an average value of the dielectric constant of each target. This explains the reason of our mathematical modeling here by a 1D hyperbolic PDE rather than by its 3D analog. This mathematical model, along with its frequency domain analog, was considered in all above cited previous publications of our research group, which were working with the experimental data of this paper by a variety of globally convergent inversion techniques.

An estimate of the dielectric constant of a target cannot differentiate with a good certainty between an explosive and a non explosive. However, we hope that these estimates, combined with current targets classification procedures, might potentially result in a decrease of the false alarm rate for those targets. We believe therefore that our results for experimental data might help to address an important application to the standoff detection and identification of antipersonnel land mines and IEDs.

Given a numerical method for a CIP, we call this method globally convergent if:

1. 1.

   There is a theorem claiming that this method delivers at least one point in a sufficiently small neighborhood of the true solution of that CIP without relying on a good initial guess about this solution.

2. 2.

   This theorem is confirmed numerically.

CIPs are both highly nonlinear and ill-posed. As a result, conventional least squares cost functionals for them are non convex and typically have many local minima and ravines, see, e.g. [49] for a numerical example of this phenomenon. Conventional numerical methods for CIPs are based on minimizations of those functionals, see, e.g. [8, 11, 12]. The goal of the convexification is to avoid local minima and ravines and to provide accurate and reliable solutions of CIPs.

In the convexification, a CWF is a part of a certain least squares cost functional $J_{\lambda}$, where $\lambda>1$ is the parameter of the CWF. The presence of the CWF ensures that, with a proper choice of the parameter $\lambda,$ the functional $J_{\lambda}$ is strictly convex on a certain bounded convex set $S$ of an arbitrary diameter $d>0$ in a Hilbert space. It was proven later in [1] that the existence and uniqueness of the minimizer of $J_{\lambda}$ on $S$ are also guaranteed. Furthermore, the gradient projection method of the minimization of $J_{\lambda}$ on $S$ converges globally to that minimizer if starting at an arbitrary point of $S$. In addition, if the level of noise $\delta>0$ in the data tends to zero, then the gradient projection method delivers a sequence, which converges to the true solution. Thus, since $d>0$ is an arbitrary number, then this is the global convergence, as long as numerical studies confirm this theory. Furthermore, it was recently established in [33, 40] that the gradient descent method of the minimization of $J_{\lambda}$ on $S$ also converges globally to the true solution.

First, we come up with the same boundary value problem (BVP) for a nonlinear and nonlocal hyperbolic PDE as we did in the convexification method for the same 1D CIP in [30]. The solution of this BVP immediately provides the target unknown coefficient. However, when solving this BVP, rather than constructing a globally strictly convex cost functional for this BVP, as it was done in [30] and other works on the convexification, we construct a sequence of linearized BVPs. Just as the original BVP, each BVP of this sequence involves a hyperbolic operator with non local terms and overdetermined boundary conditions.

Because of non local terms and the overdetermined boundary conditions, we solve each BVP of this sequence by the new version of the QRM, the Carleman Quasi-Reversibility Method (CQRM). Indeed, the classical QRM solves ill-posed/overdetermined BVPs for linear PDEs, see [37] for the originating work as well as, e.g. [19, 31] for updates. Convergence of the QRM-solution to the true solution is proven via an appropriate Carleman estimates [19, 31]. However, a CWF is not used in the optimizing quadratic functional of the conventional QRM of [19, 31]. Unlike this, a CWF is involved in CQRM. It is exactly the presence of the CWF, which enables us to derive the above mentioned convergence estimate, which is an analog of the one of the contraction principle. This allows us, in turn to establish convergence rate of the sequence of CQRM-solutions to the true solution of that BVP. The latter result ensures that our method is a globally convergent one.

Various versions of the second generation of the convexification method involving CQRM were previously published in [2, 3, 45, 44, 39]. In these publications, globally convergent numerical methods for some nonlinear inverse problems are presented. In [2, 3] CIPs for hyperbolic PDEs in $\mathbb{R}^{n}$ are considered. The main difference between [2, 3] and our work is that it is assumed in [2, 3] that one of initial conditions is not vanishing everywhere in the closed domain of interest. In other words, papers [2, 3] work exactly in the framework of the Bukhgeim-Klibanov method, see [6] for the originating work on this method and, e.g. [4, 5, 17, 22, 18, 31, 53] for some follow up publications. On the other hand, in all above cited works on the convexification, including this one, either the initial condition is the $\delta-$function, or a similar requirement holds for the Helmholtz equation. The only exception is the paper [29], which also works within the framework of the method of [6]. We refer to papers [10, 34] for some numerical methods for CIPs with the Dirichlet-to-Neumann map data. In these works the number $m$ of free variables in the data exceeds the number $n$ of free variables in the unknown coefficient, $m>n.$ In our paper $m=n=1.$ Also, $m=n$ in all other above cited works on the convexification.

There is a classical Gelfand-Levitan method [41] for solutions of 1D CIPs for some hyperbolic PDEs. This method does not rely on the optimization, and, therefore, avoids the phenomenon of local minima. It reduces the original CIP to a linear integral equation of the second kind. This equation is called “Gelfand-Levitan equation” (GL). However, the questions about uniqueness and stability of the solution of GL for the case of noisy data is open, see, e.g. Lemma 2.4 in the book [48, Chapter 2]. This lemma is valid only in the case of noiseless data. However, the realistic data are always noisy. In addition, it was demonstrated numerically in [13] that GL cannot work well for exactly the same experimental data as the ones we use in the current paper. On the other hand, it was shown in [24, 25, 51] that the convexification works well with these data. The same is true for the method of the current paper.

Uniqueness and Lipschitz stability theorems of the CIP considered here are well known. Indeed, it was shown in, e.g. [51] that, using the same change of variables as the one considered in section 2 below, one can reduce our CIP to a similar CIP for the equation $v_{tt}=v_{yy}+r\left(y\right)v,y\in\mathbb{R}$ with the unknown coefficient $r\left(y\right).$ We refer to [48, Theorem 2.6 of Chapter 2] for the Lipschitz stability estimate for the latter CIP. In addition, both uniqueness and Lipschitz stability results for our CIP actually follow from Theorem 8.1 below as well as from the convergence analysis in two other works of our research group on the convexification method for this CIP [30, 51].

This paper is arranged as follows. In section 2 we state both forward and inverse problems. In section 3 we derive a nonlinear BVP with non local terms. In section 4 we describe iterations of our method to solve that BVP. In section 5 we formulate the Carleman estimate for the principal part of the PDE operator of that BVP. In section 6 we prove the strong convexity of a functional of section 5 on an appropriate bounded set. In section 7 we formulate two methods for finding the unique minimizer of that functional: gradient descent method and gradient projection method and prove the global convergence to the minimizer for each of them. In section 8 we establish the contraction mapping property and prove the global convergence of the method of section 4. In section 9 we formulate two more global convergence theorems, which follow from results of sections 7 and 8. Numerical results with simulated and experimental data are presented in Section 10 and 11 respectively. Concluding remarks are given in section 12.

Below all functions are real valued ones. Let $b>1$ be a known number, $x\in\mathbb{R}$ be the spatial variable and the function $c(x)\in C^{3}(\mathbb{R})$, represents the spatially distributed dielectric constant. We assume that

	$\displaystyle c(x)$	$\displaystyle\in[1,b],x\in\mathbb{R},$		(2.1)
	$\displaystyle c(x)=1$	$\displaystyle\mbox{ if }x\in(-\infty,\varepsilon]\cup[1,\infty),$		(2.2)

where $\varepsilon\in\left(0,1\right)$ is a small number. Let $T$ be a positive number. We consider the following Cauchy problem for a 1D hyperbolic PDE with a variable coefficient in the principal part of the operator:

$$c(x)u_{tt}=u_{xx},\text{ }x\in\mathbb{R},t\in\left(0,T\right),$$

(2.3)

$$u\left(x,0\right)=0,u_{t}\left(x,0\right)=\delta\left(x\right).$$

(2.4)

The problem of finding the function $u(x,t)$ from conditions (2.3), (2.4) is our Forward Problem.

Let $\tau\left(x\right)$ be the travel time needed for the wave to travel from the point source at $\left\{0\right\}$ to the point $x$,

$$\tau\left(x\right)=\mathop{\displaystyle\int}\limits_{0}^{x}\sqrt{c(s)}ds$$

(2.5)

By (2.5) the following 1D analog of the eikonal equation is valid:

$$\tau^{\prime}\left(x\right)=\sqrt{c\left(x\right)}.$$

(2.6)

Let $H(z),z\in\mathbb{R}$ be the Heaviside function,

$$H(z)=\left\{\begin{array}[]{c}1,\quad z>0,\\ 0,\quad z<0.\end{array}\right.$$

Lemma 2.1 [30, Lemma 2.1]. For $x\geq 0$, the function $u(x,t)$ has the form

$$u(x,t)=\frac{H(t-\tau\left(x\right))}{2c^{1/4}(x)}+\hat{u}(x,t),$$

(2.7)

where the function $\hat{u}\in C^{2}(t\geq\tau\left(x\right)),\hat{u}(x,\tau\left(x\right))=0.$ Furthermore,

$$\lim_{t\rightarrow\tau^{\div}\left(x\right)}u(x,t)=\frac{1}{2c^{1/4}(x)}.$$

(2.8)

Lemma 2.2 [Absorbing boundary condition] [30, 50]. Let $b>\varepsilon$ be the number in (2.1). Let $x_{1}\geq b$ and $x_{2}\leq\varepsilon$ be two arbitrary numbers. Then the solution $u(x,t)$ of Forward Problem (2.3), (2.4) satisfies the absorbing boundary condition at $x=x_{1}$ and $x=x_{2}$, i.e.

$$u_{x}(x_{1},t)+u_{t}(x_{1},t)=0\mbox{ for }t\in(0,T),$$

(2.9)

$$u_{x}(x_{2},t)-u_{t}(x_{2},t)=0\mbox{ for }t\in(0,T).$$

(2.10)

We are interested in the following inverse problem:

Coefficient Inverse Problem (CIP). Suppose that the following two functions $g_{0}(t)$, $g_{1}(t)$ are known:

$$u\left(\varepsilon,t\right)=g_{0}\left(t\right),\quad u_{x}\left(\varepsilon,t\right)=g_{1}\left(t\right),\quad t\in\left(0,T\right).$$

(2.11)

Determine the function $c(x)$ for $x\in(\varepsilon,1)$, assuming that the number $b>\varepsilon$ in (2.1) is known.

Remark 2.1. Note that only the function $g_{0}\left(t\right)$ can be measured. As to the function $g_{1}\left(t\right),$ it follows from (2.10) that

$$g_{1}\left(t\right)=g_{0}^{\prime}\left(t\right).$$

(2.12)

We differentiate noisy function using the Tikhonov regularization method [52]. Since this method is well known, we do not describe it here.

We now introduce a change of variable

$$q(x,t)=u(x,t+\tau(x)).$$

(3.1)

We will consider the function $q(x,t)$ only for $t\geq 0.$ Using (2.6) and (3.1), we obtain

$$q_{xx}-2q_{xt}\tau^{\prime}-q_{t}\tau^{\prime\prime}=0.$$

(3.2)

Furthermore, it follows from (2.8)

$$q(x,0)=\frac{1}{2c^{1/4}(x)}\neq 0.$$

(3.3)

By (2.5) and (3.3)

$$\tau^{\prime\prime}(x)=-\frac{q_{x}(x,0)}{2q^{3}(x,0)}.$$

(3.4)

Substituting (2.6), (3.3) and (3.4) in (3.2), we obtain

$$L(q)=q_{xx}-q_{xt}\frac{1}{2q^{2}(x,0)}+q_{t}\frac{q_{x}(x,0)}{2q^{3}(x,0)}=0.$$

(3.5)

Equation, (3.5) is a nonlinear PDE with respect to the function $q(x,t)$ with nonlocal terms $q_{x}(x,0)$ and $q(x,0)$. We now need to obtain boundary conditions for the function $q.$ By (2.2) and (2.5) $\tau\left(x\right)=x$ for $x\in\left[0,\varepsilon\right].$ Hence, (2.11), (2.12) and (3.1) lead to

$$q\left(\varepsilon,t\right)=g_{0}\left(t+\varepsilon\right),q_{x}\left(\varepsilon,t\right)=2g_{0}^{\prime}\left(t+\varepsilon\right),\text{ }t\in\left(0,T\right),$$

(3.6)

We will solve equation (3.5) in the rectangle

$$\Omega=\left\{{(x,t):x\in(\varepsilon,b),t\in(0,T)}\right\}.$$

(3.7)

By (3.1)

$$q_{x}(x,t)=u_{x}(x,t+\tau(x))+u_{t}(x,t+\tau(x))\tau^{\prime}(x)$$

(3.8)

By (2.1), (2.2) and (2.5) $\tau^{\prime}(b)=1$. Hence, using (2.9) and (3.8), we obtain

$$q_{x}(b,t)=0.$$

(3.9)

It follows from (2.1) and (3.3) that

$$\frac{1}{2b^{1/4}}\leq q(x,0)\leq\frac{1}{2},\quad x\in[\varepsilon,b].$$

(3.10)

In addition, we need $q\in C^{2}(\overline{\Omega})$ and we also need to bound the norm $\left\|q\right\|_{C^{2}(\overline{\Omega})}$ from the above. Let $R>0$ be an arbitrary number. We define the set $B(R,g_{0})$ as

$$B(R,g_{0})=\left\{\begin{array}[]{c}q\in H^{4}(\Omega):\left\|q\right\|_{H^{4}(\Omega)}<R,\\ q(\varepsilon,t+\varepsilon)=g_{0}(t+\varepsilon),q_{x}(\varepsilon,t)=2g_{0}^{\prime}\left(t+\varepsilon\right),\\ q_{x}(b,t)=0,\\ \frac{1}{2b^{1/4}}\leq q\left(x,0\right)\leq 1/2,x\in[\varepsilon,b].\end{array}\right.$$

(3.11)

We assume below that

$$B(R,g_{0})\neq\varnothing.$$

(3.12)

By embedding theorem $\overline{B(R,g_{0})}\subset C^{2}(\overline{\Omega})$ and

$$\left\|q\right\|_{C^{2}(\overline{\Omega})}\leq KR,\text{ }\forall q\in\overline{B(R,g_{0})}.$$

(3.13)

where the constant $K=K\left(\Omega\right)>0$ depends only on the domain $\Omega.$ Using (3.10), we define the function $q^{0}\left(x,0\right)$ as:

	$\displaystyle q^{0}\left(x,0\right)$	$\displaystyle=$	$\displaystyle\left\{\begin{array}[]{c}q\left(x,0\right)\text{ if }q\left(x,0\right)\in\left[1/\left(2b^{1/4}\right),1/2\right],\\ 1/\left(2b^{1/4}\right),\text{ if }q\left(x,0\right)<1/\left(2b^{1/4}\right),\\ 1/2\text{ if }q\left(x,0\right)>1/2.\end{array}\right.\text{ }\forall q\in\overline{B(R,g_{0})},$		(3.17)
	$\displaystyle\forall x$	$\displaystyle\in$	$\displaystyle\left[\varepsilon,b\right].$

Then the function $q^{0}\left(x,0\right)$ is piecewise continuously differentiable in $\left[\varepsilon,b\right]$ and by (3.13) and (3.17)

$$\left\|q^{0}\left(x,0\right)\right\|_{C\left[\varepsilon,b\right]},\max_{\left[\varepsilon,b\right]}\left|q_{x}^{0}\left(x,0\right)\right|\leq KR,\text{ }\forall q\in B(R,g_{0}),$$

(3.18)

$$\frac{1}{2b^{1/4}}\leq q^{0}(x,0)\leq\frac{1}{2},\quad x\in[\varepsilon,b].$$

(3.19)

Thus, (3.5), (3.6), (3.9), (3.17) and (3.19) result in the following BVP for a nonlinear PDE with non-local terms:

$$q_{xx}(x,t)-q_{xt}(x,t)\frac{1}{2\left(q^{0}(x,0)\right)^{2}}+q_{t}(x,t)\frac{q_{x}(x,0)}{2\left(q^{0}(x,0)\right)^{3}}=0\text{ in }\Omega,$$

(3.20)

$$q(\varepsilon,t)=g_{0}(t+\varepsilon),\text{ }q_{x}(\varepsilon,t)=2g_{0}^{\prime}(t+\varepsilon),\text{ }q_{x}(b,t)=0.$$

(3.21)

Thus, we focus below on the numerical solution of the following problem:

Problem 3.1. Find a function $q\in B(R,g_{0})$ satisfying conditions (3.20), (3.21), where the function $q^{0}(x,0)$ is defined in (3.17).

Suppose that we have solved Problem 3.1. Then, using (3.3) and (3.17), we set

$$c\left(x\right)=\frac{1}{\left(2q^{0}\left(x,0\right)\right)^{4}}.$$

(3.22)

In this section we formulate the Carleman estimate, which is the main tool of our construction. This estimate follows from Theorem 3.1 of [30] as well as from (3.18). Let $q\left(x,t\right)\in\overline{B(R,g_{0})}$ be an arbitrary function and let the function $q^{0}\left(x,0\right)$ be constructed from the function $q\left(x,t\right)$ as in (3.17). Consider the operator $L_{0}:H^{2}(\Omega)\rightarrow L_{2}(\Omega),$

$$L_{0}u=u_{xx}(x,t)-u_{xt}(x,t)\frac{1}{2\left(q^{0}(x,0)\right)^{2}},\quad\mbox{for all }(x,t)\in\Omega.$$

Theorem 5.1 (Carleman estimate [30]). There exists a number $\alpha_{0}=\alpha_{0}\left(R,\Omega\right)>0$ depending only on $R,\Omega$ such that for any $\alpha\in(0,\alpha_{0})$ there exists a sufficiently large number $\lambda_{0}=\lambda_{0}\left(R,\Omega,\alpha\right)>1$ depending only on $R$,$\Omega,\alpha$ such that for all $\lambda\geq\lambda_{0}$ and for all functions $v\in H^{2}(\Omega)$ the following Carleman estimate holds:

$$\mathop{\displaystyle\int}\limits_{\Omega}\left(L_{0}v\right)^{2}e^{2\lambda\varphi}dxdt$$

$$\geq C\mathop{\displaystyle\int}\limits_{\Omega}\big{(}\lambda\left(v_{x}^{2}+v_{t}^{2}\right)+\lambda^{3}v^{2}\big{)}e^{2\lambda\varphi}dxdt+C\mathop{\displaystyle\int}\limits_{\varepsilon}^{b}\big{(}\lambda v_{x}^{2}(x,0)+\lambda^{3}v^{2}(x,0)\big{)}e^{-2\lambda x}dx$$

$$-C\mathop{\displaystyle\int}\limits_{0}^{T}\left(\lambda\left(v_{x}^{2}(\varepsilon,t)+v_{t}^{2}(\varepsilon,t)\right)+\lambda^{3}v^{2}(\varepsilon,t)\right)e^{-2\lambda(\varepsilon+\alpha t)}dt$$

(5.1)

$$-C\mathop{\displaystyle\int}\limits_{\varepsilon}^{b}\left(\lambda v_{x}^{2}(x,T)+\lambda^{3}v^{2}(x,T)\right)e^{-2\lambda(x+\alpha T)}dx.$$

 

Remark 5.1. Here and everywhere below $C=C\left(R,\Omega,\alpha\right)>0$ denotes different constants depending only on listed parameters.

Functional (4.8) is quadratic. We prove in this section that it is strictly convex on the set $\overline{B(R,g_{0})}.$ In addition, we prove existence and uniqueness of its minimizer on this set. Although similar results were proven in many of the above cited publications on the convexification, see, e.g. [30] for the closest one, there are some peculiarities here, which are important for our convergence analysis, see Remarks 6.1, 6.2 below.

Introduce the subspace $H_{0}^{4}\left(\Omega\right)\subset H^{4}\left(\Omega\right)$ as:

$$H_{0}^{4}\left(\Omega\right)=\left\{v\in H^{4}\left(\Omega\right):v\left(\varepsilon,t\right)=v_{x}\left(\varepsilon,t\right)=v_{x}\left(b,t\right)=0\right\}.$$

(6.1)

Denote $\left[,\right]$ the scalar product in the space $H^{4}\left(\Omega\right).$ Also, denote

$$A_{n}\left(q\right)\left(x,t\right)=q_{xx}(x,t)-q_{xt}(x,t)\frac{1}{2\left(q_{n-1,\min}^{0}(x,0)\right)^{2}},$$

(6.2)

$$B_{n}\left(x,t\right)=\partial_{t}q_{\left(n-1\right),\min}\left(x,t\right)\frac{\partial_{x}q_{\left(n-1\right),\min}(x,0)}{2\left(q_{n-1,\min}^{0}(x,0)\right)^{3}}.$$

(6.3)

Theorem 6.1. 1. For any set of parameters $\lambda,\beta$ and for any $q\in\overline{B(R,g_{0})}$ the functional $J_{\lambda,\beta}^{\left(n\right)}$ has the Fréchet derivative $\left(J_{\lambda,\beta}^{\left(n\right)}(q)\right)^{\prime}\in H_{0}^{4}\left(\Omega\right).$ The formula for $\left(J_{\lambda,\beta}^{\left(n\right)}(q)\right)^{\prime}$ is:

$$\left(J_{\lambda,\beta}^{\left(n\right)}(q)\right)^{\prime}\left(h\right)=2\mathop{\displaystyle\int}\limits_{\Omega}\left(A_{n}\left(q\right)\left(x,t\right)+B_{n}\left(x,t\right)\right)A_{n}\left(h\right)\left(x,t\right)e^{2\lambda\varphi}dxdt+2\beta\left[q,h\right],$$

(6.4)

for all $h\in H_{0}^{4}\left(\Omega\right).$

2. This derivative is Lipschitz continuous in $\overline{B(R,g_{0})}$, i.e. there exists a constant $D>0$ such that

$$\left\|\left(J_{\lambda,\beta}^{\left(n\right)}(q_{2})\right)^{\prime}-\left(J_{\lambda,\beta}^{\left(n\right)}(q_{1})\right)^{\prime}\right\|_{H^{4}\left(\Omega\right)}\leq D\left\|q_{2}-q_{1}\right\|_{H^{4}\left(\Omega\right)},\quad\mbox{for all }q_{1},q_{2}\in\overline{B(R,g_{0})}.$$

(6.5)

3. Let $\alpha_{0}=\alpha_{0}\left(R,\Omega\right)>0,\alpha\in\left(0,\alpha_{0}\right)$ and $\lambda_{0}=\lambda_{0}\left(R,\Omega,\alpha\right)\geq 1$ be the numbers of Theorem 5.1. Then there exists a sufficiently large constant

$$\lambda_{1}=\lambda_{1}\left(R,\Omega,\alpha\right)\geq\lambda_{0}>1$$

(6.6)

depending only on listed parameters such that for all $\lambda\geq\lambda_{1}$ and for all $\beta\in\left[2e^{-\lambda\alpha T},1\right)$ the functional $J_{\lambda,\beta}^{\left(n\right)}(q)$ is strictly convex on the set $\overline{B\left(R,g\right)}.$ More precisely, let $q\in\overline{B(R,g)}$ be an arbitrary function and also let the function $q+h\in\overline{B(R,g)}.$ Then the following inequality holds:

$$J_{\lambda,\beta}^{\left(n\right)}(q+h)-J_{\lambda,\beta}^{\left(n\right)}(q)-\left(J_{\lambda,\beta}^{\left(n\right)}(q)\right)^{\prime}(h)\geq C\mathop{\displaystyle\int}\limits_{\Omega}\Big{[}\lambda\left(h_{x}^{2}+h_{t}^{2}\right)+\lambda^{3}h^{2}\Big{]}e^{2\lambda\varphi}dxdt$$

$$+C\mathop{\displaystyle\int}\limits_{\varepsilon}^{b}\big{(}\lambda h_{x}^{2}(x,0)+\lambda^{3}h^{2}(x,0)\big{)}e^{-2\lambda x}dx+\frac{\beta}{2}\left\|h\right\|_{H^{4}\left(\Omega\right)}^{2},\text{ }\forall\lambda\geq\lambda_{1}.$$

(6.7)

4. For any $\lambda\geq\lambda_{1}$ there exists unique minimizer

$$q_{n,\min}\in\overline{B(R,g_{0})}$$

(6.8)

 of the functional $J_{\lambda,\beta}^{\left(n\right)}(q)$ on the set $\overline{B(R,g_{0})}$ and the following inequality holds:

$$\left[\left(J_{\lambda,\beta}^{\left(n\right)}(q)\right)^{\prime},q-q_{n,\min}\right]\geq 0,\text{ }\forall q\in\overline{B(R,g_{0})},$$

(6.9)

where $\left[,\right]$ denotes the scalar product in $H^{4}\left(\Omega\right).$

Remark 6.1. Even though the expression in the right hand side of (6.4) is linear with respect to the function $q,$ we cannot use Riesz theorem here to prove existence and uniqueness of the minimizer $q_{n,\min},$ at which $\left(J_{\lambda,\beta}^{\left(n\right)}(q_{n,\min})\right)^{\prime}=0.$ Rather, all what we can prove is (6.9). This is because we need to ensure that the function $q_{n,\min}\in\overline{B(R,g_{0})}$.

Proof of Theorem 6.1. Since both function $q,q+h\in\overline{B(R,g)}$ satisfy the same boundary conditions, then

$$h\in H_{0}^{4}\left(\Omega\right).$$

(6.10)

By (4.8) and (6.10)

$$J_{\lambda,\beta}^{(n)}(q+h)-J_{\lambda,\beta}^{(n)}(q)=2\mathop{\displaystyle\int}\limits_{\Omega}\left(A_{n}\left(q\right)\left(x,t\right)+B_{n}\left(x,t\right)\right)A_{n}\left(h\right)\left(x,t\right)e^{2\lambda\varphi}dxdt+2\beta\left[q,h\right]$$

$$+\mathop{\displaystyle\int}\limits_{\Omega}\left[A_{n}\left(h\right)\left(x,t\right)\right]^{2}e^{2\lambda\varphi}dxdt+\beta\left\|h\right\|_{H^{4}\left(\Omega\right)}^{2},\text{ }\forall h\in H_{0}^{4}\left(\Omega\right).$$

(6.11)

The expression in the first line of (6.11) coincides with expression in the right hand side of (6.4). In fact, this is a bounded linear functional mapping $H_{0}^{4}\left(\Omega\right)$ in $\mathbb{R}$. Therefore, by Riesz theorem, there exists unique function $\widetilde{J}_{n}\left(q\right)\in H_{0}^{4}\left(\Omega\right)$ such that

$$\left[\widetilde{J}_{n}\left(q\right),h\right]=2\mathop{\displaystyle\int}\limits_{\Omega}\left(A_{n}\left(q\right)\left(x,t\right)+B_{n}\left(x,t\right)\right)A_{n}\left(h\right)\left(x,t\right)e^{2\lambda\varphi}dxdt+2\beta\left[q,h\right],\text{ }\forall h\in H_{0}^{4}\left(\Omega\right).$$

In addition, it is clear from (6.11) and (LABEL:5.9) that

$$\lim_{\left\|h\right\|_{H^{4}\left(\Omega\right)}\rightarrow 0}\left\{\frac{1}{\left\|h\right\|_{H^{4}\left(\Omega\right)}}\left[J_{\lambda,\beta}^{(n)}(q+h)-J_{\lambda,\beta}^{(n)}(q)-\left[\widetilde{J}_{n}\left(q\right),h\right]\right]\right\}=0.$$

Therefore,

$$\widetilde{J}_{n}\left(q\right)=\left(J_{\lambda,\beta}^{\left(n\right)}(q)\right)^{\prime}\in H_{0}^{4}\left(\Omega\right)$$

(6.12)

is the Fréchet derivative of the functional $J_{\lambda,\beta}^{\left(n\right)}(q):\overline{B(R,g_{0})}\rightarrow\mathbb{R}$ at the point $q\in\overline{B(R,g_{0})},$ and the right hand side of (6.4) indeed represents $\left(J_{\lambda,\beta}^{\left(n\right)}(q)\right)^{\prime}\left(h\right).$ Estimate (6.5) obviously follows from (6.4).

We now prove strict convexity estimate (6.7). To do this, we apply Carleman estimate (5.1) to the third line of (6.11). We obtain

$$J_{\lambda,\beta}^{(n)}(q+h)-J_{\lambda,\beta}^{(n)}(q)-\left(J_{\lambda,\beta}^{\left(n\right)}(q)\right)^{\prime}\left(h\right)\\ \geq C\mathop{\displaystyle\int}\limits_{\Omega}\big{(}\lambda\left(h_{x}^{2}+h_{t}^{2}\right)+\lambda^{3}h^{2}\big{)}e^{2\lambda\varphi}dxdt+C\mathop{\displaystyle\int}\limits_{\varepsilon}^{b}\big{(}\lambda h_{x}^{2}(x,0)+\lambda^{3}h^{2}(x,0)\big{)}e^{-2\lambda x}dx\\ -C\mathop{\displaystyle\int}\limits_{\varepsilon}^{b}\left(\lambda h_{x}^{2}(x,T)+\lambda^{3}h^{2}(x,T)\right)e^{-2\lambda(x+\alpha T)}dx+\beta\left\|h\right\|_{H^{4}\left(\Omega\right)}^{2},\text{ }\forall h\in H_{0}^{4}\left(\Omega\right).$$

(6.13)

By trace theorem there exists a constant $C_{1}=C_{1}\left(\Omega\right)>0$ depending only on the domain $\Omega$ such that

$$\left\|v\right\|_{H^{4}\left(\Omega\right)}^{2}\geq C_{1}\left\|v\left(x,T\right)\right\|_{H^{1}\left(\Omega\right)}^{2},\text{ }\forall v\in H^{4}\left(\Omega\right).$$

Since the regularization parameter $\beta\in\left[2e^{-\lambda\alpha T},1\right),$ then we can choose $\lambda_{1}$ so large that

$$\frac{\beta}{2}C_{1}\geq C_{1}e^{-\lambda\alpha T}>Ce^{-2\lambda(x+\alpha T)},\text{ }\forall\lambda\geq\lambda_{1},\forall x\in\left[\varepsilon,b\right].$$

Hence, for these values of $\lambda,$ the expression in the last line of (6.13) can be estimated from the below as:

$$-C\mathop{\displaystyle\int}\limits_{\varepsilon}^{b}\left(\lambda h_{x}^{2}(x,T)+\lambda^{3}h^{2}(x,T)\right)e^{-2\lambda(x+\alpha T)}dx+\beta\left\|h\right\|_{H^{4}\left(\Omega\right)}^{2}\geq\frac{\beta}{2}\left\|h\right\|_{H^{4}\left(\Omega\right)}^{2},\text{ }\forall h\in H_{0}^{4}\left(\Omega\right).$$

(6.14)

Hence, (6.13) and (6.14) imply

$$J_{\lambda,\beta}^{(n)}(q+h)-J_{\lambda,\beta}^{(n)}(q)-\left(J_{\lambda,\beta}^{\left(n\right)}(q)\right)^{\prime}\left(h\right)\geq C\mathop{\displaystyle\int}\limits_{\Omega}\big{(}\lambda\left(h_{x}^{2}+h_{t}^{2}\right)+\lambda^{3}h^{2}\big{)}e^{2\lambda\varphi}dxdt\\ +C\mathop{\displaystyle\int}\limits_{\varepsilon}^{b}\big{(}\lambda h_{x}^{2}(x,0)+\lambda^{3}h^{2}(x,0)\big{)}e^{-2\lambda x}dx+\frac{\beta}{2}\left\|h\right\|_{H^{4}\left(\Omega\right)}^{2},\forall\lambda\geq\lambda_{1}.$$

This proves (6.7). The existence and uniqueness of the minimizer $q_{\min,n}\in\overline{B(R,g_{0})}$ and inequality (6.9) follow from (6.7) as well as from a combination of Lemma 2.1 and Theorem 2.1 of [1], also see [42, Chapter 10, section 3]. $\square$