# A Carleman-Picard approach for reconstructing zero-order coefficients in
parabolic equations with limited data

Ray Abney Department of Mathematics and Statistics, University of North
Carolina at Charlotte, Charlotte, NC 28223, USA<EMAIL_ADDRESS>Thuy T.
Le Department of Mathematics, NC State University, Raleigh, NC 27695, USA,
<EMAIL_ADDRESS>Loc H. Nguyen Department of Mathematics and Statistics,
University of North Carolina at Charlotte, Charlotte, NC, 28223, USA,
<EMAIL_ADDRESS>Cam Peters Department of Mathematics and Statistics,
University of North Carolina at Charlotte, Charlotte, NC, 28223, USA,
<EMAIL_ADDRESS>
###### Abstract

We propose a globally convergent computational technique for the nonlinear
inverse problem of reconstructing the zero-order coefficient in a parabolic
equation using partial boundary data. This technique is called the “reduced
dimensional method”. Initially, we use the polynomial-exponential basis to
approximate the inverse problem as a system of 1D nonlinear equations. We then
employ a Picard iteration based on the quasi-reversibility method and a
Carleman weight function. We will rigorously prove that the sequence derived
from this iteration converges to the accurate solution for that 1D system
without requesting a good initial guess of the true solution. The key tool for
the proof is a Carleman estimate. We will also show some numerical examples.

Keywords: time reduction, Carleman Picard iteration, nonlinear, parabolic.

AMS subject classification:

##  1 Introduction

Let $d\geq 2$ be the spatial dimension. This paper aims to solve a coefficient
inverse problem for the following initial value problem

$\left\\{\begin{array}[]{rcll}u_{t}({\bf x},t)&=&\Delta u({\bf x},t)+c({\bf
x})u({\bf x},t)&({\bf x},t)\in\mathbb{R}^{d}\times(0,\infty),\\\ u({\bf
x},0)&=&p({\bf x})&{\bf x}\in\mathbb{R}^{d}.\end{array}\right.$ (1.1)

More precisely, we propose a globally convergent method to solve the following
inverse problem.

######  Problem 1.1.

Let $R$ and $T$ be two positive numbers. Define $\Omega=(-R,R)^{d}$, and

$\Gamma=\\{{\bf x}=({\bf x}^{\prime},z):{\bf
x}^{\prime}=(x_{1},\dots,x_{d-1})\mbox{ and },|x_{i}|<R,i=1,\dots,d-1,z=\pm
R\\}\subset\partial\Omega.$ (1.2)

Assume that $|p|>0$ in $\overline{\Omega}$. Given the boundary measurements

$f({\bf x},t)=u({\bf x},t)\mbox{ and }g({\bf x},t)=u_{z}({\bf x},t)$ (1.3)

for all $({\bf x},t)\in\Gamma\times(0,T)$, compute the coefficient $c({\bf
x}),$ for ${\bf x}\in\Omega.$

Problem 1.1 boasts countless real-world applications. Consider a scenario
where the internal points of the medium $\Omega$ remain inaccessible. By
recording partial boundary data of the function $u$, specifically the heat and
heat flux as discussed in this paper, over a designated time frame and by
resolving Problem 1.1, one can identify the coefficient $c({\bf x})$, ${\bf
x}\in\Omega$. This allows the examination of the medium without causing any
damage to it. An important example can be drawn from bioheat transfer, where
the coefficient $c({\bf x})$ signifies blood perfusion. Understanding this
coefficient is vital for determining the temperature of blood coursing through
tissue, as highlighted in [11]. However, the uniqueness of Problem 1.1,
especially when data collection is limited to a specific subset of
$\partial\Omega$, remains an open area and is explored within the reduced
dimensional framework of this paper. Variations of Problem 1.1, with some
internal data assumed to be known, have been addressed in [4, 9, 45].
Additionally, the uniqueness can be found in [14] when provided with the
Dirichlet to Neuman map. In this paper, the uniqueness of Problem 1.1 is
assumed. Another topic of interest is the inverse challenge of retrieving
other coefficients, such as diffusion or initial conditions, based on the
final time measurements or boundary measurements for parabolic equations. This
is an intriguing and critical issue, with theoretical findings and
computational methods elaborated in [1, 27, 30, 32, 36, 37, 46, 48].

Inverse problems of computing the coefficients for parabolic equations have
been extensively explored. To the authors’ knowledge, the widely-used
technique for resolving such issues is the optimal control approach; see the
important works [5, 10, 11, 16, 49] and other cited references. The
researchers in [5] employed the optimal control method with a preconditioner
to achieve high-quality numerical calculations of thermal conductivity.
However, a significant limitation of this technique is the necessity for a
reliable initial estimation of the true solution, which is not consistently
accessible. We would like to particularly highlight the convexification
method, as described in [2, 22, 25, 29]. This approach addresses the challenge
of obtaining an initial guess. The studies in[2, 22, 25, 29] suggested to
minimize some Carleman weighted strictly convex functionals. When minimized,
the minimizers of these functionals produce the solution to the problem at
hand. Other worthy mentions are [39] and [44], which respectively present
alternative approaches to address Problem 1.1 by iteratively solving a Picard-
like approximation and its linearization. The approaches above consider the
full boundary observation. Unlike this, our contribution is introducing a
fresh technique that does not rely on prior insights into the actual
coefficient and requests only partial observation.

Our approach to addressing Problem 1.1 is split into two phases. In the
initial phase, drawing inspiration from [39, 44], we eliminate the unknown
coefficient $c({\bf x})$ from (1.1). By this, we obtain a partial differential
equation. The equation that emerges from this phase is a complicated one,
which involves both nonlocal and nonlinear terms. On the other hand, the
boundary condition of the solution is only provided on
$\Gamma\subsetneq\partial\Omega$. As of now, there is no established numerical
method to address it. During the subsequent phase, we transform this equation
into a system of nonlinear ordinary differential equations. This
transformation is guided by truncating the Fourier series with respect to a
special basis introduced in paper [43]. This basis is named the polynomial-
exponential basis. It is the high-dimensional version of the 1D polynomial-
exponential basis originally introduced in [23]. We then deploy a predictor-
corrector strategy to solve this nonlinear system. Within this framework, the
preliminary approximation of the true solution is derived without any prior
understanding. Subsequently, the resolution to Problem 1.1 is achieved. The
corrector stage in this procedure is executed using the quasi-reversibility
method and a Carleman weight function. The quasi-reversibility method was
first introduced by Lattès and Lions in [28] for numerical solutions of ill-
posed problems for partial differential equations. It has been studied
intensively since then, see e.g., [3, 6, 7, 8, 12, 13, 15, 26, 20, 38, 42]. A
survey on this method can be found in [21]. A question arises immediately
whether or not the iteration led by the predictor-corrector procedure above
converges. In this paper, we will rigorously prove this important result. The
proof is motivated by the one in [30, 32, 40]. However, its advantage is that
we can relax a technical condition in those papers about the smoothness of the
noise. That means the noise model in this paper is more realistic than in the
earlier publications.

The paper is organized as follows. In Section 2, we introduce our
approximation dimensional model that leads to the dimensional reduction
approach. In Section 3, we establish a 1D Carleman estimate. Section 4 is for
the algorithm and the proof of its convergence. In Section 5, we present our
numerical study. Section 6 is for concluding remarks.

##  2 The reduced dimension model

Define

$v({\bf x},t)=u_{t}({\bf x},t)\quad\mbox{for all }({\bf
x},t)\in\Omega\times(0,T).$ (2.1)

Then, by differentiating both sides of the differential equation in (1.1) with
respect to $t$, we obtain

$v_{t}({\bf x},t)=\Delta v({\bf x},t)+c({\bf x})v({\bf x},t)\quad\mbox{for all
}({\bf x},t)\in\Omega\times(0,T).$ (2.2)

Since $u({\bf x},0)=p({\bf x})$, we have

$v({\bf x},0)=u_{t}({\bf x},0)=\Delta u({\bf x},0)+c({\bf x})u({\bf
x},0)=\Delta p({\bf x})+c({\bf x})p({\bf x})\quad\mbox{for all }{\bf
x}\in\Omega.$ (2.3)

Recall the assumption that $|p({\bf x})|>0$ for ${\bf x}\in\overline{\Omega}$.
Due to (2.3)

$c({\bf x})=\frac{v({\bf x},0)-\Delta p({\bf x})}{p({\bf x})}\quad\mbox{for
all }{\bf x}\in\Omega.$ (2.4)

Plugging $c({\bf x})$, computed in (2.4), into (2.2) gives

$v_{t}({\bf x},t)=\Delta v({\bf x},t)+\frac{v({\bf x},0)-\Delta p({\bf
x})}{p({\bf x})}v({\bf x},t)\quad\mbox{for all }({\bf
x},t)\in\Omega\times(0,T).$ (2.5)

Equation (2.5) is nonlinear and nonlocal. A theory to solve it is not yet
available. We propose the following dimensional reduction approach to solve
it.

######  Remark 2.1.

The change of variable in (2.1) and the elimination of $c$ to derive equation
(2.5) were adopted in [39, 44].

For each $n\in\mathbb{N}$, define $\phi_{n}(t)=t^{n-1}e^{t}$ for all
$t\in(0,T)$ and $\Phi_{n}(x)=x^{n-1}e^{x}$ for all $x\in(-R,R)$. The sets
$\\{\phi_{n}\\}_{n\geq 1}$ and $\\{\Phi_{n}\\}_{n\geq 1}$ are complete in
$L^{2}(0,T)$ and $L^{2}(-R,R)$ respectively. Applying the Gram-Schmidt
orthonormalization process on these two sets, we obtain orthonormal bases
$\\{\psi_{n}\\}$ and $\\{\Psi_{n}\\}$ of $L^{2}(0,T)$ and $L^{2}(-R,R)$
respectively. For each multi-index ${\bf
n}=(n_{1},\dots,n_{d-1},n_{t})\in\mathbb{N}^{d}$, define the $d$-dimensional
tensor-valued function ${\bf P}_{{\bf n}}$ as

${\bf P}_{\bf n}({\bf
x}^{\prime},t)=\Psi_{n_{1}}(x_{1})\dots\Psi_{n_{d-1}}(x_{d-1})\psi_{n_{t}}(t)$

for all ${\bf x}^{\prime}=(x_{1},\dots,x_{d-1})\in(-R,R)^{d-1},t\in(0,T)$. It
is obvious that the set $\\{{\bf P}_{\bf n}:\bf n\in\mathbb{N}^{d}\\}$ is an
orthonormal basis of the space $L^{2}(\Omega^{\prime}\times(0,T))$. We name
this basis the polynomial-exponential basis. The 1D version of the polynomial-
exponential basis was introduced in [23], and the higher dimension version was
defined in [43].

######  Remark 2.2.

The one-dimensional version of the polynomial-exponential basis was employed
to solve a list of inverse problems, including nonlinear inverse problems with
simulated and experimental data [17, 18, 19, 31, 33, 39, 44] and linear
inverse problems [34, 41, 44, 47].

From now on, for all ${\bf x}\in\Omega$ we write ${\bf x}=({\bf
x}^{\prime},z)\in\Omega^{\prime}\times\mathbb{R}$ where ${\bf
x}^{\prime}\in\Omega^{\prime}:=(-R,R)^{d-1}$ consists of the first $d-1$
coordinates and $z\in(-R,R)$ is the last coordinate of ${\bf x}$. Then, by
expanding the function $v({\bf x}^{\prime},z,t)$ using the basis
$\big{\\{}{\bf P}_{{\bf n}}:{\bf n}\in\mathbb{N}^{d}\big{\\}}$, we can
approximate the function $v({\bf x}^{\prime},z,t)$ as follows

$\displaystyle v({\bf x}^{\prime},z,t)$ $\displaystyle=\sum_{{\bf
n}\in\mathbb{N}^{d}}v_{{\bf n}}(z){\bf P}_{{\bf n}}({\bf
x}^{\prime},t)\simeq\sum_{{\bf n}\leqq{\bf N}}v_{{\bf n}}(z){\bf P}_{{\bf
n}}({\bf x}^{\prime},t)$
$\displaystyle=\sum_{n_{1}=1}^{N_{1}}\dots\sum_{n_{d-1}}^{N_{d-1}}\sum_{n_{t}=1}^{N_{t}}v_{(n_{1},\dots,n_{d-1},n_{t})}(z)\Psi_{n_{1}}(x_{1})\dots\Psi_{n_{d-1}}(x_{d-1})\psi_{n_{t}}(t)$
(2.6)

for $({\bf x}^{\prime},z,t)\in\Omega^{\prime}\times(-R,R)\times(0,T)$, where
${\bf N}=(N_{1},\dots,N_{d-1},N_{t})$ represents a cut-off vectors and

$v_{{\bf n}}(z)=\int_{\Omega^{\prime}\times(0,T)}v({\bf x}^{\prime},z,t){\bf
P}_{{\bf n}}({\bf x}^{\prime},t)d{\bf x}^{\prime}dt.$ (2.7)

The values of the cut-off numbers $N_{1},\dots,N_{d-1},$ and $N_{t}$ will be
chosen based on the given data in (1.3). See Section 5.2 and Figure 1 for an
illustration of a suitable choice of these numbers. In (2.6) and hereafter, we
understand ${\bf n}\leqq{\bf N}$ by the statement that

${\bf n}=(n_{1},\dots,n_{d-1},n_{t})\leqq{\bf
N}=(N_{1},\dots,N_{d-1},N_{t})\mbox{ is equivalent to }\\\ 1\leq n_{1}\leq
N_{1},\dots,1\leq n_{d-1}\leq N_{d-1},\mbox{ and }1\leq n_{t}\leq N_{t}.$
(2.8)

We assume that the approximation (2.6) is valid. Plugging (2.6) into (2.5)
gives

$\sum_{{\bf n}\leqq{\bf N}}v_{{\bf n}}(z)\frac{\partial{\bf P}_{{\bf n}}({\bf
x}^{\prime},t)}{\partial t}=\sum_{{\bf n}\leqq{\bf N}}v_{{\bf
n}}^{\prime\prime}(z){\bf P}_{{\bf n}}({\bf x}^{\prime},t)+\sum_{{\bf
n}\leqq{\bf N}}v_{{\bf n}}(z)\Delta_{{\bf x}^{\prime}}{\bf P}_{{\bf n}}({\bf
x}^{\prime},t)\\\ +\frac{1}{p({\bf x})}\Big{(}\sum_{{\bf n}\leqq{\bf
N}}v_{{\bf n}}(z){\bf P}_{{\bf n}}({\bf x}^{\prime},0)-\Delta p({\bf
x})\Big{)}\sum_{{\bf n}\leqq{\bf N}}v_{{\bf n}}(z){\bf P}_{{\bf n}}({\bf
x}^{\prime},t)$ (2.9)

for all $({\bf x}^{\prime},z,t)\in\Omega^{\prime}\times(-R,R)\times(0,T).$ For
each multi-index ${\bf m}\leqq{\bf N}$, we multiply ${\bf P}_{{\bf m}}({\bf
x}^{\prime},t)$ to both sides of (2.9), and then integrate the resulting
equation over $\Omega^{\prime}\times(0,T)$ to get

$\sum_{{\bf n}\leqq{\bf N}}v_{{\bf
n}}(z)\int_{\Omega^{\prime}\times(0,T)}\frac{\partial{\bf P}_{{\bf n}}({\bf
x}^{\prime},t)}{\partial t}{\bf P}_{{\bf m}}({\bf x}^{\prime},t)d{\bf
x}^{\prime}dt\\\ =\sum_{{\bf n}\leqq{\bf N}}v_{{\bf
n}}^{\prime\prime}(z)\int_{\Omega^{\prime}\times(0,T)}{\bf P}_{{\bf n}}({\bf
x}^{\prime},t){\bf P}_{{\bf m}}({\bf x}^{\prime},t)d{\bf
x}^{\prime}dt+\sum_{{\bf n}\leqq{\bf N}}v_{{\bf
n}}(z)\int_{\Omega^{\prime}\times(0,T)}\Delta_{{\bf x}^{\prime}}{\bf P}_{{\bf
n}}({\bf x}^{\prime},t){\bf P}_{{\bf m}}({\bf x}^{\prime},t)d{\bf
x}^{\prime}dt\\\ +\frac{1}{p({\bf x})}\Big{(}\sum_{{\bf n}\leqq{\bf N}}v_{{\bf
n}}(z){\bf P}_{{\bf n}}({\bf x}^{\prime},0)-\Delta p({\bf
x})\Big{)}\int_{\Omega^{\prime}\times(0,T)}\sum_{{\bf n}\leqq{\bf N}}v_{{\bf
n}}(z){\bf P}_{{\bf n}}({\bf x}^{\prime},t)bP_{{\bf m}}({\bf
x}^{\prime},t)d{\bf x}^{\prime}dt$ (2.10)

for all $z\in(-R,R).$ Defining

$s_{{\bf m}{\bf n}}=-\int_{\Omega^{\prime}\times(0,T)}\frac{\partial{\bf
P}_{{\bf n}}({\bf x}^{\prime},t)}{\partial t}{\bf P}_{{\bf m}}({\bf
x}^{\prime},t)d{\bf
x}^{\prime}dt+\int_{\Omega^{\prime}\times(0,T)}\Delta_{{\bf x}^{\prime}}{\bf
P}_{{\bf n}}({\bf x}^{\prime},t){\bf P}_{{\bf m}}({\bf x}^{\prime},t)d{\bf
x}^{\prime}dt$

and

${\bf F_{{\bf m}}}([v_{{\bf n}}(z)]_{{\bf n}\leqq{\bf N}})=\frac{1}{p({\bf
x})}\Big{(}\sum_{{\bf n}\leqq{\bf N}}v_{{\bf n}}(z){\bf P}_{{\bf n}}({\bf
x}^{\prime},0)-\Delta p({\bf x})\Big{)}v_{{\bf m}}(z),$ (2.11)

we obtain from (2.10) that

$v_{{\bf m}}^{\prime\prime}(z)+\sum_{{\bf n}\leqq{\bf N}}s_{{\bf m}{\bf
n}}v_{{\bf n}}(z)+{\bf F_{{\bf m}}}([v_{{\bf n}}(z)]_{{\bf n}\leqq{\bf N}})=0$
(2.12)

for all $z\in(-R,R)$ and for all ${\bf m}\leqq{\bf N}$, see (2.8) for the
definition of $\leqq$. Coupling all equations (2.12) for ${\bf m}\leqq{\bf N}$
forms a system second-order ordinary equations for the $d$-dimensional valued
tensor ${\bf v}(z)=[v_{{\bf n}}(z)]_{{\bf n}\leqq{\bf N}},$ $z\in(-R,R)$. The
Cauchy boundary conditions for the tensor ${\bf v}$ can be derived from (1.3)
and (2.7), read as

${\bf v}(\pm R)=[v_{{\bf n}}(\pm R)]_{{\bf n}\leqq{\bf
N}}=\Big{[}\int_{\Omega^{\prime}\times(0,T)}f_{t}({\bf x}^{\prime},\pm
R,t){\bf P}_{{\bf n}}({\bf x}^{\prime},t)dt\Big{]}_{{\bf n}\leqq{\bf N}}$
(2.13)

and

${\bf v}^{\prime}(\pm R)=[v_{{\bf n}}^{\prime}(\pm R)]_{{\bf n}\leqq{\bf
N}}=\Big{[}\int_{\Omega^{\prime}\times(0,T)}g_{t}({\bf x}^{\prime},\pm
R,t){\bf P}_{{\bf n}}({\bf x}^{\prime},t)dt\Big{]}_{{\bf n}\leqq{\bf N}}.$
(2.14)

Combining (2.12), (2.13), and (2.14), we obtain a system of Cauchy problem for
${\bf v}=[v_{{\bf n}}(z)]_{{\bf n}\leqq{\bf N}}$

$\left\\{\begin{array}[]{ll}v_{{\bf m}}^{\prime\prime}(z)+\sum_{{\bf
n}\leqq{\bf N}}s_{{\bf m}{\bf n}}v_{{\bf n}}(z)+{\bf F_{{\bf m}}}({\bf
v})=0&z\in(-R,R),\\\ v_{{\bf m}}(z)=P_{{\bf m}}(z)\displaystyle&z=\pm R,\\\
v_{{\bf m}}^{\prime}(z)=Q_{{\bf m}}(z)&z=\pm
R,\end{array}\right.\quad\mbox{for }{\bf m}\leqq{\bf N}.$ (2.15)

where

$\displaystyle P_{{\bf m}}(z)$
$\displaystyle=\int_{\Omega^{\prime}\times(0,T)}f_{t}({\bf
x}^{\prime},z,t){\bf P}_{{\bf m}}({\bf x}^{\prime},t)dt,$ (2.16)
$\displaystyle Q_{{\bf m}}(z)$
$\displaystyle=\displaystyle\int_{\Omega^{\prime}\times(0,T)}g_{t}({\bf
x}^{\prime},z,t){\bf P}_{{\bf m}}({\bf x}^{\prime},t)dt.$ (2.17)

Introduce the “tensor multiplication” operator

$\mathcal{S}::{\bf v}=\big{[}\sum_{{\bf n}\leqq{\bf N}}s_{{\bf m}{\bf
n}}v_{{\bf n}}\big{]}_{{\bf m}\leqq{\bf N}}$

and the notations

$\displaystyle\mathcal{F}({\bf v})$ $\displaystyle=[{\bf F}_{{\bf m}}({\bf
v})]_{{\bf m}\leqq{\bf N}},$ $\displaystyle\mathcal{P}(z)$
$\displaystyle=[P_{{\bf m}}(z)]_{{\bf m}\leqq{\bf N}},\quad z=\pm R,$
$\displaystyle\mathcal{Q}(z)$ $\displaystyle=[Q_{{\bf m}}(z)]_{{\bf
m}\leqq{\bf N}},\quad z=\pm R.$

We shorten the coupling system in (2.15) as

$\left\\{\begin{array}[]{ll}{\bf v}^{\prime\prime}(z)+\mathcal{S}::{\bf
v}(z)+\mathcal{F}({\bf v}(z))=0&z\in(-R,R),\\\ {\bf v}(z)=\mathcal{P}(z)&z=\pm
R,\\\ {\bf v}^{\prime}(z)=\mathcal{Q}(z)&z=\pm R.\end{array}\right.$ (2.18)

######  Remark 2.3.

Computing the values of ${\bf v}$ and ${\bf v}^{\prime}$ at $z=\pm R$ in
(2.13), (2.14), and (2.18) requires us to differentiate the given data $f$ and
$g$ with respect to the time $t$. This task is not trivial, especially when
the data are corrupted by noise. In this paper, we employ the new
differentiating technique in [43], in which we approximate the data by
eliminating their high-frequency terms from the Fourier expansion of the given
data with respect to the polynomial-exponential basis before differentiating.
It was numerically shown in [43] that computing derivatives using this new
technique is more accurate than the conventional ones; say the finite
difference, the cubic spline, and the Tikhonov optimization methods.

######  Remark 2.4.

The first key point of our dimension reduction approach lies in the derivation
of the approximation model (2.18), a system of first-order ODEs along the
$z-$axis. The approximation model (2.18) involves

$|{\bf N}|=|(N_{1},\dots,N_{d-1},N_{t})|=N_{t}\prod_{i=1}^{d-1}N_{i}$ (2.19)

equations versus the same numbers of unknown entries of ${\bf v}=[v_{{\bf
m}}]_{{\bf m}\leqq{\bf N}}$. This allows for the computation of the tensor-
valued function ${\bf v}(z)$ for $z\in(-R,R)$, and subsequently the function
$v({\bf x}^{\prime},z,t)$ for all $({\bf
x}^{\prime},z,t)\in\Omega^{\prime}\times(-R,R)\times(0,T)$. The solution
$c({\bf x}),$ ${\bf x}\in\Omega$, to Problem 1.1 can be computed via the
knowledge of $v$ and the reconstruction formula (2.4). However, this
convenience comes with a trade-off. The truncation in (2.5) makes system
(2.18) not exact. It should be considered as an approximation context for
Problem 1.1. Studying the behavior of (2.18) when all cut-off numbers
$N_{2},\dots,N_{d},N_{t}$ tend towards $\infty$ presents a significant
challenge. This paper does not cover this complex topic, which prioritizes
computational aspects. In exchange, we will show that our dimension reduction
method is acceptable in numerics. It can quickly deliver reliable solutions
since we have transferred a high dimensional problem into a problem along the
$z-$axis, which is a 1D problem.

As noted in Remark 2.4, once the system of ODEs in (2.18) with Cauchy boundary
data is solved, the computed solution to Problem 1.1 follows. However, this
task is challenging since (2.18) is nonlinear. There are several methods to
solve nonlinear systems of ODEs. The conventional approach is based on
optimization. For example, one can solve (2.18) by minimizing the least
squares cost functional

$J_{\rm lsq}({\bf v})=\int_{-R}^{R}\big{|}{\bf
v}^{\prime\prime}(z)+\mathcal{S}::{\bf v}+\mathcal{F}({\bf
v})\big{|}^{2}dz+\mbox{a regularization term}$ (2.20)

subject to the endpoint condition in (2.18) and then accepting the minimizer
as the computed solution. This method is effective when a good initial guess
of (2.18) is given because $J_{\rm lsq}$ might have multiple local minima. The
challenge is that such an initial guess is not always available in practical
applications. Consequently, the optimization approach is not deemed suitable
for solving (2.18). There are three approaches to solve (2.18) without
requesting a good initial guess, all based on Carleman convexification.

  1. 1.

The Carleman convexification method. The key of the Carleman convexification
method is to include a Carleman weight function; e.g.,
$W_{\lambda}(z)=e^{-\lambda z}$, $\lambda>1$, to the least squares cost
functional in (2.20). That means one can minimize the Carleman weighted
functional

$J_{\rm conv}({\bf v})=\int_{-R}^{R}W_{\lambda}(z)\big{|}{\bf
v}^{\prime\prime}(z)+\mathcal{S}::{\bf v}+\mathcal{F}({\bf
v})\big{|}^{2}dz+\mbox{a regularization term},$

subject to the boundary conditions in (2.18) where ${\bf v}=[v_{{\bf
m}}]_{{\bf m}\leqq{\bf N}}$. One can prove that $J_{\rm conv}$ is uniformly
convex in any bounded subset of the functional space containing the desired
solution provided that $\lambda$ is sufficiently large. Also, the unique
minimizer is close to the true solution to (2.18). The original
convexification method was first introduced in [24], with subsequent results
found in [2, 19, 25, 33]. Despite its efficacy in producing reliable numerical
solutions, the convexification method has a high computational cost.

  2. 2.

The Carleman contraction method. The contraction method for solving (2.18)
primarily starts with an initial function ${\bf v}^{(0)}=[v_{{\bf
m}}^{(0)}]_{{\bf m}\leqq{\bf N}}$. Note that ${\bf v}^{(0)}$ might be far away
from the true solution to (2.18). From this point, given that ${\bf
v}^{(k)}=[v_{{\bf m}}^{(k)}]_{{\bf m}\leqq{\bf N}}$, $k\geq 0$, is known, we
compute ${\bf v}^{(k+1)}=[v_{{\bf m}}^{(k+1)}]_{{\bf m}\leqq{\bf N}}$ as the
“Carleman-regularized” solution to

$\left\\{\begin{array}[]{ll}{{\bf
v}^{(k+1)}}^{\prime\prime}(z)+\mathcal{S}::{\bf v}^{(k+1)}(z)+\mathcal{F}({\bf
v}^{(k)}(z))=0&z\in(-R,R),\\\ {\bf v}(z)=\mathcal{P}(z)&z=\pm R,\\\ {\bf
v}^{\prime}(z)=\mathcal{Q}(z)&z=\pm R.\end{array}\right.$ (2.21)

By Carleman-regularized solution, we mean ${\bf v}^{(k+1)}$ is the minimizer
of

$J^{(k)}({\bf v})=\int_{-R}^{R}W_{\lambda}(z)\big{|}{\bf
v}^{\prime\prime}(z)+\mathcal{S}::{\bf v}+\mathcal{F}({\bf
v}^{(k)})])\big{|}^{2}dz+\mbox{a regularization term}$

subject to the boundary conditions in (2.18). The choice of $\lambda$,
$W_{\lambda}(z)$, and the regularization term will be specified later. The
procedure to compute ${\bf v}^{(k+1)}$ above involves the combination of the
quasi-reversibility method [28] and an appropriate Carleman estimated, as in
[30, 32, 40]. Thanks to the presence of the Carleman weight function
$W_{\lambda}(z)$, one can follow the arguments in [30, 32, 40] to prove the
convergence of the constructed sequence $\big{\\{}{\bf
v}^{(k)}\big{\\}}_{k\geq 0}$ to the true solution to (2.18).

  3. 3.

The Carleman-Newton method.  The Carleman-Newton method is similar to the
Carleman contraction method. Given an initial solution ${\bf v}^{(0)}$ that
can be chosen arbitrary, we find ${\bf v}^{(1)}$ as the Carleman regularized
solution to the linearization of (2.18) about ${\bf v}^{(0)}$. We refer the
reader to [1, 35] for details and the rigorous proof of the convergence due to
the Carleman-Newton method.

Among the three methods mentioned above, we will choose the second approach;
i.e., we will establish a 1D analog of the Carleman contraction method to
solve (2.18). This choice is appropriate due to the global convergence, the
rapid rate of convergence, and the simplicity of the computational
implementation. In the previous two sentences, we mentioned “analog” because
$\mathcal{F}$ does not satisfy the Lipschitz condition in [40], which requires
some modification in analysis.

In the next section, we establish a Carleman estimate, which plays an
important role in proving the convergence of the Carleman contraction method.

##  3 A 1D-Carleman estimate

Let $z_{0}<-R$ be a fixed number. We have the lemma.

######  Lemma 3.1.

There is a number $\lambda_{0}>1$ and a constant $C>0$ depending only on $R$
and $z_{0}$ such that for all function $w\in C^{2}([-R,R])$, we have

$\int_{-R}^{R}e^{2\lambda(z-z_{0})^{-2}}|w^{\prime\prime}(z)|^{2}dz\geq-
Ce^{2\lambda(R-z_{0})^{-2}}(\lambda^{3}|w(R)|^{2}+\lambda|w^{\prime}(R)|^{2})\\\
-Ce^{2\lambda(-R-z_{0})^{-2}}(\lambda^{3}|w(-R)|^{2}+\lambda|w^{\prime}(-R)|^{2})+C\lambda^{3}\int_{-R}^{R}e^{2\lambda(z-z_{0})^{-2}}|w(z)|^{2}dz\\\
+C\lambda\int_{-R}^{R}e^{2\lambda(z-z_{0})^{-2}}|w^{\prime}(z)|^{2}dz.$ (3.1)

###### Proof.

Step 1. Define

$y(z)=e^{\lambda(z-z_{0})^{-2}}w(z)\quad\mbox{ or equivalently }\quad
w(z)=e^{-\lambda(z-z_{0})^{-2}}y(z)$ (3.2)

for all $z\in[-R,R].$ We have

$w^{\prime}(z)=e^{-\lambda(z-z_{0})^{-2}}\Big{[}2\lambda(z-z_{0})^{-3}y(z)+y^{\prime}(z)\Big{]}$
(3.3)

and

$w^{\prime\prime}(z)=e^{-\lambda(z-z_{0})^{-2}}\Big{[}2\lambda(z-z_{0})^{-6}[3(z-z_{0})^{2}-2\lambda]y(z)+4\lambda(z-z_{0})^{-3}y^{\prime}(z)+y^{\prime\prime}(z)\Big{]}$
(3.4)

for all $z\in[-R,R]$. Thus,

$\displaystyle e^{2\lambda(z-z_{0})^{-2}}|w^{\prime\prime}(z)|^{2}$
$\displaystyle=\Big{[}2\lambda(z-z_{0})^{-6}[3(z-z_{0})^{2}-2\lambda]y(z)+4\lambda(z-z_{0})^{-3}y^{\prime}(z)+y^{\prime\prime}(z)\Big{]}^{2}$
$\displaystyle\geq
16\lambda^{2}(z-z_{0})^{-9}[3(z-z_{0})^{2}-2\lambda]y(z)y^{\prime}(z)+8\lambda(z-z_{0})^{-3}y^{\prime}(z)y^{\prime\prime}(z)$
$\displaystyle=8\lambda^{2}(z-z_{0})^{-9}[3(z-z_{0})^{2}-2\lambda]\frac{d}{dz}|y(z)|^{2}+4\lambda(z-z_{0})^{-3}\frac{d}{dz}|y^{\prime}(z)|^{2}$

for all $z\in[-R,R].$ Here, we have used the inequality $(a+b+c)^{2}\geq
2ab+2bc.$ Thus,

$(z-z_{0})^{10}e^{2\lambda(z-z_{0})^{-2}}|w^{\prime\prime}(z)|^{2}\geq
8\lambda^{2}(z-z_{0})[3(z-z_{0})^{2}-2\lambda]\frac{d}{dz}|y(z)|^{2}+4\lambda(z-z_{0})^{7}\frac{d}{dz}|y^{\prime}(z)|^{2}$
(3.5)

for all $z\in[-R,R].$ By the product rule in differentiation
$ab^{\prime}=(ab)^{\prime}-a^{\prime}b$, we have

$(z-z_{0})^{10}e^{2\lambda(z-z_{0})^{-2}}|w^{\prime\prime}(z)|^{2}\geq\frac{d}{dz}\Big{[}8\lambda^{2}(z-z_{0})[3(z-z_{0})^{2}-2\lambda]|y(z)|^{2}\Big{]}\\\
-8\lambda^{2}|y(z)|^{2}\frac{d}{dz}\Big{[}(z-z_{0})[3(z-z_{0})^{2}-2\lambda]\Big{]}+\frac{d}{dz}\Big{[}4\lambda(z-z_{0})^{7}|y^{\prime}(z)|^{2}\Big{]}\\\
-|y^{\prime}(z)|^{2}\frac{d}{dz}\Big{[}4\lambda(z-z_{0})^{7}\Big{]}$ (3.6)

for all $z\in[-R,R].$ Rearranging terms in (3.6) and simplifying the resulting
inequality, we get

$(z-z_{0})^{10}e^{2\lambda(z-z_{0})^{-2}}|w^{\prime\prime}(z)|^{2}\geq\frac{d}{dz}\Big{[}8\lambda^{2}(z-z_{0})[3(z-z_{0})^{2}-2\lambda]|y(z)|^{2}+\frac{d}{dz}\Big{[}4\lambda(z-z_{0})^{7}|y^{\prime}(z)|^{2}\Big{]}\Big{]}\\\
-8\lambda^{2}|y(z)|^{2}\Big{[}9(z-z_{0})^{2}-2\lambda\Big{]}-28|y^{\prime}(z)|^{2}\Big{[}\lambda(z-z_{0})^{6}\Big{]}$

for all $z\in[-R,R].$ Therefore,

$(z-z_{0})^{10}e^{2\lambda(z-z_{0})^{-2}}|w^{\prime\prime}(z)|^{2}\geq\frac{d}{dz}\Big{[}8\lambda^{2}(z-z_{0})[3(z-z_{0})^{2}-2\lambda]|y(z)|^{2}+\Big{[}4\lambda(z-z_{0})^{7}|y^{\prime}(z)|^{2}\Big{]}\Big{]}\\\
+16\lambda^{3}|y(z)|^{2}-72\lambda^{2}(z-z_{0})^{2}|y(z)|^{2}-28\lambda|y^{\prime}(z)|^{2}(z-z_{0})^{6}$
(3.7)

for all $z\in[-R,R].$ Integrating (3.7) over $[-R,R]$ and noting that
$\lambda^{3}\gg\lambda^{2}\gg\lambda$ as $\lambda$ large, we can find a number
$\lambda_{0}>1$ and a generic constant $C>0$, both of which depend only on
$z_{0}$ and $R$, such that

$\int_{-R}^{R}e^{2\lambda(z-z_{0})^{-2}}|w^{\prime\prime}(z)|^{2}dz\geq-C(\lambda^{3}|y(R)|^{2}+\lambda|y^{\prime}(R)|^{2})-C(\lambda^{3}|y(-R)|^{2}+\lambda|y^{\prime}(-R)|^{2})\\\
+C\lambda^{3}\int_{-R}^{R}|y(z)|^{2}dz-C\lambda\int_{-R}^{R}|y^{\prime}(z)|^{2}dz.$
(3.8)

for all $\lambda\geq\lambda_{0}.$

Step 2. Recall from (3.2) that $y(z)=e^{\lambda(z-z_{0})^{-2}}w(z).$ We have

$y^{\prime}(z)=e^{\lambda(z-z_{0})^{-2}}[-2\lambda(z-z_{0})^{-3}w(z)+w^{\prime}(z)].$

Thus, by the inequality $-(a+b)^{2}\geq-2(a^{2}+b^{2})$, we have

$-|y^{\prime}(z)|^{2}\geq-2e^{2\lambda(z-z_{0})^{-2}}\Big{[}4\lambda^{2}(z-z_{0})^{-6}|w(z)|^{2}+|w^{\prime}(z)|^{2}\Big{]}$
(3.9)

for all $z\in[-R,R].$ Combining (3.8) and (3.9) and recalling that $C$ is a
generic constant depending only on $z_{0}$ and $R$, we have

$\int_{-R}^{R}e^{2\lambda(z-z_{0})^{-2}}|w^{\prime\prime}(z)|^{2}dz\geq-
Ce^{2\lambda(R-z_{0})^{-2}}(\lambda^{3}|w(R)|^{2}+\lambda|w^{\prime}(R)|^{2})\\\
-Ce^{2\lambda(-R-z_{0})^{-2}}(\lambda^{3}|w(-R)|^{2}+\lambda|w^{\prime}(-R)|^{2})\\\
+C\lambda^{3}\int_{-R}^{R}e^{2\lambda(z-z_{0})^{-2}}|w(z)|^{2}dz-C\lambda\int_{-R}^{R}e^{2\lambda(z-z_{0})^{-2}}|w^{\prime}(z)|^{2}dz.$
(3.10)

Step 3. Using the inequality $2ab\leq a^{2}+b^{2}$, we have

$\displaystyle\int_{-R}^{R}e^{2\lambda(z-z_{0})^{-2}}\frac{d}{dz}|w^{\prime}(z)|^{2}dz$
$\displaystyle=2\int_{-R}^{R}e^{2\lambda(z-z_{0})^{-2}}w^{\prime}(z)w^{\prime\prime}(z)dz$
$\displaystyle\leq\int_{-R}^{R}e^{2\lambda(z-z_{0})^{-2}}|w^{\prime}(z)|^{2}dz+\int_{-R}^{R}e^{2\lambda(z-z_{0})^{-2}}|w^{\prime\prime}(z)|^{2}dz.$

Therefore,

$\displaystyle\int_{-R}^{R}e^{2\lambda(z-z_{0})^{-2}}|w^{\prime\prime}(z)|^{2}dz$
$\displaystyle\geq\int_{-R}^{R}e^{2\lambda(z-z_{0})^{-2}}\frac{d}{dz}|w^{\prime}(z)|^{2}dz-\int_{-R}^{R}e^{2\lambda(z-z_{0})^{-2}}|w^{\prime}(z)|^{2}dz$
$\displaystyle=\int_{-R}^{R}\frac{d}{dz}\Big{[}e^{2\lambda(z-z_{0})^{-2}}|w^{\prime}(z)|^{2}\Big{]}dz-\int_{-R}^{R}|w^{\prime}(z)|^{2}\frac{d}{dz}e^{2\lambda(z-z_{0})^{-2}}dz$
$\displaystyle\hskip
184.9429pt-\int_{-R}^{R}e^{2\lambda(z-z_{0})^{-2}}|w^{\prime}(z)|^{2}dz.$

As a result,

$\int_{-R}^{R}e^{2\lambda(z-z_{0})^{-2}}|w^{\prime\prime}(z)|^{2}dz\geq-C\Big{[}e^{2\lambda(R-z_{0})^{-2}}|w^{\prime}(R)|^{2}+e^{2\lambda(-R-z_{0})^{-2}}|w^{\prime}(-R)|^{2}\Big{]}\\\
+2\lambda\int_{-R}^{R}(z-z_{0})^{-3}e^{2\lambda(z-z_{0})^{-2}}|w^{\prime}(z)|^{2}dz-\int_{-R}^{R}e^{2\lambda(z-z_{0})^{-2}}|w^{\prime}(z)|^{2}dz.$
(3.11)

Adding (3.10) and (3.11) and recalling that $\lambda\geq\lambda_{0}\gg 1$, we
obtain (3.1).

∎

##  4 A Picard-like iteration to solve (2.18)

In this section, we employ the Carleman estimate in Lemma 3.1 to construct a
sequence that converges to the solution to (2.18), provided that this true
solution exists. We consider the circumstance that the boundary data
$\mathcal{P}$ and $\mathcal{Q}$ of (2.18) contain noise. Let $\mathcal{P}^{*}$
and $\mathcal{Q}^{*}$ be the unknown exact values of the boundary data
$\mathcal{P}$ and $\mathcal{Q}$, respectively. Let ${\bf v}^{*}$ be the
solution to (2.18) with $\mathcal{P}$ and $\mathcal{Q}$ being replaced by
$\mathcal{P}^{*}$ and $\mathcal{Q}^{*}$, respectively. That means, ${\bf
v}^{*}$ solves

$\left\\{\begin{array}[]{ll}{{\bf v}^{*}}^{\prime\prime}(z)+\mathcal{S}::{\bf
v}^{*}(z)+\mathcal{F}({\bf v}^{*}(z))=0&z\in(-R,R),\\\ {\bf
v}^{*}(z)=\mathcal{P}^{*}(z)&z=\pm R,\\\ {{\bf
v}^{*}}^{\prime}(z)=\mathcal{Q}^{*}(z)&z=\pm R.\end{array}\right.$ (4.1)

In this section, we assume the existence of the solution ${\bf v}^{*}$ to
(4.1). We now consider the case when noise is introduced to the data. Let
$\delta>0$ be the noise level. That means,

$\max_{z\in\\{-R,R\\}}\big{\\{}|\mathcal{P}(z)-\mathcal{P}^{*}(z),|\mathcal{Q}(z)-\mathcal{Q}^{*}(z)|\big{\\}}<\delta.$
(4.2)

######  Remark 4.1 (Noise model).

In this section, for simplicity, we assume that noise is introduced into the
indirect data $\mathcal{P}(\pm R)$ and $\mathcal{Q}(\pm R)$ as in (4.2) rather
than to the direct data, $f$ and $g$. This assumption serves theoretical
purposes only. In our computational study, we study the more realistic case
where the direct data $f^{*}$ and $g^{*}$ are impacted by noise as in (5.4)
and (5.5). Recall that the the entries $P_{{\bf m}}(\pm R)$ and $Q_{{\bf
m}}(\pm R)$ of indirect data $\mathcal{P}(\pm R)$ and $\mathcal{Q}(\pm R)$ are
computed by the knowledge of the derivatives of $f$ and $g$ via (2.16) and
(2.17). Given that differentiating noisy data presents significant challenges
and can greatly amplify errors, even minor noise in $f$ and $g$ can lead to
substantial inaccuracies in $P_{{\bf m}}(\pm R)$ and $Q_{{\bf m}}(\pm R)$. To
address this issue, we employ a novel differentiation approach as presented in
[43]. In [43], this method has been demonstrated to have superior stability
compared to traditional methods like finite difference, cubic splines, or the
Tikhonov regularization technique.

Consider the space of admissible solutions

$H=\Big{\\{}\bm{\varphi}=[\varphi_{{\bf m}}]_{{\bf m}\leqq{\bf
N}}:\varphi_{{\bf m}}\in H^{2}(-R,R)\mbox{ for all }{\bf m}\leqq{\bf
N}\Big{\\}}=H^{2}(-R,R)^{|{\bf N}|}$

where $|{\bf N}|$ is as in (2.19). Fix an arbitrary large number $M>0$, define
the close ball in $H$

$B_{M}=\\{\bm{\varphi}:\|\bm{\varphi}\|_{H^{2}(-R,R)^{|{\bf N}|}}\leq M\\}.$
(4.3)

For each ${\bf v}\in B_{M}$, define the functional

$J_{\bf
v}^{(\lambda,\epsilon)}(\bm{\varphi})=\int_{-R}^{R}e^{2\lambda(z-z_{0})^{-2}}\Big{|}\bm{\varphi}^{\prime\prime}+\mathcal{S}::\bm{\varphi}+\mathcal{F}({\bf
v})\Big{|}^{2}dz+\lambda^{4}e^{2\lambda(R-z_{0})^{-2}}\left|\bm{\varphi}(R)-\mathcal{P}(R)\right|^{2}\\\
+\lambda^{4}e^{2\lambda(-R-z_{0})^{-2}}\left|\bm{\varphi}(-R)-\mathcal{P}(-R)\right|^{2}+\lambda^{4}e^{2\lambda(R-z_{0})^{-2}}\left|\bm{\varphi}^{\prime}(R)-\mathcal{Q}(R)\right|^{2}\\\
+\lambda^{4}e^{2\lambda(-R-z_{0})^{-2}}\left|\bm{\varphi}^{\prime}(-R)-\mathcal{Q}(-R)\right|^{2}+\epsilon\|\bm{\varphi}\|_{H^{2}(-R,R)^{|{\bf
N}|}}^{2}$ (4.4)

for all $\bm{\varphi}$ where $\lambda\geq\lambda_{0}$ and $\epsilon>0$ will be
chosen later. For all ${\bf v}\in B_{M}$, the functional $J_{\bf
v}^{(\lambda,\epsilon)}$ is uniformly convex in the close and convex set
$B_{M}$ of $H$. It has a unique minimizer. We define the map
$\Phi^{(\lambda,\epsilon)}:B_{M}\to B_{M}$ that sends ${\bf v}$ to such a
minimizer. More precisely,

$\Phi^{(\lambda,\epsilon)}({\bf v})=\underset{\bm{\varphi}\in B_{M}}{{\rm
argmin}}J_{\bf v}^{(\lambda,\epsilon)}(\bm{\varphi})\quad\mbox{for all }{\bf
v}\in B_{M}.$ (4.5)

We define the sequence $\\{{\bf v}^{(k)}\\}_{k\geq 0}$ as follows:

$\left\\{\begin{array}[]{ll}{\bf v}^{(0)}\mbox{ chosen arbitrarily in
}B_{M}\\\ {\bf v}^{(k+1)}=\Phi^{(\lambda,\epsilon)}({\bf v}^{(k)})\quad k\geq
0.\end{array}\right.$ (4.6)

The following theorem guarantees the convergence of the sequence $\\{{\bf
v}^{(k)}\\}_{k\geq 0}$ to ${\bf v}^{*}$.

######  Theorem 4.1.

Let $M$ be a large number such that both ${\bf v}^{*}$ and ${\bf v}^{(0)}$ are
in $B_{M}$. Let $\lambda_{0}$ be the number in Lemma 3.1. Then, there exists
$\lambda_{1}>\lambda_{0}$ depending only on $M,$ $p$, $N,$ $R$, and $[{\bf
P}_{{\bf m}\leqq{\bf N}}]$ such that

$\int_{-R}^{R}e^{2\lambda(z-z_{0})^{-2}}|{\bf v}^{(k+1)}-{\bf
v}^{*}|^{2}dz\leq\Big{(}\frac{C}{\lambda^{3}}\Big{)}^{k+1}\int_{-R}^{R}e^{2\lambda(z-z_{0})^{-2}}\big{|}{\bf
v}^{(0)}-{\bf v}^{*}\big{|}^{2}dz\\\
+\frac{C/\lambda^{3}}{1-C/\lambda^{3}}\Big{[}\lambda^{4}e^{2\lambda(R-z_{0})^{-2}}|\mathcal{P}(R)-\mathcal{P}^{*}(R))|^{2}+\lambda^{4}e^{2\lambda(-R-z_{0})^{-2}}|\mathcal{P}(-R)-\mathcal{P}^{*}(-R)|^{2}\\\
+\lambda^{4}e^{2\lambda(R-z_{0})^{-2}}|\mathcal{Q}(R)-\mathcal{Q}^{*}(R)|^{2}+\lambda^{4}e^{2\lambda(-R-z_{0})^{-2}}|\mathcal{Q}(-R)-\mathcal{Q}^{*}(-R)|^{2}+\epsilon\|{\bf
v}^{*}\|_{H^{2}(-R,R)^{|{\bf N}|}}^{2}\Big{]}$ (4.7)

for all $k\geq 0.$

######  Remark 4.2.

Theorem 4.1 and its proof are stated and proved using similar arguments in
[30, 32, 40]. However, we still need some important modifications:

  1. 1.

The nonlinearity $\mathcal{F}$ in [30, 32, 40] needs to satisfy the Lipschitz
condition. However, the function $\mathcal{F}$ in the current work does not
meet this requirement. To address this issue, it is necessary to confine the
computational domain to a bounded set $B_{M}$ for an arbitrarily large number
$M$. Within this bounded domain, the Lipschitz condition is automatically
satisfied.

  2. 2.

In [30, 32], the analysis of noise was not explored, whereas it was somewhat
examined in [40]. By “somewhat,” it means that in [40], a technical condition
had to be imposed. The noise in the Dirichlet observations and the noise in
the Neumann measurements are not independent. Specifically, it was assumed
that the noise in the Dirichlet observation is the trace of a function, and
the noise in the Neumann measurement needs to be the trace of that function’s
normal derivative. Given that this circumstance is somewhat impractical, we
opt to relax it in the present paper.

###### Proof of Theorem 4.1.

In the proof, we will employ the dot product

$\bm{\varphi}\cdot{\bf h}=\sum_{{\bf m}\leqq{\bf N}}\varphi_{{\bf m}}h_{{\bf
m}}$

for all $\bm{\varphi}=[\varphi_{{\bf m}}]_{{\bf m}\leq{\bf N}}$ and ${\bf
h}=[h_{{\bf m}}]_{{\bf m}\leq{\bf N}}$ in $H$. Fix $k\geq 0$. Set

${\bf h}={\bf v}^{(k+1)}-{\bf v}^{*}.$

Since ${\bf v}^{(k+1)}$ is the minimizer of $J^{\lambda,\epsilon}_{{\bf
v}^{(k)}}$ in $B_{M}$ and ${\bf v}^{*}$ is in the interior of $B_{M}$, we have

$\int_{-R}^{R}e^{2\lambda(z-z_{0})^{-2}}\big{(}{{\bf
v}^{(k+1)}}^{\prime\prime}(z)+\mathcal{S}::{\bf v}^{(k+1)}(z)+\mathcal{F}({\bf
v}^{(k)})\big{)}\cdot\big{(}{\bf h}^{\prime\prime}(z)+\mathcal{S}::{\bf
h}(z)\big{)}dz\\\ +\lambda^{4}e^{2\lambda(R-z_{0})^{-2}}({\bf
v}^{(k+1)}(R)-\mathcal{P}(R))\cdot{\bf
h}(R)+\lambda^{4}e^{2\lambda(-R-z_{0})^{-2}}({\bf
v}^{(k+1)}(-R)-\mathcal{P}(-R))\cdot{\bf h}(-R)\\\
\lambda^{4}e^{2\lambda(R-z_{0})^{-2}}({{\bf
v}^{(k+1)}}^{\prime}(R)-\mathcal{Q}(R))\cdot{\bf
h}^{\prime}(R)+\lambda^{4}e^{2\lambda(-R-z_{0})^{-2}}({{\bf
v}^{(k+1)}}^{\prime}(-R)-\mathcal{Q}(-R))\cdot{\bf h}^{\prime}(-R)\\\
+\epsilon\langle{\bf v}^{(k+1)},{\bf h}\rangle_{H^{2}(-R,R)^{|{\bf N}|}}\leq
0.$ (4.8)

On the other hand, since ${\bf v}^{*}$ is the true solution to (4.1), we have

$\int_{-R}^{R}e^{2\lambda(z-z_{0})^{-2}}\big{(}{{\bf
v}^{*}}^{\prime\prime}(z)+\mathcal{S}::{\bf v}^{*}(z)+\mathcal{F}({\bf
v}^{*})\big{)}\cdot\big{(}{\bf h}^{\prime\prime}(z)+\mathcal{S}::{\bf
h}(z)\big{)}dz\\\ +\lambda^{4}e^{2\lambda(R-z_{0})^{-2}}({\bf
v}^{*}(R)-\mathcal{P}^{*}(R))\cdot{\bf
h}(R)+\lambda^{4}e^{2\lambda(-R-z_{0})^{-2}}({\bf
v}^{*}(-R)-\mathcal{P}^{*}(-R))\cdot{\bf h}(-R)\\\
\lambda^{4}e^{2\lambda(R-z_{0})^{-2}}({{\bf
v}^{*}}^{\prime}(R)-\mathcal{Q}^{*}(R))\cdot{\bf
h}^{\prime}(R)+\lambda^{4}e^{2\lambda(-R-z_{0})^{-2}}({{\bf
v}^{*}}^{\prime}(-R)-\mathcal{Q}^{*}(-R))\cdot{\bf h}^{\prime}(-R)\\\
+\epsilon\langle{\bf v}^{*},{\bf h}\rangle_{H^{2}(-R,R)^{|{\bf
N}|}}=\epsilon\langle{\bf v}^{*},{\bf h}\rangle_{H^{2}(-R,R)^{|{\bf N}|}}.$
(4.9)

Subtracting (4.8) from (4.9) gives

$\int_{-R}^{R}e^{2\lambda(z-z_{0})^{-2}}\big{(}{{\bf
v}^{(k+1)}}^{\prime\prime}(z)-{{\bf
v}^{*}}^{\prime\prime}(z)+\mathcal{S}::({\bf v}^{(k+1)}(z)-{{\bf
v}^{*}}^{\prime\prime}(z))+\mathcal{F}({\bf v}^{(k)})-\mathcal{F}({\bf
v}^{*})\big{)}\cdot\big{(}{\bf h}^{\prime\prime}(z)+\mathcal{S}::{\bf
h}(z)\big{)}dz\\\ +\lambda^{4}e^{2\lambda(R-z_{0})^{-2}}({\bf
v}^{(k+1)}(R)-{\bf v}^{*}(R))\cdot{\bf
h}(R)+\lambda^{4}e^{2\lambda(-R-z_{0})^{-2}}({\bf v}^{(k+1)}(-R)-{\bf
v}^{*}(-R))\cdot{\bf h}(-R)\\\ +\lambda^{4}e^{2\lambda(R-z_{0})^{-2}}({{\bf
v}^{(k+1)}}^{\prime}(R)-{{\bf v}^{*}}^{\prime}(R))\cdot{\bf
h}^{\prime}(R)+\lambda^{4}e^{2\lambda(-R-z_{0})^{-2}}({{\bf
v}^{(k+1)}}^{\prime}(-R)-{{\bf v}^{*}}^{\prime}(-R))\cdot{\bf
h}^{\prime}(-R)\\\ +\epsilon\langle{\bf v}^{(k+1)}-{\bf v}^{*},{\bf
h}\rangle_{H^{2}(-R,R)^{|{\bf N}|}}\\\
\leq\lambda^{4}e^{2\lambda(R-z_{0})^{-2}}(\mathcal{P}(R)-\mathcal{P}^{*}(R))\cdot{\bf
h}(R)+\lambda^{4}e^{2\lambda(-R-z_{0})^{-2}}(\mathcal{P}(-R)-\mathcal{P}^{*}(-R))\cdot{\bf
h}(-R)\\\
+\lambda^{4}e^{2\lambda(R-z_{0})^{-2}}(\mathcal{Q}(R)-\mathcal{Q}^{*}(R))\cdot{\bf
h}^{\prime}(R)+\lambda^{4}e^{2\lambda(-R-z_{0})^{-2}}(\mathcal{Q}(-R)-\mathcal{Q}^{*}(-R))\cdot{\bf
h}^{\prime}(-R)\\\ -\epsilon\langle{\bf v}^{*},{\bf
h}\rangle_{H^{2}(-R,R)^{|{\bf N}|}}.$ (4.10)

Recalling ${\bf h}={\bf v}^{(k+1)}-{\bf v}^{*}$, we have

$\int_{-R}^{R}e^{2\lambda(z-z_{0})^{-2}}\big{|}{\bf
h}^{\prime\prime}(z)+\mathcal{S}::{\bf
h}(z)\big{|}^{2}dz+\lambda^{4}e^{2\lambda(R-z_{0})^{-2}}\ |{\bf
h}(R)|^{2}+\lambda^{4}e^{2\lambda(-R-z_{0})^{-2}}|{\bf h}(-R)|^{2}\\\
+\lambda^{4}e^{2\lambda(R-z_{0})^{-2}}|{\bf
h}^{\prime}(R)|^{2}+\lambda^{4}e^{2\lambda(-R-z_{0})^{-2}}|{\bf
h}^{\prime}(-R)|^{2}+\epsilon\|{\bf h}\|_{H^{2}(-R,R)^{|{\bf N}|}}^{2}\\\
\leq-\int_{-R}^{R}e^{2\lambda(z-z_{0})^{-2}}\big{(}\mathcal{F}({\bf
v}^{(k)})-\mathcal{F}({\bf v}^{*})\big{)}\cdot\big{(}{\bf
h}^{\prime\prime}(z)+\mathcal{S}::{\bf
h}(z)\big{)}dz+\lambda^{4}e^{2\lambda(R-z_{0})^{-2}}(\mathcal{P}(R)-\mathcal{P}^{*}(R))\cdot{\bf
h}(R)\\\
+\lambda^{4}e^{2\lambda(-R-z_{0})^{-2}}(\mathcal{P}(-R)-\mathcal{P}^{*}(-R))\cdot{\bf
h}(-R)+\lambda^{4}e^{2\lambda(R-z_{0})^{-2}}(\mathcal{Q}(R)-\mathcal{Q}^{*}(R))\cdot{\bf
h}^{\prime}(R)\\\
+\lambda^{4}e^{2\lambda(-R-z_{0})^{-2}}(\mathcal{Q}(-R)-\mathcal{Q}^{*}(-R))\cdot{\bf
h}^{\prime}(-R)-\epsilon\langle{\bf v}^{*},{\bf h}\rangle_{H^{2}(-R,R)^{|{\bf
N}|}}.$ (4.11)

Rearranging terms (4.11) and using the inequality
$|ab|\leq\frac{1}{2}(a^{2}+b^{2})$, we have

$\int_{-R}^{R}e^{2\lambda(z-z_{0})^{-2}}\big{|}{\bf
h}^{\prime\prime}(z)+\mathcal{S}::{\bf
h}(z)\big{|}^{2}dz+\lambda^{4}e^{2\lambda(R-z_{0})^{-2}}\ |{\bf
h}(R)|^{2}+\lambda^{4}e^{2\lambda(-R-z_{0})^{-2}}|{\bf h}(-R)|^{2}\\\
+\lambda^{4}e^{2\lambda(R-z_{0})^{-2}}|{\bf
h}^{\prime}(R)|^{2}+\lambda^{4}e^{2\lambda(-R-z_{0})^{-2}}|{\bf
h}^{\prime}(-R)|^{2}+\epsilon\|{\bf h}\|_{H^{2}(-R,R)^{|{\bf N}|}}^{2}\\\
\leq-\int_{-R}^{R}e^{2\lambda(z-z_{0})^{-2}}\big{|}\mathcal{F}({\bf
v}^{(k)})-\mathcal{F}({\bf
v}^{*})\big{|}^{2}dz+\lambda^{4}e^{2\lambda(R-z_{0})^{-2}}|\mathcal{P}(R)-\mathcal{P}^{*}(R))|^{2}\\\
+\lambda^{4}e^{2\lambda(-R-z_{0})^{-2}}|\mathcal{P}(-R)-\mathcal{P}^{*}(-R)|^{2}+\lambda^{4}e^{2\lambda(R-z_{0})^{-2}}|\mathcal{Q}(R)-\mathcal{Q}^{*}(R)|^{2}\\\
+\lambda^{4}e^{2\lambda(-R-z_{0})^{-2}}|\mathcal{Q}(-R)-\mathcal{Q}^{*}(-R)|^{2}+\epsilon\|{\bf
v}^{*}\|_{H^{2}(-R,R)^{|{\bf N}|}}^{2}.$ (4.12)

Apply the inequality $(a+b)^{2}\geq\frac{1}{2}a^{2}-b^{2}$ for the first term
of (4.12). We have

$\frac{1}{2}\int_{-R}^{R}e^{2\lambda(z-z_{0})^{-2}}\big{|}{\bf
h}^{\prime\prime}(z)\big{|}^{2}dz+\lambda^{4}e^{2\lambda(R-z_{0})^{-2}}\ |{\bf
h}(R)|^{2}+\lambda^{4}e^{2\lambda(-R-z_{0})^{-2}}|{\bf h}(-R)|^{2}\\\
+\lambda^{4}e^{2\lambda(R-z_{0})^{-2}}|{\bf
h}^{\prime}(R)|^{2}+\lambda^{4}e^{2\lambda(-R-z_{0})^{-2}}|{\bf
h}^{\prime}(-R)|^{2}+\epsilon\|{\bf h}\|_{H}^{2}\leq
2\int_{-R}^{R}e^{2\lambda(z-z_{0})^{-2}}|\mathcal{S}::{\bf h}|^{2}dz\\\
-\int_{-R}^{R}e^{2\lambda(z-z_{0})^{-2}}\big{|}\mathcal{F}({\bf
v}^{(k)})-\mathcal{F}({\bf
v}^{*})\big{|}^{2}dz+\lambda^{4}e^{2\lambda(R-z_{0})^{-2}}|\mathcal{P}(R)-\mathcal{P}^{*}(R))|^{2}\\\
+\lambda^{4}e^{2\lambda(-R-z_{0})^{-2}}|\mathcal{P}(-R)-\mathcal{P}^{*}(-R)|^{2}+\lambda^{4}e^{2\lambda(R-z_{0})^{-2}}|\mathcal{Q}(R)-\mathcal{Q}^{*}(R)|^{2}\\\
+\lambda^{4}e^{2\lambda(-R-z_{0})^{-2}}|\mathcal{Q}(-R)-\mathcal{Q}^{*}(-R)|^{2}+\epsilon\|{\bf
v}^{*}\|_{H^{2}(-R,R)^{|{\bf N}|}}^{2}.$ (4.13)

Since $B_{M}$ is bounded in $H$, by the Sobolev embedding theorem in 1D,
$B_{M}$ is bounded in $C([-R,R])^{|{\bf N}|}$. We can find a number $C$
depending only on $M$, $\mathcal{F}$ (and hence $p,$ $N$, $R$, $[{\bf P}_{{\bf
m}}]_{{\bf m}\leqq{\bf N}}$) such that

$|\mathcal{F}({\bf v}^{(k)}(z))-\mathcal{F}(\mathcal{{\bf v}^{*}})(z)|\leq
C|{\bf v}^{(k)}(z)-\mathcal{{\bf v}^{*}}(z)|\quad\mbox{for all }z\in(-R,R).$
(4.14)

Combining (4.13) and (4.14), we can find a constant $C$ depending only $M$,
$p$, $N$, $R$, $[{\bf P}_{{\bf m}}]_{{\bf m}\leqq{\bf N}}$ such that

$\int_{-R}^{R}e^{2\lambda(z-z_{0})^{-2}}\big{|}{\bf
h}^{\prime\prime}(z)\big{|}^{2}dz+\lambda^{4}e^{2\lambda(R-z_{0})^{-2}}\ |{\bf
h}(R)|^{2}+\lambda^{4}e^{2\lambda(-R-z_{0})^{-2}}|{\bf h}(-R)|^{2}\\\
+\lambda^{4}e^{2\lambda(R-z_{0})^{-2}}|{\bf
h}^{\prime}(R)|^{2}+\lambda^{4}e^{2\lambda(-R-z_{0})^{-2}}|{\bf
h}^{\prime}(-R)|^{2}+\epsilon\|{\bf h}\|_{H^{2}(-R,R)^{|{\bf N}|}}^{2}\\\ \leq
C\Big{[}\int_{-R}^{R}e^{2\lambda(z-z_{0})^{-2}}\big{|}{\bf
h}\big{|}^{2}dz+\int_{-R}^{R}e^{2\lambda(z-z_{0})^{-2}}\big{|}{\bf
v}^{(k)}-{\bf
v}^{*}\big{|}^{2}dz+\lambda^{4}e^{2\lambda(R-z_{0})^{-2}}|\mathcal{P}(R)-\mathcal{P}^{*}(R))|^{2}\\\
+\lambda^{4}e^{2\lambda(-R-z_{0})^{-2}}|\mathcal{P}(-R)-\mathcal{P}^{*}(-R)|^{2}+\lambda^{4}e^{2\lambda(R-z_{0})^{-2}}|\mathcal{Q}(R)-\mathcal{Q}^{*}(R)|^{2}+\lambda^{4}e^{2\lambda(-R-z_{0})^{-2}}|\mathcal{Q}(-R)-\mathcal{Q}^{*}(-R)|^{2}\\\
+\epsilon\|{\bf v}^{*}\|_{H^{2}(-R,R)^{|{\bf N}|}}^{2}\Big{]}.$ (4.15)

Applying the Carleman estimate in Lemma 3.1 for each entry of ${\bf h}$, we
have

$\int_{-R}^{R}e^{2\lambda(z-z_{0})^{-2}}|{\bf h}^{\prime\prime}(z)|^{2}dz\geq-
Ce^{2\lambda(R-z_{0})^{-2}}(\lambda^{3}|{\bf h}(R)|^{2}+\lambda|{\bf
h}^{\prime}(R)|^{2})\\\ -Ce^{2\lambda(-R-z_{0})^{-2}}(\lambda^{3}|{\bf
h}(-R)|^{2}+\lambda|{\bf
h}^{\prime}(-R)|^{2})+C\lambda^{3}\int_{-R}^{R}e^{2\lambda(z-z_{0})^{-2}}|{\bf
h}(z)|^{2}dz\\\ +C\lambda\int_{-R}^{R}e^{2\lambda(z-z_{0})^{-2}}|{\bf
h}^{\prime}(z)|^{2}dz.$ (4.16)

Combining (4.15) and (4.16) and noting that $\lambda^{4}\gg\lambda^{3}$ , we
have

$\lambda^{3}\int_{-R}^{R}e^{2\lambda(z-z_{0})^{-2}}|{\bf
h}(z)|^{2}dz+\lambda\int_{-R}^{R}e^{2\lambda(z-z_{0})^{-2}}|{\bf
h}^{\prime}(z)|^{2}dz+\lambda^{4}e^{2\lambda(R-z_{0})^{-2}}\ |{\bf
h}(R)|^{2}\\\ +\lambda^{4}e^{2\lambda(-R-z_{0})^{-2}}|{\bf
h}(-R)|^{2}+\lambda^{4}e^{2\lambda(R-z_{0})^{-2}}|{\bf
h}^{\prime}(R)|^{2}+\lambda^{4}e^{2\lambda(-R-z_{0})^{-2}}|{\bf
h}^{\prime}(-R)|^{2}\\\ \leq
C\Big{[}\int_{-R}^{R}e^{2\lambda(z-z_{0})^{-2}}\big{|}{\bf
h}\big{|}^{2}dz+\int_{-R}^{R}e^{2\lambda(z-z_{0})^{-2}}\big{|}{\bf
v}^{(k)}-{\bf
v}^{*}\big{|}^{2}dz+\lambda^{4}e^{2\lambda(R-z_{0})^{-2}}|\mathcal{P}(R)-\mathcal{P}^{*}(R))|^{2}\\\
+\lambda^{4}e^{2\lambda(-R-z_{0})^{-2}}|\mathcal{P}(-R)-\mathcal{P}^{*}(-R)|^{2}+\lambda^{4}e^{2\lambda(R-z_{0})^{-2}}|\mathcal{Q}(R)-\mathcal{Q}^{*}(R)|^{2}\\\
+\lambda^{4}e^{2\lambda(-R-z_{0})^{-2}}|\mathcal{Q}(-R)-\mathcal{Q}^{*}(-R)|^{2}+\epsilon\|{\bf
v}^{*}\|_{H^{2}(-R,R)^{|{\bf N}|}}^{2}\Big{]}.$ (4.17)

It follows from (4.17) and the fact ${\bf h}={\bf v}^{(k+1)}-{\bf v}^{*}$ that

$\int_{-R}^{R}e^{2\lambda(z-z_{0})^{-2}}|{\bf v}^{(k+1)}-{\bf
v}^{*}|^{2}dz\leq\frac{C}{\lambda^{3}}\Big{[}\int_{-R}^{R}e^{2\lambda(z-z_{0})^{-2}}\big{|}{\bf
v}^{(k)}-{\bf
v}^{*}\big{|}^{2}dz+\lambda^{4}e^{2\lambda(R-z_{0})^{-2}}|\mathcal{P}(R)-\mathcal{P}^{*}(R))|^{2}\\\
+\lambda^{4}e^{2\lambda(-R-z_{0})^{-2}}|\mathcal{P}(-R)-\mathcal{P}^{*}(-R)|^{2}+\lambda^{4}e^{2\lambda(R-z_{0})^{-2}}|\mathcal{Q}(R)-\mathcal{Q}^{*}(R)|^{2}\\\
+\lambda^{4}e^{2\lambda(-R-z_{0})^{-2}}|\mathcal{Q}(-R)-\mathcal{Q}^{*}(-R)|^{2}+\epsilon\|{\bf
v}^{*}\|_{H^{2}(-R,R)^{|{\bf N}|}}^{2}\Big{]}.$ (4.18)

Recall that ${\bf h}={\bf v}^{(k+1)}-{\bf v}^{*}$. Applying (4.18) when $k+1$
is replaced by $k$ and combining the resulting estimate with (4.18), we have

$\int_{-R}^{R}e^{2\lambda(z-z_{0})^{-2}}|{\bf v}^{(k+1)}-{\bf
v}^{*}|^{2}dz\leq\frac{C}{\lambda^{3}}\Big{[}\frac{C}{\lambda^{3}}\Big{[}\int_{-R}^{R}e^{2\lambda(z-z_{0})^{-2}}\big{|}{\bf
v}^{(k-1)}-{\bf v}^{*}\big{|}^{2}dz\\\
+\lambda^{4}e^{2\lambda(R-z_{0})^{-2}}|\mathcal{P}(R)-\mathcal{P}^{*}(R))|^{2}+\lambda^{4}e^{2\lambda(-R-z_{0})^{-2}}|\mathcal{P}(-R)-\mathcal{P}^{*}(-R)|^{2}\\\
+\lambda^{4}e^{2\lambda(R-z_{0})^{-2}}|\mathcal{Q}(R)-\mathcal{Q}^{*}(R)|^{2}+\lambda^{4}e^{2\lambda(-R-z_{0})^{-2}}|\mathcal{Q}(-R)-\mathcal{Q}^{*}(-R)|^{2}\\\
+\epsilon\|{\bf v}^{*}\|_{H^{2}(-R,R)^{|{\bf
N}|}}^{2}\Big{]}+\lambda^{4}e^{2\lambda(R-z_{0})^{-2}}|\mathcal{P}(R)-\mathcal{P}^{*}(R))|^{2}+\lambda^{4}e^{2\lambda(-R-z_{0})^{-2}}|\mathcal{P}(-R)-\mathcal{P}^{*}(-R)|^{2}\\\
+\lambda^{4}e^{2\lambda(R-z_{0})^{-2}}|\mathcal{Q}(R)-\mathcal{Q}^{*}(R)|^{2}+\lambda^{4}e^{2\lambda(-R-z_{0})^{-2}}|\mathcal{Q}(-R)-\mathcal{Q}^{*}(-R)|^{2}+\epsilon\|{\bf
v}^{*}\|_{H^{2}(-R,R)^{|{\bf N}|}}^{2}\Big{]}.$

Continuing the procedure, we obtain

$\int_{-R}^{R}e^{2\lambda(z-z_{0})^{-2}}|{\bf v}^{(k+1)}-{\bf
v}^{*}|^{2}dz\leq\Big{(}\frac{C}{\lambda^{3}}\Big{)}^{k+1}\int_{-R}^{R}e^{2\lambda(z-z_{0})^{-2}}\big{|}{\bf
v}^{(0)}-{\bf v}^{*}\big{|}^{2}dz\\\
+\sum_{i=1}^{k}\Big{(}\frac{C}{\lambda^{3}}\Big{)}^{k}\Big{[}\lambda^{4}e^{2\lambda(R-z_{0})^{-2}}|\mathcal{P}(R)-\mathcal{P}^{*}(R))|^{2}+\lambda^{4}e^{2\lambda(-R-z_{0})^{-2}}|\mathcal{P}(-R)-\mathcal{P}^{*}(-R)|^{2}\\\
+\lambda^{4}e^{2\lambda(R-z_{0})^{-2}}|\mathcal{Q}(R)-\mathcal{Q}^{*}(R)|^{2}+\lambda^{4}e^{2\lambda(-R-z_{0})^{-2}}|\mathcal{Q}(-R)-\mathcal{Q}^{*}(-R)|^{2}+\epsilon\|{\bf
v}^{*}\|_{H^{2}(-R,R)^{|{\bf N}|}}^{2}\Big{]}.$ (4.19)

Choose $\lambda>\lambda_{1}>\lambda_{0}$ such that $C/\lambda^{3}\in(0,1)$ for
some $\lambda_{1}$ depending only on $C$ and therefore only on $M,$ $p$, $N,$
$R$, and $[{\bf P}_{{\bf m}\leqq{\bf N}}]$. Estimate (4.7) is a direct
consequence of (4.19). ∎

######  Remark 4.3.

Estimate (4.7) is interesting in the sense that when the data has noise,
although the over-determined problem (2.18) might have no solution, we are
still able to provide a resonably accurate numerical solution. In fact, the
sequence $\\{{\bf v}^{(k)}\\}_{k\geq 0}$ is well-defined regardless whether or
not (2.18) is solvable. This sequence is defined via minimizing a strictly
convex cost functional as in (4.5) when ${\bf v}$ is replaced by ${\bf
v}^{(k)}$, $k\geq 0$. Due to (4.2) and the result of Theorem 4.1 in (4.9),
when $\lambda>\lambda_{1}$ is fixed so that $C/\lambda^{3}\in(0,1)$, we have

$\|{\bf v}^{(k+1)}-{\bf v}^{*}\|^{2}_{L^{2}(-R,R)}\leq
C_{1}\Big{[}\Big{(}\frac{C}{\lambda^{3}}\Big{)}^{k+1}\|{\bf v}^{(0)}-{\bf
v}^{*}\|^{2}_{L^{2}(-R,R)}+\frac{\delta^{2}C/\lambda^{3}}{1-C/\lambda^{3}}+\epsilon\|{\bf
v}^{*}\|_{H^{2}(-R,R)^{|{\bf N}|}}^{2}\Big{]}$ (4.20)

where $C$ is the constant in Theorem 4.1 and $C_{1}$ is a constant depending
only on $\lambda,$ $z_{0}$, and $R$. On the other hand, (4.20) guarantees the
global convergence. Although ${\bf v}^{(0)}$ is not a good initial guess of
${\bf v}^{*}$, the approximating sequence converges to a numerical solution
with the error $O(\delta+\sqrt{\epsilon}\|{\bf v}^{*}\|_{H^{2}(-R,R)^{|{\bf
N}|}}).$

Motivated from (2.6), define

$v^{(k)}({\bf
x}^{\prime},z,t)=\sum_{n_{1}=1}^{N_{1}}\dots\sum_{n_{d-1}}^{N_{d-1}}\sum_{n_{t}=1}^{N_{t}}v_{(n_{1},\dots,n_{d-1},n_{t})}(z)\Psi_{n_{1}}(x_{1})\dots\Psi_{n_{d-1}}(x_{d-1})\psi_{n_{t}}(t)$

for all $({\bf x},t)=({\bf x}^{\prime},z,t)\in\Omega\times(0,T)$. Due to
(2.4), set

$c^{(k)}({\bf x})=\frac{v^{(k)}({\bf x},0)-\Delta p({\bf x})}{p({\bf x})}$
(4.21)

for all ${\bf x}\in\Omega.$ By Remark 4.3 and Theorem 4.1, we can find
$\lambda_{1}$ as in Theorem 4.1 such that for all $\lambda>\lambda_{1}$, we
have

$\|c^{(k+1)}-c^{*}\|^{2}_{L^{2}(-R,R)}\leq
C_{1}\Big{[}\Big{(}\frac{C}{\lambda^{3}}\Big{)}^{k+1}\|{\bf v}^{(0)}-{\bf
v}^{*}\|^{2}_{L^{2}(-R,R)}+\frac{\delta^{2}C/\lambda^{3}}{1-C/\lambda^{3}}+\epsilon\|{\bf
v}^{*}\|_{H^{2}(-R,R)^{|{\bf N}|}}^{2}\Big{]}$ (4.22)

where $C_{1}$ is a constant depending only on $\lambda,$ $z_{0}$, $R$, ${\bf
N}$, and $[{\bf P}_{{\bf n}}]_{{\bf n}\leqq{\bf N}}.$ Here $c^{*}$ is the true
coefficient defined by

$c^{*}({\bf x})=\frac{v^{*}({\bf x},0)-\Delta p({\bf x})}{p({\bf x})}$

for all ${\bf x}\in\Omega$.

The analysis above leads to Algorithm 1 to solve the inverse problem under
consideration. Having the data $f$ and $g$ in hand, we can follow Algorithm 1
to compute a numerical solution to the inverse problem under consideration.

Algorithm 1 The procedure to compute the numerical solution to (2.18)

1: Choose cut-off numbers $N_{1},\dots,N_{d-1}$, and $N_{t}$. Set ${\bf
N}=(N_{1},N_{2},\dots,N_{d-1},N_{t})$

2: Choose Carleman parameters $z_{0},$ and $\lambda$, and a regularization
parameter $\epsilon$. Choose a large number $M$.

3: Set $n=0$. Choose an arbitrary initial solution ${\bf v}^{(0)}\in B_{M}.$

4: Compute ${\bf v}_{(k+1)}=\Phi^{(\lambda,\epsilon)}({\bf v}^{(k)})$ where
$\Phi^{(\lambda,\epsilon)}$ is defined in (4.5).

5:if $\|{\bf v}_{(k+1)}-{\bf v}_{(k)}\|_{L^{2}(-R,R)}>\kappa_{0}$ (for some
fixed number $\kappa_{0}>0$) then

6: Replace $k$ by $k+1.$

7: Go back to Step 4.

8:else

9: Set the computed solution ${\bf v}^{\rm comp}={\bf v}_{(k+1)}.$

10:end if

11: Write ${\bf v}^{\rm comp}=[v_{{\bf m}}^{(k+1)}]_{{\bf m}\leqq{\bf N}}$ and
set the desired solution as in (2.18).

12:Due to (2.6), compute  $v^{\rm comp}({\bf x},t)=v({\bf
x}^{\prime},z,t)=\sum_{n_{1}=1}^{N_{1}}\dots\sum_{n_{d-1}}^{N_{d-1}}\sum_{n_{t}=1}^{N_{t}}v_{(n_{1},\dots,n_{d-1},n_{t})}^{\rm
comp}(z)\Psi_{n_{1}}(x_{1})\dots\Psi_{n_{d-1}}(x_{d-1})\psi_{n_{t}}(t)$ for
${\bf x}\in\Omega$, $t\in(0,T).$

13: By (2.4), we obtain the reconstructed coefficient $c$ as  $c^{\rm
comp}({\bf x})=\frac{v^{\rm comp}({\bf x},0)-\Delta p({\bf x})}{p({\bf x})}$
(4.23) for all ${\bf x}\in\Omega$.

##  5 Numerical study

To illustrate our method, we numerically study the inverse problem, Problem
1.1, in 2D. That means, for simplicity, we set $d=2$. Let $\Omega=(-R,R)^{2}$
and therefore $\Omega=(-R,R)$. Rather than using the notation ${\bf
x}^{\prime}$, we write $x$ for an element of $\Omega^{\prime}.$ In this
section, we set $R=1$.

###  5.1 The forward problem

To generate simulated data, we numerically compute the solution to (1.1).
Since we only need the data on $\Gamma,$ a part of $\partial\Omega$ and
$t\in(0,T)$, it is not necessary to solve (1.1) on the whole unbounded domain
$\mathbb{R}^{d}\times(0,\infty)$. Rather, we choose a domain
$G=(-R_{1},R_{1})^{2}$ for some $R_{1}>R$, in which $\Omega$ is compactly
contained, and a positive number $T$. We solve

$\left\\{\begin{array}[]{rcll}u_{t}({\bf x},t)&=&\Delta u({\bf x},t)+c({\bf
x})u({\bf x},t)&({\bf x},t)\in G\times(0,T),\\\ u({\bf x},0)&=&p({\bf x})&{\bf
x}\in G\end{array}\right.$ (5.1)

by the explicit method in the finite difference scheme. More precisely, we
discretize $G$ by arranging a grid of points

$\mathcal{G}=\big{\\{}{\bf x}_{i,j}=(x_{i}=-R_{1}+(i-1)d_{\bf
x},z_{j}=-R_{1}+(j-1)d_{\bf x}):i,j=1,\dots,N_{\bf
x}^{1}\big{\\}}\subset\overline{G}.$

where $N_{\bf x}^{1}$ is an integer and $d_{\bf x}=2R_{1}/(N_{\bf x}^{1}-1)$.
On the time domain, we arrange the partition of $[0,T]$ as

$\mathcal{T}=\big{\\{}t_{l}=(l-1)d_{t}:l=1,\dots,N_{t}\big{\\}}$

where $N_{t}$ is an integer and $d_{t}=T/(N_{t}-1).$ We set $R=1$, $R_{1}=3$,
$T=0.5$, $N_{\bf x}^{1}=241$, and $N_{t}=4001$. The finite difference version
of (5.1) is read as

$\frac{u({\bf x}_{i,j},t_{l+1})-u({\bf x}_{i,j},t_{l})}{d_{t}}=\frac{u({\bf
x}_{i+1,j},t_{l})+u({\bf x}_{i-1,j},t_{l})+u({\bf x}_{i,j+1},t_{l})+u({\bf
x}_{i,j-1},t_{l})-4u({\bf x}_{i,j},t_{l})}{d_{\bf x}^{2}}\\\ +c({\bf
x}_{i,j})u({\bf x}_{i,j},t_{l})$ (5.2)

for all ${\bf x}_{i,j}\in\mathcal{G}$, $i,j\in\\{1,\dots,N_{\bf x}^{1}\\}$,
and $t_{l}\in\mathcal{T},$ $l\in\\{1,\dots,N_{t}\\}.$ It follows from (5.2)
that

$u({\bf x}_{i,j},t_{l+1})=u({\bf x}_{i,j},t_{l})+d_{t}\Big{[}\frac{u({\bf
x}_{i+1,j},t_{l})+u({\bf x}_{i-1,j},t_{l})+u({\bf x}_{i,j+1},t_{l})+u({\bf
x}_{i,j-1},t_{l})-4u({\bf x}_{i,j},t_{l})}{d_{\bf x}^{2}}\\\ +c({\bf
x}_{i,j})u({\bf x}_{i,j},t_{l})\Big{]}.$ (5.3)

So, given $u({\bf x}_{i,j},t_{1})=p({\bf x}_{i,j})$ for
$i,j\in\\{1,\dots,N_{\bf x}^{1}\\}$, we can compute $u({\bf x}_{i,j},t_{2})$
via (5.3), and then continue get $u({\bf x}_{i,j},t_{N_{t}})$. In our
computation, we set $p({\bf x})=2$ for all ${\bf x}\in G.$ We then easily
extract the noiseless data $f^{*}$ and $g^{*}$ on $\Gamma$. We pretend not to
know $f^{*}$ and $g^{*}$. We only know the noisy version of $f^{*}$ and
$g^{*}$ as

$f=f^{*}(1+\mbox{random numbers in }[-\delta,\delta])$ (5.4)

and

$g=g^{*}(1+\mbox{random numbers in }[-\delta,\delta])$ (5.5)

where $\delta$ is the noise level. In our computational program, we choose the
initial function $u({\bf x},0)=p({\bf x})=2$, for all ${\bf x}\in\Omega$,
$R_{1}=3$ and $T=0.5$.

###  5.2 The implementation of Algorithm 1

Step 1 of Algorithm 1. In 2D, we only need we determine $N_{1}$ and $N_{t}$.
These numbers are chosen by a data-driven procedure. Recall that the data is
measured on $\Gamma$ defined as in (1.2). The measurement set consists of two
parts, namely $\Gamma=\Gamma^{+}\cup\Gamma^{-}$ where

$\Gamma^{+}=\\{{\bf x}=(x,z):x\in(-R,R)\mbox{and }z=R\\}$

and

$\Gamma^{-}=\\{{\bf x}=(x,z):x\in(-R,R)\mbox{and }z=-R\\}.$

One can examine the approximation formula

$u(x,z,t)\simeq\sum_{n_{1}=1}^{N_{1}}\sum_{n_{t}=1}^{N_{t}}u_{(n_{1},n_{t})}(z){\bf
P}_{(n_{1},n_{t})}(x,t)$ (5.6)

for the data $u({\bf x},t)=u(x,z,t)=f({\bf x},t)$ on either $\Gamma^{+}$ or
$\Gamma^{-}$. Note that (5.6) is an analog of (2.6) where $v$ is replaced by
the function $u$. We do so on $\Gamma^{-}.$ Define the mismatch function

$\varphi_{(N_{1},N_{t})}(x,t)=\Big{|}u(x,-R,t)-\sum_{n_{1}=1}^{N_{1}}\sum_{n_{t}=1}^{N_{t}}u_{(n_{1},n_{t})}(-R){\bf
P}_{(n_{1},n_{t})}(x,t)\Big{|}$

for $(x,t)\in\Gamma^{-}\times(0,T).$ Then, we test the smallness of
$\varphi(N_{1},N_{t})$ by manually and gradually increasing $N_{1}$ and
$N_{t}$ until $\varphi(N_{1},N_{t})$ is sufficiently small. For example, we
take the data in the numerical test 1 below. We increase those two numbers so
that $\|\varphi_{(N_{1},N_{t})}\|_{L^{\infty}(\Gamma^{-}\times(0,T))}<5\times
10^{-4}$. By this, we find $N_{1}=15$ and $N_{t}=10$. We chose these two
numbers for all of our numerical tests. In Figure 1, we display the graphs of
the function $\varphi_{(N_{1},N_{t})}$ for some values of $N_{1}$ and $N_{t}.$

(a) $\varphi_{(5,5)}$

(b) $\varphi_{(15,8)}$

(c) $\varphi_{(15,10)}$

Figure 1:  The difference $\varphi_{(N_{1},N_{t})}$ of the data $u(x,-R,t)$
and its truncated Fourier approximation
$\sum_{n_{1}=1}^{N_{1}}\sum_{n_{t}=1}^{N_{t}}u_{(n_{1},n_{t})}(-R){\bf
P}_{(n_{1},n_{t})}(x,t)$, for some values of $(N_{1},N_{t})$. The data is
taken in Test 1 below. Evidently, for Step 1 in Algorithm 1 during the
calculation of the desired coefficient in Test 1, we can choose $N_{1}=15$ and
$N_{t}=10.$

Step 2 of Algorithm 1. In Step 2, we select the parameters through a trial-
error procedure involving experimentation and adjustments. To do so, we begin
with a reference test where the accurate solution is known. Using this
reference, we adjust the parameters until Algorithm 1 produces a satisfactory
numerical outcome with data free of noise, i.e., $\delta=0$. These chosen
parameters are then applied across all subsequent tests and various noise
levels $\delta$. The reference for our adjustments is Test 1, mentioned below.
In all our numerical analyses, these values are set as $\epsilon=10^{-6.5}$,
$\kappa_{0}=10^{-3}$, $z_{0}=-10$, and $\lambda=10$.

Step 3 of Algorithm 1. We need to choose a function ${\bf v}^{(0)}inH$ . A
convenient method to compute this function is solving the linear problem
obtained by excluding the nonlinearity $\mathcal{F}$. More precisely, we solve
the following problem

$\left\\{\begin{array}[]{ll}{{\bf
v}^{(0)}}^{\prime\prime}(z)+\mathcal{S}::{\bf v}^{(0)}(z)=0&z\in(-R,R),\\\
{\bf v}^{(0)}(z)=\mathcal{P}(z)&z=\pm R,\\\ {{\bf
v}^{(0)}}^{\prime}(z)=\mathcal{Q}(z)&z=\pm R\end{array}\right.$ (5.7)

for the function ${\bf v}^{(0)}.$ We can compute the numerical solution to
(5.7) by the quasi-reversibility method. We do not present the numerical
implementation to solve linear PDEs using the Carleman quasi-reversibility
method in this paper. The reader can find the details about this in [32, 39,
44].

Step 4 of Algorithm 1. In Step 4, we minimize $J_{{\bf v}^{(k)}}$ and set
${\bf v}^{(k+1)}$ as its minimizer. The obtained minimizer ${\bf v}^{(k+1)}$
is actually the solution to

$\left\\{\begin{array}[]{ll}{{\bf
v}^{(k+1)}}^{\prime\prime}(z)+\mathcal{S}::{\bf v}^{(k+1)}(z)+\mathcal{F}({\bf
v}^{(k)}(z))=0&z\in(-R,R),\\\ {\bf v}^{(k+1)}(z)=\mathcal{P}(z)&z=\pm R,\\\
{{\bf v}^{(k+1)}}^{\prime}(z)=\mathcal{Q}(z)&z=\pm R.\end{array}\right.$ (5.8)

The details in implementation to compute the regularized solution ${\bf
v}^{(k+1)}$ to (5.8) were presented in [32, 39, 44], in which we employ the
optimization package already built in Matlab. We do not repeat it here.

Steps 5–13. Executing these steps is direct and not complicated. For brevity,
we do not present them here.

###  5.3 Numerical examples

We show in this section three (3) numerical examples computed by Algorithm 1.

Test 1. In test 1, the true coefficient $c$ is given by

$c^{\rm
true}(x,z)=\left\\{\begin{array}[]{ll}e^{\frac{0.35x^{2}+(z-0.4)^{2}}{0.55^{2}-(0.35x^{2}+(z-0.4)^{2})}}&0.35x^{2}+(z-0.4)^{2}<0.55^{2},\\\
0&\mbox{otherwise}.\end{array}\right.$

for $(x,z)\in\Omega.$ It is characterized by an “ellipse” inclusion. The true
and computed coefficient $c$ are displayed in Figure 2.

(a) The true coefficient $c^{\rm true}$

(b) The computed coefficient $c^{\rm comp}$

(c) The consecutive relative error $\displaystyle\frac{\|{\bf v}^{(k)}-{\bf
v}^{(k-1)}\|_{L^{\infty}}}{\|{\bf v}^{(k)}\|_{L^{\infty}}}$.

Figure 2:  (a) The true coefficient, (b) the computed one from data on
$\Gamma$ corrupted with $5\%$ of noise, and (c) the relative difference of the
reconstructed of the solution ${\bf v}$ to (2.18) at the $k^{\rm th}$
iteration. Although the data are missing in $\partial\Omega\setminus\Gamma$,
the reconstruction is acceptable. The convergence of our method is numerically
confirmed.

It is evident that Algorithm 1 generates a satisfactory numerical solution. It
is evident that the “ellipse inclusion” was successfully detected. The maximum
value of the function $c$ inside the inclusion is $1$. The constructed value
is $0.952$ (relative error 4.8%). Due to Figure 2c, the stopping criterium of
Algorithm 1 is met after only seven iterations.

Test 2. We test the case when the true coefficient is characterized by two
horizontal inclusions. More precisely, we set

$c^{\rm
true}(x,z)=\left\\{\begin{array}[]{ll}1&\max\\{0.25|x|,4|z-0.6|\\}<0.8,\\\
1&\max\\{0.25|x|,4|z+0.6|\\}<0.8,\\\
0&\mbox{otherwise},\end{array}\right.\quad\mbox{for all }(x,z)\in\Omega.$

The true and computed solutions to Problem 1.1 are shown in Figure 3

(a) The true coefficient $c^{\rm true}$

(b) The computed coefficient $c^{\rm comp}$

(c) The consecutive relative error $\displaystyle\frac{\|{\bf v}^{(k)}-{\bf
v}^{(k-1)}\|_{L^{\infty}}}{\|{\bf v}^{(k)}\|_{L^{\infty}}}$.

Figure 3:  (a) The true coefficient, (b) the computed one from data on
$\Gamma$ corrupted with $5\%$ of noise, and (c) the relative difference of the
reconstructed of the solution ${\bf v}$ to (2.18) at the $k^{\rm th}$
iteration. ”Despite the absence of data in $\partial\Omega\setminus\Gamma$,
the reconstruction remains satisfactory. Our method’s fast convergence has
been verified numerically.

As in Test 1, we can see that Algorithm 1 provides a satisfactory numerical
solution. It is evident that both horizontal inclusions were successfully
identified. The maximum value of the function $c$ inside each inclusion is
$1$. The constructed value of $c$ inside the upper inclusion is $0.963$
(relative error 3.7%), and the one inside the lower inclusion is $0.875$
(relative error 12.4%). Due to Figure 3c, Algorithm 1 stops at only nine
iterations.

Test 3. Test 3 checks the case when the true coefficient $c$ has a $T$
inclusion. That means the function $c^{\rm true}$ takes the value 1 inside a
letter $T$ and $0$ otherwise. We refer the reader to Figure 4 for the image of
the true and computed coefficients.

(a) The true coefficient $c^{\rm true}$

(b) The computed coefficient $c^{\rm comp}$

(c) The consecutive relative error $\displaystyle\frac{\|{\bf v}^{(k)}-{\bf
v}^{(k-1)}\|_{L^{\infty}}}{\|{\bf v}^{(k)}\|_{L^{\infty}}}$.

Figure 4:  (a) The true coefficient, (b) the computed one from data on
$\Gamma$ corrupted with $5\%$ of noise, and (c) the relative difference of the
reconstructed of the solution ${\bf v}$ to (2.18) at the $k^{\rm th}$
iteration. ”Despite the absence of data in $\partial\Omega\setminus\Gamma$,
the reconstruction of the letter $T$ is acceptable. Our method’s fast
convergence has been verified numerically.

As in Test 1 and Test 2, we can see that the $T$ inclusion was successfully
found. The maximum value of the function $c$ inside each inclusion is $1$. The
constructed value of $c$ inside the inclusion is $0.98$ (relative error 2%).
Due to Figure 4c, Algorithm 1 stops at only 13 iterations.

######  Remark 5.1.

The computational cost of Algorithm 1 is not expensive. In fact, we only need
to solve several 1D linear over-determined problems in Steps 5-7. We have used
a personal computer iMac with a 3.2GHz Intel Core i5 Processor and memory of
24GB, not a professional workstation, to compute the numerical solutions
above. It took 46.93, 56.56, and 77.12 seconds to complete all computational
tasks for the inverse problem in Test 1, Test 2, and Test 3, respectively.

######  Remark 5.2.

It is worth mentioning that the images of the inclusions in the tests above
are not perfect. One can find some artifacts occurring in the numerical
results. However, these errors are acceptable because we solve the inverse
problem when the data are given only on the part $\Gamma$ of $\partial\Omega.$

##  6 Concluding remarks

In this study, we provide a numerical method to solve a nonlinear coefficient
inverse problem for parabolic equations. This inverse problem has numerous
practical applications. To obtain these solutions, we employ the polynomial-
differential basis, converting the inverse problem into a set of 1D nonlinear
equations. Next, we introduce a method to address this nonlinear system. Our
method is based on a combination of the Picard iteration, the quasi-
reversibility method, and a Carleman estimate. We provide analytic proof of
the method’s convergence and demonstrate its effectiveness through some
numerical examples.

## Acknowledgement

The works of RA, TTL, and LHN were partially supported by National Science
Foundation grant DMS-2208159, by funds provided by the Faculty Research Grant
program at UNC Charlotte Fund No. 111272, and by the CLAS small grant provided
by the College of Liberal Arts & Sciences, UNC Charlotte. The works of CP were
supported in part by National Science Foundation grant DMS-2150179.

## References

  * [1] A. Abhishek, T. T. Le, L. H. Nguyen, and T. Khan.  The Carleman-Newton method to globally reconstruct a source term for nonlinear parabolic equation.  preprint, arXiv:2209.08011, 2022.
  * [2] 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.
  * [3] E. Bécache, L. Bourgeois, L. Franceschini, and J. Dardé.  Application of mixed formulations of quasi-reversibility to solve ill-posed problems for heat and wave equations: The 1d case.  Inverse Problems & Imaging, 9(4):971–1002, 2015.
  * [4] L. Beilina and M. V. Klibanov.  Approximate Global Convergence and Adaptivity for Coefficient Inverse Problems.  Springer, New York, 2012.
  * [5] L. Borcea, V. Druskin, A. V. Mamonov, and M. Zaslavsky.  A model reduction approach to numerical inversion for a parabolic partial differential equation.  Inverse Problems, 30:125011, 2014.
  * [6] L. Bourgeois.  Convergence rates for the quasi-reversibility method to solve the Cauchy problem for Laplace’s equation.  Inverse Problems, 22:413–430, 2006.
  * [7] L. Bourgeois and J. Dardé.  A duality-based method of quasi-reversibility to solve the Cauchy problem in the presence of noisy data.  Inverse Problems, 26:095016, 2010.
  * [8] L. Bourgeois, D. Ponomarev, and J. Dardé.  An inverse obstacle problem for the wave equation in a finite time domain.  Inverse Probl. Imaging, 13(2):377–400, 2019.
  * [9] 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.
  * [10] K. Cao and D. Lesnic.  Determination of space-dependent coefficients from temperature measurements using the conjugate gradient method.  Numer Methods Partial Differential Eq., 34:1370–1400, 2018.
  * [11] K. Cao and D. Lesnic.  Simultaneous reconstruction of the perfusion coefficient and initial temperature from time-average integral temperature measurements.  Applied Mathematical Modelling, 68:523–539, 2019.
  * [12] C. Clason and M. V. Klibanov.  The quasi-reversibility method for thermoacoustic tomography in a heterogeneous medium.  SIAM J. Sci. Comput., 30:1–23, 2007.
  * [13] J. Dardé.  Iterated quasi-reversibility method applied to elliptic and parabolic data completion problems.  Inverse Problems and Imaging, 10:379–407, 2016.
  * [14] V. Isakov.  Some inverse problems for the diffusion equation.  Inverse Problems, 15(1):3–10, 1999.
  * [15] B. Kaltenbacher and W. Rundell.  Regularization of a backwards parabolic equation by fractional operators.  Inverse Probl. Imaging, 13(2):401–430, 2019.
  * [16] Y. L. Keung and J. Zou.  Numerical identifications of parameters in parabolic systems.  Inverse Problems, 14:83–100, 1998.
  * [17] 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.
  * [18] 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\.
  * [19] 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.
  * [20] 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.
  * [21] M. V. Klibanov.  Carleman estimates for the regularization of ill-posed Cauchy problems.  Applied Numerical Mathematics, 94:46–74, 2015.
  * [22] M. V. Klibanov.  Carleman weight functions for solving ill-posed Cauchy problems for quasilinear PDEs.  Inverse Problems, 31:125007, 2015.
  * [23] M. V. Klibanov.  Convexification of restricted Dirichlet to Neumann map.  J. Inverse and Ill-Posed Problems, 25(5):669–685, 2017.
  * [24] 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.
  * [25] M. V. Klibanov, L. H. Nguyen, and H. V. Tran.  Numerical viscosity solutions to Hamilton-Jacobi equations via a Carleman estimate and the convexification method.  Journal of Computational Physics, 451:110828, 2022.
  * [26] M. V. Klibanov and F. Santosa.  A computational quasi-reversibility method for Cauchy problems for Laplace’s equation.  SIAM J. Appl. Math., 51:1653–1675, 1991.
  * [27] M. V. Klibanov and A. G. Yagola.  Convergent numerical methods for parabolic equations with reversed time via a new Carleman estimate.  preprint, 2019.
  * [28] R. Lattès and J. L. Lions.  The Method of Quasireversibility: Applications to Partial Differential Equations.  Elsevier, New York, 1969.
  * [29] P. N. H. Le, T. T. Le, and L. H. Nguyen.  The Carleman convexification method for Hamilton-Jacobi equations on the whole space.  preprint, arXiv:2206.09824, 2022.
  * [30] T. T. Le.  Global reconstruction of initial conditions of nonlinear parabolic equations via the Carleman-contraction method.  In D-L. Nguyen, L. H. Nguyen, and T-P. Nguyen, editors, Advances in Inverse problems for Partial Differential Equations, volume 784 of Contemporary Mathematics, pages 23–42. American Mathematical Society, 2023.
  * [31] T. T. Le, V. A. Khoa, M. V. Klibanov, L. H. Nguyen, G. W. Bidney, and V. N. Astratov.  Numerical verification of the convexification method for a frequency-dependent inverse scattering problem with experimental data.  to appear in Journal of Applied and Industrial Mathematics, preprint arXiv:2306.00761, 2023.
  * [32] 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,, 30(2):265–286, 2022\.
  * [33] 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.  Journal of Scientific Computing, 91(3):74, 2022.
  * [34] T. T. Le, L. H. Nguyen, T-P. Nguyen, and W. Powell.  The quasi-reversibility method to numerically solve an inverse source problem for hyperbolic equations.  Journal of Scientific Computing, 87:90, 2021.
  * [35] T. T. Le, L. H. Nguyen, and H. V. Tran.  A Carleman-based numerical method for quasilinear elliptic equations with over-determined boundary data and applications.  Computers and Mathematics with Applications, 125:13–24, 2022.
  * [36] Q. Li and L. H. Nguyen.  Recovering the initial condition of parabolic equations from lateral Cauchy data via the quasi-reversibility method.  Inverse Problems in Science and Engineering, 28:580–598, 2020.
  * [37] H. T. Nguyen, V. A. Khoa, and V. A. Vo.  Analysis of a quasi-reversibility method for a terminal value quasi-linear parabolic problem with measurements.  SIAM Journal on Mathematical Analysis, 51:60–85, 2019.
  * [38] L. H. Nguyen.  An inverse space-dependent source problem for hyperbolic equations and the Lipschitz-like convergence of the quasi-reversibility method.  Inverse Problems, 35:035007, 2019.
  * [39] 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\.
  * [40] L. H. Nguyen.  The Carleman contraction mapping method for quasilinear elliptic equations with over-determined boundary data.  Acta Mathematica Vietnamica, DOI: https://doi.org/10.1007/s40306-023-00500-w, 2023.
  * [41] L. H. Nguyen and M. V. Klibanov.  Carleman estimates and the contraction principle for an inverse source problem for nonlinear hyperbolic equations.  Inverse Problems, 38:035009, 2022.
  * [42] L. H. Nguyen, Q. Li, and M. V. Klibanov.  A convergent numerical method for a multi-frequency inverse source problem in inhomogenous media.  Inverse Problems and Imaging, 13:1067–1094, 2019.
  * [43] P. M. Nguyen, T. T. Le, L. H. Nguyen, and M. V. Klibanov.  Numerical differentiation by the polynomial-exponential basis.  to appear in Journal of Applied and Industrial Mathematics, preprint arXiv:2304.05909, 2023.
  * [44] 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.
  * [45] A. I. Prilepko and A. B. Kostin.  On certain inverse problems for parabolic equations with final and integral observation.  Russ. Acad. Sci. Sb. Math., 75:473–490, 1993.
  * [46] A. I. Prilepko, D. G. Orlovsky, and I. A. Vasin.  Methods for solving inverse problems in mathematical physics, volume 321.  Pure and Applied Mathematics, Marcel Dekker, New Youk, 2000.
  * [47] A. V. Smirnov, M. V. Klibanov, and L. H. Nguyen.  On an inverse source problem for the full radiative transfer equation with incomplete data.  SIAM Journal on Scientific Computing, 41:B929–B952, 2019.
  * [48] N. H. Tuan, V. V. Au, V. A. Khoa, and D. Lesnic.  Identification of the population density of a species model with nonlocal diffusion and nonlinear reaction.  Inverse Problems, 33:055019, 2017.
  * [49] L. Yang, J-N. Yu, and Y-C. Deng.  An inverse problem of identifying the coefficient of parabolic equation.  Applied Mathematical Modelling, 32:1984–1995, 2008.