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

Thuy T. Le11footnotemark: 1 Micheal V. Klibanov Department of Mathematics and
Statistics, University of North Carolina at Charlotte, Charlotte, NC, 28223,
USA<EMAIL_ADDRESS><EMAIL_ADDRESS>(corresponding author)
<EMAIL_ADDRESS>Loc H. Nguyen11footnotemark: 1 Anders Sullivan US Army
Research Laboratory, 2800 Powder Mill Road, Adelphi, MD 20783-1197, USA,
<EMAIL_ADDRESS><EMAIL_ADDRESS>Lam
Nguyen22footnotemark: 2

###### Abstract

It is shown that the contraction mapping principle with the involvement of a
Carleman Weight Function works for a Coefficient Inverse Problem for a 1D
hyperbolic equation. Using a Carleman estimate, the global convergence of the
corresponding numerical method is established. Numerical studies for both
computationally simulated and experimentally collected data are presented. The
experimental part is concerned with the problem of computing dielectric
constants of explosive-like targets in the standoff mode using severely
underdetermined data.

Key words: one-dimensional wave equation, Carleman estimate, iterative method,
contraction principle, numerical solution, experimental data, global
convergence.

AMS subject classification: 35R25, 35R30

##  1 Introduction

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.

##  2 Statements of Forward and Inverse Problems

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._

##  3 A Boundary Value Problem for a Nonlinear PDE With Non-Local Terms

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)

##  4 Numerical Method for Problem 3.1

###  4.1 The function $q_{0}\left(x,t\right)$

We now find the first approximation $q_{0}(x,t)$ for the function $q(x,t).$
Using (2.2), we choose $c\left(x\right)\equiv 1$ as the first guess for the
function $c\left(x\right).$ Hence, by (3.3),

$q_{0}(x,0)\equiv\frac{1}{2}.$ (4.1)

We now need to find the function $q_{0}(x,t).$ To do this, drop the nonlinear
third term in the left hand side of equation (3.20) and, using (4.1) and
(3.17)_,_ set $1/\left(2q^{0}(x,0)\right)^{2}:=2$. Then (3.20), (3.21)__
become:

$q_{0xx}(x,t)-2q_{0xt}(x,t)=0\text{ in }\Omega,$ (4.2)
$q_{0}(\varepsilon,t)=g_{0}(t+\varepsilon),\text{
}q_{0x}(\varepsilon,t)=2g_{0}^{\prime}(t+\varepsilon),q_{0x}(b,t)=0.$ (4.3)

BVP (4.2), (4.3) has overdetermined boundary conditions (4.3). Typically, QRM
works well for BVPs for PDEs with overdetermined boundary conditions [19, 31].
Therefore, we solve BVP (4.2), (4.3) via CQRM. This means that we consider the
following minimization problem:

Minimization Problem Number $0$. _Assuming ( 3.12), minimize the functional
_$J_{\lambda,\beta}^{(0)}:\overline{B(R,g_{0})}\rightarrow\mathbb{R}$_on the
set_ ,__

$J_{\lambda,\beta}^{(0)}(q_{0})=\mathop{\displaystyle\int}\limits_{\Omega}\left(q_{0xx}-2q_{0xt}\right)^{2}e^{2\lambda\varphi}dxdt+\beta\|q_{0}\|_{H^{4}(\Omega)}^{2},$
(4.4)

_where_ $e^{2\lambda\varphi}$_is the Carleman Weight Function for the
operator_ $\partial_{x}-2\partial_{t}$ [50, 51]

$e^{2\lambda\varphi}=e^{-2\lambda\left(x+\alpha t\right)},$ (4.5)

_where_ $\alpha\in\left(0,1/2\right)$_is the parameter, and_
$\beta\in\left(0,1\right)$ _is the regularization parameter. Both parameters
will be chosen later._

Theorem 6.1 guarantees that, for appropriate values of parameters
$\lambda,\beta,$ there exists unique minimizer $q_{0,\min}$__
$\overline{B(R,g_{0})}$ of the functional $J_{\lambda,\beta}^{(0)}(q_{0}).$

###  4.2 The function $q_{n}\left(x,t\right)$ for $n\geq 1$

Assume that functionals
$J_{\lambda,\beta}^{(m)}(q_{m}):\overline{B(R,g_{0})}\rightarrow\mathbb{R}$
are defined for $m=0,...,n-1$ and their minimizers functions
$\left\\{q_{m,\min}\right\\}_{m=0}^{n-1}\subset\overline{B(R,g_{0})}$ are
constructed already, all for the same values of parameters
$\lambda,\alpha,\beta$. Replace in (3.20) $q\left(x,t\right)$ with
$q_{n}\left(x,t\right),q^{0}(x,0)$ with
$q_{n-1,\min}^{0}(x,0),q_{x}\left(x,t\right)$ with
$\partial_{x}q_{\left(n-1\right),\min}\left(x,t\right)$ and
$q_{t}\left(x,t\right)$ with
$\partial_{t}q_{\left(n-1\right),\min}\left(x,t\right).$ Then problem (3.20),
(3.21) becomes a linear one with respect to the function
$q_{n}\left(x,t\right),$

$q_{nxx}(x,t)-\frac{q_{nxt}(x,t)}{2\left(q_{n-1,\min}^{0}(x,0)\right)^{2}}+\frac{\partial_{t}q_{\left(n-1\right),\min}\left(x,t\right)\partial_{x}q_{\left(n-1\right)\min}(x,0)}{2\left(q_{n-1,\min}^{0}(x,0)\right)^{3}}=0\text{
in }\Omega,$ (4.6) $q_{n}(\varepsilon,t)=g_{0}(t+\varepsilon),\text{
}q_{nx}(\varepsilon,t)=2g_{0}^{\prime}(t+\varepsilon),q_{nx}(b,t)=0.$ (4.7)

To solve problem (4.6), (4.7), we consider the following minimization problem:

Minimization Problem Number $n$. _Assuming ( 3.12), minimize the functional
_$J_{\lambda,\beta}^{(n)}:H^{4}\left(\Omega\right)\rightarrow\mathbb{R}$_on
the set_ $B(R,g_{0}),$

$J_{\lambda,\beta}^{(n)}(q_{n})=\int_{\Omega}\Big{(}q_{nxx}(x,t)-\frac{q_{nxt}(x,t)}{2\big{(}q_{n-1,\min}^{0}(x,0)\big{)}^{2}}+\frac{\partial_{t}q_{\left(n-1\right),\min}\left(x,t\right)\partial_{x}q_{\left(n-1\right),\min}(x,0)}{2\big{(}q_{n-1,\min}^{0}(x,0)\big{)}^{3}}\Big{)}^{2}e^{2\lambda\varphi}dxdt$
$+\beta\left\|q_{n}\right\|_{H^{4}\left(\Omega\right)}^{2}.$ (4.8)

Suppose that there exists unique minimizer
$q_{n,\min}\in\overline{B(R,g_{0})}$ of the functional
$J_{\lambda,\beta}^{(n)}(q_{n}).$ Then, following (3.22), (3.17) and (3.22),
we set

$c_{n}\left(x\right)=\frac{1}{\left(2q_{n,\min}^{0}\left(x,0\right)\right)^{4}}.$
(4.9)

The rest of the analytical part of this paper is devoted to the convergence
analysis of the iterative numerical method presented in this section.

##  5 The Carleman Estimate

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._

##  6 Strict Convexity of Functional (4.8) on the Set $B(R,g_{0})$, Existence
and Uniqueness of Its Minimizer

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$

Remark 6.2. _Since the functional_ $J_{\lambda,\beta}^{\left(n\right)}(q)$_is
quadratic, then its strict convexity on the whole space_
$H^{4}\left(\Omega\right)$_follows immediately from the presence of the
regularization term_
$\beta\left\|q_{n}\right\|_{H^{4}\left(\Omega\right)}^{2}$_in it. However, in
addition to the claim of its strict convexity, we actually need in our
convergence analysis below those terms in the right hand side of the strict
convexity estimate ( 6.7), which are different from the term
_$\beta\left\|h\right\|_{H^{4}\left(\Omega\right)}^{2}/2.$_These terms are
provided by Carleman estimate ( 6.7). The condition
_$\beta\in\left[2e^{-\lambda\alpha T},1\right)$_of Theorem 6.1 is imposed to
dominate the negative term in the last line of ( 6.7)._

##  7 How to Find the Minimizer

Since we search for the minimizer $q_{n,\min}$ of functional (4.8) on the
bounded set $\overline{B(R,g_{0})}$ rather than on the whole space
$H^{4}\left(\Omega\right),$ then we cannot just use Riesz theorem to find this
minimizer, also see [1] and [42, Chapter 10, section 3] for the case of
finding a minimizer of a strictly convex functional on a bounded set. Two ways
of finding the minimizer $q_{n,\min}\in\overline{B(R,g_{0})}$ are described in
this section.

###  7.1 Gradient descent method

Keeping in mind (3.12), choose an arbitrary function $q_{0,n}\in B(R,g_{0}).$
We arrange the gradient descent method of the minimization of functional (4.8)
as follows:

$q_{k,n}=q_{\left(k-1\right),n}-\eta\left(J_{\lambda,\beta}^{\left(n\right)}(q_{\left(k-1\right),n})\right)^{\prime},\text{
}k=1,2,...,$ (7.1)

where $\eta\in\left(0,1\right)$ is a small number, which is chosen later. It
is important to note that since functions
$J_{\lambda,\beta}^{\left(n\right)}(q_{\left(k-1\right),n})\in
H_{0}^{4}\left(\Omega\right),$ then boundary conditions (4.7) are kept the
same for all functions $q_{k,n}\left(x,t\right).$ Also, using (3.17) and
(4.9), we set

$c_{k,n}\left(x\right)=\frac{1}{\left(2q_{k,n}^{0}\left(x,0\right)\right)^{4}},\text{
}x\in\left[\varepsilon,b\right],$ (7.2)
$c_{n,\min}\left(x\right)=\frac{1}{\left(2q_{n,\min}^{0}\left(x,0\right)\right)^{4}},\text{
}x\in\left[\varepsilon,b\right].$ (7.3)

Theorem 7.1 claims the global convergence of the gradient descent method (7.1)
in the case when $q_{\min,n}\in B(R/3,g_{0}).$

Theorem 7.1. _Let the number_
$\lambda_{1}=\lambda_{1}\left(R,\Omega,\alpha\right)\geq\lambda_{0}>1$_be the
one defined in ( 6.6). Let _$\lambda\geq\lambda_{1}$_. For this value of_
$\lambda,$_let_ $q_{n,\min}\in\overline{B\left(R,g_{0}\right)}$_be the unique
minimizer of the functional_ $J_{\lambda,\beta}^{\left(n\right)}(q_{n})$_on
the set_ $\overline{B\left(R,g_{0}\right)}$_with the regularization parameter_
$\beta\in\left[2e^{-\lambda\alpha T},1\right)$_(Theorem 6.1)_. _Assume that
the function_ $q_{n,\min}\in B(R/3,g_{0}).$_For each_ $n,$ _choose the
starting point of the gradient descent method ( 7.1) as _$q_{0,n}\in
B(R/3,g_{0}).$ _Then, there exists a number_
$\eta_{0}\in\left(0,1\right)$_such that for any_ $\eta\in(0,\eta_{0})$_all
functions_ $q_{k,n}\in B(R,g_{0}).$_Furthermore, there exists a number_
$\theta=\theta\left(\eta\right)\in(0,1)$_such that the following convergence
estimates are valid:_

$\left\|q_{k,n}-q_{n,\min}\right\|_{H^{4}\left(\Omega\right)}\leq\theta^{k}\left\|q_{0,n}-q_{n,\min}\right\|_{H^{4}\left(\Omega\right)},k=1,...,$
(7.4) $\left\|c_{k,n}-c_{n,\min}\right\|_{H^{3}\left(\varepsilon,b\right)}\leq
C\theta^{k}\left\|q_{0,n}-q_{n,\min}\right\|_{H^{4}\left(\Omega\right)},\text{
}k=1,...$ (7.5)

Proof. Estimate (7.4) follows immediately from [33, Theorem 4.6] combined with
“corrections” of functions
$q_{k,n}\left(x,0\right),q_{n,\min}\left(x,0\right)$ via (3.17). Estimate
(7.5) follows immediately from trace theorem, (3.17) and (7.2)-(7.4).
$\square$

###  7.2 Gradient projection method

Suppose now that there is no information on whether or not the function __
$q_{n,\min}\in B(R/3,g_{0}).$ In this case we construct the gradient
projection method. We introduce the function $F\left(x,t\right)$ since it is
easy to construct the projection operator on a ball with the center at
$\left\\{0\right\\}.$

Suppose that there exists a function $F\in H^{4}\left(\Omega\right)$ such that

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

This function exists, for example, in the case when the function $g_{0}(t)\in
H^{5}\left(0,T+\varepsilon\right).$ Indeed, as one of examples of this
function (among many), one can take, e.g.

$F\left(x,t\right)=\chi\left(x\right)\left(g_{0}(t+\varepsilon)+2xg_{0}^{\prime}(t+\varepsilon)\right),$

where the function $\chi\left(x\right)$ is such that

$\chi\left(x\right)\in
C^{4}\left[\varepsilon,b\right],\chi\left(x\right)=\left\\{\begin{array}[]{c}1,x\in\left[\varepsilon,b/4\right],\\\
0,x\in\left[b/2,b\right],\\\ \text{between }0\text{ and }1\text{ for
}x\in\left(b/4,b/2\right).\end{array}\right.$

The existence of such functions $\chi\left(x\right)$ is well known from the
Analysis course. Denote

$p_{n}\left(x,t\right)=q_{n}\left(x,t\right)-F(x,t).$ (7.7)

Then (4.6), (4.7) become:

$p_{nxx}(x,t)-p_{nxt}(x,t)\frac{1}{2\left(q_{n-1,\min}^{0}(x,0)\right)^{2}}$
(7.8)
$+F_{xx}(x,t)-F_{xt}(x,t)\frac{1}{2\left(q_{n-1,\min}^{0}(x,0)\right)^{2}}+\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}}=0\text{
in }\Omega,$ $p_{n}(\varepsilon,t)=0,\text{
}p_{nx}(\varepsilon,t)=0,p_{nx}(b,t)=0.$ (7.9)

Assume that

$\left\|F\right\|_{H^{4}\left(\Omega\right)}<R.$ (7.10)

By (7.7), (7.9), (7.10) and triangle inequality

$p_{n}\in B_{0}\left(2R\right)=\left\\{p\in
H_{0}^{4}\left(\Omega\right):\left\|p\right\|_{H_{0}^{4}\left(\Omega\right)}<2R\right\\}.$
(7.11)

To find the function $p_{n}\in B_{0}\left(2R\right)$ satisfying conditions
(7.8), (7.9), we minimize the following functional
$I_{\lambda,\beta}^{\left(n\right)}(p_{n}):\overline{B_{0}\left(2R\right)}\rightarrow\mathbb{R}$

$I_{\lambda,\beta}^{\left(n\right)}(p_{n})=J_{\lambda,\beta}^{\left(n\right)}(p_{n}+F),\text{
}p_{n}\in B_{0}\left(2R\right).$ (7.12)

Remark 7.1. _An obvious analog of Theorem 6.1 is valid of course for
functional ( 7.12). But in this case we should have instead of (6.6)
_$\lambda\geq\widetilde{\lambda}_{1}=\lambda_{1}\left(2R,\Omega,\alpha\right)\geq\lambda_{0}>1.$
_In particular, there exists unique minimizer_
$p_{n,\min}\in\overline{B_{0}\left(2R\right)}$_of this functional on the
closed ball_ $\overline{B_{0}\left(2R\right)}.$_We omit the formulation of
this theorem, since it is an obvious reformulation of Theorem 6.1._

Let
$P_{\overline{B_{0}\left(2R\right)}}:H_{0}^{4}\left(\Omega\right)\rightarrow\overline{B_{0}\left(2R\right)}$
be the projection operator of the space $H_{0}^{4}\left(\Omega\right)$ on the
closed ball $\overline{B_{0}\left(2R\right)}\subset
H_{0}^{4}\left(\Omega\right).$ Then this operator can be easily constructed:

$P_{\overline{B_{0}\left(2R\right)}}\left(p\right)=\left\\{\begin{array}[]{c}p\text{
if }p\in\overline{B_{0}\left(2R\right)},\\\
\frac{2R}{\left\|p\right\|_{H^{4}\left(\Omega\right)}}p\text{ if
}p\notin\overline{B_{0}\left(2R\right)}.\end{array}\right.$

We now construct the gradient projection method of the minimization of the
functional $I_{\lambda,\beta}^{\left(n\right)}(p_{n})$ on the set
$\overline{B_{0}\left(2R\right)}.$ Let $p_{0,n}\in B_{0}\left(2R\right)$ be an
arbitrary function. Then the sequence of the gradient projection method is:

$p_{k,n}=P_{\overline{B_{0}\left(2R\right)}}\left(p_{\left(k-1\right),n}-\eta\left(I_{\lambda,\beta}^{\left(n\right)}(p_{\left(k-1\right),n})\right)^{\prime}\right),\text{
}k=1,2,...,$ (7.13)

where $\eta\in\left(0,1\right)$ is a small number, which is chosen later.
Using (7.2), (7.3) and (7.7), we set

$\widetilde{c}_{k,n}\left(x\right)=\frac{1}{\left(2\widetilde{q}_{k,n}^{0}\left(x,0\right)\right)^{4}},\text{
}x\in\left[\varepsilon,b\right],$ (7.14)
$\widetilde{c}_{n,\min}\left(x\right)=\frac{1}{\left(2\widetilde{q}_{n,\min}^{0}\left(x,0\right)\right)^{4}},\text{
}x\in\left[\varepsilon,b\right].$ (7.15)

Here the function $\widetilde{q}_{n,k}^{0}\left(x,0\right)$ is obtained as
follows: First, we consider the function
$\left(p_{n,k}+F\right)\left(x,0\right).$ Next, we apply procedure (3.17) to
this function. Similarly for $\widetilde{q}_{n,\min}^{0}\left(x,0\right).$

Denote $p_{n,\min}\in\overline{B_{0}\left(2R\right)}$ the unique minimizer of
functional (7.12) on the set $\overline{B_{0}\left(2R\right)}$ (Remark 7.1).
Following (7.7), denote

$\widetilde{q}_{n,\min}=p_{n,\min}+F,\text{ }\widetilde{q}_{k,n}=p_{k,n}+F.$
(7.16)

We omit the proof of Theorem 7.2 since it is very similar with the proof of
Theorem 7.1. The only difference is that instead of Theorem 4.6 of [33] one
should use Theorem 4.1 of [30].

Theorem 7.2. _Let ( 7.6), (7.10) hold. Let the number
_$\lambda_{1}=\lambda_{1}\left(R,\Omega,\alpha\right)\geq\lambda_{0}>1$_be the
one defined in ( 6.6) and let_

$\lambda\geq\widetilde{\lambda}_{1}=\lambda_{1}\left(2R,\Omega,\alpha\right)\geq\lambda_{1}\left(R,\Omega,\alpha\right).$

_For this value of_ $\lambda$ _and for_ $\beta\in\left[2e^{-\lambda\alpha
T},1\right)$_let_ $p_{n,\min}\in\overline{B_{0}\left(2R\right)}$_be the unique
minimizer of functional ( 7.12) on the set _$\overline{B_{0}\left(2R\right)}$
_(Remark 7.1). Let notations ( 7.14) hold. For each _$n,$ _choose the starting
point of the gradient projection method ( 7.1) as _$p_{0,n}\in B_{0}(2R).$
_Then there exists a number_ $\eta_{0}\in\left(0,1\right)$_such that for any_
$\eta\in(0,\eta_{0})$_there exists a number_
$\theta=\theta\left(\eta\right)\in(0,1)$_such that the following convergence
estimates are valid for the iterative process ( 7.13):_

$\left\|\widetilde{q}_{k,n}-\widetilde{q}_{n,\min}\right\|_{H^{4}\left(\Omega\right)}\leq\theta^{k}\left\|\widetilde{q}_{0,n}-\widetilde{q}_{n,\min}\right\|_{H^{4}\left(\Omega\right)},$
$\left\|\widetilde{c}_{k,n}-\widetilde{c}_{n,\min}\right\|_{H^{4}\left(\Omega\right)}\leq
C\theta^{k}\left\|\widetilde{q}_{0,n}-\widetilde{q}_{n,\min}\right\|_{H^{4}\left(\Omega\right)},$

_where functions_
$\widetilde{c}_{k,n},\widetilde{c}_{n,\min},\widetilde{q}_{k,n},\widetilde{q}_{n,\min}$_are
defined in ( 7.14)-(7.16)._

Remark 7.2. _Both theorems 7.1 and 7.2 claim the global convergence of
corresponding minimization procedures. This is because the starting function
in both cases is an arbitrary one either in_ $B\left(R/3,g_{0}\right)$_or in_
$B_{0}\left(2R\right)$_and a smallness condition is not imposed on the number_
$R$_._

##  8 Contraction Mapping and Global Convergence

In this section, we prove the global convergence of the numerical method of
section 4 for solving Problem 3.1. To do this, we introduce first the exact
solution of our CIP. Recall that the concept of the existence of the exact
solution is one of the main concepts of the theory of ill-posed problems [4,
52]. In particular, an estimate in our global convergence theorem is very
similar to the one in contraction mapping.

Suppose that there exists a function $c^{\ast}\left(x\right)$ satisfying
conditions (2.1), (2.2). Let $u^{\ast}\left(x,t\right)$ be the solution of
problem (2.3), (2.4) with $c:=c^{\ast}.$ We assume that the corresponding data
$g_{0}^{\ast}\left(t\right),\partial_{t}g_{0}^{\ast}\left(t\right)$ for the
CIP are noiseless, see (2.11), (2.12). Let $q^{\ast}\left(x,t\right)$ be the
function $q^{\ast}\left(x,t\right)$ which is constructed from the function
$u^{\ast}\left(x,t\right)$ as in (3.1). Following (3.11), we assume that

$q^{\ast}\in B(R,g_{0}^{\ast}),$ (8.1)
$B(R,g_{0}^{\ast})=\left\\{\begin{array}[]{c}q\in
H^{4}(\Omega):\left\|q\right\|_{H^{4}(\Omega)}<R,\\\
q(\varepsilon,t+\varepsilon)=g_{0}^{\ast}(t+\varepsilon),q_{x}(\varepsilon,t)=2\left(g_{0}^{\ast}\right)^{\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\\}.$ (8.2)

By (3.20)-(3.21)

$q_{xx}^{\ast}(x,t)-q_{xt}^{\ast}(x,t)\frac{1}{2\left(q^{\ast}(x,0)\right)^{2}}+q_{t}^{\ast}(x,t)\frac{q_{x}^{\ast}(x,0)}{2\left(q^{\ast}(x,0)\right)^{3}}=0\text{
in }\Omega,$ (8.3) $q^{\ast}(\varepsilon,t)=g_{0}^{\ast}(t+\varepsilon),\text{
}q_{x}^{\ast}(\varepsilon,t)=2\partial_{t}g_{0}^{\ast}\left(t+\varepsilon\right),q_{x}^{\ast}(b,t)=0.$
(8.4)

By (3.3)

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

It is important in the formulation of Theorem 6.1 that both functions $q$ and
$q+h$ should have the same boundary conditions as prescribed in $B(R,g_{0}).$
However, boundary conditions for functions $q_{n}$ and $q^{\ast}$ are
different. Hence, similarly to subsection 7.2, we consider a function
$F^{\ast}\in H^{4}\left(\Omega\right)$ such that (see (7.6))

$F^{\ast}(\varepsilon,t)=g_{0}^{\ast}(t+\varepsilon),\text{
}F_{x}^{\ast}(\varepsilon,t)=2\partial_{t}g_{0}^{\ast}(t+\varepsilon),\text{
}F_{x}^{\ast}(b,t)=0.$ (8.6)

We assume similarly to (7.10) that

$\left\|F^{\ast}\right\|_{H^{4}\left(\Omega\right)}<R.$ (8.7)

Also, similarly to (7.7), we introduce the function $p^{\ast}\left(x,t\right)$
as:

$p^{\ast}\left(x,t\right)=q^{\ast}\left(x,t\right)-F^{\ast}(x,t).$ (8.8)

Let $\lambda_{1}$ be the number defined in (6.6), let $\lambda\geq\lambda_{1}$
and the function $q_{n,\min}\in\overline{B\left(R,g_{0}\right)}$ (see (6.8))
be the unique minimizer of the functional
$J_{\lambda,\beta}^{\left(n\right)}(q_{n})$ on the set
$\overline{B\left(R,g_{0}\right)},$ the existence of which is established in
Theorem 6.1. Following (7.7), denote

$p_{n,\min}\left(x,t\right)=q_{n,\min}\left(x,t\right)-F\left(x,t\right).$
(8.9)

By (6.8), (7.11), (8.1), (8.2) and (8.7)-(8.9)

$p_{n,\min},p^{\ast}\in\overline{B_{0}\left(2R\right)}.$ (8.10)

Also, it follows from embedding theorem, (3.13), (6.8), (7.10), (8.7) and
(8.10) that

$\left\|q^{\ast}\right\|_{C^{2}\left(\overline{\Omega}\right)},\left\|p^{\ast}\right\|_{C^{2}\left(\overline{\Omega}\right)},\left\|q_{n,\min}\right\|_{{}_{C^{2}\left(\overline{\Omega}\right)}},\left\|p_{n,\min}\right\|_{{}_{C^{2}\left(\overline{\Omega}\right)}}\leq
C.$ (8.11)

We assume that the data $g_{0},g_{0}^{\prime}$ for our CIP are given with
noise of the level $\delta,$ where the number $\delta>0$ is sufficiently
small. More precisely, we assume that

$\left\|F-F^{\ast}\right\|_{H^{4}\left(\Omega\right)}<\delta.$ (8.12)

Observe that (3.17), (8.1) and (8.2) imply that

$\left|q_{n-1,\min}^{0}\left(x,0\right)-q^{\ast}\left(x,0\right)\right|\leq\left|q_{n-1,\min}\left(x,0\right)-q^{\ast}\left(x,0\right)\right|,\text{
}x\in\left[\varepsilon,b\right].$ (8.13)

By (8.3), (8.4), (8.6) and (8.8)

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

in $\Omega\times[0,T]$ and

$p^{\ast}(\varepsilon,t)=p_{x}^{\ast}(\varepsilon,t)=p_{x}^{\ast}(b,t)=0.$
(8.15)

Theorem 8.1 (contraction mapping and global convergence of the method of
section 3). _Let functions_ $F$_,_ $F^{\ast}\in
H^{4}\left(\Omega\right)$_satisfy conditions ( 7.6), (7.10), (8.6), (8.7) and
(8.12). Let a sufficiently large number
_$\lambda_{1}=\lambda_{1}\left(R,\Omega,\alpha\right)\geq\lambda_{0}>1$_be the
one defined in ( 6.6). Let_

$\lambda\geq\widetilde{\lambda}_{1}=\lambda_{1}\left(2R,\Omega,\alpha\right)\geq\lambda_{1}\left(R,\Omega,\alpha\right).$
(8.16)

_For this value of_ $\lambda,$_let_
$q_{n,\min}\in\overline{B\left(R,g_{0}\right)}$_be the unique minimizer of the
functional_ $J_{\lambda,\beta}^{\left(n\right)}(q_{n})$_on the set_
$\overline{B\left(R,g_{0}\right)}$_with the regularization parameter_
$\beta\in\left[2e^{-\lambda\alpha T},1\right)$_(Theorem 6.1)_._Let_

$\overline{q}_{n}=q_{n,\min}-q^{\ast},\text{
}\overline{c}_{n}=c_{n,\min}-c^{\ast},\emph{\ }$ (8.17)

_where_ $c_{n,\min}$_is defined in ( 7.3). Then the following convergence
estimate holds_

$\mathop{\displaystyle\int}\limits_{\Omega}\left(\overline{q}_{nx}^{2}+\overline{q}_{nt}^{2}+\overline{q}_{n}^{2}\right)\left(x,t\right)e^{2\lambda\varphi}dxdt+\mathop{\displaystyle\int}\limits_{\varepsilon}^{b}\left(\overline{q}_{nx}^{2}+\overline{q}_{n}^{2}\right)\left(x,0\right)e^{-2\lambda
x}dx$
$\leq\frac{C}{\lambda}\mathop{\displaystyle\int}\limits_{\Omega}\left(\overline{q}_{\left(n-1\right)x}^{2}+\overline{q}_{\left(n-1\right)t}^{2}+\overline{q}_{n-1}^{2}\right)\left(x,t\right)e^{2\lambda\varphi}dxdt$
(8.18)
$+\frac{C}{\lambda}\mathop{\displaystyle\int}\limits_{\varepsilon}^{b}\left(\overline{q}_{\left(n-1\right)x}^{2}\left(x,0\right)+\overline{q}_{n-1}^{2}\right)\left(x,0\right)e^{-2\lambda
x}dx+C\left(\delta^{2}+\beta\right),$

_which leads to:_

$\mathop{\displaystyle\int}\limits_{\Omega}\left(\overline{q}_{nx}^{2}+\overline{q}_{nt}^{2}+\overline{q}_{n}^{2}\right)\left(x,t\right)e^{2\lambda\varphi}dxdt$
(8.19)
$\leq\frac{C^{n}}{\lambda^{n}}\mathop{\displaystyle\int}\limits_{\Omega}\left(\overline{q}_{0x}^{2}+\overline{q}_{0t}^{2}+\overline{q}_{0}^{2}\right)\left(x,t\right)e^{2\lambda\varphi}dxdt+C\left(\delta^{2}+\beta\right).$

_In addition,_

$\left\|\overline{c}_{n}\right\|_{H^{1}\left(\varepsilon,b\right)}^{2}\leq\frac{C^{n}}{\lambda^{n}}\mathop{\displaystyle\int}\limits_{\Omega}\left(\overline{q}_{0x}^{2}+\overline{q}_{0t}^{2}+\overline{q}_{0}^{2}\right)\left(x,t\right)e^{2\lambda\varphi}dxdt+C\left(\delta^{2}+\beta\right).$
(8.20)

Remarks 8.1:

1\. _Due to the presence of the term_ $C/\lambda$_with a sufficiently large_
$\lambda$_, estimate ( 8.18) is similar to the one in the classic contraction
mapping principle, although we do not claim here the existence of the fixed
point. This explains the title of our paper. _

_2\. When computing the unique minimizer_ $q_{0,\min}$_of functional ( 4.4) on
the set _$\overline{B(R,g_{0})},$_we do not impose a smallness condition on
the number_ $R$_. Therefore, Theorem 8.1 claims the global convergence of the
method of section 4, see Introduction for our definition of the global
convergence. The same is true for Theorems 9.1 and 9.2 in the next section._

Proof Theorem 8.1. Denote $h_{n}=p^{\ast}-p_{n,\min}.$ By (8.9)
$h_{n}=-\overline{q}_{n}+\left(F-F^{\ast}\right).$ Hence, (8.12) and embedding
theorem imply:

$h_{n}^{2}\geq\frac{1}{2}\overline{q}_{n}^{2}-C\delta^{2},h_{nx}^{2}\geq\frac{1}{2}\overline{q}_{nx}^{2}-C\delta^{2},h_{nt}^{2}\geq\frac{1}{2}\overline{q}_{nt}^{2}-C\delta^{2}\text{
in }\overline{\Omega},$ (8.21) $h_{n}^{2}+h_{nx}^{2}+h_{nt}^{2}\leq
C\left(\overline{q}_{n}^{2}+\overline{q}_{nx}^{2}+\overline{q}_{nt}^{2}+\delta^{2}\right)\text{
in }\overline{\Omega}.$ (8.22)

Consider the functional
$I_{\lambda,\beta}^{\left(n\right)}(p^{\ast})=J_{\lambda,\beta}^{\left(n\right)}(p^{\ast}+F).$
Since both functions $p_{n,\min}$ and $p^{\ast}$ have the same zero boundary
conditions (7.9) and since by (8.10) both of them belong to the set
$\overline{B_{0}\left(2R\right)}$, then the analog of Theorem 6.1, which is
mentioned in Remark 7.1, implies (see (6.7))

$I_{\lambda,\beta}^{\left(n\right)}(p^{\ast})-I_{\lambda,\beta}^{\left(n\right)}(p_{n,\min})-\left(I_{\lambda,\beta}^{\left(n\right)}(p_{n,\min})\right)^{\prime}(h_{n})$
$\geq
C\mathop{\displaystyle\int}\limits_{\Omega}\Big{[}\lambda\left(h_{nx}^{2}+h_{nt}^{2}\right)+\lambda^{3}h_{n}^{2}\Big{]}e^{2\lambda\varphi}dxdt$
(8.23) $+C\mathop{\displaystyle\int}\limits_{\varepsilon}^{b}\big{(}\lambda
h_{nx}^{2}(x,0)+\lambda^{3}h_{n}^{2}(x,0)\big{)}e^{-2\lambda
x}dx+\frac{\beta}{2}\left\|h_{n}\right\|_{H^{4}\left(\Omega\right)}^{2},\text{
}\forall\lambda\geq\widetilde{\lambda}_{1}.$

By (6.9)

$-\left(I_{\lambda,\beta}^{\left(n\right)}(p_{n,\min})\right)^{\prime}(h_{n})\leq
0.$

Hence, the left hand side of (8.23) can be estimated as:

$I_{\lambda,\beta}^{\left(n\right)}(p^{\ast})-I_{\lambda,\beta}^{\left(n\right)}(p_{n,\min})-\left(I_{\lambda,\beta}^{\left(n\right)}(p_{n,\min})\right)^{\prime}(h)\leq
I_{\lambda,\beta}^{\left(n\right)}(p^{\ast}).$ (8.24)

We now estimate the right hand side of (8.24) from the above. It follows from
(4.8), (7.12) and (8.3) that

$I_{\lambda,\beta}^{(n)}(p^{\ast})=\mathop{\displaystyle\int}\limits_{\Omega}G_{n}^{2}e^{2\lambda\varphi}dxdt+\beta\left\|p^{\ast}+F\right\|_{H^{4}\left(\Omega\right)}^{2},$
(8.25)

where

$G_{n}=p_{xx}^{\ast}(x,t)-p_{xt}^{\ast}(x,t)\frac{1}{2\left(q_{\left(n-1\right),\min}^{0}(x,0)\right)^{2}}+\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_{\left(n-1\right),\min}^{0}(x,0)\right)^{3}}$
$+F_{xx}-F_{xt}(x,t)\frac{1}{2\left(q_{\left(n-1\right),\min}^{0}(x,0)\right)^{2}}$
$=q_{xx}^{\ast}(x,t)-q_{xt}^{\ast}(x,t)\frac{1}{2\left(q^{\ast}(x,0)\right)^{2}}+q_{t}^{\ast}(x,t)\frac{q_{x}^{\ast}(x,0)}{2\left(q^{\ast}(x,0)\right)^{3}}$
(8.26)
$+\left(F-F^{\ast}\right)_{xx}-\left(F-F^{\ast}\right)_{xt}\frac{1}{2\left(q_{\left(n-1\right),\min}^{0}(x,0)\right)^{2}}$
$-q_{xt}^{\ast}\left(\frac{1}{2\left(q_{\left(n-1\right),\min}^{0}(x,0)\right)^{2}}-\frac{1}{2\left(q^{\ast}(x,0)\right)^{2}}\right)$
$+\left(\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_{\left(n-1\right),\min}^{0}(x,0)\right)^{3}}-q_{t}^{\ast}(x,t)\frac{q_{x}^{\ast}(x,0)}{2\left(q^{\ast}(x,0)\right)^{3}}\right).$

By (8.3), the third line of (8.26) equals zero. Hence, (8.26) becomes

$G_{n}=\left(F-F^{\ast}\right)_{xx}+\left(F-F^{\ast}\right)_{xt}\frac{1}{2\left(q_{\left(n-1\right),\min}^{0}(x,0)\right)^{2}}$
$-q_{xt}^{\ast}\left(\frac{1}{2\left(q_{\left(n-1\right),\min}^{0}(x,0)\right)^{2}}-\frac{1}{2\left(q^{\ast}(x,0)\right)^{2}}\right)$
$+\left(\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_{\left(n-1\right),\min}^{0}(x,0)\right)^{3}}-q_{t}^{\ast}(x,t)\frac{q_{x}^{\ast}(x,0)}{2\left(q^{\ast}(x,0)\right)^{3}}\right).$

Hence, by (3.17), (8.12) and embedding theorem

$\left|G_{n}\right|\leq
C\delta+C\left|\frac{1}{2\left(q_{\left(n-1\right),\min}^{0}(x,0)\right)^{2}}-\frac{1}{2\left(q^{\ast}(x,0)\right)^{2}}\right|$
$+\left|\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_{\left(n-1\right),\min}^{0}(x,0)\right)^{3}}-q_{t}^{\ast}(x,t)\frac{q_{x}^{\ast}(x,0)}{2\left(q^{\ast}(x,0)\right)^{3}}\right|.$
(8.27)

Using (8.12) and (8.13), we obtain

$\left|\frac{1}{2\left(q_{\left(n-1\right),\min}^{0}(x,0)\right)^{2}}-\frac{1}{2\left(q^{\ast}(x,0)\right)^{2}}\right|$
$=\frac{\left|\left(q_{\left(n-1\right),\min}^{0}\left(x,0\right)-q^{\ast}(x,0)\right)-\left(F-F^{\ast}\right)\left(x,0\right)\right|\left|q_{\left(n-1\right),\min}^{0}(x,0)+q^{\ast}(x,0)\right|}{2\left(q_{\left(n-1\right),\min}^{0}(x,0)\right)^{2}\left(q^{\ast}(x,0)\right)^{2}.}$
$\leq C\delta+C\left|h_{n-1}\left(x,0\right)\right|.$

Combining this with (8.27), we obtain

$\left|G_{n}\right|\leq C\delta+C\left|h_{n-1}\left(x,0\right)\right|$
$+\left|\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_{\left(n-1\right),\min}^{0}(x,0)\right)^{3}}-q_{t}^{\ast}(x,t)\frac{q_{x}^{\ast}(x,0)}{2\left(q^{\ast}(x,0)\right)^{3}}\right|.$
(8.28)

Next,

$\frac{1}{2\left(q_{\left(n-1\right),\min}^{0}(x,0)\right)^{3}}=\frac{1}{2\left(q^{\ast}(x,0)\right)^{3}}$
$+\left(\frac{1}{2\left(q_{\left(n-1\right),\min}^{0}\left(x,0\right)\right)^{3}}-\frac{1}{2\left(q^{\ast}(x,0)\right)^{3}}\right)$
$=\frac{1}{2\left(q^{\ast}(x,0)\right)^{3}}+S_{n-1}\left(x,0\right)\left[\left(q_{\left(n-1\right),\min}^{0}\left(x,0\right)-q^{\ast}(x,0)\right)-\left(F-F^{\ast}\right)\left(x,0\right)\right],$

where the function $S_{n-1}\left(x,0\right)$ can be estimated as

$\left|S_{n-1}\left(x,0\right)\right|\leq C.$ (8.29)

Hence,

$\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_{\left(n-1\right),\min}^{0}(x,0)\right)^{3}}=$
$\frac{\partial_{t}q_{\left(n-1\right),\min}\left(x,t\right)\partial_{x}q_{\left(n-1\right),\min}(x,0)}{2\left(q^{\ast}(x,0)\right)^{3}}+$
$+\left[\partial_{t}q_{\left(n-1\right),\min}\left(x,t\right)\partial_{x}q_{\left(n-1\right),\min}(x,0)\right]\times$
$\times
S_{n-1}\left(x,0\right)\left[\left(q_{\left(n-1\right),\min}^{0}\left(x,0\right)-q^{\ast}(x,0)\right)-\left(F-F^{\ast}\right)\left(x,0\right)\right].$

Hence, using (3.11), (3.17), (8.1), (8.13) and (8.29), we obtain

$\left|\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_{\left(n-1\right),\min}^{0}(x,0)\right)^{3}}-q_{t}^{\ast}(x,t)\frac{q_{x}^{\ast}(x,0)}{2\left(q^{\ast}(x,0)\right)^{3}}\right|$
$\leq\frac{1}{2\left(q^{\ast}(x,0)\right)^{3}}\left|\partial_{t}q_{\left(n-1\right),\min}\left(x,t\right)\partial_{x}q_{\left(n-1\right),\min}(x,0)-q_{t}^{\ast}(x,t)q_{x}^{\ast}(x,0)\right|$
(8.30) $+C\left(\delta+\left|h_{n-1}\left(x,0\right)\right|\right)$ $\leq
C\left(\delta+\left|h_{n-1}\left(x,0\right)\right|+\left|\partial_{x}h_{n-1,x}\left(x,0\right)\right|+\left|h_{n-1,t}\left(x,t\right)\right|\right).$

Combining (8.28) with (8.30), we obtain

$\left|G_{n}\left(x,t\right)\right|\leq
C\left(\delta+\left|h_{n-1}\left(x,0\right)\right|+\left|h_{n-1,x}\left(x,0\right)\right|+\left|h_{n-1,t}\left(x,t\right)\right|\right).$

Hence, by (8.25)

$I_{\lambda,\beta}^{(n)}(p^{\ast})\leq
C\mathop{\displaystyle\int}\limits_{\Omega}\left(\delta^{2}+h_{n-1}^{2}\left(x,0\right)+h_{n-1,x}^{2}\left(x,0\right)+h_{n-1,t}^{2}\left(x,t\right)\right)e^{2\lambda\varphi}dxdt+C\beta.$

Substituting this in (8.24) and then using (8.23), we obtain

$\mathop{\displaystyle\int}\limits_{\Omega}\left(h_{nx}^{2}+h_{nt}^{2}+h_{n}^{2}\right)e^{2\lambda\varphi}dxdt+\mathop{\displaystyle\int}\limits_{\varepsilon}^{b}\left(h_{nx}^{2}(x,0)+h_{n}^{2}(x,0)\right)e^{-2\lambda
x}dx$
$\leq\frac{C}{\lambda}\mathop{\displaystyle\int}\limits_{\Omega}\left(\delta^{2}+h_{n-1}^{2}\left(x,0\right)+h_{n-1,x}^{2}\left(x,0\right)+h_{n-1,t}^{2}\left(x,t\right)\right)e^{2\lambda\varphi}dxdt+C\beta.$
(8.31)

Obviously

$\mathop{\displaystyle\int}\limits_{\Omega}\left(h_{n-1}^{2}\left(x,0\right)+h_{n-1,x}^{2}\left(x,0\right)\right)e^{2\lambda\varphi}dxdt\leq\frac{1}{2\lambda\alpha}\mathop{\displaystyle\int}\limits_{\varepsilon}^{b}\left(h_{\left(n-1\right)x}^{2}(x,0)+h_{n-1}^{2}(x,0)\right)e^{-2\lambda
x}dx,$ (8.32)
$\mathop{\displaystyle\int}\limits_{\Omega}\delta^{2}e^{2\lambda\varphi}dxdt\leq
C\frac{\delta^{2}}{\lambda^{2}}.$ (8.33)

Denote

$y_{n}=\mathop{\displaystyle\int}\limits_{\Omega}\left(h_{nx}^{2}+h_{nt}^{2}+h_{n}^{2}\right)e^{2\lambda\varphi}dxdt+\mathop{\displaystyle\int}\limits_{\varepsilon}^{b}\left(h_{nx}^{2}(x,0)+h_{n}^{2}(x,0)\right)e^{-2\lambda
x}dx.$ (8.34)

Then (8.31)-(8.33) imply

$y_{n}\leq\frac{C}{\lambda}y_{n-1}+C\left(\frac{\delta^{2}}{\lambda^{2}}+\beta\right).$
(8.35)

Iterating (8.35) with respect to $n,$ we obtain

$\mathop{\displaystyle\int}\limits_{\Omega}\left(h_{nx}^{2}+h_{nt}^{2}+h_{n}^{2}\right)e^{2\lambda\varphi}dxdt+\mathop{\displaystyle\int}\limits_{\varepsilon}^{b}\left(h_{nx}^{2}(x,0)+h_{n}^{2}(x,0)\right)e^{-2\lambda
x}dx$
$\leq\frac{C}{\lambda^{n}}\left[\mathop{\displaystyle\int}\limits_{\Omega}\left(h_{0x}^{2}+h_{0t}^{2}+h_{0}^{2}\right)e^{2\lambda\varphi}dxdt+\mathop{\displaystyle\int}\limits_{\varepsilon}^{b}\left(h_{0x}^{2}(x,0)+h_{0}^{2}(x,0)\right)e^{-2\lambda
x}dx\right]$ (8.36) $+C\left(\frac{\delta^{2}}{\lambda^{2}}+\beta\right).$

Apply (8.21) to the left hand side of estimate (8.36). Also, apply (8.22) at
$n=0$ to the right hand side of (8.36). We obtain (8.19). Estimate (8.20)
follows from an obvious combination of (8.19) with (7.3), (8.5), (8.13) and
(LABEL:8.150). Finally, estimate (8.18) follows immediately from (8.21),
(8.22) and (8.31). $\square$

##  9 Global Convergence of the Gradient and Gradient Projection Methods to
the Exact Solution

First, we consider the gradient method of the minimization of functionals
$J_{\lambda,\beta}^{\left(n\right)}(q_{\left(k-1\right),n})$ on the set
$\overline{B(R,g_{0})},$ see (7.1). The proof of Theorem 9.1 follows
immediately from the triangle inequality combined with Theorems 7.1 and 8.1.

Theorem 9.1. _Let_ $\alpha_{0}$_and_ $\lambda_{0}$_be the numbers of Theorem
5.1._ _Let the sufficiently large number_
$\lambda_{1}=\lambda_{1}\left(R,\Omega,\alpha\right)\geq\lambda_{0}>1$_be the
one defined in ( 6.6). Let the number _$\widetilde{\lambda}_{1}$_be the same
as in ( 8.16),_

$\widetilde{\lambda}_{1}=\lambda_{1}\left(2R,\Omega,\alpha\right)\geq\lambda_{1}\left(R,\Omega,\alpha\right).$

_Let_ $\lambda\geq\lambda_{1}$_and let the regularization parameter_
$\beta\in\left[2e^{-\lambda\alpha T},1\right).$_Assume that the functions_
$q_{\min,n}\in B(R/3,g_{0})$_for all_ $n.$_For each_ $n$, _choose the starting
point of the gradient method ( 7.1) as _$q_{0,n}\in B(R/3,g_{0}).$ _Then there
exists a number_ $\eta_{0}\in\left(0,1\right)$_such that for any_
$\eta\in(0,\eta_{0})$_all functions_ $q_{k,n}\in B(R,g_{0}),k=1,...;$
$n=1,....$_Furthermore, there exists a number_
$\theta=\theta\left(\eta\right)\in(0,1)$_such that the following convergence
estimate is valid:_

$\left\|c_{k,n}-c^{\ast}\right\|_{H^{1}\left(\varepsilon,b\right)}\leq
C\theta^{k}\left\|q_{0,n}-q_{n,\min}\right\|_{H^{4}\left(\Omega\right)}$
$+\frac{C^{n/2}}{\lambda^{n/2}}\left(\mathop{\displaystyle\int}\limits_{\Omega}\left(\overline{q}_{0x}^{2}+\overline{q}_{0t}^{2}+\overline{q}_{0}^{2}\right)\left(x,t\right)e^{2\lambda\varphi}dxdt\right)^{1/2}+C\left(\delta+\sqrt{\beta}\right),$

_where the function_ $\overline{q}_{0}$_is defined in ( 8.17)._

Consider now the gradient projection method of the minimization of the
functionals
$I_{\lambda,\beta}^{\left(n\right)}(p_{n})=J_{\lambda,\beta}^{\left(n\right)}(p_{n}+F)$
in (7.12) on the set $\overline{B_{0}\left(2R\right)}$, see (7.13).__ We use
notations (7.14). Theorem 9.2 follows immediately from the triangle inequality
combined with Theorems 7.2 and 8.1.

Theorem 9.2. _Let the number_
$\lambda_{1}=\lambda_{1}\left(R,\Omega,\alpha\right)\geq\lambda_{0}>1$_be the
one defined in ( 6.6). Let the number _$\widetilde{\lambda}_{1}$_be the same
as in ( 8.16),_

$\widetilde{\lambda}_{1}=\lambda_{1}\left(2R,\Omega,\alpha\right)\geq\lambda_{1}\left(R,\Omega,\alpha\right).$

_Let_ $\lambda\geq\lambda_{1}$_and let the regularization parameter_
$\beta\in\left[2e^{-\lambda\alpha T},1\right).$_Consider the gradient
projection method ( 7.13). For each_ $n$, _choose the starting point_
$p_{0,n}$ _of this method as an arbitrary point of the ball_ $B_{0}(2R).$
_Then there exists a number_ $\eta_{0}\in\left(0,1\right)$_such that for any_
$\eta\in(0,\eta_{0})$_there exists a number_
$\theta=\theta\left(\eta\right)\in(0,1)$_such that the following convergence
estimate is valid:_

$\left\|\widetilde{c}_{k,n}-c^{\ast}\right\|_{H^{1}\left(\varepsilon,b\right)}\leq
C\theta^{k}\left\|\widetilde{q}_{0,n}-\widetilde{q}_{n,\min}\right\|_{H^{4}\left(\Omega\right)}$
$+\frac{C^{n/2}}{\lambda^{n/2}}\left(\mathop{\displaystyle\int}\limits_{\Omega}\left(\overline{q}_{0x}^{2}+\overline{q}_{0t}^{2}+\overline{q}_{0}^{2}\right)\left(x,t\right)e^{2\lambda\varphi}dxdt\right)^{1/2}+C\left(\delta+\sqrt{\beta}\right),$

_where functions_
$\widetilde{c}_{k,n},\widetilde{q}_{0,n},\widetilde{q}_{n,\min}$_and_
$\overline{q}_{0}$_are defined in ( 7.14), (7.16) and (8.17) respectively. _

##  10 Numerical Studies

###  10.1 Numerical implementation

To generate the simulated data, by Lemma 2.2, we solve problem (2.3), (2.4)
for the case when the whole real line is replaced by a large interval $(-a,a)$
with the absorbing boundary conditions (2.9)–(2.10). More precisely, just as
in Section 6.1 of [30], we use the implicit scheme to numerically solve

$\left\\{\begin{array}[]{rcll}c(x)u_{tt}(x,t)&=&u_{xx}(x,t)&(x,t)\in(-a,a)\times(0,T),\\\
u(-a,t)-u_{x}(-a,t)&=&0&t\in(0,T),\\\ u(a,t)+u_{x}(a,t)&=&0&t\in(0,T),\\\
u(x,0)&=&0&x\in\mathbb{R},\\\
u_{t}(x,0)&=&\widetilde{\delta}(x)&x\in\mathbb{R},\end{array}\right.$ (10.1)

where $a=5$, $T=6$ and

$\widetilde{\delta}(x)=\frac{30}{\sqrt{2\pi}}e^{-\frac{(30x)^{2}}{2}}$

is a smooth approximation of the Dirac function $\delta\left(x\right)$. We
solve problem (10.1) by the implicit finite difference method. In the finite
difference scheme, we arrange a uniform partition for the interval $[-a,a]$ as
$\\{y_{0}=-a,y_{1},\dots,y_{N}=a\\}\subset[-a,a]$ with $y_{i}=a+2ia/N_{x}$,
$i=0,\dots,N_{x}$, where $N_{x}$ is a large number. In the time domain, we
split the interval $[0,T]$ into $N_{t}+1$ uniform sub-intervals
$[t_{j},t_{j+1}]$, $j=0,\dots,N_{t}$, with $t_{j}=jT/N_{t},$ where $N_{t}$ is
a large number. In our computational setting, $N_{x}=3001$ and $N_{t}=301$.
These numbers are the same as in [30].

We observe a computational error for the function $u$ near $(x=0,t=0)$. This
is due to the fact that the function $\widetilde{\delta}\left(x\right)$ is not
exactly equal to the Dirac function. We correct the error as follows. It
follows from (2.2), (2.7) and (2.8) that $u(x,t)=1/2$ in a neighborhood of the
point $\left(x,t\right)=\left(0,0\right).$ We, therefore, replace the data
$u(\epsilon,t)$ by $1/2$ when $|t|$ is small. In our computation, we set
$u(x,t)=1/2$ for $(x,t)\in[0,0.0067]\times[0,0.26]$. This data correction step
is exactly the same as in Section 6.1 of [30] and is illustrated by Figure 2
in that publication. Then, we can extract the noiseless data $g_{0}^{\ast}$
easily. We next add the noise into the data via the formula

$g_{0}=g_{0}^{\ast}(1+\delta\mbox{rand})$ (10.2)

where $\delta$ is the noise level and rand is the function that generates
uniformly distributed random numbers in the range $[-1,1].$ In all numerical
tests with simulated data below, the noise level $\delta=0.05,$ i.e. $5\%.$
Due to (2.12), the function $g_{1}=g_{0}^{\prime}.$ Due to the presence of
noise, see (10.2), we cannot compute $g_{1}=g_{0}^{\prime}$ by the finite
difference method. Hence, the function $g_{0}^{\prime}$ is computed by the
Tikhonov regularization method. The version of the Tikhonov regularization
method for this problem is well-known. Hence, we do not describe this step
here.

Having the data for the function $q$ in hand, we proceed as in Algorithm 1.

Algorithm 1 A numerical method to solve Problem 3.1

1: Choose a set of parameters $\lambda$, $\alpha$ and $\beta$.

2: Compute the function $q_{0}$ by minimizing the functional
$J^{(0)}_{\lambda,\beta}$ defined in (4.4). Due to (3.3), the initial
reconstruction is given by
$c_{\mathrm{i}nit}(x)=\frac{1}{(2q_{0}(x,0))^{4}}\quad\mbox{for all
}x\in[\epsilon,b].$

3: Assume that the function $q_{n-1}$ is known. We compute the function
$q_{n}$ by minimizing the function $J^{(n)}_{\lambda,\beta}$ defined in (4.8).

4: Set $q_{\mathrm{c}omp}=q_{n}$ when $n=n^{*}$ is large enough.

5: Due to (3.3), the function $c_{\mathrm{c}omp}$ is set to be
$c_{\mathrm{c}omp}(x)=\frac{1}{(2q_{\mathrm{c}omp}(x,0))^{4}}\quad\mbox{for
all }x\in[\epsilon,b].$

In step 1 of Algorithm 1, we choose $\lambda=2$, $\alpha=0.3$ and
$\beta=10^{-11}.$ These parameters were chosen by a trial-error process that
is similar to the one in [50]. Just as in [50], we choose a reference
numerical test in which we know the true solution. In fact, test 1 of section
10.2 was our reference test. We tried several values of $\lambda,$ $\alpha,$
and $\beta$ until the numerical solution to that reference test became
acceptable. Then, we have used the same values of these parameters for all
other tests, including all five (5) available cases of experimental data.

We next implement Steps 2 and 3 of Algorithm 1. We write differential
operators in the functionals $J_{\lambda,\beta}^{(0)}$ and
$J_{\lambda,\beta}^{(n)}$ in the finite differences with the step size in
space $\Delta x=0.0033$ and the step size in time $\Delta t=0.02$ and minimize
resulting functionals with respect to values of corresponding functions at
grid points. Since the integrand in the definitions of the functional
$J_{\lambda,\beta}^{(n)},n=0,1,...$ is the square of linear functions, then
its minimizer is its critical point. In finite differences, we can write a
linear system whose solution is the desired critical point. We solve this
system by the command “lsqlin” of Matlab. The details in implementation by the
finite difference method including the discretization, the derivation of the
linear system for the critical point, and the use of “lsqlin” are very similar
to the scheme in [43, 47]. Recall that in our theory, in the definition of the
functional $J_{\lambda,\beta}^{(n)}$ acting on $q_{n}$, see (4.8), we replaced
$q_{n-1}$ with its analog $q_{n-1,\mathrm{min}}$ which belongs to the bounded
set $\overline{B(R,g_{0})},$ and also replaced $q_{n-1,\mathrm{min}}(x,0)$
with $q_{n-1,\mathrm{min}}^{0}(x,0)$. These replacements are only for the
theoretical purpose to avoid the case when $q_{n-1}$ blows up. However, in the
numerical studies, these steps can be relaxed. This means that in Step 3, we
have minimized the finite difference analog of the functional

$q_{n}\mapsto\int_{\Omega}\Big{(}q_{nxx}(x,t)-\frac{q_{nxt}(x,t)}{2\big{(}q_{n-1}(x,0)\big{)}^{2}}+\frac{\partial_{t}q_{n-1}\left(x,t\right)\partial_{x}q_{n-1}(x,0)}{2\big{(}q_{n-1}(x,0)\big{)}^{3}}\Big{)}^{2}e^{2\lambda\varphi}dxdt+\beta\left\|q_{n}\right\|_{H^{2}\left(\Omega\right)}^{2},$
(10.3)

subject to the boundary conditions in lines 2 and 3 of (3.11). Another
numerical simplification is that rather than using the $H^{4}-$norm in the
regularization term, we use the $H^{2}-$ norm in (10.3). Although the
theoretical analysis supporting the above simplifications is missing, we did
not experience any difficulties in computing the numerical solutions to
Problem 3.1. All of our numerical results are satisfactory.

###  10.2 Numerical results for computationally simulated data

To test Algorithm 1, we present four (4) numerical examples.

Test 1 (the reference test). We first test the case of one inclusion with a
high inclusion/background contrast. The true dielectric constant function
$c(x)$ has the following form

$c_{\mathrm{true}}(x)=\left\\{\begin{array}[]{ll}1+14e^{\frac{(x-1)^{2}}{(x-1)^{2}-0.2^{2}}}&\mbox{if
}|x-1|<0.2,\\\ 1&\mbox{otherwise}.\end{array}\right.$ (10.4)

(a)

(b)

Figure 1: Test 1. The true and reconstructed functions $c(x),$ where
$c_{\mathrm{true}}$ is given in (10.4). (a) The functions $c_{\mathrm{init}}$
and $c_{\mathrm{comp}}$ are obtained by Step 2 and Step 5 of Algorithm 1
respectively. (b) The consecutive relative error
$\|c_{n}-c_{n-1}\|_{L^{\infty}(\epsilon,M)}/\|c_{n}\|_{L^{\infty}(\epsilon,M)}$,
$n=1,\dots,10.$ The data is with $\delta=5\%$ noise.

Thus, the inclusion/background contrast as 15:1. It is evident from Figure 1
that we can successfully detect an object. The diameter of this object is 0.4
and the distance between this object and the source is 1. The true
inclusion/background contrast is 15/1. The maximal value of the computed
dielectric constant is 15.28. The relative error in this maximal value is
1.89% while the noise level in the data is 5%. Although the contrast is high,
our method provides good numerical solution without requiring any knowledge of
the true solution. Our method converges fast. Although the initial
reconstruction $c_{\mathrm{init}}$ computed in Step 2 of Algorithm 1 is not
very good, see Figure 1a, one can see in Figure 1b, that the convergence
occurs after only five (5) iterations. This fact verifies numerically the
“contraction property” of Theorem 8.1 including the key estimate (8.20).

Test 2. The true function $c$ in this test has two “inclusions”,

$c_{\mathrm{true}}(x)=\left\\{\begin{array}[]{ll}1+5e^{\frac{(x-0.6)^{2}}{(x-0.6)^{2}-0.2^{2}}}&\mbox{
if }\left|x-0.6\right|<0.2,\\\
1+8e^{\frac{(x-1.4)^{2}}{(x-1.4)^{2}-0.3^{2}}}&\mbox{ if
}\left|x-1.4\right|<0.3,\\\ 1&\mbox{ otherwise. }\end{array}\right.$ (10.5)

Numerical results of this test are displayed in Figure 2.

(a)

(b)

Figure 2: Test 2. The true and reconstructed functions $c(x),$ where
$c_{\mathrm{true}}$ is given in (10.5). (a) The functions $c_{\mathrm{init}}$
and $c_{\mathrm{comp}}$ are obtained by Step 2 and Step 5 of Algorithm 1
respectively. (b) The consecutive relative error
$\|c_{n}-c_{n-1}\|_{L^{\infty}(\epsilon,M)}/\|c_{n}\|_{L^{\infty}(\epsilon,M)}$,
$n=1,\dots,10.$ The data is with $\delta=5\%$ noise.

Test 2 is more complicated than Test 1. However, we still obtain good
numerical results. It is evident from Figure 2a that we can precisely detect
the two inclusions without any initial guess. The true maximal values of the
dielectric constant of the inclusions in the left and the right are 6 and 9
respectively. The reconstructed maximal values in those inclusions are
acceptable. They are 5.31 (relative error 11.5%) and 7.8 (relative error
13.3%). As in Test 1, the initial reconstruction $c_{\mathrm{init}}$ computed
in Step 2 of Algorithm 1 is far from $c_{\mathrm{true}}(x).$ Still, our
iterative procedure converges after 7 iterations, see Figure 2b.

Test 3 We test the case when the function $c_{\mathrm{true}}(x)$ is
discontinuous. It is a piecewise constant function given by

$c_{\mathrm{true}}(x)=\left\\{\begin{array}[]{ll}10&\mbox{ if
}\left|x-1\right|<0.15,\\\ 1&\mbox{ otherwise. }\end{array}\right.$ (10.6)

The numerical solution for this test is presented in Figure 3.

(a)

(b)

Figure 3: Test 3. The true and reconstructed functions $c(x),$ where
$c_{\mathrm{true}}$ is given in (10.6). (a) The functions $c_{\mathrm{init}}$
and $c_{\mathrm{comp}}$ are obtained by Step 2 and Step 5 of Algorithm 1
respectively. (b) The consecutive relative error
$\|c_{n}-c_{n-1}\|_{L^{\infty}(\epsilon,M)}/\|c_{n}\|_{L^{\infty}(\epsilon,M)}$,
$n=1,\dots,10.$ The data is with $\delta=5\%$ noise.

Although the function $c_{\text{true}}$ is not smooth and actually not even
continuous, Algorithm 1 works and provides a reliable numerical solution. The
computed maximal value of the dielectric constant of the object is 9.25
(relative error 7.5%), which is acceptable. The location of the object is
detected precisely, see Figure 3a. As in the previous two tests, Algorithm 1
converges fast. After the fifth iteration, the reconstructed function $c_{n}$
becomes unchanged. Again, this fact numerically confirms both Theorem 8.1 and
the robustness of our method.

Test 4 In this test, the function $c_{\mathrm{true}}(x)$ has the following
form:

$c_{\mathrm{true}}(x)=\left\\{\begin{array}[]{ll}3.5+0.3\ast\sin(\pi(x-1.35))&\mbox{
if }\left|x-0.9\right|<0.5,\\\ 8&\mbox{ if }\left|x-2\right|<0.37,\\\ 1&\mbox{
otherwise. }\end{array}\right.$ (10.7)

This test is interesting since the graph of the function (10.7) consists of a
curve and an inclusion. The numerical solution for this case is presented in
Figure 4.

(a)

(b)

Figure 4: Test 4. The true and reconstructed functions $c(x),$ where
$c_{\mathrm{true}}$ is given in (10.7). (a) The functions $c_{\mathrm{init}}$
and $c_{\mathrm{comp}}$ are obtained by Step 2 and Step 5 of Algorithm 1
respectively. (b) The consecutive relative error
$\|c_{n}-c_{n-1}\|_{L^{\infty}(\epsilon,M)}/\|c_{n}\|_{L^{\infty}(\epsilon,M)}$,
$n=1,\dots,10.$ The data is with $\delta=5\%$ noise.

One can observe from Figure 4a that our method to compute the initial
reconstruction in Step 2 of Algorithm 1 is not very effective. However, after
only 6 iterations, good numerical results are obtained. The curve in the first
inclusion locally coincides with the true one and the maximal value of the
computed dielectric constant within inclusion is quite accurate: it is 7.83
(relative error 2.12%). Our method converges at the iteration number 6.

######  Remark 10.1.

It follows from all above tests that Algorithm 1 is robust in solving a highly
nonlinear and severely ill-posed Problem 3.1. It provides satisfactory
numerical results with a fast convergence. For each test, the computational
time to compute the numerical solution is about 29 seconds on a MacBook Pro
2019 with 2.6 GHz processor and 6 Intel i7 cores. This is almost a real time
computation.

##  11 Numerical Results for Experimental Data

The experimental data we use to test Algorithm 1 were collected by the Forward
Looking Radar built in the US Army Research Laboratory [46]. The data were
initially collected to detect and identify targets mimicking shallow anti-
personnel land mines and IEDs. Five mimicking devices were: a bush, a wood
stake, a metal box, a metal cylinder, and a plastic cylinder. The bush and the
wood stake were placed in the air, while the other three objects were buried
at a few centimeters depth in the ground. We refer the reader to [46] for
details about the experimental set up and the data collection. This was
recently repeated in [30]. We, therefore, do not describe here the device and
the data collection procedure again. Since the location of the targets can be
detected by the Ground Position System (GPS), we are only interested in
computing the values of the dielectric constants of those targets. We do so by
our new method written in Algorithm 1.

As in the earlier works using these data [13, 30, 23, 24, 25, 35, 36, 50], we
compute the function $c_{\mathrm{rel}}(x),$ where

$c_{\mathrm{rel}}(x)=\left\\{\begin{array}[]{ll}\frac{c_{\mathrm{target}}}{c_{\mathrm{bckgr}}}\left(x\right)\text{
if
}\max\frac{c_{\mathrm{target}}}{c_{\mathrm{bckgr}}}\left(x\right)>1&\text{and
}x\in D,\\\ 1&\text{otherwise},\end{array}\right.$ (11.1)
$c_{\mathrm{rel}}(x)=\left\\{\begin{array}[]{ll}\min\frac{c_{\mathrm{target}}}{c_{\mathrm{bckgr}}}\left(x\right)\text{
if }\max\frac{c_{\mathrm{target}}}{c_{\mathrm{bckgr}}}\left(x\right)\leq
1&\text{and }x\in D,\\\ 1&\text{otherwise},\end{array}\right.$ (11.2)

where $D$ is a sub interval of the interval $[\epsilon,M],$ which is occupied
by the target. Next, we define the computed value of $c_{\mathrm{target}}$ as
[23]:

$c_{\mathrm{comp}}=c_{\mathrm{bckgr}}\times\left\\{\begin{array}[]{c}\max
c_{\mathrm{rel}}(x)\text{ if }\max c_{\mathrm{rel}}(x)>1,\\\ \min
c_{\mathrm{rel}}(x)\text{ if }\max c_{\mathrm{rel}}(x)<1.\end{array}\right.$
(11.3)

As in the above cited publications, we have to preprocess the raw data of [46]
before importing them into our solver. The data preprocessing is exactly the
same as in [30, Section 7.1]. First, we observe that the $L_{\infty}-$norm of
the experimental data far exceeds the $L_{\infty}-$norm of the computationally
simulated data. Therefore, we need to scale the experimental data by dividing
it by a calibration factor $\mu>0,$ i.e. we replace the raw experimental data
$f_{\text{raw}}\left(t\right)$ with
$f_{\text{scale}}\left(t\right)=f_{\text{raw}}\left(t\right)/\mu$. Then we
work only with $f_{\text{scale}}\left(t\right).$ To find the calibration
factor, we use a trial-and-error process. First, we select a reference target.
We then try many values of $\mu$ such that the reconstruction of the reference
target is satisfactory, i.e. the computed dielectric constant is in its
published range, see below in this section about the published range. Then
this calibration factor is used is all other tests. When the object is in the
air, our reference target is bush. In this case, the calibration factor
$\mu_{\mathrm{air}}=534,592.$ When the object is buried under the ground, our
reference target is the metal box and the calibration factor was
$\mu_{\mathrm{ground}}=265,223$.

Next, we preprocess the data $f_{\text{scale}}\left(t\right)$ as follows.
First, we first use a lower envelop (built in Matlab) to bound the data. We
then truncate the data that is not in a small interval centered at the global
minimizer of the data, see [30, Section 7.1] for the choice of this small
interval. But in the case of the plastic cylinder we use the upper envelop.
The choice of the upper or lower envelopes is as follows. We look at the
function $f_{\text{scale}}\left(t\right)$ and find the three extrema with
largest absolute values. If the middle extremal value among these three is a
minimum, then we bound the data by a lower envelop. Otherwise, we use an upper
envelop. See [30, Section 7.1] for more details and the reason of this choice.
In particular, Figures 4a, 4b, 4d, 5a, 5b, 5d, 5e, 5g, 5h of [30] provide
illustrations. Likewise, our Figures 5 displays computed functions
$c_{\text{target}}\left(x\right)$ for our five targets. The computed
dielectric constants $c_{\mathrm{comp}}$ defined in (11.3) by Algorithm 1 is
listed in Table 1.

(a) $c_{\mathrm{t}arget}$ for bush

(b) $c_{\mathrm{t}arget}$ for wood stake

(c) $c_{\mathrm{t}arget}$ for metal box

(d) $c_{\mathrm{t}arget}$ for metal cylinder

(e) $c_{\mathrm{t}arget}$ for plastic cylinder

Figure 5: Computed functions $c_{\mathrm{target}}(x,y)$ for our five targets, also see (11.1)- (11.3) and Table 1. Target | $c_{\mathrm{bckgr}}$ | computed $c_{\mathrm{rel}}$ | $c_{\mathrm{bckgr}}$ | computed $c_{\text{target}}$ | True $c_{\text{target}}$
---|---|---|---|---|---
Bush | 1 | 7.62 | 1 | 10.99 | $[3,20]$
Wood stake | 1 | 2.01 | 1 | 2.26 | $[2,6]$
Metal box | 4 | 4.00 | $[3,5]$ | $[12.00,20.00]$ | $[10,30]$
Metal cylinder | 4 | 4.01 | $[3,5]$ | $[12.3,20.5]$ | $[10,30]$
Plastic cylinder | 4 | 0.59 | $[3,5]$ | $[1.6,2.95]$ | $\left[1.1,3.2\right]$
Table 1: Computed dielectric constants of five targets

The true values of dielectric constants of our targets were not measured in
the experiment. Therefore, we compare our computed values with the published
ones. The published values of the dielectric constants of our targets are
listed in the last column of Table 1. They can be found on the website of
Honeywell (Table of dielectric constants, https://goo.gl/kAxtzB). Also, see
[9] for the experimentally measured range of the dielectric constants of
vegetation samples, which we assume have the same range as the dielectric
constant of bush. In the table of dielectric constants of Honeywell as well as
in [9], any dielectric constant is not a number. Rather, each dielectric
constant of these references is given within a certain interval. As to the
metallic targets, it was established in [35] that they have the so-called
“apparent” dielectric constant whose values are in the interval
$\left[10,30\right].$

Conclusion. It is clear from Table 1 that our computed dielectric constants
are consistent with the intervals of published ones. Therefore, our results
for experimental data are satisfactory.

######  Remark 11.1 (Speed of computations).

Our experimental data are sparse. The size of the data in time is $N_{t}=80$,
which is a lot smaller than that in the dense simulated data ($N_{t}=300$).
Therefore, the speed of computations is much faster than for the case of
simulated data of section 10. Thus, all results of this section were computed
in real time on the same computer (MacBook Pro 2019 with 2.6 GHz processor and
6 Intel i7 cores).

##  12 Concluding Remarks

We have developed a new globally convergent numerical method for a 1-D
Coefficient Inverse Problem with backscattering data for the wave-like PDE
(2.3). This is the second generation of the above cited convexification method
of our research group. The main novelty here is that, rather than minimizing a
globally strictly convex weighted cost functional arising in the
convexification, we solve on each iterative step a linear boundary value
problem. This is done using the so-called Carleman Quasi-Reversibility Method.
Just like in the convexification, the key element of the technique of this
paper, on which our convergence analysis is based, is the presence of the
Carleman Weight Function in each quadratic functional which we minimize. The
convergence estimate is similar to the key estimate of the classical
contraction mapping principle. The latter explains the title of this paper. We
have proven a global convergence theorem of our method. Our numerical results
for computationally simulated data demonstrate high reconstruction accuracies
in the presence of 5% random noise in the data.

Furthermore, our numerical results for experimentally collected data have
satisfactory accuracy via providing values of computed dielectric constants of
explosive-like targets within their published ranges.

A practically important observation here is that our computations for
experimental data were performed in _real time_.

## Acknowledgement

The effort of TTL, MVK and LHN was supported by US Army Research Laboratory
and US Army Research Office grant W911NF-19-1-0044.

## References

  * [1] A. B. Bakushinskii, M. V. Klibanov, and N. A. Koshev.  Carleman weight functions for a globally convergent numerical method for ill-posed Cauchy problems for some quasilinear PDEs.  Nonlinear Anal. Real World Appl., 34:201–224, 2017.
  * [2] L. Baudouin, M. de Buhan, and S. Ervedoza.  Convergent algorithm based on Carleman estimates for the recovery of a potential in the wave equation.  SIAM J. Nummer. Anal., 55:1578–1613, 2017.
  * [3] L. Baudouin, M. de Buhan, S. Ervedoza, and A. Osses.  Carleman-based reconstruction algorithm for the waves.  _SIAM Journal on Numerical Analysis_ , 59(2):998–1039, 2021.
  * [4] L. Beilina and M. V. Klibanov.  Approximate Global Convergence and Adaptivity for Coefficient Inverse Problems.  Springer, New York, 2012.
  * [5] M. Bellassoued and M. Yamamoto.  Carleman Estimates and Applications to Inverse Problems for Hyperbolic Systems.  Springer, Japan, 2017.
  * [6] A. L. Bukhgeim and M. V. Klibanov.  Uniqueness in the large of a class of multidimensional inverse problems.  _Soviet Math. Doklady_ , 17:244–247, 1981.
  * [7] T. Carleman.  Sur un probleme d’unicit e pur les systeé mes d’equations aux dériveés partielles à deux varibales indépendantes.  Ark. Mat. Astr. Fys., 26B:1–9, 1939.
  * [8] G. Chavent.  Nonlinear Least Squares for Inverse Problems: Theoretical Foundations and Step-by-Step Guide for Applications, Scientic Computation.  Springer, New York, 2009.
  * [9] H. T. Chuah, K. Y. Lee, and T. W. Lau.  Dielectric constants of rubber and oil palm leaf samples at X-band.  IEEE Trans. Geosci. Remote Sens, 33, 221–223, 1995.
  * [10] M. de Hoop, P. Kepley, and L. Oksanen.  Recovery of a smooth metric via wave field and coordinate transformation reconstruction.  SIAM J. Appl. Math., 78:1931–1953, 2018.
  * [11] A. V. Goncharsky and S. Y. Romanov.  Iterative methods for solving coefficient inverse problems of wave tomography in models with attenuation.  Inverse Problems, 33:025003, 2017.
  * [12] V. Goncharsky and S. Y. Romanov.  A method of solving the coefficient inverse problems of wave tomography.  Comput. Math. Appl., 77:967–980, 2019.
  * [13] A.L. Karchevsky, M.V. Klibanov, L. Nguyen, N. Pantong, and A. Sullivan.  The Krein method and the globally convergent method for experimental data.  Applied Numerical Mathematics, 74:111–127, 2013.
  * [14] V. A. Khoa, G. W. Bidney, M. V. Klibanov, L. H. Nguyen, L. Nguyen, A. Sullivan, and V. N. Astratov.  Convexification and experimental data for a 3D inverse scattering problem with the moving point source.  Inverse Problems, 36:085007, 2020.
  * [15] V. A. Khoa, G. W. Bidney, M. V. Klibanov, L. H. Nguyen, L. Nguyen, A. Sullivan, and V. N. Astratov.  An inverse problem of a simultaneous reconstruction of the dielectric constant and conductivity from experimental backscattering data.  Inverse Problems in Science and Engineering, 29(5):712–735, 2021.
  * [16] V. A. Khoa, M. V. Klibanov, and L. H. Nguyen.  Convexification for a 3D inverse scattering problem with the moving point source.  SIAM J. Imaging Sci., 13(2):871–904, 2020.
  * [17] M. V. Klibanov.  Inverse problems and Carleman estimates.  Inverse Problems, 8:575–596, 1992.
  * [18] M. V. Klibanov.  Carleman estimates for global uniqueness, stability and numerical methods for coefficient inverse problems.  J. Inverse and Ill-Posed Problems, 21:477–560, 2013\.
  * [19] M. V. Klibanov.  Carleman estimates for the regularization of ill-posed Cauchy problems.  Applied Numerical Mathematics, 94:46–74, 2015.
  * [20] M. V. Klibanov and O. V. Ioussoupova.  Uniform strict convexity of a cost functional for three-dimensional inverse scattering problem.  SIAM J. Math. Anal., 26:147–179, 1995.
  * [21] M. V. Klibanov.  Global convexity in a three-dimensional inverse acoustic problem.  SIAM J. Math. Anal., 28:1371–1388, 1997.
  * [22] M. V. Klibanov and A. Timonov.  Carleman Estimates for Coefficient Inverse Problems and Numerical Applications.  Inverse and Ill-Posed Problems Series. VSP, Utrecht, 2004.
  * [23] M. V. Klibanov, L. H. Nguyen, A. Sullivan, and L. Nguyen.  A globally convergent numerical method for a 1-d inverse medium problem with experimental data.  Inverse Problems and Imaging, 10:1057–1085, 2016.
  * [24] M. V. Klibanov, A. E. Kolesov, L. Nguyen, and A. Sullivan.  Globally strictly convex cost functional for a 1-D inverse medium scattering problem with experimental data.  SIAM J. Appl. Math., 77:1733–1755, 2017.
  * [25] M. V. Klibanov, A. E. Kolesov, L. Nguyen, and A. Sullivan.  A new version of the convexification method for a 1-D coefficient inverse problem with experimental data.  Inverse Problems, 34:35005, 2018.
  * [26] M. V. Klibanov, A. E. Kolesov, and D-L. Nguyen.  Convexification method for an inverse scattering problem and its performance for experimental backscatter data for buried targets.  _SIAM J. Imaging Sci._ , 12:576–603, 2019.
  * [27] M. V. Klibanov, J. Li, and W. Zhang.  Convexification of electrical impedance tomography with restricted Dirichlet-to-Neumann map data.  Inverse Problems, 35:035005, 2019\.
  * [28] M. V. Klibanov, Z. Li, and W. Zhang.  Convexification for the inversion of a time dependent wave front in a heterogeneous medium.  SIAM J. Appl. Math., 79:1722–1747, 2019\.
  * [29] M. V. Klibanov, J. Li, and W. Zhang.  Convexification for an inverse parabolic problem.  Inverse Problems, 36:085008, 2020.
  * [30] M. V. Klibanov, T. T. Le, L. H. Nguyen, A. Sullivan, and L. Nguyen.  Convexification-based globally convergent numerical method for a 1D coefficient inverse problem with experimental data.  preprint, Arxiv:2104.11392, 2021.
  * [31] M. V. Klibanov and J. Li.  Inverse Problems and Carleman Estimates: Global Uniqueness, Global Convergence and Experimental Data.  , De Gruyter, 2021.
  * [32] M. V. Klibanov, A. V. Smirnov, V. A. Khoa, A. Sullivan, and L. Nguyen.  Through-the-wall nonlinear SAR imaging.  _IEEE Transactions on Geoscience and Remote Sensing,_ 59: 7474-7486, 2021\.
  * [33] M. V. Klibanov, V. A. Khoa, A. V. Smirnov, L. H. Nguyen, G. W. Bidney, L. Nguyen, A. Sullivan, and V. N. Astratov.  Convexification inversion method for nonlinear SAR imaging with experimentally collected data.  to appear in J. Applied and Industrial Mathematics, see also Arxiv:2103.10431, 2021.
  * [34] J. Korpela, M. Lassas, and L. Oksanen.  Regularization strategy for an inverse problem for a $1+1$ dimensional wave equation.  Inverse Problems, 32:065001, 2016.
  * [35] A. Kuzhuget, L. Beilina, M.V. Klibanov, A. Sullivan, Lam Nguyen, and M. A. Fiddy.  Blind backscattering experimental data collected in the field and an approximately globally convergent inverse algorithm.  Inverse Problems, 28:095007, 2012.
  * [36] A. V. Kuzhuget, L. Beilina, M. V. Klibanov, A. Sullivan, L. Nguyen, and M. A. Fiddy.  Quantitative image recovery from measured blind backscattered data using a globally convergent inverse method.  IEEE Trans. Geosci. Remote Sens, 51:2937–2948, 2013.
  * [37] R. Lattès and J. L. Lions.  The Method of Quasireversibility: Applications to Partial Differential Equations.  Elsevier, New York, 1969.
  * [38] M. M. Lavrent’ev, V. G. Romanov, and S. P. Shishat$\cdot$skiĭ.  Ill-Posed Problems of Mathematical Physics and Analysis.  Translations of Mathematical Monographs. AMS, Providence: RI, 1986.
  * [39] T. T. Le and L. H. Nguyen.  A convergent numerical method to recover the initial condition of nonlinear parabolic equations from lateral Cauchy data.  Journal of Inverse and Ill-posed Problems, DOI: https://doi.org/10.1515/jiip-2020-0028, 2020.
  * [40] T. T. Le and L. H. Nguyen.  The gradient descent method for the convexification to solve boundary value problems of quasi-linear PDEs and a coefficient inverse problem.  preprint Arxiv:2103.04159, 2021.
  * [41] B. M. Levitan.  Inverse Sturm–Liouville Problems.  VSP, Utrecht, 1987.
  * [42] M. Minoux.  Mathematical Programming: Theory and Algorithms.  John Wiley & Sons, New York, 1986.
  * [43] L. H. Nguyen.  A new algorithm to determine the creation or depletion term of parabolic equations from boundary measurements.  Computers and Mathematics with Applications, 80:2135–2149, 2020.
  * [44] L. H. Nguyen and M. V. Klibanov.  Carleman estimates and the contraction principle for an inverse source problem for nonlinear hyperbolic equations.  preprint arXiv:2108.03500, 2021.
  * [45] L. H. Nguyen and M. V. Klibanov.  Carleman estimates and the contraction principle for an inversesource problem for nonlinear hyperbolic equations.  preprint arXix:2108.03500v3, 2021.
  * [46] N. Nguyen, D. Wong, M. Ressler, F. Koenig, Stanton B., G. Smith, J. Sichina, and K. Kappra.  Obstacle avolidance and concealed target detection using the Army Research Lab ultra-wideband synchronous impulse Reconstruction (UWB SIRE) forward imaging radar.  Proc. SPIE, 6553:65530H (1)–65530H (8), 2007.
  * [47] P. M. Nguyen and L. H. Nguyen.  A numerical method for an inverse source problem for parabolic equations and its application to a coefficient inverse problem.  Journal of Inverse and Ill-posed Problems, 38:232–339, 2020.
  * [48] V.G. Romanov. _Inverse Problems of Mathematical Physics_ , VNU Press, 1986.
  * [49] J. A. Scales, M. L. Smith, and T. L. Fischer.  Global optimization methods for multimodal inverse problems.  J. Computational Physics, 103:258–268, 1992.
  * [50] A. V. Smirnov, M. V. Klibanov, and L. H. Nguyen.  Convexification for a 1D hyperbolic coefficient inverse problem with single measurement data.  Inverse Probl. Imaging, 14(5):913–938, 2020.
  * [51] A. V. Smirnov, M. V. Klibanov, A. Sullivan, and L. Nguyen.  Convexifcation for an inverse problem for a 1D wave equation with experimental data.  Inverse Problems, 36:095008, 2020.
  * [52] A. N. Tikhonov, A. Goncharsky, V. V. Stepanov, and A. G. Yagola.  Numerical Methods for the Solution of Ill-Posed Problems.  Kluwer Academic Publishers Group, Dordrecht, 1995.
  * [53] M. Yamamoto.  Carleman estimates for parabolic equations. Topical Review.  Inverse Problems, 25:123013, 2009\.