# Convexification Numerical Method for a Coefficient Inverse Problem for the
Radiative Transport Equation

Michael V. Klibanov  Department of Mathematics and Statistics, University of
North Carolina at Charlotte, Charlotte, NC, 28223, USA<EMAIL_ADDRESS>,
Jingzhi Li  Department of Mathematics & National Center for Applied
Mathematics Shenzhen & SUSTech International Center for Mathematics, Southern
University of Science and Technology, Shenzhen 518055, P. R. China
<EMAIL_ADDRESS>,  Loc H. Nguyen  Department of Mathematics and
Statistics, University of North Carolina at Charlotte, Charlotte, NC, 28223,
USA<EMAIL_ADDRESS>and  Zhipeng Yang  Department of Mathematics,
Southern University of Science and Technology, Shenzhen 518055, P. R. China
<EMAIL_ADDRESS>
###### Abstract.

An $\left(n+1\right)-$D coefficient inverse problem for the radiative
stationary transport equation is considered for the first time. A globally
convergent so-called convexification numerical method is developed and its
convergence analysis is provided. The analysis is based on a Carleman
estimate. In particular, convergence analysis implies a certain uniqueness
theorem. Extensive numerical studies in the 2-D case are presented.

###### Key words and phrases:

coefficient inverse problem, radiative transport equation, convexification
method, a special orthonormal basis, global convergence analysis, numerical
experiments

###### 2020 Mathematics Subject Classification:

Primary 35R30; Secondary 65M32

MK and LN were supported by US Army Research Laboratory and US Army Research
Office Grant W911NF-19-1-0044. LN was supported in part by NSF Grant
#DMS-2208159 and by funds provided by the Faculty Research Grant program at
UNC Charlotte Fund # 111272. JL was partially supported by the NSF of China
No. 11971221, Guangdong NSF Major Fund No. 2021ZDZX1001, the Shenzhen Sci-Tech
Fund No. RCJC20200714114556020, JCYJ20200109115422828 and
JCYJ20190809150413261.

##  1\. Introduction

The stationary radiative transfer equation (RTE) governs the propagation of
the radiation field in media with absorbing, emitting and scattering
radiation. RTE has broad applications in optics, including diffuse optical
tomography [17], astrophysics, atmospheric science and other applied
disciplines. For example, in the single particle emission tomography the
coefficient we reconstruct is the emission coefficient [29, formula (2.1)].
The RTE we consider here is a more general one since we introduce an integral
operator in it.

For the first time, we develop in this paper a globally convergent numerical
method for a Coefficient Inverse Problem (CIP) of the recovery of a spatially
distributed coefficient of RTE in the $\left(n+1\right)-$D case, $n\geq 1$.
All CIPs are both nonlinear and ill-posed. Numerical methods for inverse
source problems for RTE were developed in [11, 12, 13, 31]. In the case of
single particle emission tomography, i.e. when the kernel of the integral
operator in RTE is the identical zero, an inversion formula for the inverse
source problem was first derived in [30] and tested numerically in [16].
Inverse problems of [11, 12, 13, 16, 30, 31] are linear ones.

Our CIP is formally determined, i.e. the number of free variables in the data
equals the number of free variables in the unknown coefficient. Our data are
incomplete, i.e. the source runs along an interval of a straight line and the
data are measured only at a part of the boundary $\partial\Omega$ of the
domain of interest $\Omega$. This is unlike the classical case of X-ray
tomography when the source runs all around $\Omega$ and the data are measured
on the whole boundary $\partial\Omega$.

A significant attention has been paid to the questions of uniqueness and
stability of CIPs for RTE, see, e.g. [2, 14, 21, 28, 32] and the references
cited therein. We refer here only to a limited number of works on this topic
since our interest in this paper is in a numerical development.

Our goal is to construct a globally convergent numerical method for our CIP,
to conduct its convergence analysis and to confirm our method via numerical
experiments. To achieve this goal, we develop here a new version of the so-
called convexification method [26]. The convexification constructs a least
squares cost functional, which is strictly convex on a convex set in an
appropriate Hilbert space. The diameter of this set is fixed and is an
arbitrary one. We prove existence and uniqueness of the minimizer of that
functional on that set and estimate convergence rate of minimizers to the true
solution depending on the level of the noise in the data. As a by-product, we
obtain a certain uniqueness result for our CIP, see Theorem 4.3 below. Also,
we establish the global convergence of the gradient descent method of the
minimization of our functional. Recall that only local convergence of this
method can be proven in the non-convex case. In addition, we present numerical
experiments in the 2-D case.

We call a numerical method for a CIP  _globally convergent_ if a theorem is
proven, which claims that this method delivers at least one point in a
sufficiently small neighborhood of the solution of this CIP without an
advanced knowledge of that neighborhood, also, see [26, Definiion 1.4.2] for a
similar statement. In other words, a good approximation for the true solution
can be obtained without a knowledge of a good first guess about this solution.

Conventional numerical methods for CIPs are based on the minimizations of
least squares cost functionals, see, e.g. [9, 15]. However, these functionals
are usually non-convex. This means, in turn that they typically have multiple
local minima and ravines. On the other hand, since any gradient-like
optimization method can stop at any point of a local minimum, then this
technique needs to have a good first guess about the solution. Besides, there
is no rigorous guarantee of neither the existence of the global minimizer nor
that such a minimizer, if it exists, is indeed close to the true solution.
These considerations have prompted the first author to work on the
developments of the convexification technique with the first publications [24,
20].

The key element of the convexification is the presence of the so-called
Carleman Weight Function (CWF) in the above mentioned functional. The presence
of the CWF ensures the strict convexity property of that functional on that
bounded set. The CWF is the function, which is used as the weight in the
Carleman estimate for the corresponding differential operator, see, e.g. books
[5, 26] for Carleman estimates. In particular, we demonstrate in our numerical
experiments that the absence of the CWF leads to a significant deterioration
of the numerical solution.

The convexification uses the idea of the so-called Bukhgeim-Klibanov method
(BK), which is based on Carleman estimates. BK was originally introduced in
the field of CIPs in 1981 [8] only for the proofs of uniqueness theorems for
multidimensional CIPs. The work [8] has generated many publications of many
authors since then, see, e.g. books [5, 26], papers [14, 21, 22, 28] for a few
of such publications as well as references cited therein. The majority of
currently known publications about BK is also dedicated to the issues of
uniqueness and stability of CIPs. Unlike these, it was originally proposed in
[24, 20] to use the idea of BK for the construction of the above mentioned
globally strictly convex cost functionals for CIPs.

While initially only analytical results for the convexification were derived,
an active exploration of this method for numerical studies of CIPs has started
from the work [1], which has removed some obstacles for real computations,
see, e.g., publications [19, 27, 26] for a few examples out of many. In
particular, [19] and [26, Chapter 10] demonstrate the performance of the
convexification for experimentally collected data for backscatter microwaves.
An updated version of the convexification, which combines Carleman estimates
with the contraction principle, can be found in [3, 4, 7].

In our approach, we use a certain approximate mathematical model. Our model
amounts to the truncation of a Fourier-like series with respect to a special
orthonormal basis in $L_{2}\left(a,b\right)$ as well as to the so-called
“partial finite differences”, see subsection 3.3 for details. This basis was
originally introduced in [23], also see [26, Chapter 6]. We do not know how to
prove convergence of our method in the case when the number of terms of this
series $N\rightarrow\infty.$ Therefore, the _only way_ to verify the validity
of this model is via numerical simulations. In this regard, we note that
truncations of Fourier-like series with respect to the same basis were done
for a variety of CIPs in [19, 27, 31] and in [26, Chapters 7,10,12]. In each
of these cases, the validity of the corresponding approximate mathematical
model was verified numerically well. Similar cases of truncated Fourier series
without proofs of convergence at $N\rightarrow\infty$ can be observed for
inverse problems considered by some other authors, see, e.g. [16, 18]. And
successful numerical verifications also took place in these references.

Philosophically, the situation of our approximate mathematical model being
verified in numerical studies is somewhat similar with the well known
situation of the Huygens-Fresnel diffraction theory in optics. This theory is
not yet derived rigorously from the Maxwell’s system. This is why it is stated
in section 8.1 of the classical textbook of Born and Wolf [6] that “ _because
of mathematical difficulties, approximate mathematical models must be used in
most cases of practical interest. Of these the theory of Huygens and Fresnel
is by far most powerful and adequate for the treatment of the majority of
problems encountered in experimental optics._ ” The conclusion we draw from
this is that in the case of significant challenges for a certain applied
mathematical problem, it is reasonable to introduce an approximate
mathematical model. However, this model must be verified numerically.

All functions considered below are real valued ones. Let $B$ be a Banach space
of functions and $s\geq 2$ be an integer. Below $B_{s}=\underbrace{B\times
B\times...\times B},$ $s$ terms, is the space of $s-$dimensional vector
functions generated by $B$ and with the obvious extension of the $B-$norm.

In section 2 we pose forward and inverse problems and prove existence and
uniqueness theorem for the solution of the forward problem. In section 3 we
describe our transformation procedure. In section 4 we introduce our
convexification functional and formulate five theorems. These theorems are
proven in section 5. In section 6 we present our numerical results.

##  2\. Statements of Forward and Inverse Problems

For $n\geq 1,$ points in $\mathbb{R}^{n+1}$ are denoted below as
$\mathbf{x}=(x_{1},x_{2},...,x_{n-1},y)\in\mathbb{R}^{n+1}.$ Let numbers
$A,a,b,d>0$, where

(2.1) $1<a<b.$

Define the rectangular prism $\Omega\subset\mathbb{R}^{n+1}$ and parts
$\partial_{1}\Omega,\partial_{2}\Omega,\partial_{3}\Omega$ of its boundary
$\partial\Omega$ as:

(2.2) $\Omega=\\{\mathbf{x}:-A<x_{1},...,x_{n}<A,a<y<b\\},$ (2.3)
$\partial_{1}\Omega=\left\\{\mathbf{x}:-A<x_{1},...,x_{n}<A,y=a\right\\},$
(2.4)
$\partial_{2}\Omega=\left\\{\mathbf{x:}-A<x_{1},...,x_{n}<A,y=b\right\\},$
(2.5) $\partial_{3}\Omega=\left\\{x_{i}=\pm
A,y\in\left(a,b\right),i=1,...,n,\right\\}.$

Let $\Gamma_{d}$ be the line where the external sources are,

(2.6) $\Gamma_{d}=\\{\mathbf{x}_{\alpha}=(\alpha,0,...,0):\alpha\in[-d,d]\\}.$

Hence, $\Gamma_{d}$ is a part of the $x_{1}-$axis. It follows from (2.1) and
(2.2) that $\Gamma_{d}\cap\overline{\Omega}=\varnothing.$

Let the points of external sources $\mathbf{x}_{\alpha}$ run along
$\Gamma_{d},$ $\mathbf{x}_{\alpha}\in\Gamma_{d}$. Let $\varepsilon>0$ be a
sufficiently small number. To avoid dealing with singularities, we model the
$\delta\left(\mathbf{x}\right)-$function as:

(2.7)
$f\left(\mathbf{x}\right)=C_{\varepsilon}\left\\{\begin{array}[]{c}\exp\left(\frac{\left|\mathbf{x}\right|^{2}}{\varepsilon^{2}-\left|\mathbf{x}\right|^{2}}\right),\left|\mathbf{x}\right|<\varepsilon,\\\
0,\left|\mathbf{x}\right|\geq\varepsilon\end{array}\right.,$

where the constant $C_{\varepsilon}$ is chosen such that

(2.8)
$C_{\varepsilon}\mathop{\displaystyle\int}\limits_{\left|\mathbf{x}\right|<\varepsilon}\exp\left(\frac{\left|\mathbf{x}\right|^{2}}{\varepsilon^{2}-\left|\mathbf{x}\right|^{2}}\right)d\mathbf{x}=1.$

Hence, the function
$f\left(\mathbf{x}-\mathbf{x}_{\alpha}\right)=f\left(x_{1}-\alpha,x_{2},...,x_{n},y\right)\in
C^{\infty}\left(\mathbb{R}^{n+1}\right)$ plays the role of the source function
for the source $\mathbf{x}_{\alpha}.$ We choose $\varepsilon$ so small that

(2.9)
$f\left(\mathbf{x}-\mathbf{x}_{\alpha}\right)=0,\forall\mathbf{x}\in\overline{\Omega},\forall\mathbf{x}_{\alpha}\in\Gamma_{d}.$

Let $u(\mathbf{x},\alpha)$ denotes the steady-state radiance at the point
$\mathbf{x}$ generated by the source function
$f\left(\mathbf{x}-\mathbf{x}_{\alpha}\right).$ We assume that the function
$u(\mathbf{x},\alpha)$ is governed by the RTE of the following form [17]:

$\nu(\mathbf{x,}\alpha)\cdot\nabla_{\mathbf{x}}u(\mathbf{x},\alpha)+a\left(\mathbf{x}\right)u(\mathbf{x},\alpha)$
(2.10)
$=\mu_{s}(\mathbf{x})\int_{\Gamma_{d}}K(\mathbf{x},\alpha,\beta)u(\mathbf{x},\beta)d\beta+f\left(\mathbf{x}-\mathbf{x}_{\alpha}\right),\text{
}\mathbf{x}\in\mathbb{R}^{n+1},\mathbf{x}_{\alpha}\in\Gamma_{d}.$

The kernel $K(\mathbf{x},\alpha,\beta)$ of the integral operator in (2.10) is
called the “phase function”,

(2.11) $K(\mathbf{x},\alpha,\beta)\geq 0,\mathbf{x\in}\overline{\Omega};\text{
}\alpha,\beta\in\left[-d,d\right],\text{ }$

see [17]. In addition, we assume that

(2.12) $K(\mathbf{x},\alpha,\beta)\in
C^{2}\left(\overline{\Omega}\times\left[-d,d\right]^{2}\right).$

In equation (2.10),

(2.13)
$a\left(\mathbf{x}\right)=\mu_{a}\left(\mathbf{x}\right)+\mu_{s}(\mathbf{x}),$

where $\mu_{a}\left(\mathbf{x}\right)$ and $\mu_{s}(\mathbf{x})$ are the
absorption and scattering coefficients respectively. As stated in section 1,
$a\left(\mathbf{x}\right)$ is the emission coefficient [29, formula (2.1)]. We
assume that

(2.14) $\mu_{a}\left(\mathbf{x}\right),\mu_{s}(\mathbf{x})\geq
0,\mu_{a}\left(\mathbf{x}\right)=\mu_{s}(\mathbf{x})=0\text{,
}\mathbf{x}\in\mathbb{\mathbb{R}}^{n+1}\setminus\Omega,$ (2.15)
$\mu_{a}\left(\mathbf{x}\right),\mu_{s}(\mathbf{x})\in
C^{2}\left(\mathbb{\mathbb{R}}^{n+1}\right).$

For two arbitrary points $\mathbf{x},\mathbf{z}\in\mathbb{R}^{n+1}$ let
$L\left(\mathbf{x},\mathbf{z}\right)$ be the line segment connecting these
points and let $ds$ be the element of the euclidean length on
$L\left(\mathbf{x},\mathbf{z}\right).$ In (2.10) $\nu(\mathbf{x,}\alpha)$
denotes the unit vector, which is parallel to
$L\left(\mathbf{x},\mathbf{x}_{\alpha}\right),$

(2.16)
$\nu(\mathbf{x,}\alpha)=\frac{\mathbf{x}-\mathbf{x}_{\alpha}}{\left|\mathbf{x}-\mathbf{x}_{\alpha}\right|}=\frac{\left(x_{1}-\alpha,x_{2},...,x_{n},y\right)}{\sqrt{\left(x_{1}-\alpha\right)^{2}+x_{2}^{2}+...+x_{n}^{2}+y^{2}}}.$

Denote

(2.17)
$D^{n+1}=\left\\{\left(x_{1},...,x_{n},y\right)\in\mathbb{R}^{n}\times\left[0,b\right]\right\\},\text{
}D_{a}^{n+1}=D^{n+1}\cap\left\\{y>a\right\\}.$

Forward Problem. _Let ( 2.1)-(2.16) hold. Find the function _

$u\left(\mathbf{x},\alpha\right)\in
C^{1}\left(\overline{D^{n+1}}\times\left[-d,d\right]\right)$_satisfying
equation ( 2.10) and the initial condition_

(2.18) $u(\mathbf{x}_{\alpha},\mathbf{x}_{\alpha})=0\text{ for
}\mathbf{x}_{\alpha}\in\Gamma_{d}.$

Coefficient Inverse Problem (CIP). _Let ( 2.1)-(2.16) hold. Let the function
_$u\left(\mathbf{x},\alpha\right)\in
C^{1}\left(\overline{D^{n+1}}\times\left[-d,d\right]\right)$_be the solution
of the Forward Problem. Assume that the coefficient_
$a\left(\mathbf{x}\right)$_of equation ( 2.10) is unknown. Determine the
function _$a\left(\mathbf{x}\right),$_assuming that the following function_
$g\left(\mathbf{x},\alpha\right)$_is known:_

(2.19)
$g\left(\mathbf{x},\alpha\right)=u\left(\mathbf{x},\alpha\right),\forall\mathbf{x}\in\partial\Omega\diagdown\partial_{1}\Omega,\forall\alpha\in\left(-d,d\right).$

First, we formulate and prove an existence and uniqueness theorem for the
solution of the forward problem. We refer to [31, 32] for some other existence
and uniqueness results for the forward problem for equation (2.10). Unlike
Theorem 2.1, the positivity of the function $u$ was not discussed in these
references. On the other hand, we need this property of $u$ for our numerical
method.

Theorem 2.1. _Assume that ( 2.1), (2.2), (2.6)-(2.8) and (2.11)-(2.16) hold.
Then_ _there exists unique solution_ $u\left(\mathbf{x},\alpha\right)\in
C^{1}\left(\overline{D^{n+1}}\times\left[-d,d\right]\right)$_of equation (
2.10) with the initial condition (2.18). In addition, the following inequality
holds:_

(2.20) $u\left(\mathbf{x},\alpha\right)\geq m>0\text{ for
}\left(\mathbf{x},\alpha\right)\in\overline{\Omega}\times\left[-d,d\right],$
(2.21)
$m=\min_{\left(\mathbf{x},\alpha\right)\in\partial_{1}\Omega\times\left[-d.d\right]}\left[\mathop{\displaystyle\int}\limits_{L(\mathbf{x},\mathbf{x}_{\alpha})}f\left(\mathbf{x}\left(s\right)-\mathbf{x}_{\alpha}\right)ds\right]>0.$

_In addition, there exists a sufficiently large number_ $X>A$_such that_

(2.22) $u\left(\mathbf{x},\alpha\right)=0\text{ for
}\left|x_{1}\right|,...,\left|x_{n}\right|>X,\alpha\in\left(-d,d\right).$

Proof. Equation (2.10) can be rewritten as

(2.23)
$D_{\nu}u(\mathbf{x},\alpha)+a(\mathbf{x})u(\mathbf{x},\alpha)=\mu_{s}(\mathbf{x})\mathop{\displaystyle\int}\limits_{\Gamma_{d}}K(\mathbf{x},\alpha,\beta)u(\mathbf{x},\beta)d\beta+f\left(\mathbf{x}-\mathbf{x}_{\alpha}\right),$

where $D_{\nu}$ is the directional derivative of $u$ in the direction of the
vector $\nu.$ The first line of (2.23) can be treated as the first order
linear ordinary differential operator along the line segment
$L\left(\mathbf{x},\mathbf{x}_{\alpha}\right)$. Denote

(2.24)
$p\left(\mathbf{x},\alpha\right)=\mathop{\displaystyle\int}\limits_{L(\mathbf{x,x}_{\alpha})}a(\mathbf{x}\left(s\right))ds.$

Obviously,
$L(\mathbf{x},\mathbf{x}_{\alpha})=\left\\{\mathbf{z}\left(t\right)=\left(\alpha+t\left(x_{1}-\alpha\right),tx_{2},...,tx_{n},ty\right),t\in\left(0,1\right)\right\\}.$
Hence, $ds=\left|\mathbf{x}-\mathbf{x}_{\alpha}\right|dt$ in (2.24). We obtain

(2.25)
$p\left(\mathbf{x},\alpha\right)=\left|\mathbf{x}-\mathbf{x}_{\alpha}\right|\mathop{\displaystyle\int}\limits_{0}^{1}a\left(\alpha+t\left(x_{1}-\alpha\right),tx_{2},...,tx_{n},ty\right)dt.$

Since by (2.13) and (2.15) $a\in C^{1}\left(\mathbb{R}^{n+1}\right),$ then
(2.16), (2.25) and elementary calculations imply that

(2.26) $D_{\nu}p\left(\mathbf{x},\alpha\right)=a\left(\mathbf{x}\right),$

where $D_{\nu}$ is the operator of the directional derivative in the direction
of the vector $\nu\left(\mathbf{x},\alpha\right).$ Consider now the
integration factor $c(\mathbf{x},\alpha),$

(2.27) $c(\mathbf{x},\alpha)=e^{p\left(\mathbf{x},\alpha\right)}.$

Then (2.26) and (2.27) imply

(2.28)
$D_{\nu}c(\mathbf{x},\alpha)=a\left(\mathbf{x}\right)c(\mathbf{x},\alpha).$

Multiply both sides of (2.23) by $c(\mathbf{x},\alpha)$. We obtain

(2.29)
$c(\mathbf{x,}\alpha)D_{\nu}u(\mathbf{x},\alpha)+c(\mathbf{x},\alpha)a(\mathbf{x})u(\mathbf{x},\alpha)=$
$=c(\mathbf{x,}\alpha)\mu_{s}(\mathbf{x})\mathop{\displaystyle\int}\limits_{\Gamma_{d}}K(\mathbf{x},\alpha,\beta)u(\mathbf{x},\beta)d\beta+c(\mathbf{x,}\alpha)f\left(\mathbf{x}-\mathbf{x}_{\alpha}\right).$

Using (2.28), we obtain

(2.30)
$cD_{\nu}u+cau=D_{\nu}\left(cu\right)-uD_{\nu}c+cau=D_{\nu}\left(cu\right)-cau+cau=D_{\nu}\left(cu\right).$

Since by (2.7)-(2.9), (2.13) and (2.14) $c(\mathbf{x},\alpha)=1$ for all
points $\mathbf{x}$ where $f\left(\mathbf{x}-\mathbf{x}_{\alpha}\right)\neq
0,$ then
$c(\mathbf{x},\alpha)f\left(\mathbf{x}-\mathbf{x}_{\alpha}\right)=f\left(\mathbf{x}-\mathbf{x}_{\alpha}\right).$
Hence, (2.6)-(2.9), (2.14) and (2.30) imply that problem equation (2.29) is
equivalent with

(2.31)
$D_{\nu}(\left(cu\right)(\mathbf{x,}\alpha))=c(\mathbf{x,}\alpha)\mu_{s}(\mathbf{x})\mathop{\displaystyle\int}\limits_{\Gamma_{d}}K(\mathbf{x},\alpha,\beta)u(\mathbf{x},\beta)d\beta+f\left(\mathbf{x}-\mathbf{x}_{\alpha}\right).$

Integrating (2.31) over the line $L(\mathbf{x},\mathbf{x}_{\alpha})$ and using
the initial condition (2.18), we obtain

$u(\mathbf{x,}\alpha)=\frac{1}{c(\mathbf{x,}\alpha)}\mathop{\displaystyle\int}\limits_{L(\mathbf{x},\mathbf{x}_{\alpha})}c(\mathbf{x}\left(s\right),\alpha)\mu_{s}(\mathbf{x}\left(s\right))\left(\mathop{\displaystyle\int}\limits_{\Gamma_{d}}K(\mathbf{x}\left(s\right),\alpha,\beta)u(\mathbf{x}\left(s\right),\beta)d\beta\right)ds+$
(2.32)
$+\frac{1}{c(\mathbf{x,}\alpha)}\mathop{\displaystyle\int}\limits_{L(\mathbf{x},\mathbf{x}_{\alpha})}f\left(\mathbf{x}\left(s\right)-\mathbf{x}_{\alpha}\right)ds.$

Assume now that in (2.32)
$\mathbf{x}\in\left\\{y\in\left(0,a\right)\right\\}.$ Hence, by (2.2)
$\mathbf{x\notin}\overline{\Omega}.$ Hence, (2.14), (2.24), (2.27) and (2.32)
imply

(2.33)
$u(\mathbf{x,}\alpha)=\mathop{\displaystyle\int}\limits_{L(\mathbf{x},\mathbf{x}_{\alpha})}f\left(\mathbf{x}\left(s\right)-\mathbf{x}_{\alpha}\right)ds,\text{
}\mathbf{x}\in\left\\{y\in\left(0,a\right)\right\\}.$

This and (2.3) imply

(2.34)
$u(\mathbf{x,}\alpha)=u_{0}(\mathbf{x,}\alpha)=\mathop{\displaystyle\int}\limits_{L(\mathbf{x},\mathbf{x}_{\alpha})}f\left(\mathbf{x}\left(s\right)-\mathbf{x}_{\alpha}\right)ds,\text{
}(\mathbf{x,}\alpha)\in\partial_{1}\Omega\times\left(-d,d\right),$ (2.35)
$u_{0}(\mathbf{x,}\alpha)\geq m,$

where the number $m$ is defined in (2.21). In addition, (2.7), (2.14) and
(2.32) imply that (2.22) holds for a sufficiently large number $X>A.$ It
follows from the above construction that problem (2.32)- (2.34) is equivalent
with problem (2.18), (2.29).

Let now $\mathbf{x}\in D_{a}^{n+1},$ where the set $D_{a}^{n+1}$ is defined in
(2.17). Hence, similarly with (2.24) and (2.25), for any appropriate function
$\varphi\left(\mathbf{z}\right)$

$\mathop{\displaystyle\int}\limits_{L(\mathbf{x},\mathbf{x}_{\alpha})}\varphi\left(\mathbf{z}\left(s\right)\right)ds=\left|\mathbf{x}-\mathbf{x}_{\alpha}\right|\mathop{\displaystyle\int}\limits_{0}^{1}\varphi\left(\alpha+t\left(x_{1}-\alpha\right),tx_{2},...,tx_{n},ty\right)dt=$
$=\frac{\left|\mathbf{x}-\mathbf{x}_{\alpha}\right|}{y}\mathop{\displaystyle\int}\limits_{0}^{y}\varphi\left(\alpha+\frac{\left(x_{1}-\alpha\right)}{y}z,\frac{x_{2}}{y}z,...,\frac{x_{n}}{y}z,z\right)dz$
$=\frac{\left|\mathbf{x}-\mathbf{x}_{\alpha}\right|}{y}\mathop{\displaystyle\int}\limits_{0}^{a}\varphi\left(\alpha+\frac{\left(x_{1}-\alpha\right)}{y}z,\frac{x_{2}}{y}z,...,\frac{x_{n}}{y}z,z\right)dz$
$+\frac{\left|\mathbf{x}-\mathbf{x}_{\alpha}\right|}{y}\mathop{\displaystyle\int}\limits_{a}^{y}\varphi\left(\alpha+\frac{\left(x_{1}-\alpha\right)}{y}z,\frac{x_{2}}{y}z,...,\frac{x_{n}}{y}z,z\right)dz.$

Hence, (2.32) and (2.34) imply

$u(x_{1},x_{2},...,y\mathbf{,}\alpha)=$ (2.36)
$=\frac{\left|\mathbf{x}-\mathbf{x}_{\alpha}\right|}{yc(\mathbf{x})}\mathop{\displaystyle\int}\limits_{a}^{y}\left(c\mu_{s}\right)(\mathbf{x}\left(z\right),\alpha)\left(\mathop{\displaystyle\int}\limits_{\Gamma_{d}}K(\mathbf{x}\left(z\right),\alpha,\beta)u(\mathbf{x}\left(z\right),\beta)d\beta\right)dz+u_{0}(\mathbf{x,}\alpha),$
(2.37)
$\mathbf{x}\left(z\right)=\left(\alpha+\frac{\left(x_{1}-\alpha\right)}{y}z,\frac{x_{2}}{y}z,...,\frac{x_{n}}{y}z,z\right),\mathbf{x}\in
D_{a}^{n+1},\alpha\in\left(-d,d\right).$

Thus, by (2.22), we have obtained the $\alpha-$dependent family of integral
equations (2.22), (2.36), (2.37) of the Volterra type in the bounded domain

$D_{a,b,X}^{n+1}=\left\\{\mathbf{x}=\left(x_{1},...,x_{2},y\right):\left|x_{1}\right|,\left|x_{2}\right|,...\left|x_{n}\right|<X,y\in\left(a,b\right)\right\\}.$

Furthermore, any solution of (2.36), (2.37) satisfies (2.22). Since problem
(2.33), (2.36), (2.37) is equivalent with equation (2.29) with the initial
condition (2.18), then it is sufficient to solve problem (2.33), (2.36),
(2.37). It follows from the well known classical results about Volterra
equations that there exists unique function
$u\left(\mathbf{x},\alpha\right)\in
C\left(\overline{D_{a,b,X}^{n+1}}\times\left[-d,d\right]\right)$ satisfying
(2.36), (2.37), and this function can be obtained via the following iterative
process:

$u_{n}\left(\mathbf{x},\alpha\right)=$ (2.38)
$=\frac{\left|\mathbf{x}-\mathbf{x}_{\alpha}\right|}{yc(\mathbf{x,}\alpha)}\mathop{\displaystyle\int}\limits_{a}^{y}\left(c\mu_{s}\right)(\mathbf{x}\left(z\right),\alpha)\left(\mathop{\displaystyle\int}\limits_{\Gamma_{d}}K(\mathbf{x}\left(z\right),\alpha,\beta)u_{n-1}(\mathbf{x}\left(z\right),\beta)d\beta\right)dz+u_{0}(\mathbf{x,}\alpha),$

where $n=1,2,...$ It follows from (2.38) that for $n=0,1...$

(2.39)
$\left|u_{n}\left(\mathbf{x},\alpha\right)\right|\leq\mathop{\displaystyle\sum}\limits_{k=0}^{n}\frac{\left(M_{1}\left(y-a\right)\right)^{k}}{k!},\text{
}\left(\mathbf{x},\alpha\right)\in\overline{D_{a,b,X}^{n+1}}\times\left[-d,d\right],$

where the number
$M_{1}=M_{1}\left(a,b,d,X,\left\|a\left(x\right)\right\|_{C\left(\overline{\Omega}\right)},\left\|K\right\|_{C\left(\overline{\Omega}\times\left[-d,d\right]^{2}\right)}\right)>0$
depends only on listed parameters. Next, (2.11), (2.34), (2.35), (2.38) and
(2.39) imply that estimate (2.34) holds. Finally, it follows from (2.12),
(2.15) and (2.38) that functions $u_{n}\left(\mathbf{x},\alpha\right)$ can be
differentiated once with respect to $x_{1},...,x_{n},y,\alpha$ and, similarly
with (2.39), the following estimates hold for

(2.40)
$\left|\nabla_{\mathbf{x}}u_{n}\left(\mathbf{x},\alpha\right)\right|,\left|\partial_{\alpha}u_{n}\left(\mathbf{x},\alpha\right)\right|\leq\mathop{\displaystyle\sum}\limits_{k=0}^{n}\frac{\left(M_{2}\left(y-a\right)\right)^{k}}{k!}$

for
$\left(\mathbf{x},\alpha\right)\in\overline{D_{a,b,X}^{n+1}}\times\left[-d,d\right]$
and $n=0,1...$ Here the number

$M_{2}=M_{2}\left(a,b,d,X,\left\|a\left(x\right)\right\|_{C\left(\overline{\Omega}\right)},\left\|K\right\|_{C\left(\overline{\Omega}\times\left[-d,d\right]^{2}\right)},\left\|f\right\|_{C^{1}\left(\left|\mathbf{x}\right|\leq\varepsilon\right)}\right)>0$
depends only on listed parameters. Estimates (2.39) and (2.40) together with
(2.22) imply that we have found the unique function
$u\left(\mathbf{x},\alpha\right)\in
C^{1}\left(\overline{D^{n+1}}\times\left[-d,d\right]\right)$ satisfying
equation (2.29) with the initial condition (2.18). Furthermore, it follows
from (2.11), (2.35), (2.38) and (2.39) that (2.20) holds. $\ \square$

##  3\. Transformation

###  3.1. An integral differential equation without the unknown coefficient
$a\left(\mathbf{x}\right)$

The first step of the convexification is a transformation of the above CIP to
a boundary value problem for a certain PDE of the second order, in which the
unknown coefficient $a\left(\mathbf{x}\right)$ is not involved. By (2.20) we
can introduce a new function $w\left(\mathbf{x},\alpha\right),$

(3.1) $w\left(\mathbf{x},\alpha\right)=\ln
u\left(\mathbf{x},\alpha\right),\text{
}\left(\mathbf{x},\alpha\right)\in\overline{\Omega}\times\left[-d,d\right].$

Substituting (3.1) in (2.10) and (2.19) and using (2.34), we obtain

(3.2)
$\nu(\mathbf{x,}\alpha)\cdot\nabla_{\mathbf{x}}w(\mathbf{x},\alpha)+a\left(\mathbf{x}\right)$
$=e^{-w\left(\mathbf{x},\alpha\right)}\mu_{s}(\mathbf{x})\mathop{\displaystyle\int}\limits_{\Gamma_{d}}K(\mathbf{x},\alpha,\beta)e^{w\left(\mathbf{x},\beta\right)}d\beta,\mathbf{x}\in\Omega,\alpha\in\left(-d,d\right),$
(3.3) $w(\mathbf{x},\alpha)\mid_{\partial\Omega}=\ln
g_{1}\left(\mathbf{x},\alpha\right),$ (3.4)
$g_{1}\left(\mathbf{x},\alpha\right)=\left\\{\begin{array}[]{c}g\left(\mathbf{x},\alpha\right),\mathbf{x}\in\partial\Omega\diagdown\partial_{1}\Omega,\alpha\in\left(-d,d\right),\\\
u_{0}\left(\mathbf{x},\alpha\right),\mathbf{x}\in\partial_{1}\Omega,\alpha\in\left(-d,d\right).\end{array}\right.$

Differentiate both sides of (3.3) with respect to $\alpha$ and use
$\partial_{\alpha}a\left(\mathbf{x}\right)\equiv 0.$ We obtain an integral
differential equation with the derivatives up to the second order,

$\nu(\mathbf{x,}\alpha)\cdot\nabla_{\mathbf{x}}w_{\alpha}(\mathbf{x},\alpha)+\partial_{\alpha}\nu(\mathbf{x,}\alpha)\cdot\nabla_{\mathbf{x}}w(\mathbf{x},\alpha)$
(3.5)
$=\mu_{s}(\mathbf{x})\frac{\partial}{\partial\alpha}\left[e^{-w\left(\mathbf{x},\alpha\right)}\int_{\Gamma_{d}}K(\mathbf{x},\alpha,\beta)e^{w\left(\mathbf{x},\beta\right)}d\beta\right],\mathbf{x}\in\Omega,\alpha\in\left(-d,d\right).$

We have
$\nu(\mathbf{x,}\alpha)=\left(\nu_{1}(\mathbf{x,}\alpha),\nu_{2}(\mathbf{x,}\alpha),...,\nu_{n+1}(\mathbf{x,}\alpha)\right),$
where by (2.16)

$\nu_{1}(\mathbf{x,}\alpha)=\frac{x_{1}-\alpha}{\sqrt{\left(x_{1}-\alpha\right)^{2}+x_{2}^{2}+...+x_{n}^{2}+y^{2}}},$
$\nu_{k}(\mathbf{x,}\alpha)=\frac{x_{k}}{\sqrt{\left(x_{1}-\alpha\right)^{2}+x_{2}^{2}+...+x_{n}^{2}+y^{2}}},\text{
}k=2,...,n,$
$\nu_{n+1}(\mathbf{x,}\alpha)=\frac{y}{\sqrt{\left(x_{1}-\alpha\right)^{2}+x_{2}^{2}+...+x_{n}^{2}+y^{2}}}.$

Hence, (2.1)-(2.6) imply that with a constant $M_{3}=M_{3}\left(A,d\right)>0$
depending only on numbers $A,d$ the following estimates are valid:

(3.6)
$\left|\frac{\partial_{\alpha}\nu_{k}(\mathbf{x,}\alpha)}{\nu_{n+1}(\mathbf{x,}\alpha)}\right|,\left|\frac{\nu_{k}(\mathbf{x,}\alpha)}{\nu_{n+1}(\mathbf{x,}\alpha)}\right|\leq\frac{M_{3}}{a},k=1,...,n,\mathbf{x\in}\overline{\Omega},\alpha\in\left[-d,d\right].$

Dividing (3.5) by $\nu_{n+1}(\mathbf{x,}\alpha),$ we obtain

$\partial_{\alpha}w_{y}(\mathbf{x,}\alpha)+\frac{\partial_{\alpha}\nu_{n+1}(\mathbf{x,}\alpha)}{\nu_{n+1}(\mathbf{x,}\alpha)}w_{y}(\mathbf{x,}\alpha)$
(3.7)
$+\mathop{\displaystyle\sum}\limits_{k=1}^{n}\frac{\nu_{k}(\mathbf{x,}\alpha)}{\nu_{n+1}(\mathbf{x,}\alpha)}\partial_{\alpha}w_{x_{k}}(\mathbf{x,}\alpha)+\mathop{\displaystyle\sum}\limits_{k=1}^{n}\frac{\partial_{\alpha}\nu_{k}(\mathbf{x,}\alpha)}{\nu_{n+1}(\mathbf{x,}\alpha)}w_{x_{k}}(\mathbf{x,}\alpha)$
$-\frac{\mu_{s}(\mathbf{x})}{\nu_{n+1}(\mathbf{x,}\alpha)}\frac{\partial}{\partial\alpha}\left[e^{-w\left(\mathbf{x},\alpha\right)}\int_{\Gamma_{d}}K(\mathbf{x},\alpha,\beta)e^{w\left(\mathbf{x},\beta\right)}d\beta\right]=0,\mathbf{x}\in\Omega,\alpha\in\left(-d,d\right),$
(3.8) $w(\mathbf{x},\alpha)\mid_{\partial\Omega}=\ln
g_{1}\left(\mathbf{x},\alpha\right).$

The new equation (3.7) does not contain the unknown coefficient
$a\left(\mathbf{x}\right).$ We need now to solve problem (3.7), (3.8).

###  3.2. An orthonormal basis in $L_{2}(-d,d)$

We now describe the orthonormal basis in $L_{2}(-d,d)$ of [23] and [26,
section 6.2.3], which was mentioned in section 1. Consider the set of linearly
independent functions $\\{\alpha^{s}e^{\alpha}\\}_{s=0}^{\infty}$. This set is
complete in the space $L_{2}(-d,d)$. Applying the Gram-Schmidt
orthonormalization procedure to this set, we obtain the orthonormal basis
$\\{\Psi_{s}\left(\alpha\right)\\}_{s=0}^{\infty}$ in $L_{2}(-d,d)$. The
function $\Psi_{s}(\alpha)$ has the form

(3.9) $\Psi_{s}(\alpha)=P_{s}(\alpha)e^{\alpha},\forall s\geq 0,$

where $P_{s}(\alpha)$ is a polynomial of the degree $s$. Even though the Gram-
Schmidt orthonormalization procedure is unstable for the infinite number of
elements, it is still stable for a non large number of elements, and this was
observed in numerical studies of, e.g. [26, Chapters 7,10,12], [27] and
references cited therein. Consider the $N\times N$ matrix
$M_{N}=\left(a_{s,k}\right)_{\left(s,k\right)=\left(0,0\right)}^{\left(N-1,N-1\right)}$,

(3.10)
$a_{s,k}=\mathop{\displaystyle\int}\limits_{-d}^{d}\Psi_{s}^{\prime}\left(\alpha\right)\Psi_{k}\left(\alpha\right)d\alpha.$

It was proven in [23], [26, section 6.2.3] that

(3.11) $a_{s,k}=\left\\{\begin{array}[]{c}1\text{ if }s=k,\\\ 0\text{ if
}s<k.\end{array}\right.$

By (3.11) $\det M_{N}=1.$ Hence, the matrix $M_{N}$ is invertible.

Remark 3.1. _The invertibility of the matrix_ $M_{N}$_is the single most
important property of the orthonormal basis_
$\\{\Psi_{s}\left(\alpha\right)\\}_{s=0}^{\infty}.$_Consider, for example any
basis of either classical orthonormal polynomials or orthonormal trigonometric
functions in_ $L_{2}(-d,d)$_. Since the first function of such a basis is a
constant, then it follows from ( 3.10) that the first raw of the analog
_$\widetilde{M}_{N}$ _of the matrix_ $M_{N}$_consists of only zeros. Hence,
the matrix_ $\widetilde{M}_{N}$_is not invertible._

###  3.3. A coupled system of nonlinear integral differential equations

Represent the function $w(\mathbf{x,}\alpha)$ as truncated Fourier series with
respect to the above basis $\left\\{\Psi_{s}(\alpha)\right\\}_{s=0}^{\infty}$
and also assume that the derivative $w_{\alpha}(\mathbf{x,}\alpha)$ can be
represented as the sum of the term-by-term derivatives of that series

(3.12)
$w(\mathbf{x,}\alpha)=\mathop{\displaystyle\sum}\limits_{s=0}^{N-1}w_{s}(\mathbf{x})\Psi_{s}(\alpha),\text{
}w_{\alpha}(\mathbf{x,}\alpha)=\mathop{\displaystyle\sum}\limits_{n=0}^{N-1}w_{s}(\mathbf{x})\Psi_{n}^{\prime}(\alpha).$

Therefore, we now need to find the vector function $W(\mathbf{x})$ of
coefficients $w_{s}(\mathbf{x}),$ where

(3.13) $W(\mathbf{x})=\left(w_{0},...,w_{N-1}\right)^{T}(\mathbf{x}).$

Substituting (3.12)__ in (3.7), we obtain

$\mathop{\displaystyle\sum}\limits_{s=0}^{N-1}\partial_{y}w_{s}(\mathbf{x})\Psi_{s}^{\prime}(\alpha)+\frac{\partial_{\alpha}\nu_{n+1}(\mathbf{x,}\alpha)}{\nu_{n+1}(\mathbf{x,}\alpha)}\mathop{\displaystyle\sum}\limits_{s=0}^{N-1}\partial_{y}w_{sy}(\mathbf{x})\Psi_{n}(\alpha)$
$+\mathop{\displaystyle\sum}\limits_{s=0}^{N-1}\mathop{\displaystyle\sum}\limits_{i=1}^{n}\frac{\nu_{i}(\mathbf{x,}\alpha)}{\nu_{n+1}(\mathbf{x,}\alpha)}\partial_{x_{i}}w_{s}(\mathbf{x})\Psi_{s}^{\prime}(\alpha)+\mathop{\displaystyle\sum}\limits_{s=0}^{N-1}\mathop{\displaystyle\sum}\limits_{i=1}^{n}\frac{\partial_{\alpha}\nu_{i}(\mathbf{x,}\alpha)}{\nu_{n+1}(\mathbf{x,}\alpha)}\partial_{x_{i}}w_{s}(\mathbf{x})\Psi_{s}(\alpha)$
(3.14) $-\frac{\mu_{s}(\mathbf{x})}{\nu_{n+1}(\mathbf{x,}\alpha)}\times$
$\times\frac{\partial}{\partial\alpha}\left[\exp\left(-\mathop{\displaystyle\sum}\limits_{s=0}^{N-1}w_{s}(\mathbf{x})\Psi_{s}(\alpha)\right)\mathop{\displaystyle\int}\limits_{\Gamma_{d}}K(\mathbf{x},\alpha,\beta)\exp\left(\mathop{\displaystyle\sum}\limits_{n=0}^{N-1}w_{s}(\mathbf{x})\Psi_{s}(\beta)\right)d\beta\right]$
$=0,\mathbf{x}\in\Omega,\alpha\in\left(-d,d\right).$

Using (3.8), (3.12) and (3.13), we obtain the boundary condition for the
vector function $W\left(\mathbf{x}\right),$

(3.15)
$W\left(\mathbf{x}\right)\mid_{\partial\Omega}=P\left(\mathbf{x}\right)=\left(p_{0},...,p_{N-1}\right)^{T}(\mathbf{x}),$
(3.16)
$p_{s}(\mathbf{x})=\mathop{\displaystyle\int}\limits_{-d}^{d}\ln\left[g_{1}\left(\mathbf{x},\alpha\right)\right]\Psi_{s}\left(\alpha\right)d\alpha,n=0,...,N-1.$

Multiply sequentially equation (3.14) by functions
$\Psi_{k}\left(\alpha\right),k=0,...,N-1$ and integrate with respect to
$\alpha\in\left(-d,d\right).$ We obtain

(3.17)
$\left(M_{N}+A_{n+1}\left(\mathbf{x}\right)\right)W_{y}\left(\mathbf{x}\right)+\mathop{\displaystyle\sum}\limits_{i=1}^{n}A_{i}\left(\mathbf{x}\right)W_{x_{i}}\left(\mathbf{x}\right)+F\left(W\left(\mathbf{x}\right),\mathbf{x}\right)=0,\text{
}\mathbf{x}\in\Omega,$

where $A_{n+1}\in C_{N^{2}}\left(\overline{\Omega}\right)$ and $A_{i}\in
C_{N^{2}}\left(\overline{\Omega}\right),i=1,...,n$ are $N\times N$ matrices,
the $N-$D vector function

(3.18) $F\left(s,x\right)\in C^{2}\left(\mathbb{R}^{N+n+1}\right)$

is generated by the operator in the fourth line of (3.14). Here, (3.18)
follows from (2.12), (2.13), (2.15), (3.5) and (3.14). Obviously the vector
function $F\left(W\left(\mathbf{x}\right),\mathbf{x}\right)$ is nonlinear with
respect to $W\left(\mathbf{x}\right).$ By (3.6) and (LABEL:3.7)

(3.19)
$\left\|A_{i}\left(\mathbf{x}\right)\right\|_{C_{N^{2}}\left(\overline{\Omega}\right)}\leq\frac{C_{1}}{a},$

where the number $C_{1}=C_{1}\left(A,d,N,b\right)>0$ depends only on listed
parameters.

It follows from (2.1) and (3.19) that there exists such a number $a_{0}>1$
that the matrix

(3.20)
$D_{N}\left(\mathbf{x}\right)=\left(M_{N}+A_{n+1}\left(\mathbf{x}\right)\right)=M_{N}\left(I+M_{N}^{-1}A_{1}\left(\mathbf{x}\right)\right)$

is invertible with the inverse matrix

(3.21) $D_{N}^{-1}\left(\mathbf{x}\right),\forall a\geq
a_{0}=a_{0}\left(A,d,N\right)>1,\text{
}\forall\mathbf{x}\in\overline{\Omega},$

where the number $a_{0}\left(A,d,N\right)$ depends only on listed parameters.
Thus, the matrix $D_{N}^{-1}\left(\mathbf{x}\right)$ exists if the domain
$\Omega$ is located sufficiently far from the line of sources $\Gamma_{d}.$

Remarks 3.2:

  1. (1)

_To work with the convexification method, we need to assume below that the
derivatives_ $W_{x_{i}}$ _in ( 3.17) are written in finite differences, unlike
the _$y-$_derivative, and the grid step size_ $h\geq h_{0}>0$_with a fixed
constant_ $h_{0}$_. The latter assumption is a practical one since one does
not allow the grid step sizes tend to zero in practical computations._

  2. (2)

_The assumptions that the sign “_ $\approx"$_is replaced with the sign “_
$="$_in ( 3.12), (3.14) and (3.17) as well as the assumption of item 1 form
our approximate mathematical model, see section 1 for a discussion of the
issue of approximate mathematical models. Our specific model is verified
computationally in the 2D case in section 6._

###  3.4. Partial finite differences

Let $m>1$ be an integer. Consider $n$ partitions of the interval
$\left(-A,A\right),$ see (2.2):

(3.22)
$-A=x_{i,0}<x_{i,1}<...<x_{i,m}=A,x_{i,j+1}-x_{i,j}=h,j=0,...,m-1,i=1,...,n.$

We assume that

(3.23) $h\geq h_{0}=const.>0.$

Define the discrete subset $\Omega^{h}\hskip 5.0pt$of the domain $\Omega$ as:

(3.24)
$\Omega_{1}^{h}=\left\\{x_{i,j}\right\\}_{\left(i,j\right)=\left(1,0\right)}^{\left(i,j\right)=\left(n,m\right)},\text{
}$ (3.25)
$\Omega^{h}=\Omega_{1}^{h}\times\left(a,b\right)=\left\\{\left(x_{i,j},y\right):\left\\{x_{i,j}\right\\}_{\left(i,j\right)=\left(1,0\right)}^{\left(i,j\right)=\left(n,m\right)}\in\Omega_{1}^{h},y\in\left(a,b\right)\right\\}.$

We denote
$\mathbf{x}^{h}=\left\\{\left(x_{i,j},y\right):x_{i,j}\in\Omega_{1}^{h},y\in\left(a,b\right)\right\\}.$
By (2.3)-(2.5) and (3.24), (3.25) the boundary $\partial\Omega^{h}$ of the
domain $\Omega^{h}$ is:

$\partial\Omega^{h}=\partial_{1}\Omega^{h}\cup\partial_{2}\Omega^{h}\cup\partial_{3}\Omega^{h},$
$\partial_{1}\Omega^{h}=\Omega_{1}^{h}\times\left\\{y=a\right\\},\text{
}\partial_{2}\Omega^{h}=\Omega_{1}^{h}\times\left\\{y=b\right\\},$
$\partial_{3}\Omega^{h}=\left\\{\left(x_{i,0},y\right),\left(x_{i,m},y\right):y\in\left(a,b\right),i=1,...,n\right\\}.$

Let the vector function $Q(\mathbf{x})\in C_{N}^{1}(\overline{\Omega})$.
Denote

$Q_{ij}^{h}\left(y\right)=Q(x_{i,j},y),\hskip
5.0pti=1,...,n;j=0,...,m;y\in\left(a,b\right),$
$Q^{h}(\mathbf{x}^{h})=\left\\{Q^{h}\left(x_{1,j},x_{2,j},...,x_{n,j},y\right)\right\\},y\in\left(a,b\right).$

Thus, $Q^{h}(\mathbf{x}^{h})$ is an $N-D$ vector function of discrete
variables $x_{i,j}\in\Omega^{h}$ and continuous variable
$y\in\left(a,b\right).$ Note that the boundary terms at
$\partial_{3}\Omega^{h}$ of this vector function, which correspond to
$Q(\mathbf{x})\mid_{\partial_{3}\Omega^{h}},$are
$\left\\{Q_{i,0}^{h}\left(y\right)\right\\}\cup\left\\{Q_{i,m}^{h}\left(y\right)\right\\},i=1,...,n.$
For two vector functions
$Q^{\left(1\right)}(\mathbf{x})=\left(Q_{1}^{\left(1\right)}(\mathbf{x}),...,Q_{N}^{\left(1\right)}(\mathbf{x})\right)^{T}$
and
$Q^{\left(2\right)}(\mathbf{x})=\left(Q_{1}^{\left(2\right)}(\mathbf{x}),...,Q_{N}^{\left(2\right)}(\mathbf{x})\right)^{T}$
their scalar product $Q^{\left(1\right)}(\mathbf{x})\cdot
Q^{\left(2\right)}(\mathbf{x})$ is defined as the scalar product in
$\mathbb{R}^{N}$ and
$\left(Q^{\left(1\right)}(\mathbf{x})\right)^{2}=Q^{\left(1\right)}(\mathbf{x})\cdot
Q^{\left(1\right)}(\mathbf{x}).$ Respectively

(3.26) $Q^{\left(1\right)h}(\mathbf{x}^{h})\cdot
Q^{\left(2\right)h}(\mathbf{x}^{h})=\mathop{\displaystyle\sum}\limits_{k=1}^{N}\mathop{\displaystyle\sum}\limits_{\left(i,j\right)=\left(1,1\right)}^{\left(i,j\right)=\left(1,m-1\right)}Q_{k}^{\left(1\right)h}(x_{i,j},y)Q_{k}^{\left(2\right)h}(x_{i,j},y),$
(3.27) $\left(Q^{h}(\mathbf{x}^{h})\right)^{2}=Q^{h}(\mathbf{x}^{h})\cdot
Q^{h}(\mathbf{x}^{h}),\text{
}\left|Q^{h}(\mathbf{x}^{h})\right|=\sqrt{Q^{h}(\mathbf{x}^{h})\cdot
Q^{h}(\mathbf{x}^{h})}.$

We will use formulas (3.26), (3.27) everywhere below without further
mentioning. We exclude $j=0$ and $j=m$ here since we work below with finite
difference derivatives as defined in the next paragraph.

We define finite difference derivatives of $Q^{h}(\mathbf{x}^{h})$ with
respect to $x_{1},...,x_{n}$ only at interior points of the domain
$\Omega^{h}$ with as:

(3.28)
$\partial_{x_{k}}Q^{h}(\mathbf{x}^{h})=Q_{x_{k}}^{h}(\mathbf{x}^{h})=\left\\{Q^{h}\left(x_{1j},x_{2j},...,x_{nj},y\right)\right\\}_{x_{k}}=$
$=\left\\{\frac{Q^{h}\left(x_{1j},...,x_{k-1,j,}x_{k,j+1},x_{k+1,j}...,x_{n,j},y\right)-Q^{h}\left(x_{1j},...,x_{k-1,j,}x_{k,j-1},x_{k+1,j}...,x_{n,j},y\right)}{2h}\right\\}.$

The second line of (3.28) should be adjusted in an obvious fashion for $k=1$
and $k=n.$ Also, in that line $j=1,...,m-1$. We now define discrete analogs of
spaces $C_{N}\left(\overline{\Omega}\right)$, $H_{N}^{1}\left(\Omega\right)$
and $L_{N}^{2}\left(\Omega\right)$ as:

$C_{N}^{h}\left(\overline{\Omega^{h}}\right)=\left\\{Q^{h}(\mathbf{x}^{h}):\left\|Q^{h}(\mathbf{x}^{h})\right\|_{C^{h}\left(\overline{\Omega^{h}}\right)}=\max_{y\in\left[a,b\right]}\left(\max_{i=\overline{1,n},j=\overline{0,m}}\left|Q_{ij}^{h}\left(y\right)\right|\right)<\infty\right\\},$
$H_{N}^{1,h}\left(\Omega^{h}\right)=\left\\{\begin{array}[]{c}Q^{h}(\mathbf{x}^{h}):\left\|Q^{h}(\mathbf{x}^{h})\right\|_{H_{N}^{1,h}\left(\Omega^{h}\right)}^{2}=\\\
\mathop{\displaystyle\sum}\limits_{\left(i,j.k\right)=\left(1,1,1\right)}^{\left(i,j,k\right)=\left(n,m-1,n\right)}\mathop{\displaystyle\int}\limits_{a}^{b}\left[\left(Q_{ij}^{h}\left(y\right)\right)^{2}+\left(Q_{x_{k}}^{h}(\mathbf{x}^{h})\right)^{2}+\left(\partial_{y}Q_{ij}^{h}\left(y\right)\right)^{2}\right]dy<\infty\end{array}\right\\},$
$H_{N,0}^{1,h}\left(\Omega^{h}\right)=\left\\{Q^{h}(\mathbf{x}^{h})\in
H_{N}^{1,h}\left(\Omega^{h}\right):Q^{h}(\mathbf{x}^{h})\mid_{\partial\Omega^{h}}=0\right\\},$
$L_{N}^{2,h}\left(\Omega^{h}\right)=\left\\{\begin{array}[]{c}Q^{h}(\mathbf{x}^{h}):\left\|Q^{h}(\mathbf{x}^{h})\right\|_{L_{N}^{2,h}\left(\Omega^{h}\right)}^{2}=\\\
=\mathop{\displaystyle\sum}\limits_{\left(i,j\right)=\left(1,1\right)}^{\left(i,j\right)=\left(n,m-1\right)}\left\\{\mathop{\displaystyle\int}\limits_{a}^{b}\left(Q_{ij}^{h}\left(y\right)\right)^{2}dy\right\\}<\infty\end{array}\right\\}.$

By embedding theorem $H_{N}^{1,h}\left(\Omega^{h}\right)\subset
C_{N}^{h}\left(\overline{\Omega^{h}}\right)$ and

(3.29)
$\left\|Q^{h}(\mathbf{x}^{h})\right\|_{C_{N}^{h}\left(\overline{\Omega^{h}}\right)}\leq
C_{2}\left\|Q^{h}(\mathbf{x}^{h})\right\|_{H_{N}^{1,h}\left(\overline{\Omega^{h}}\right)},\forall
Q^{h}\in H_{N}^{1,h}\left(\Omega^{h}\right),$

where the number $C_{2}=C_{2}\left(h_{0},A,a,b,\Omega\right)>0$ depends only
on listed parameters, and the number $h_{0}$ is defined in (3.23).

Remark 3.3. _Since we work with finite difference derivatives ( 3.28) only at
interior points of the discrete domain _$\Omega^{h}$_, then in any
differential operator below we use_
$W^{h}\left(\mathbf{x}^{h}\right)=\left\\{W_{ij}^{h}\left(y\right)\right\\}_{\left(i,j\right)=\left(1,1\right)}^{\left(i,j\right)=\left(n,m-1\right)},y\in\left(a,b\right),$_i.e.
boundary points_ $x_{i,0}$_and_ $x_{i,m}$_are involved only in finite
difference derivatives_ $W_{i}^{h}\left(\mathbf{x}^{h}\right)$_at_
$x_{i,1}$_and_ $x_{i,m-1}$_. Points_ $x_{i,0}$_and_ $x_{i,m}$_are not included
in ( 3.26) for the same reason._

Using (3.22)-(3.28), we now rewrite problem (3.15)-(3.17) in the form of
finite differences with respect to $x_{1},...,x_{n}$ as:

(3.30)
$D_{N}^{h}\left(\mathbf{x}^{h}\right)W_{y}^{h}\left(\mathbf{x}^{h}\right)+\mathop{\displaystyle\sum}\limits_{i=1}^{n}A_{i}^{h}\left(\mathbf{x}^{h}\right)W_{x_{i}}^{h}\left(\mathbf{x}^{h}\right)+F^{h}\left(W^{h}\left(\mathbf{x}^{h}\right),\mathbf{x}^{h}\right)=0,\mathbf{x}^{h}\in\Omega^{h},$
(3.31)
$W^{h}\left(\mathbf{x}^{h}\right)\mid_{\partial\Omega^{h}}=P^{h}\left(\mathbf{x}^{h}\right).$

In (3.30) $N\times N$ matrices $A_{i}^{h}\left(\mathbf{x}^{h}\right)\in
C_{N}\left(\overline{\Omega^{h}}\right)$, the matrix
$D_{N}^{h}\left(\mathbf{x}^{h}\right)\in C\left(\overline{\Omega^{h}}\right)$
is the discrete analog of the matrix $D_{N}\left(\mathbf{x}\right)$ defined in
(3.20) By (3.19) and (3.20)

(3.32)
$\left\|D_{N}^{h}\right\|_{C_{N^{2}}^{h}\left(\overline{\Omega^{h}}\right)}\leq
C_{3},\forall a\geq a_{0}=a_{0}\left(A,d,N\right)>1.$

The vector function
$F^{h}\left(W^{h}\left(\mathbf{x}^{h}\right),\mathbf{x}^{h}\right)$ the
discrete analog of the vector function $F$ in (3.17), and by (3.18)
$F^{h}\left(W^{h}\left(\mathbf{x}^{h}\right),\mathbf{x}^{h}\right)$ it is
twice continuously differentiable with respect to components of
$W^{h}\left(\mathbf{x}^{h}\right)$. By (3.21) there exists the inverse matrix
$\left(D_{N}^{h}\right)^{-1}\left(\mathbf{x}^{h}\right)$ and

(3.33)
$\left\|\left(D_{N}^{h}\right)^{-1}\right\|_{C_{N^{2}}^{h}\left(\overline{\Omega^{h}}\right)}\leq
C_{3},\forall a\geq a_{0}=a_{0}\left(A,d,N\right)>1.$

Furthermore, (3.19) and (3.33) imply for $i=1,...,n$

(3.34)
$\left\|\left(\left(D_{N}^{h}\right)^{-1}A_{i}^{h}\right)\left(\mathbf{x}^{h}\right)\right\|_{C_{N^{2}}^{h}\left(\overline{\Omega^{h}}\right)}\leq
C_{3},\forall a\geq a_{0}=a_{0}\left(A,d,N\right)>1,$

Here and everywhere below
$C_{3}=C_{3}\left(A,d,N,a,b,h_{0},\left\|K\right\|_{C^{1}\left(\overline{\Omega}\times\left[-d,d\right]^{2}\right)}\right)>0$
denotes different constants depending only on listed parameters.

The following formulas are discrete analogs of (3.12) and (3.13):

(3.35)
$w^{h}(\mathbf{x}^{h}\mathbf{,}\alpha)=\mathop{\displaystyle\sum}\limits_{n=0}^{N-1}w_{n}^{h}(\mathbf{x}^{h})\Psi_{n}(\alpha),\text{
}\partial_{\alpha}w^{h}(\mathbf{x}^{h}\mathbf{,}\alpha)=\mathop{\displaystyle\sum}\limits_{n=0}^{N-1}w_{n}^{h}(\mathbf{x}^{h})\Psi_{n}^{\prime}(\alpha),$
(3.36)
$W^{h}(\mathbf{x}^{h})=\left(w_{0}^{h},...,w_{N-1}^{h}\right)^{T}(\mathbf{x}^{h}).$

Suppose that we have found the vector function $W^{h}(\mathbf{x}^{h})$ in
(3.36). Then, to find the discrete analog $a^{h}\left(\mathbf{x}^{h}\right)$
of the unknown coefficient $a\left(\mathbf{x}\right),$ we use (3.2) and (3.35)
as:

(3.37)
$a^{h}\left(\mathbf{x}^{h}\right)=-\frac{1}{2d}\mathop{\displaystyle\int}\limits_{-d}^{d}\nu(\mathbf{x}^{h}\mathbf{,}\alpha)\cdot\nabla_{\mathbf{x}^{h}}w^{h}(\mathbf{x}^{h},\alpha)d\alpha+$
$+\frac{1}{2d}\mathop{\displaystyle\int}\limits_{-d}^{d}\left(e^{-w^{h}\left(\mathbf{x}^{h},\alpha\right)}\mu_{s}(\mathbf{x}^{h})\mathop{\displaystyle\int}\limits_{-d}^{d}K(\mathbf{x}^{h},\alpha,\beta)e^{w^{h}\left(\mathbf{x}^{h},\beta\right)}d\beta\right)d\alpha,\mathbf{x}^{h}\in\Omega^{h}.$

##  4\. Convexification Method for Problem (3.30), (3.31)

Let $R>0$ be an arbitrary number. Define the set $B\left(R,P^{h}\right)$ as

(4.1) $B\left(R,P^{h}\right)=\left\\{W^{h}\in
H_{N}^{1,h}\left(\Omega^{h}\right):W^{h}\left(\mathbf{x}^{h}\right)\mid_{\partial\Omega^{h}}=P^{h}\left(\mathbf{x}^{h}\right),\left\|W^{h}\right\|_{H_{N}^{1,h}\left(\Omega^{h}\right)}<R\right\\},$

where $P^{h}\left(\mathbf{x}^{h}\right)$ is the boundary condition in (3.31).
Consider matrices $G_{i}^{h}\left(\mathbf{x}^{h}\right)$ and the vector
function
$\Phi^{h}\left(W^{h}\left(\mathbf{x}^{h}\right),\mathbf{x}^{h}\right),$

(4.2)
$G_{i}^{h}\left(\mathbf{x}^{h}\right)=\left(\left(D_{N}^{h}\right)^{-1}A_{i}^{h}\right)\left(\mathbf{x}^{h}\right),\text{
}$ (4.3)
$\Phi^{h}\left(W^{h}\left(\mathbf{x}^{h}\right),\mathbf{x}^{h}\right)=\left(D_{N}^{h}\left(\mathbf{x}^{h}\right)\right)^{-1}F^{h}\left(W^{h}\left(\mathbf{x}^{h}\right),\mathbf{x}^{h}\right).$

By (3.18), (3.33), (3.34), (4.1)-(4.3) and the multidimensional analog of
Taylor formula [34]

(4.4)
$\left\|G_{i}^{h}\left(\mathbf{x}^{h}\right)\right\|_{C_{N^{2}}^{h}\left(\overline{\Omega^{h}}\right)}\leq
C_{3},\text{ }\forall a\geq a_{0}>1,\text{ }i=1,...,n,$ (4.5)
$\left\|\Phi^{h}\left(W^{h}\left(\mathbf{x}^{h}\right),\mathbf{x}^{h}\right)\right\|_{C_{N}^{h}\left(\overline{\Omega^{h}}\right)}\leq
C_{3},\text{ }\forall W^{h}\in\overline{B\left(R,P^{h}\right)},$ (4.6)
$\Phi^{h}\left(W_{2}^{h}\left(\mathbf{x}^{h}\right),\mathbf{x}^{h}\right)=\Phi^{h}\left(W_{1}^{h}\left(\mathbf{x}^{h}\right),\mathbf{x}^{h}\right)+\Phi_{1}\left(W_{1}^{h}\left(\mathbf{x}^{h}\right),\mathbf{x}^{h}\right)\left(W_{2}^{h}-W_{1}^{h}\right)\left(\mathbf{x}^{h}\right)$
$+\Phi_{2}^{h}\left(W_{1}^{h}\left(\mathbf{x}^{h}\right),W_{2}^{h}\left(\mathbf{x}^{h}\right),\mathbf{x}^{h}\right),\text{
}\forall
W_{1}^{h},W_{2}^{h}\in\overline{B\left(R,P^{h}\right)},\forall\mathbf{x}^{h}\in\overline{\Omega^{h}},$

where the vector function
$\Phi_{1}\left(W_{1}^{h}\left(\mathbf{x}^{h}\right),\mathbf{x}^{h}\right)$ is
independent on $W_{2}^{h}\left(\mathbf{x}^{h}\right),$ the vector function
$\Phi_{2}^{h}\left(W_{1}^{h}\left(\mathbf{x}^{h}\right),W_{2}^{h}\left(\mathbf{x}^{h}\right),\mathbf{x}^{h}\right)$
is nonlinear with respect to
$\left(W_{2}^{h}-W_{1}^{h}\right)\left(\mathbf{x}^{h}\right),$ both these
vector functions are continuous with respect to their variables for
$\mathbf{x}^{h}\in\Omega^{h}$ and the following estimates hold for all
$W_{1}^{h},W_{2}^{h}\in\overline{B\left(R,P^{h}\right)}$ and for all
$\mathbf{x}^{h}\in\Omega^{h}$

(4.7)
$\left|\Phi_{1}\left(W_{1}^{h}\left(\mathbf{x}^{h}\right),\mathbf{x}^{h}\right)\right|\leq
C_{3},\text{ }\forall
W_{1}^{h}\in\overline{B\left(R,P^{h}\right)},\forall\mathbf{x}^{h}\in\Omega^{h},$
(4.8)
$\left|\Phi_{2}^{h}\left(W_{1}^{h}\left(\mathbf{x}^{h}\right),W_{2}^{h}\left(\mathbf{x}^{h}\right),\mathbf{x}^{h}\right)\right|\leq
C_{3}\left(W_{2}^{h}-W_{1}^{h}\right)^{2}\left(\mathbf{x}^{h}\right),\text{
}\forall\mathbf{x}^{h}\in\Omega^{h},$

also, see (3.27).

Lemma 4.1. _Let_ $A$_be a_ $k\times k$_matrix which has the inverse_
$A^{-1}.$_Then there exists a number_ $\mu=\mu\left(A\right)>0$_such that_
$\left\|Ax\right\|^{2}\geq\mu\left\|x\right\|^{2},\forall x\in\mathbb{R}^{k},$
_where_ $\left\|\cdot\right\|$_is_ _the euclidean norm._

We omit the proof since this lemma is well known.

Corollary 4.1. _The following inequality holds:_

$\left(D_{N}^{h}\left(\mathbf{x}^{h}\right)W_{y}^{h}\left(\mathbf{x}^{h}\right)+\mathop{\displaystyle\sum}\limits_{i=1}^{n}A_{i}^{h}\left(\mathbf{x}^{h}\right)W_{x_{i}}^{h}\left(\mathbf{x}^{h}\right)\right)^{2}\geq$
$\geq
C_{3}\left(W_{y}^{h}\left(\mathbf{x}^{h}\right)+\mathop{\displaystyle\sum}\limits_{i=1}^{n}G_{i}^{h}\left(\mathbf{x}^{h}\right)W_{x_{i}}^{h}\left(\mathbf{x}^{h}\right)\right)^{2},\forall\mathbf{x}^{h}\in\Omega^{h}.$

Proof. Denote

(4.9)
$Y\left(\mathbf{x}^{h}\right)=W_{y}^{h}\left(\mathbf{x}^{h}\right)+\mathop{\displaystyle\sum}\limits_{i=1}^{n}G_{i}^{h}\left(\mathbf{x}^{h}\right)W_{x_{i}}^{h}\left(\mathbf{x}^{h}\right)+\Phi^{h}\left(W^{h}\left(\mathbf{x}^{h}\right),\mathbf{x}^{h}\right).$

We have

(4.10)
$\left(D_{N}^{h}\left(\mathbf{x}^{h}\right)W_{y}^{h}\left(\mathbf{x}^{h}\right)+\mathop{\displaystyle\sum}\limits_{i=1}^{n}A_{i}^{h}\left(\mathbf{x}^{h}\right)W_{x_{i}}^{h}\left(\mathbf{x}^{h}\right)\right)^{2}=\left(D_{N}^{h}\left(\mathbf{x}^{h}\right)Y\left(\mathbf{x}^{h}\right)\right)^{2}.$

The rest of the proof follows immediately from (4.2), (4.3) and Lemma 4.1.
$\square$

Introduce the following weighted cost functional
$J_{\lambda}\left(W^{h}\right)$:

(4.11)
$J_{\lambda}\left(W^{h}\right)=\left\|\left(D_{N}^{h}W_{y}^{h}+\mathop{\displaystyle\sum}\limits_{i=1}^{n}A_{i}^{h}W_{x_{i}}^{h}+F^{h}\left(W^{h}\left(\mathbf{x}^{h}\right),\mathbf{x}^{h}\right)\right)e^{\lambda
y}\right\|_{L_{N}^{2,h}\left(\Omega^{h}\right)}^{2}.$

Minimization Problem. _Minimize functional ( 4.11) on the set_
$\overline{B\left(R,P^{h}\right)}.$

Theorems 4.1-4.5 are our analytical results about the functional
$J_{\lambda,\gamma}\left(W^{h}\right).$

Theorem 4.1 (Carleman estimate). _Assume that the number_ $a\geq a_{0},$_as in
( 3.33). Then there exists a sufficiently large number
_$\lambda_{0}=\lambda_{0}\left(A,d,N,a,b,h_{0}\right)\geq 1$_depending only on
listed parameters such that the following Carleman estimate holds_

(4.12)
$\left\|\left(W_{y}^{h}\left(\mathbf{x}^{h}\right)+\mathop{\displaystyle\sum}\limits_{i=1}^{n}G_{i}^{h}\left(\mathbf{x}^{h}\right)W_{x_{i}}^{h}\left(\mathbf{x}^{h}\right)\right)e^{\lambda
y}\right\|_{L_{N}^{2,h}\left(\Omega^{h}\right)}^{2}\geq$
$\geq\left\|W_{y}^{h}e^{\lambda
y}\right\|_{L_{N}^{2,h}\left(\Omega^{h}\right)}^{2}+\lambda^{2}\left\|W^{h}e^{\lambda
y}\right\|_{L_{N}^{2,h}\left(\Omega^{h}\right)}^{2},\text{ }\forall W^{h}\in
H_{N,0}^{1,h}\left(\Omega^{h}\right),\forall\lambda\geq\lambda_{0}.$

__

Theorem 4.2 (the central analytical result). _The following three assertions
hold:_

1\. _The functional_ $J_{\lambda}\left(W^{h}\right)$_in ( 4.11) has the
Fréchet derivative _$J_{\lambda}^{\prime}\left(W^{h}\right)\in
H_{N,0}^{1,h}\left(\Omega^{h}\right)$_at any point_
$W^{h}\in\overline{B\left(R,P^{h}\right)}$_and for any value of the parameter_
$\lambda\geq 0.$_This function satisfies the Lipschitz condition_

(4.13)
$\left\|J_{\lambda}^{\prime}\left(W_{1}^{h}\right)-J_{\lambda}^{\prime}\left(W_{2}^{h}\right)\right\|_{H_{N}^{1,h}\left(\Omega^{h}\right)}\leq
C_{4}\left\|W_{1}^{h}-W_{2}^{h}\right\|_{H_{N}^{1,h}\left(\Omega^{h}\right)},\forall
W_{1}^{h},W_{2}^{h}\in\overline{B\left(R,P^{h}\right)}$

_for all_ $\lambda\geq 0,$ _where the number_ $C_{4}>0$_depends on the same
parameters as ones in_ $C_{3}$_as well as on_ $\lambda.$

_Assume that the number_ $a\geq a_{0},$_as in ( 3.33). Then:_

_2\. There exists a sufficiently large number_ $\lambda_{1}$__

(4.14)
$\lambda_{1}=\lambda_{1}\left(R,A,d,N,a,b,h_{0}\right)\geq\lambda_{0}\geq 1$

_depending only on listed parameters such that the functional_
$J_{\lambda}\left(W^{h}\right)$_in ( 4.11) is strictly convex on the set
_$\overline{B\left(R,P^{h}\right)},$_i.e. the following inequality holds:_

(4.15)
$J_{\lambda}\left(W_{2}^{h}\right)-J_{\lambda}\left(W_{1}^{h}\right)-J_{\lambda}^{\prime}\left(W_{1}^{h}\right)\left(W_{2}^{h}-W_{1}^{h}\right)\geq
C_{3}e^{2\lambda
a}\left\|W_{2}^{h}-W_{1}^{h}\right\|_{H_{N}^{1,h}\left(\Omega^{h}\right)}^{2},$
(4.16) $\forall\lambda\geq\lambda_{1},\forall
W_{1}^{h},W_{2}^{h}\in\overline{B\left(R,P^{h}\right)}.$

_3\. For each_ $\lambda\geq\lambda_{1}$_there exists unique minimizer_
$W_{\min,\lambda}^{h}\in\overline{B\left(R,P^{h}\right)}$_of the functional_
$J_{\lambda}\left(W^{h}\right)$_on the set_
$\overline{B\left(R,P^{h}\right)}$_. Furthermore, the following inequality
holds:_

(4.17)
$J_{\lambda}^{\prime}\left(W_{\min,\lambda}^{h}\right)\left(W^{h}-W_{\min,\lambda}^{h}\right)\geq
0,\text{ }\forall W^{h}\in\overline{B\left(R,P^{h}\right)}.$

Theorem 4.3 follows immediately from (4.15) and (4.16). This is a certain
uniqueness result for our CIP, which is obtained as a by-product. A further
discussion of the uniqueness issue is outside of the scope of this paper.

Theorem 4.3. _Assume that the number_ $a\geq a_{0},$_as in ( 3.33). Then there
exists at most one pair of functions _$\left(W^{h},a^{h}\right)\in
H_{N}^{1,h}\left(\Omega^{h}\right)\times L_{N}^{2,h}\left(\Omega^{h}\right)$
_satisfying conditions ( 3.30), (3.37)._

We now estimate the accuracy of the minimizer $W_{\min,\lambda}^{h}$ depending
on the level of the noise $\delta>0$ in the data. Following the concept of
Tikhonov for ill-posed problems [33], we assume the existence of the exact
solution

(4.18) $W^{h\ast}\in B\left(R,P^{h\ast}\right)$

of problem (3.30)-(3.31) with the exact, i.e. noiseless data $P^{h\ast},$ i.e.
for $\mathbf{x}^{h}\in\Omega^{h}$

(4.19)
$D_{N}^{h}\left(\mathbf{x}^{h}\right)W_{y}^{h\ast}\left(\mathbf{x}^{h}\right)+\mathop{\displaystyle\sum}\limits_{i=1}^{n}A_{i}^{h}\left(\mathbf{x}^{h}\right)W_{x_{i}}^{h\ast}\left(\mathbf{x}^{h}\right)+F^{h}\left(W^{h\ast}\left(\mathbf{x}^{h}\right),\mathbf{x}^{h}\right)=0,$
(4.20)
$W^{h\ast}\left(\mathbf{x}^{h}\right)\mid_{\partial\Omega^{h}}=P^{h\ast}\left(\mathbf{x}^{h}\right).$

Suppose that there exists a vector function $S^{h}\in
H_{N}^{1,h}\left(\Omega^{h}\right)$ such that

(4.21)
$S^{h}\left(\mathbf{x}^{h}\right)\mid_{\partial\Omega^{h}}=P^{h}\left(\mathbf{x}^{h}\right),\text{
}\left\|S^{h}\right\|_{H_{N}^{1,h}\left(\Omega^{h}\right)}\,<R.$

Let $S^{h\ast}\in H_{N}^{1,h}\left(\Omega^{h}\right)$ be such a vector
function that

(4.22)
$S^{h\ast}\left(\mathbf{x}^{h}\right)\mid_{\partial\Omega^{h}}=P^{h\ast}\left(\mathbf{x}^{h}\right),\text{
}\left\|S^{h\ast}\right\|_{H_{N}^{1,h}\left(\Omega^{h}\right)}\,<R.$

We assume that

(4.23)
$\left\|S^{h}-S^{h\ast}\right\|_{H_{N}^{1,h}\left(\Omega^{h}\right)}<\delta.$

Theorem 4.4. _Assume that the number_ $a\geq a_{0},$_as in ( 3.33). Suppose
that conditions (4.19)-(4.23) hold. Also, consider the number
_$\lambda_{2},$__

(4.24) $\lambda_{2}=\lambda_{1}\left(2R,A,d,N,a,b,h_{0}\right),$

_where_ $\lambda_{1}\left(R,A,d,N,a,b,h_{0}\right)$_is the number in ( 4.14).
Let _$W_{\min,\lambda_{2}}^{h}$_be the minimizer of functional ( 4.11) on the
set _$\overline{B\left(R,P^{h}\right)},$_which was found in Theorem 4.2. Let_
$\alpha\in\left(0,R\right)$ _be a number. Suppose that ( 4.18) is replaced
with_

(4.25) $W^{h\ast}\in B\left(R-\alpha,P^{h\ast}\right)\text{ and
}C_{3}\delta<\alpha.$

_Then the vector function_ $W_{\min,\lambda_{2}}^{h}$_belongs to the open set_
$B\left(R,P^{h}\right)$_and the following accuracy estimate holds:_

(4.26)
$\left\|W_{\min,\lambda_{2}}^{h}-W^{h\ast}\right\|_{H_{N}^{1,h}\left(\Omega^{h}\right)}\leq
C_{3}\delta.$

__

Consider now the gradient descent method of the minimization of functional _(
4.11) on the set _$\overline{B\left(R,P^{h}\right)}.$ Let $W_{0}^{h}\in
B\left(R/3,P^{h}\right)$ be an arbitrary point of this set. We treat it as the
starting point of the latter method. The sequence of this method is:

(4.27) $W_{n}^{h}=W_{n-1}^{h}-\gamma
J_{\lambda_{2}}^{\prime}\left(W_{n-1}^{h}\right),n=1,2,...,$

where $\gamma>0$ is a small number and $\lambda_{2}$ is the same as in (4.24).
Note that since by Theorem 4.2 functions
$J_{\lambda_{2}}^{\prime}\left(W_{n-1}^{h}\right)\in
H_{N,0}^{1,h}\left(\Omega^{h}\right)$, then all vector functions $W_{n}^{h}$
have the same boundary conditions $P^{h}.$

Theorem 4.5. _Let conditions of Theorem 4.4 hold, except that ( 4.25) is
replaced with_

(4.28) $W^{h\ast}\in B\left(\left(R-\alpha\right)/3,P^{h\ast}\right)\text{ and
}C_{3}\delta<\alpha/3.$

_Then there exists a sufficiently small number_ $\gamma>0$_and a number_
$\theta=\theta\left(\gamma\right)\in\left(0,1\right)$_such that in ( 4.27) all
functions _$W_{n}^{h}\in B\left(R,P^{h}\right)$_and the following convergence
estimates hold:_

(4.29)
$\left\|W_{n}^{h}-W_{\min,\lambda_{2}}^{h}\right\|_{H_{N}^{1,h}\left(\Omega^{h}\right)}\leq\theta^{n}\left\|W_{0}^{h}-W_{\min,\lambda_{2}}^{h}\right\|_{H_{N}^{1,h}\left(\Omega^{h}\right)},$
(4.30)
$\left\|W_{n}^{h}-W^{h\ast}\right\|_{H_{N}^{1,h}\left(\Omega^{h}\right)}\leq
C_{3}\delta+\theta^{n}\left\|W_{0}^{h}-W_{\min,\lambda_{2}}^{h}\right\|_{H_{N}^{1,h}\left(\Omega^{h}\right)},$
(4.31)
$\left\|a_{n}^{h}-a^{h\ast}\right\|_{L_{N}^{2,h}\left(\Omega^{h}\right)}\leq
C_{3}\delta+\theta^{n}\left\|W_{0}^{h}-W_{\min,\lambda_{2}}^{h}\right\|_{H_{N}^{1,h}\left(\Omega^{h}\right)},$

_where_ $a_{n}^{h}\left(\mathbf{x}^{h}\right)$_and_
$a_{n}^{h\ast}\left(\mathbf{x}^{h}\right)$_are functions which are obtained
from_ $W_{n}^{h}$_and_ $W^{h\ast}$_respectively via ( 3.35) and (3.37)._

Remarks 4.1:

  1. (1)

_Estimates ( 4.29)-(4.31) guarantee the global convergence of the gradient
descent method (4.27) since _$R>0$_is an arbitrary number and the starting
point_ $W_{0}^{h}$_is an arbitrary point of the set_
$B\left(R/3,P^{h}\right)$_, see section 1 for our definition of the global
convergence. Note that any gradient-like method converges only locally for a
non-convex functional, i.e. it needs a good first guess about the correct
solution._

  2. (2)

_Although the above results are valid only for sufficiently large values of
the parameter_ $\lambda,$_we have established computationally in section 6
that_ $\lambda=5$_works quite well. It was computationally established in a
number of works on the convexification for a variety of CIPs that values_
$\lambda\in\left[1,3\right]$_work well numerically, see, e.g.[26, Chapters
7-10], [27] and references cited therein._

  3. (3)

_An analogy to the issue of item 2 is that any asymptotic theory basically
states that if a parameter_ $X$_is sufficiently large/small, then a certain
‘nice’ formula provides a good approximation for something. In practice,
however, given a specific problem with specific values of its parameters, only
numerical experiments can clarify which exactly values of_ $X$_are so
sufficiently large/small that this ‘nice’ formula indeed provides that good
approximation._

##  5\. Proofs

We use in this section (3.26)-(3.28) and Remark 3.3 without further
mentioning.

###  5.1. Proof of Theorem 4.1

In this proof, $W^{h}\in H_{N,0}^{1,h}\left(\Omega^{h}\right)$ is an arbitrary
function. It follows from (3.23) and (3.28) that

(5.1) $\left\|W_{x_{i}}^{h}e^{\lambda
y}\right\|_{L_{N}^{2,h}\left(\Omega^{h}\right)}^{2}\leq
C_{3}\left\|W^{h}e^{\lambda
y}\right\|_{L_{N}^{2,h}\left(\Omega^{h}\right)}^{2},\text{ }\forall W^{h}\in
H_{N,0}^{1,h}\left(\Omega^{h}\right).$

By (4.4) and (5.1)

(5.2)
$\left\|\left(\mathop{\displaystyle\sum}\limits_{i=1}^{n}G_{i}^{h}\left(\mathbf{x}^{h}\right)W_{x_{i}}^{h}\left(\mathbf{x}^{h}\right)\right)e^{\lambda
y}\right\|_{L_{N}^{2,h}\left(\Omega^{h}\right)}^{2}\leq
C_{3}\left\|W^{h}e^{\lambda y}\right\|.$

Using Cauchy-Schwarz inequality and (5.2), we obtain

$\left\|\left(W_{y}^{h}\left(\mathbf{x}^{h}\right)+\mathop{\displaystyle\sum}\limits_{i=1}^{n}G_{i}^{h}\left(\mathbf{x}^{h}\right)W_{x_{i}}^{h}\left(\mathbf{x}^{h}\right)\right)e^{\lambda
y}\right\|_{L_{N}^{2,h}\left(\Omega^{h}\right)}^{2}\geq$ (5.3) $\geq
C_{3,1}\left\|W_{y}^{h}e^{\lambda
y}\right\|_{L_{N}^{2,h}\left(\Omega^{h}\right)}^{2}-C_{3,2}\left\|W^{h}e^{\lambda
y}\right\|_{L_{N}^{2,h}\left(\Omega^{h}\right)}^{2}.$

It makes sense to use two numbers $C_{3,1},C_{3,2}>0$ here, both depend on the
same parameters as ones in $C_{3}.$ Consider now the first term in the second
line of (5.3). Introduce a new vector function
$V^{h}\left(\mathbf{x}^{h}\right)=W^{h}\left(\mathbf{x}^{h}\right)e^{\lambda
y}.$ Then

$W^{h}\left(\mathbf{x}^{h}\right)=V^{h}\left(\mathbf{x}^{h}\right)e^{-\lambda
y},W_{y}^{h}\left(\mathbf{x}^{h}\right)=\left(V_{y}^{h}\left(\mathbf{x}^{h}\right)-\lambda
V^{h}\left(\mathbf{x}^{h}\right)\right)e^{-\lambda y},$
$\left(W_{y}^{h}\left(\mathbf{x}^{h}\right)\right)^{2}e^{2\lambda
y}=\left(V_{y}^{h}\left(\mathbf{x}^{h}\right)-\lambda
V^{h}\left(\mathbf{x}^{h}\right)\right)^{2}\geq-2\lambda
V_{y}^{h}\left(\mathbf{x}^{h}\right)V^{h}\left(\mathbf{x}^{h}\right)+\lambda^{2}\left(V^{h}\left(\mathbf{x}^{h}\right)\right)^{2}$
$=\partial_{y}\left(-\lambda
V^{h}\left(\mathbf{x}^{h}\right)\right)^{2}+\lambda^{2}\left(V^{h}\left(\mathbf{x}^{h}\right)\right)^{2}.$

Hence,

$\left\|W_{y}^{h}e^{\lambda
y}\right\|_{L_{N}^{2,h}\left(\Omega^{h}\right)}^{2}\geq\mathop{\displaystyle\sum}\limits_{\left(i,j\right)=\left(1,1\right)}^{\left(i,j\right)=\left(n,m-1\right)}\left[-\lambda\left(W^{h}\left(\mathbf{x}^{h}\right)\right)^{2}\mid_{y=b}+\lambda\left(W^{h}\left(\mathbf{x}^{h}\right)\right)^{2}\mid_{y=a}\right]+$
(5.4) $+\lambda^{2}\left\|W^{h}e^{\lambda
y}\right\|_{L_{N}^{2,h}\left(\Omega^{h}\right)}^{2}.$

Since $W^{h}\left(\mathbf{x}^{h}\right)\mid_{y=b}=0,$ $\forall W^{h}\in
H_{N,0}^{1,h}\left(\Omega^{h}\right),$ then (5.4) implies

(5.5) $\left\|W_{y}^{h}e^{\lambda
y}\right\|_{L_{N}^{2,h}\left(\Omega^{h}\right)}^{2}\geq\lambda^{2}\left\|W^{h}e^{\lambda
y}\right\|_{L_{N}^{2,h}\left(\Omega^{h}\right)}^{2}.$

Adding the term $\left\|W_{y}^{h}e^{\lambda
y}\right\|_{L_{N}^{2,h}\left(\Omega^{h}\right)}^{2}$ to both sides of (5.5)
and then dividing by “2”, we obtain for all $\lambda>0$

(5.6) $\left\|W_{y}^{h}e^{\lambda
y}\right\|_{L_{N}^{2,h}\left(\Omega^{h}\right)}^{2}\geq\frac{1}{2}\left\|W_{y}^{h}e^{\lambda
y}\right\|_{L_{N}^{2,h}\left(\Omega^{h}\right)}^{2}+\frac{\lambda^{2}}{2}\left\|W^{h}e^{\lambda
y}\right\|_{L_{N}^{2,h}\left(\Omega^{h}\right)}^{2}.$

Hence, taking $\lambda_{0}^{2}\geq 4C_{3,2}/C_{3,1},\lambda\geq\lambda_{0}$
and using (5.3) and (5.6), we obtain the target estimate (4.12) of this
theorem. $\square$

###  5.2. Proof of Theorem 4.2

Let
$W_{1}^{h}\left(\mathbf{x}^{h}\right),W_{2}^{h}\left(\mathbf{x}^{h}\right)\in\overline{B\left(R,P^{h}\right)}$
be two arbitrary points. Denote
$v^{h}\left(\mathbf{x}^{h}\right)=W_{2}^{h}\left(\mathbf{x}^{h}\right)-W_{1}^{h}\left(\mathbf{x}^{h}\right).$
By (4.1)

(5.7) $v^{h}\in H_{N,0}^{1,h}\left(\Omega^{h}\right).$

By (4.11)

(5.8)
$J_{\lambda}\left(W_{2}^{h}\right)=J_{\lambda}\left(W_{1}^{h}+v^{h}\right)=\mathop{\displaystyle\sum}\limits_{\left(i,j\right)=\left(1,1\right)}^{\left(i,j\right)=\left(n,m-1\right)}\times$
$\times\mathop{\displaystyle\int}\limits_{a}^{b}\left(D_{N}^{h}\left(W_{1y}^{h}+v_{y}^{h}\right)+\mathop{\displaystyle\sum}\limits_{i=1}^{n}A_{i}^{h}\left(W_{1x_{i}}^{h}+v_{x_{i}}\right)+F^{h}\left(\left(W^{h}+v^{h}\right),\mathbf{x}^{h}\right)\right)^{2}e^{2\lambda
y}dy.$

Consider the integrand in the second line of (5.8) without, however, the term
$e^{2\lambda y}.$ By (4.2), (4.3), (4.9) and (4.10)

(5.9)
$\left(D_{N}^{h}\left(W_{1y}^{h}+v_{y}^{h}\right)+\mathop{\displaystyle\sum}\limits_{i=1}^{n}A_{i}^{h}\left(\mathbf{x}^{h}\right)\left(W_{1x_{i}}^{h}+v_{x_{i}}^{h}\right)+F^{h}\left(W^{h}+v^{h},\mathbf{x}^{h}\right)\right)^{2}$
$=\left[D_{N}^{h}\left(\left(W_{1y}^{h}+v_{y}^{h}\right)+\mathop{\displaystyle\sum}\limits_{i=1}^{n}G_{i}^{h}\left(W_{1x_{i}}^{h}+v_{x_{i}}^{h}\right)+\Phi^{h}\left(W_{1}^{h}+v^{h},\mathbf{x}^{h}\right)\right)\right]^{2}.$

Consider the term $\Phi^{h}\left(W_{1}^{h}+v^{h},\mathbf{x}^{h}\right).$ By
(4.6)-(4.8) we have for all $\mathbf{x}^{h}\in\Omega^{h}$

(5.10)
$\Phi^{h}\left(W_{1}^{h}+v^{h},\mathbf{x}^{h}\right)=\Phi^{h}\left(W_{1}^{h},\mathbf{x}^{h}\right)+\Phi_{1}\left(W_{1}^{h},\mathbf{x}^{h}\right)v^{h}\left(\mathbf{x}^{h}\right)+\Phi_{2}^{h}\left(W_{1}^{h},W_{1}^{h}\mathbf{+}v^{h},\mathbf{x}^{h}\right),$
(5.11) $\left|\Phi_{1}\left(W_{1}^{h},\mathbf{x}^{h}\right)\right|\leq
C_{3},\text{
}\left|\Phi_{2}^{h}\left(W_{1}^{h},W_{1}^{h}\mathbf{+}v^{h},\mathbf{x}^{h}\right)\right|\leq
C_{3}\left(v^{h}\left(\mathbf{x}^{h}\right)\right)^{2}\text{.}$

It follows from (5.10) and (5.11) that the second line of (5.9) is:

$\left[D_{N}^{h}\left(\left(W_{1y}^{h}+v_{y}^{h}\right)+\mathop{\displaystyle\sum}\limits_{i=1}^{n}G_{i}^{h}\left(W_{1x_{i}}^{h}+v_{x_{i}}^{h}\right)+\Phi^{h}\left(W_{1}^{h}+v^{h},\mathbf{x}^{h}\right)\right)\right]^{2}$
$=\left[D_{N}^{h}\left(W_{1y}^{h}+\mathop{\displaystyle\sum}\limits_{i=1}^{n}G_{i}^{h}W_{1x_{i}}^{h}+\Phi^{h}\left(W_{1}^{h},\mathbf{x}^{h}\right)\right)\right]^{2}$
$+2D_{N}^{h}\left(W_{1y}^{h}+\mathop{\displaystyle\sum}\limits_{i=1}^{n}G_{i}^{h}W_{1x_{i}}^{h}+\Phi^{h}\left(W_{1}^{h},\mathbf{x}^{h}\right)\right)\cdot$
(5.12) $\cdot
D_{N}^{h}\left(v_{y}^{h}+\mathop{\displaystyle\sum}\limits_{i=1}^{n}G_{i}^{h}v_{x_{i}}^{h}+\Phi_{1}\left(W_{1}^{h},\mathbf{x}^{h}\right)v^{h}\left(\mathbf{x}^{h}\right)\right)$
$+2D_{N}^{h}\left(W_{1y}^{h}+\mathop{\displaystyle\sum}\limits_{i=1}^{n}G_{i}^{h}W_{1x_{i}}^{h}+\Phi^{h}\left(W_{1}^{h},\mathbf{x}^{h}\right)\right)\cdot\left(D_{N}^{h}\Phi_{2}^{h}\left(W_{1}^{h},W_{1}^{h}\mathbf{+}v^{h},\mathbf{x}^{h}\right)\right)$
$+\left[D_{N}^{h}\left(v_{y}^{h}+\mathop{\displaystyle\sum}\limits_{i=1}^{n}G_{i}^{h}v_{x_{i}}^{h}+\Phi_{1}\left(W_{1}^{h},\mathbf{x}^{h}\right)v^{h}\left(\mathbf{x}^{h}\right)+\Phi_{2}^{h}\left(W_{1}^{h},W_{1}^{h}\mathbf{+}v^{h},\mathbf{x}^{h}\right)\right)\right]^{2}.$

The linear with respect to $v^{h}$ term in (5.12) is the scalar product of the
third and fourth lines of this equality. We denote the latter as
$Lin\left(v^{h}\left(\mathbf{x}^{h}\right),\mathbf{x}^{h}\right).$ The sum of
fifth and sixth lines of (5.12) contains only nonlinear terms with respect to
$v^{h}.$ We denote this sum as
$Nonlin\left(v^{h}\left(\mathbf{x}^{h}\right),\mathbf{x}^{h}\right).$ It
follows from Corollary 4.1, (3.32), (4.2)-(4.8), (5.11) and (5.12) and Cauchy-
Schwarz inequality that

(5.13)
$\left|Nonlin\left(v^{h}\left(\mathbf{x}^{h}\right),\mathbf{x}^{h}\right)\right|\geq
C_{3,3}\left(v_{y}^{h}+\mathop{\displaystyle\sum}\limits_{i=1}^{n}G_{i}^{h}v_{x_{i}}^{h}\right)^{2}-C_{3,4}\left(v^{h}\left(\mathbf{x}^{h}\right)\right)^{2},\forall\mathbf{x}^{h}\in\overline{\Omega^{h}}.$

Just as in the proof of Theorem 4.1, it again makes sense to use two positive
constants $C_{3,3},C_{3,4}>0$ here, both depend on the same parameters as ones
involved in $C_{3}.$ Taking into account (5.8), (5.9) and the second line of
(5.12), we obtain

(5.14)
$J_{\lambda}\left(W_{2}^{h}\right)=J_{\lambda}\left(W_{1}^{h}+v^{h}\right)=J_{\lambda}\left(W_{1}^{h}\right)+$
$+\mathop{\displaystyle\sum}\limits_{\left(i,j\right)=\left(1,1\right)}^{\left(i,j\right)=\left(n,m-1\right)}\mathop{\displaystyle\int}\limits_{a}^{b}Lin\left(v^{h}\left(\mathbf{x}^{h}\right),\mathbf{x}^{h}\right)e^{2\lambda
y}dy+\mathop{\displaystyle\sum}\limits_{\left(i,j\right)=\left(1,1\right)}^{\left(i,j\right)=\left(n,m-1\right)}\mathop{\displaystyle\int}\limits_{a}^{b}Nonlin\left(v^{h}\left(\mathbf{x}^{h}\right),\mathbf{x}^{h}\right)e^{2\lambda
y}dy.$

It follows from (4.4)-(4.6), (5.7) and from the third and fourth lines of
(5.12), which form
$Lin\left(v^{h}\left(\mathbf{x}^{h}\right),\mathbf{x}^{h}\right),$ that the
first term in the second line of (5.14) is a bounded linear functional with
respect to $v^{h}$ mapping $H_{N,0}^{1,h}\left(\Omega^{h}\right)$ in
$\mathbb{R}.$ We denote this functional as
$\widetilde{J}_{\lambda}\left(W_{1}^{h}\right)\left(v^{h}\right).$ By Riesz
theorem there exists unique function
$J_{\lambda}^{\prime}\left(W_{1}^{h}\right)\in
H_{N,0}^{1,h}\left(\Omega^{h}\right)$ such that

(5.15)
$\widetilde{J}_{\lambda}\left(W_{1}^{h}\right)\left(v^{h}\right)=\left(J_{\lambda}^{\prime}\left(W_{1}^{h}\right),v^{h}\right),\forall
v^{h}\in H_{N,0}^{1,h}\left(\Omega^{h}\right),$

where $\left(,\right)$ is the scalar product in
$H_{N,0}^{1,h}\left(\Omega^{h}\right).$ It follows from (3.33), (3.34),
(4.4)-(4.7) and the fifth and sixth lines of (5.12) combined with (5.11) that
for $\left\|v^{h}\right\|_{H_{N}^{1,h}\left(\Omega^{h}\right)}<1$

$\left|\mathop{\displaystyle\sum}\limits_{\left(i,j\right)=\left(1,1\right)}^{\left(i,j\right)=\left(n,m-1\right)}\mathop{\displaystyle\int}\limits_{a}^{b}Nonlin\left(v^{h}\left(\mathbf{x}^{h}\right),\mathbf{x}^{h}\right)e^{2\lambda
y}dy\right|\leq C_{3}e^{2\lambda
b}\left\|v^{h}\right\|_{H_{N}^{1,h}\left(\Omega^{h}\right)}^{2}.\text{ }$

Hence, (5.14) implies that

$\frac{J_{\lambda}\left(W_{1}^{h}+v^{h}\right)-J_{\lambda}\left(W_{1}^{h}\right)-\left(J_{\lambda}^{\prime}\left(W_{1}^{h}\right),v^{h}\right)}{\left\|v^{h}\right\|_{H_{N}^{1,h}\left(\Omega^{h}\right)}}=O\left(\left\|v^{h}\right\|_{H_{N}^{1,h}\left(\Omega^{h}\right)}\right),\text{
as }\left\|v^{h}\right\|_{H_{N}^{1,h}\left(\Omega^{h}\right)}\rightarrow 0.$

Thus,
$J_{\lambda}^{\prime}\left(W_{1}^{h}\right):H_{N,0}^{1,h}\left(\Omega^{h}\right)\rightarrow\mathbb{R}$
is the Fréchet derivative of the functional $J_{\lambda}\left(W^{h}\right)$ at
the point $W_{1}^{h}.$ We omit the proof of estimate (4.13) since this proof
is completely similar with the proof of Theorem 3.1 of [1].

We now prove (4.15), (4.16). Using (5.13)-(5.15), we obtain

$J_{\lambda}\left(W_{1}^{h}+v^{h}\right)-J_{\lambda}\left(W_{1}^{h}\right)-\left(J_{\lambda}^{\prime}\left(W_{1}^{h}\right),v^{h}\right)\geq$
(5.16) $\geq
C_{3,3}\left\|\left(v_{y}^{h}+\mathop{\displaystyle\sum}\limits_{i=1}^{n}G_{i}^{h}v_{x_{i}}^{h}\right)e^{\lambda
y}\right\|_{L_{N}^{2,h}\left(\Omega^{h}\right)}^{2}-C_{3,4}\left\|ve^{\lambda
y}\right\|_{L_{N}^{2,h}\left(\Omega^{h}\right)}^{2}.$

Let the number $\lambda_{0}$ be the one chosen in Theorem 4.1. Recalling (5.7)
and using Theorem 4.1, we estimate the second line of (5.16) for all
$\lambda\geq\lambda_{0},$

(5.17)
$C_{3,3}\left\|\left(v_{y}^{h}+\mathop{\displaystyle\sum}\limits_{i=1}^{n}G_{i}^{h}v_{x_{i}}^{h}\right)e^{\lambda
y}\right\|_{L_{N}^{2,h}\left(\Omega^{h}\right)}^{2}-C_{3,4}\left\|v^{h}e^{\lambda
y}\right\|_{L_{N}^{2,h}\left(\Omega^{h}\right)}^{2}\geq$ $\geq
C_{3}\left(\left\|v_{y}^{h}e^{\lambda
y}\right\|_{L_{N}^{2,h}\left(\Omega^{h}\right)}^{2}+\lambda^{2}\left\|v^{h}e^{\lambda
y}\right\|_{L_{N}^{2,h}\left(\Omega^{h}\right)}^{2}\right)-C_{3,4}\left\|v^{h}e^{\lambda
y}\right\|_{L_{N}^{2,h}\left(\Omega^{h}\right)}^{2}.$

Choosing $\lambda_{1}\geq\lambda_{0}$ so large that
$\lambda_{1}^{2}C_{3}/2>2C_{3,4}$ and using (5.17), we obtain

(5.18)
$C_{3,3}\left\|\left(v_{y}^{h}+\mathop{\displaystyle\sum}\limits_{i=1}^{n}G_{i}^{h}v_{x_{i}}^{h}\right)e^{\lambda
y}\right\|_{L_{N}^{2,h}\left(\Omega^{h}\right)}^{2}-C_{3,4}\left\|v^{h}e^{\lambda
y}\right\|_{L_{N}^{2,h}\left(\Omega^{h}\right)}^{2}\geq$ $\geq
C_{3}\left(\left\|v_{y}e^{\lambda
y}\right\|_{L_{N}^{2,h}\left(\Omega^{h}\right)}^{2}+\lambda^{2}\left\|ve^{\lambda
y}\right\|_{L_{N}^{2,h}\left(\Omega^{h}\right)}^{2}\right).$

Combining (5.16) with (5.18), we obtain (4.15), (4.16).

Finally, given (4.15) and (4.16), existence and uniqueness of the minimizer
$W_{\min,\lambda}^{h}$ of the functional $J_{\lambda}\left(W^{h}\right)$ on
the set $\overline{B\left(R,P^{h}\right)}$ as well as estimate (4.17) follow
immediately from a combination of Lemma 2.1 with Theorem 2.1 of [1]. $\ \
\square$

###  5.3. Proof of Theorem 4.4

Denote

(5.19) $B_{0}\left(2R\right)=\left\\{V^{h}\in
H_{N,0}^{1,h}\left(\Omega^{h}\right):\left\|V^{h}\right\|_{H_{N}^{1,h}\left(\Omega^{h}\right)}<2R\right\\}.$

Consider the vector functions $V^{h},V^{h\ast}$

(5.20) $V^{h}=W^{h}-S^{h},\forall W^{h}\in B\left(R,P^{h}\right);\text{
}V^{h\ast}=W^{h\ast}-S^{h\ast}.$

By (4.18), (4.21), (4.22), (5.19) and (5.20)

(5.21) $V^{h},V^{h\ast}\in B_{0}\left(2R\right).$

Consider now the functional $I_{\lambda}\left(V^{h}\right),$

(5.22)
$I_{\lambda}\left(V^{h}\right):B_{0}\left(2R\right)\rightarrow\mathbb{R},I_{\lambda}\left(V^{h}\right)=J_{\lambda}\left(V^{h}+S^{h}\right).$

An obvious analog of Theorem 4.2 holds for $I_{\lambda}\left(V^{h}\right).$
However, it follows from (5.21) that we need to replace $R$ with $2R$ in
(4.14) in this case, i.e. we need to use now the number $\lambda_{2}$ in
(4.24). Let $V_{\min,\lambda_{2}}^{h}$ be the minimizer of
$I_{\lambda_{2}}\left(V^{h}\right)$ on the set $B_{0}\left(2R\right).$

Consider
$I_{\lambda}\left(V^{h\ast}\right)=J_{\lambda}\left(V^{h\ast}+S^{h}\right).$
By (4.15), (4.16) and (5.22)

$I_{\lambda_{2}}\left(V^{h\ast}\right)-I_{\lambda_{2}}\left(V_{\min,\lambda_{2}}^{h}\right)-I_{\lambda_{2}}^{\prime}\left(V_{\min,\lambda_{2}}^{h}\right)\left(V^{h\ast}-V_{\min,\lambda_{2}}^{h}\right)\geq
C_{3}\left\|V^{h\ast}-V_{\min,\lambda_{2}}^{h}\right\|_{H_{N}^{1,h}\left(\Omega^{h}\right)}^{2}.$

By (4.17)
$-I_{\lambda_{2}}^{\prime}\left(V_{\min,\lambda_{2}}^{h}\right)\left(V^{h\ast}-V_{\min,\lambda_{2}}^{h}\right)\leq
0.$ Since $-I_{\lambda_{2}}\left(V_{\min,\lambda_{2}}^{h}\right)\leq 0$ as
well, then the latter estimate implies

(5.23)
$\left\|V^{h\ast}-V_{\min,\lambda_{2}}^{h}\right\|_{H_{N}^{1,h}\left(\Omega^{h}\right)}^{2}\leq
C_{3}I_{\lambda_{2}}\left(V^{h\ast}\right).$

Next, by (5.20) and (5.22)

(5.24)
$I_{\lambda_{2}}\left(V^{h\ast}\right)=J_{\lambda_{2}}\left(V^{h\ast}+S^{h}\right)=J_{\lambda_{2}}\left(W^{h\ast}+\left(S^{h}-S^{h\ast}\right)\right).$

By (4.19) the right hand side of (4.11) equals zero if $W^{h}$ is replaced
with $W^{h\ast}.$ Hence, using (4.2)-(4.8), (4.23) and (5.20), we obtain
$J_{\lambda_{2}}\left(W^{h\ast}+\left(S^{h}-S^{h\ast}\right)\right)\leq
C_{3}\delta^{2}.$ Hence, (5.23) and (5.24) lead to:

(5.25)
$\left\|V^{h\ast}-V_{\min,\lambda_{2}}^{h}\right\|_{H_{N}^{1,h}\left(\Omega^{h}\right)}\leq
C_{3}\delta.$

Using again (4.23), (5.20), triangle inequality and (5.25), we obtain

(5.26)
$\left\|W^{h\ast}-\widetilde{W}_{\min,\lambda_{2}}^{h}\right\|_{H_{N}^{1,h}\left(\Omega^{h}\right)}\leq
C_{3}\delta.$

Here $\widetilde{W}_{\min,\lambda_{2}}^{h}=V_{\min,\lambda_{2}}^{h}-S^{h}.$ By
(4.25), (5.26) and the triangle inequality
$\left\|\widetilde{W}_{\min,\lambda_{2}}^{h}\right\|_{H_{N}^{1,h}\left(\Omega^{h}\right)}<R.$
Hence, (4.1), (4.21), (5.19) and (5.21) imply that

(5.27) $\widetilde{W}_{\min,\lambda_{2}}^{h}\in B\left(R,P^{h}\right).$

Now, let $W_{\min,\lambda_{2}}^{h}$ be the minimizer of the functional
$J_{\lambda_{2}}\left(W^{h}\right)$ on the set
$\overline{B\left(R,P^{h}\right)},$ which is claimed by Theorem 4.2. Let
$Q_{\min,\lambda_{2}}^{h}=W_{\min,\lambda_{2}}^{h}-S^{h}.$ Then
$Q_{\min,\lambda_{2}}^{h}\in B_{0}\left(2R\right).$ Hence,

$J_{\lambda_{2}}\left(\widetilde{W}_{\min,\lambda_{2}}^{h}\right)=J_{\lambda_{2}}\left(V_{\min,\lambda_{2}}^{h}+G\right)\leq
J_{\lambda_{2}}\left(Q_{\min,\lambda_{2}}^{h}+G\right)=\min_{\overline{B\left(R,P^{h}\right)}}J_{\lambda_{2}}\left(W^{h}\right).$

Hence, by (5.27) $\widetilde{W}_{\min,\lambda_{2}}^{h}$ is also a minimizer of
the functional $J_{\lambda_{2}}\left(W^{h}\right)$ on the set
$\overline{B\left(R,P^{h}\right)}.$ However, since such a minimizer is unique
by Theorem 4.2, then
$\widetilde{W}_{\min,\lambda_{2}}^{h}=W_{\min,\lambda_{2}}^{h}\in
B\left(R,P^{h}\right).$ Hence, estimate (5.26) holds when
$\widetilde{W}_{\min,\lambda_{2}}^{h}$ is replaced with
$W_{\min,\lambda_{2}}^{h}.$ The latter immediately implies (4.26).$\ \square$

###  5.4. Proof of Theorem 4.5

Consider again the minimizer $W_{\min,\lambda_{2}}^{h}$ of the functional
$J_{\lambda_{2}}\left(W^{h}\right)$ on the set
$\overline{B\left(R,P^{h}\right)}.$ Then (4.28) and Theorem 4.4 imply that
$W_{\min,\lambda_{2}}^{h}\in B\left(R/3,P^{h}\right).$ To prove (4.29), we now
modify the proof of Theorem 6 of [25]: to adapt it for our specific case.

We set $\lambda=\lambda_{2},$ where $\lambda=\lambda_{2}$ is defined in
(4.24). It obviously follows from (4.15) that

(5.28)
$\left(J_{\lambda_{2}}^{\prime}\left(W_{1}^{h}\right)-J_{\lambda_{2}}^{\prime}\left(W_{2}^{h}\right)\right)\left(W_{1}^{h}-W_{2}^{h}\right)\geq
C_{3}\left\|W_{2}^{h}-W_{1}^{h}\right\|_{H_{N}^{1,h}\left(\Omega^{h}\right)}^{2}$

for all $W_{1}^{h},W_{2}^{h}\in B\left(R,P^{h}\right).$ Consider the operator
$T\left(W^{h}\right):B\left(R,P^{h}\right)\rightarrow
H_{N}^{1,h}\left(\Omega^{h}\right),$

(5.29) $T\left(W^{h}\right)=W^{h}-\gamma
J_{\lambda_{2}}^{\prime}\left(W^{h}\right).$

Using (4.13) and (5.28), a contraction property of the operator
$X\left(W^{h}\right),$ it was proven in [1, Theorem 2.1] and in [26, page 97]
that there exists a sufficiently small number $\gamma>0$ and a number
$\theta=\theta\left(\gamma\right)\in\left(0,1\right)$ such that

(5.30)
$\left\|T\left(W_{1}^{h}\right)-T\left(W_{2}^{h}\right)\right\|_{H_{N}^{1,h}\left(\Omega^{h}\right)}\leq\theta\left\|W_{1}^{h}-W_{2}^{h}\right\|_{H_{N}^{1,h}\left(\Omega^{h}\right)}$

for all $W_{1}^{h},W_{2}^{h}\in B\left(R,P^{h}\right).$ However, it was not
proven in [1, 26] that the operator $T$ maps $B\left(R,P^{h}\right)$ in
itself.

Since $W_{\min,\lambda_{2}}^{h}$ is an interior point of the set
$B\left(R,P^{h}\right),$ then
$J_{\lambda_{2}}^{\prime}\left(W_{\min,\lambda_{2}}^{h}\right)=0.$ Hence, by
(5.29)

(5.31) $W_{\min,\lambda_{2}}^{h}=T\left(W_{\min,\lambda_{2}}^{h}\right).$

Consider now the vector function
$W_{1}^{h}=T\left(W_{0}^{h}\right)=W_{0}^{h}-\gamma
J_{\lambda_{2}}^{\prime}\left(W_{0}^{h}\right).$ Since by Theorem 4.2 the
function $J_{\lambda_{2}}^{\prime}\left(W_{0}^{h}\right)\in
H_{N,0}^{1,h}\left(\Omega^{h}\right),$ then
$W_{1}^{h}\mid_{\partial\Omega^{h}}=W_{0}^{h}\mid_{\partial\Omega^{h}}=P^{h}.$
Next by (5.30), (5.31) and triangle inequality

$\left\|W_{1}^{h}-W_{\min,\lambda_{2}}^{h}\right\|_{H_{N}^{1,h}\left(\Omega^{h}\right)}=\left\|T\left(W_{0}^{h}\right)-T\left(W_{\min,\lambda_{2}}^{h}\right)\right\|_{H_{N}^{1,h}\left(\Omega^{h}\right)}\leq$
(5.32)
$\leq\theta\left\|W_{0}^{h}-W_{\min,\lambda_{2}}^{h}\right\|_{H_{N}^{1,h}\left(\Omega^{h}\right)}\leq\frac{2R}{3}\theta<\frac{2R}{3}.$

Hence, (4.25), (4.26) and triangle inequality imply
$\left\|W_{1}^{h}\right\|_{H_{N}^{1,h}\left(\Omega^{h}\right)}<R.$ Hence, by
(4.1) $W_{1}^{h}\in B\left(R,P^{h}\right).$ Similarly by (5.32)

$\left\|W_{2}^{h}-W_{\min,\lambda_{2}}^{h}\right\|_{H_{N}^{1,h}\left(\Omega^{h}\right)}\leq\theta\left\|W_{1}^{h}-W_{\min,\lambda_{2}}^{h}\right\|_{H_{N}^{1,h}\left(\Omega^{h}\right)}\leq$
$\leq\theta^{2}\left\|W_{0}^{h}-W_{\min,\lambda_{2}}^{h}\right\|_{H_{N}^{1,h}\left(\Omega^{h}\right)}<\frac{2R}{3}.$

Hence, applying again triangle inequality, we obtain $W_{2}^{h}\in
B\left(R,P^{h}\right).$ Continuing this process, we obtain (4.29).

To prove (4.30), we note that by triangle inequality, (4.26) and (4.29)

$\left\|W_{n}^{h}-W^{h\ast}\right\|_{H_{N}^{1,h}\left(\Omega^{h}\right)}\leq\left\|W_{\min,\lambda_{2}}^{h}-W^{h\ast}\right\|_{H_{N}^{1,h}\left(\Omega^{h}\right)}+\left\|W_{n}^{h}-W_{\min,\lambda_{2}}^{h}\right\|_{H_{N}^{1,h}\left(\Omega^{h}\right)}\leq$
$\leq
C_{3}\delta+\theta^{n}\left\|W_{0}^{h}-W_{\min,\lambda_{2}}^{h}\right\|_{H_{N}^{1,h}\left(\Omega^{h}\right)}.$

The proof of estimate (4.31) follows immediately from (3.35)-(3.37) and
(4.30). $\square$

##  6\. Numerical Studies

In order not to introduce new and complicated notations, we slightly abuse
below some notations of the previous sections. Nevertheless, the substance is
always clear from the context presented below.

###  6.1. Numerical implementation

We have conducted our numerical studies in the 2-D case. Below
$\mathbf{x}=\left(x,y\right)$ and, according to (2.2) and (2.6)

(6.1)
$\Omega=\\{\mathbf{x}:x\in\left(-A,A\right),y\in\left(a,b\right),A=1/2,a=1,b=2,$
$\Gamma_{d}=\left\\{\mathbf{x}_{\alpha}=(\alpha,0):\alpha\in[-d,d]\right\\},d=1/2.$

As to the kernel $K(\mathbf{x},\alpha,\beta)$ of the integral operator in
(2.10), we work below with the 2D Henyey-Greenstein function [17]:

(6.2)
$K(\mathbf{x},\alpha,\beta)=H(\alpha,\beta)=\frac{1}{2d}\left[\frac{1-g^{2}}{1+g^{2}-2g\cos(\alpha-\beta)}\right],\quad
g=\frac{1}{2}.$

Here, $g=1/2$ means an anisotropic scattering, which is half ballistic with
$g=0$ an half isotropic scattering with $g=1$ [11, 12, 13]. We take the same
function $f\left(\mathbf{x}\right)$ as one in (2.7), (2.8) with
$\varepsilon=0.05.$

We assume that

(6.3)
$\mu_{s}(\mathbf{x})=\left\\{\begin{array}[]{ll}5,&\quad\mathbf{x}\in\Omega,\\\
0,&\quad\mathbf{x}\in\mathbb{R}^{2}\setminus\Omega.\end{array}\right.$

Given (6.3), we use formula (2.13) for the coefficient function
$a(\mathbf{x})$ and we take in this formula

(6.4) $\mu_{a}(\mathbf{x})=\left\\{\begin{array}[]{cc}c=const.>0,&\text{inside
the tested inclusion,}\\\ 0,&\text{outside the tested
inclusion.}\end{array}\right.$

We perform the numerical tests with a variety of the value of the parameter
$c=5,10,15,20,30$, see below. Therefore, by (6.3) and (6.4)

(6.5) $\text{inclusion/background contrast}=1+\frac{c}{5}.$

From the Physics standpoint $\mu_{s}(\mathbf{x})=5$ in the domain $\Omega$
means that an average particle scatters every 1/5 of the unit. Hence, by (6.1)
and (6.3) the maximal average number of scattering events for a particle
emitted from a point $\mathbf{x}_{\alpha}\in\Gamma_{d}$ and entering $\Omega$
is around $5$ before this particle leaves $\Omega$ [11, 12, 13]. This might
happen in optics before the true diffusion occurs, e.g. in the case of the so-
called “snake photons” [10].

In computational results below $c\neq const.$ inside of computed shapes of
letters. Hence, we set for computed inclusion/background contrasts we replace
(6.5) with:

(6.6) $\text{computed inclusion/background
contrast}=1+\frac{1}{5}\max_{\overline{\Omega}}\mu_{a}(\mathbf{x}).$

To solve the forward problem (2.10), (2.18), we have solved integral equation
(2.36) with the function $u_{0}$ taken from (2.34). To do this, we have used
the discrete form of (2.34), (2.36) and the trapezoidal rule. The
discretization steps with respect to $x,y,\alpha$ where
$h_{x}=h_{y}=h_{\alpha}=1/40.$ The discretized integral equation (2.36) was
solved as a linear system using the Matlab function ’$\backslash$’. Thus, the
solution of this forward problem has provided computationally simulated data
for the inverse problem to us.

To minimize the convexification functional $J_{\lambda}\left(W^{h}\right)$ in
(4.11), we have written in the finite differences form not only the
$x-$derivative as in section 4 but the $y-$derivative as well. Also, integrals
with respect to $\alpha$ were written in the discrete form using the
trapezoidal rule. The mesh sizes were different from ones for the forward
problem. They were:

(6.7) $h_{x}=h_{y}=h_{\alpha}=h=1/20.$

Then we have minimized the resulting functional
$J_{\lambda,dis}\left(W^{h}\right)$,

(6.8)
$J_{\lambda,dis}\left(W^{h}\right)=\left\|\left(D_{N}^{h}W_{y}^{h}+A^{h}W_{x}^{h}+F^{h}\left(W^{h}\left(\mathbf{x}^{h}\right),\mathbf{x}^{h}\right)\right)e^{\lambda
y}\right\|_{L_{N}^{2,h}\left(\Omega^{h}\right)}^{2}$

in its fully discrete form with respect to the values of the vector function
$W^{h}$ at the grid points. Vector functions and matrices in (6.8) are full
analogs of those in (4.11) with the only difference that they are fully
discrete in the above sense, rather than ‘partially’ discrete as in sections
4,5. The same is true for the norm
$\left\|\cdot\right\|_{L_{N}^{2,h}\left(\Omega^{h}\right)}.$

It follows from (6.1) and (6.7) that we had total $20\times 20\times N$
unknown parameters in our minimization procedure. To solve the minimization
problem, we have used the Matlab’s built-in function fminunc with the quasi-
newton algorithm. The iterations of the function fminunc were stopping when
the following inequality occurred at the iteration number $k$:

$\left|J_{\lambda,dis}\left(W_{k}^{h}\right)\right|<10^{-2}.$

By (3.14) our technique requires computations of first derivatives with
respect to $\alpha.$ Note, however, that (3.8) and (3.12) imply that it is not
necessary to calculate the $\alpha-$derivative of the boundary data, which is
an advantage, since boundary data are noisy. The derivatives
$\partial_{\alpha}$ of functions $\Psi_{s}\left(\alpha\right)$ and the
function $K(\mathbf{x},\alpha,\beta)$ were calculated via finite differences.

We have introduced the random noise in the boundary data
$g_{1}(\mathbf{x},\alpha)$ in (3.4) on the boundary $\partial\Omega$,

(6.9)
$g_{1}(\mathbf{x},\alpha)=g_{1}(\mathbf{x},\alpha)\left(1+\sigma\zeta_{\mathbf{x}}\right).$

Here $\zeta_{\mathbf{x}}$ is the uniformly distributed random variable in the
interval $[0,1]$ depending on the point $\mathbf{x}\in\partial\Omega$ with
$\sigma=0.03$ and $\sigma=0.05,$ which correspond respectively to $3\%$ and
$5\%$ noise level.

To solve the minimization problem, we need to provide the starting
$W_{0}^{h}(\mathbf{x}^{h})$ for iterations. Due to the global convergence
property of our method, the vector function
$W_{0}^{h}(\mathbf{x}^{h})=\left(w_{0,0}^{h}(\mathbf{x}^{h}),...,w_{N-1,0}^{h}(\mathbf{x}^{h})\right)^{T}$
should not have any information about the exact solution
$W^{h\ast}(\mathbf{x}^{h}).$ On the other hand, due to (4.1), we should have
$W_{0}^{h}\left(\mathbf{x}^{h}\right)\mid_{\partial\Omega^{h}}=P^{h}\left(\mathbf{x}^{h}\right).$
Therefore, in all numerical tests below we choose the starting point as the
discrete version of the following vector function:

$w_{s,0}(x,y)=\frac{1}{2}\left(\frac{(A-x)}{2A}w_{s}(-A,y)+\frac{(x+A)}{2A}w_{s}(A,y)\right)+$

(6.10)
$+\frac{1}{2}\left(\frac{(b-y)}{b-a}w_{s}(x,a)+\frac{(y-a)}{b-a}w_{s}(x,b)\right),s=0,...,N-1.$

Expression (6.10) represents the average of linear interpolations inside of
the square $\Omega$ with respect to $x-$direction and $y-$direction of the
boundary condition for $w_{s}\left(x,y\right)$ .

###  6.2. Numerical results

Recall that we reconstruct the coefficient
$a\left(\mathbf{x}\right)=\mu_{s}\left(\mathbf{x}\right)+\mu_{a}\left(\mathbf{x}\right),$
where functions $\mu_{s}\left(\mathbf{x}\right)$ and
$\mu_{a}\left(\mathbf{x}\right)$ are given in (6.3) and (6.4) respectively.
Our results for Tests 1-4 are for noiseless data and the results for Test 5
are for noisy data as in (6.9).

To demonstrate a good performance of our technique, we intentionally test it
for rather complicated shapes of inclusions, which are non convex and have
voids. More precisely, our inclusions are letters ‘A’, ‘$\Omega$’ and two
letters jointly ‘SZ’. ‘SZ’ stands for ‘Shenzhen’, the city where the campus of
the Southern University of Science and Technology, the work place of the third
and fourth authors, is located. Thus, in our tests the coefficients
$\mu_{a}(\mathbf{x})$ in (6.4) have shapes of those letters located inside of
the $1\times 1$ square $\Omega$ defined in (6.1).

Test 1. We test the letter ‘A’ with $c=5$ in (6.4). This is our reference
case. More precisely, we use this test to figure out optimal values of
parameters $N$ and $\lambda.$ As soon as optimal parameters are selected, we
use them then for all other tests.

First, we select $N$. To do this, we solve the forward problem (2.10)_,_(2.18)
for the case when the functions $\mu_{s}\left(\mathbf{x}\right)$ and (6.3)__
and (6.4) respectively and in $c=5.$ Hence, by (6.5) the inclusion/background
contrast is 2:1 in this case. Next, we calculate norms
$\left\|w_{s}\left(\mathbf{x}\right)\right\|_{L_{2}\left(\Omega\right)}$ and
compare them. We have observed that the $L_{2}\left(\Omega\right)-$norm of the
function $w_{s}(\mathbf{x})$ decreases very rapidly when the number $s$ is
growing. More precisely, we have obtained that

(6.11)
$\frac{\mathop{\displaystyle\sum}\limits_{s=3}^{11}\left\|w_{s}\left(\mathbf{x}\right)\right\|_{L_{2}\left(\Omega\right)}}{\mathop{\displaystyle\sum}\limits_{s=0}^{11}\left\|w_{s}\left(\mathbf{x}\right)\right\|_{L_{2}\left(\Omega\right)}}=0.0084,$

which means less than 1%. The values of norms
$\left\|w_{s}\left(\mathbf{x}\right)\right\|_{L_{2}\left(\Omega\right)}$ for
$s=0,...,11$ are displayed in Table 1 and Figure 1. One can observe that
starting from $s=3$, these norms are much less than those for $s=0,1,2.$ We
conclude therefore, that we should take in our tests

(6.12) $N=3.$

Table 1. The $L_{2}\left(\Omega\right)-$norms of functions $w_{s}\left(\mathbf{x}\right)$, $s=0,1,...,11$ for the reference Test 1 with $c=5$ in (6.4). $s$ | 0 | 1 | 2 | 3 | 4 | 5
---|---|---|---|---|---|---
$\left\|w_{s}(\mathbf{x})\right\|_{L_{2}}$ | 5.7122 | 1.6383 | 0.1630 | 0.0118 | 0.0091 | 0.0077
$s$ | 6 | 7 | 8 | 9 | 10 | 11
$\left\|w_{s}(\mathbf{x})\right\|_{L_{2}}$ | 0.0067 | 0.0061 | 0.0055 | 0.0057 | 0.0058 | 0.0054
Figure 1. The decrease with respect to $s$ of the
$L_{2}\left(\Omega\right)-$norms of functions $w_{s}\left(\mathbf{x}\right)$,
$s=0,1,...,11$ for the reference Test 1 with $c=5$ in (6.4).

Next, given the optimal value of $N=3$ in (6.12), we select the optimal value
of the parameter $\lambda$ of the Carleman Weight Function $e^{2\lambda y}$ in
(6.8). To do this, we test the same letter ‘A’ with $c=5$ inside of it for
values of the parameter $\lambda=0,1,2,3,4,5,6,7.$ Our numerical results are
presented on Figure 2. We observe that the images have a low quality for
$\lambda=0,1.$ Then the quality is improved and is stabilized at $\lambda=5.$
Thus, we treat $\lambda=5$ as the optimal value of this parameter. We use this
value in all subsequent tests. The value $\lambda=5$ tells us that even though
our theorems 4.1-4.5 require sufficiently large values of the parameter
$\lambda,$ the computational practice shows that a reasonable value of
$\lambda$ can be chosen. The same observation was made in all previous works
on the convexification of this research group, see items 2 and 3 in Remarks
4.1 in the end of section 4.

Figure 2. The reconstructed coefficient $a\left(\mathbf{x}\right)$, where the
function $\mu_{a}\left(\mathbf{x}\right)$ is given in (6.4) with $c=5$ inside
of the letter ‘A’. The goal of here is to test different values of the
parameter $\lambda=0,1,2,3,4,5,6,7$ for $N=3$ as in (6.12). The value of
$\lambda$ can be seen on the top side of each square. The images have a low
quality for $\lambda=0,1.$ Then the quality is improved and is stabilized at
$\lambda=5.$ Thus, we select $\lambda=5$ as an optimal value of this parameter
for all follow up tests.

Now we want to demonstrate numerically again that $N=3$ is indeed a good
choice of $N$ for our optimal value of $\lambda=5.$ Taking $\lambda=5,$ we
test the same letter ‘A’ as above with $c=5$ in it, but for $N=1,2,3,5,7,12.$
The results are displayed in Figure 3. One can observe that reconstructions
have a low quality for $N=1,2$. Next, the reconstructions are basically the
same for $N=3,5,7,12.$ However, the computational cost increases very rapidly
with the increase of $N$. Thus, using also Table 1, Figure 1 and (6.12), we
conclude that to balance between the reconstruction accuracy and the
computational cost, we should use $N=3$. Thus, in all subsequent computations
we use

(6.13) $N=3,\lambda=5.$

Figure 3. The reconstructed coefficient $a\left(\mathbf{x}\right)$, where the
function $\mu_{a}\left(\mathbf{x}\right)$ is given in (6.4) with $c=5$ inside
of the letter ‘A’. We took the optimal value of the parameter $\lambda=5$ (see
Figure 2) and have tested different values of the parameter $N=1,2,3,5,7,12$.
A low quality can be observed for $N=1,2$. The reconstructions are basically
the same for $N=3,5,7,12$. However, the computational cost increases very
rapidly with the increase of $N$, which is explained by (6.11), Table 1 and
Figure 1. We conclude, therefore, that to balance between the reconstruction
accuracy and the computational cost, we should use $N=3$ as in (6.12). Thus,
we use below $\lambda=5$ and $N=3$.

Test 2. We test the reconstruction of the coefficient $a(\mathbf{x})$ with the
shape of the letter ‘A’ where the function $\mu_{a}\left(\mathbf{x}\right)$ is
given in (6.4) with different values of the parameter $c=15,20,30$ inside of
the letter ‘A’. Thus, by (6.6) the inclusion/background contrasts now are
respectively $4:1$, $5:1$ and $6:1$. Our computational results for this test
are displayed on Figure 4. One can observe that the quality of these images is
good for all four cases, although it slightly deteriorates for $c=20$ and
$c=30.$ The computed inclusion/background contrast is accurate, see (6.6).

Figure 4. Test 2. Exact (left) and reconstructed (right) coefficient
$a(\mathbf{x})$ for $c=10,15,20,30$ inside of the letter ‘A’ as in (6.4) for
$N=3,\lambda=5$, see (6.13). Thus, by (6.6) the inclusion/background contrasts
now are respectively $4:1$, $5:1$ and $6:1$. The image quality remains
basically the same for all these values of the parameter $c$, although some
deterioration of this quality can be observed for $c=20$ and $c=30$. The
computed inclusion/background contrasts (6.6) are accurate.

Test 3. We test the reconstruction of the coefficient $a(\mathbf{x})$ with the
shape of the letter ‘$\Omega$’ where the function
$\mu_{a}\left(\mathbf{x}\right)$ is given in (6.4) with $c=5$ inside of the
letter ‘$\Omega$’. Results are presented on Figure 5. We again observe an
accurate reconstruction.

Figure 5. Test 3. Exact (left) and reconstructed (right) coefficient
$a(\mathbf{x})$ for the case when the function
$\mu_{a}\left(\mathbf{x}\right)$ is given in (6.4) with $c=5$ inside of the
letter ‘$\Omega$’. The reconstruction is accurate.

Test 4. We test the reconstruction of the coefficient $a(\mathbf{x})$ with the
shape of two letters ‘SZ’ where the function $\mu_{a}\left(\mathbf{x}\right)$
is given in (6.4) with $c=5$ inside of each of these two letters and
$\mu_{a}\left(\mathbf{x}\right)=0$ outside of each of these two letters. In
this test, $N=3,\lambda=5$ as in (6.13). Results are presented on Figure 6.
The image quality is lower than one for the case of the single letter
‘$\Omega$’ on Figure 5. Nevertheless, the quality is still good and the
computed inclusion/background contrasts (6.6) are accurate in both letters.

Figure 6. Test 4. Exact (left) and reconstructed (right) coefficient
$a(\mathbf{x})$ for the case when the function
$\mu_{a}\left(\mathbf{x}\right)$ is given in (6.4) with $c=5$ with the shape
of two letters ‘SZ’. In (6.4) $c=5$ inside of each of these two letters and
$\mu_{a}\left(\mathbf{x}\right)=0$ outside of each of these two letters. Here
$N=3,\lambda=5$ as in (6.13). The image quality is lower than one for the case
of the single letter $\Omega$ on Figure 5. Nevertheless, the quality is still
good and the computed inclusion/background contrasts (6.6) are accurate in
both letters.

Test 5. In this test we use noisy data as in (6.9) with $\sigma=0.03$ and
$\sigma=0.05,$ i.e. with 3% and 5% noise level. We test the reconstruction of
the coefficient $a(\mathbf{x})$ with the shape of either the letter ‘A’ or the
letter ‘$\Omega$’, where the function $\mu_{a}\left(\mathbf{x}\right)$ is
given in (6.4) with $c=5$ inside of each of these two letters. Again, $N=3$,
$\lambda=5$ as in (6.13). The results are shown in Figure 7. One can observe
accurate reconstructions in all four cases. In particular, the
inclusion/background contrasts (6.6) are reconstructed accurately.

Figure 7. Reconstructed coefficient $a(\mathbf{x})$ with the shape of letters
‘A’ (top) and ‘$\Omega$’ (bottom) with $c=5$ from noise polluted observation
data as in (6.9) with $\sigma=0.03$ (left) and $\sigma=0.05$ (right), i.e.
with 3% (left) and 5% (right) noise level. Here $N=3$ and $\lambda=5$ as in
(6.13). One can observe accurate reconstructions in all four cases. In
particular, the inclusion/background contrasts (6.6) are reconstructed
accurately.

## 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 Analysis: Real World Applications 34 (2017), 201–224.
  * [2] G. Bal and A. Jollivet, _Generalized stability estimates in inverse transport theory_ , Inverse Problems and Imaging 12 (2018), 59–90.
  * [3] 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. Numerical Analysis 55 (2017), 1578–1613.
  * [4] L. Baudouin, M. de Buhan, S. Ervedoza, and A. Osses, _Carleman-based reconstruction algorithm for the waves_ , SIAM Journal on Numerical Analysis 59 (2021), 998–1039.
  * [5] M. Bellassoued and M. Yamamoto, _Carleman estimates and applications to inverse problems for hyperbolic systems_ , Springer, Japan, 2017\.
  * [6] M. Born and E. Wolf, _Principles of optics_ , 7th ed., Cambridge University Press, 1999.
  * [7] M Boulakia, M. de Buhan, and E. Schwindt, _Numerical reconstruction based on Carleman estimates of a source term in a reaction-diffusion equation_ , ESAIM: Control, Optimisation and Calculus of Variations 27 (2021), S27, 1–34.
  * [8] A. L. Bukhgeim and M. V. Klibanov, _Uniqueness in the large of a class of multidimensional inverse problems_ , Soviet Math. Doklady 17 (1981), 244–247.
  * [9] G. Chavent, _Nonlinear Least Squares for Inverse Problems: Theoretical Foundations and Step-by-Step Guide for Applications_ , Springer Science & Business Media, Berlin, 2010.
  * [10] B. B. Das, Feng Liu and R. R. Alfano, _Time-resolved fluorescence and photon migration studies in biomedical and model random media_ , Rep. Prog. Phys. 60 (1997) 227–292.
  * [11] H. Fujiwara, K. Sadiq, and A. Tamasan, _A Fourier approach to the inverse source problem in an absorbing and anisotropic scattering medium_ , Inverse Problems 36 (2020), 015005.
  * [12] H. Fujiwara, K. Sadiq, and A. Tamasan, _Numerical reconstruction of radiative sources in an absorbing and nondiffusing scattering medium in two dimensions_ , SIAM J. Imaging Sciences 13 (2020), 535–555.
  * [13] H. Fujiwara, K. Sadiq, and A. Tamasan, _A source reconstrution method in two dimensional radiative transport using boundary data measured on an arc,_ Inverse Problems 37 (2021), 115005.
  * [14] F. Gölgeleyen and M. Yamamoto, _Stability for some inverse problems for transport equations_ , SIAM J. Mathematical Analysis 48 (2016), 2319–2344.
  * [15] A. V. Goncharsky and S. Y. Romanov, _A method of solving the coefficient inverse problems of wave tomography_ , Computers and Mathematics with Applications 77 (2019), 967–980.
  * [16] J.P. Guillement and R. G. Novikov, _Inversion of weighted Radon transforms via finite Fourier series weight approximation_ , Inverse Problems in Science and Engineering 22 (2013), 787–802.
  * [17] J. Heino, S. Arridge, J. Sikora, and E. Somersalo, _Anisotropic effects in highly scattering media_ , Physical Review E 68 (2003), 03198.
  * [18] S. I. Kabanikhin, K. K. Sabelfeld, N. S. Novikov, and M. A. Shishlenin, _Numerical solution of the multidimensional Gelfand-Levitan equation_ , J. Inverse and Ill-Posed Problems 23 (2015), 439–450.
  * [19] V. A. Khoa, G. W. Bidney, M. V. Klibanov, L. H. Nguyen, A. Sullivan, L. Nguyen, and V. N. Astratov, _Convexification and experimental data for a 3D inverse scattering problem with the moving point source_ , Inverse Problems 36 (2020), 085007.
  * [20] M. V. Klibanov, _Global convexity in a three-dimensional inverse acoustic problem_ , SIAM J. Math. Anal. 28 (1997), 1371–1388.
  * [21] M. V. Klibanov and S.E. Pamyatnykh, _Global uniqueness for a coefficient inverse problem for the non-stationary transport equation via Carleman estimate_ , J. Mathematical Analysis and Applications 343 (2008), 352–365.
  * [22] M. V. Klibanov, _Carleman estimates for global uniqueness, stability and numerical methods for coefficient inverse problems_ , J. Inverse and Ill-Posed Problems 21 (2013), 477–510.
  * [23] M. V. Klibanov, _Convexification of restricted Dirichlet to Neumann map_ , J. Inverse and Ill-Posed Problems 25 (2017), no. 5, 669–685.
  * [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 (1995), 147–179.
  * [25] 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_ , Journal of Applied and Industrial Mathematics 15 (2021), 413–436.
  * [26] M. V. Klibanov and J. Li, _Inverse Problems and Carleman Estimates: Global Uniqueness, Global Convergence and Experimental Data_ , De Gruyter, 2021.
  * [27] M. V. Klibanov, J. Li, and W. Zhang, _Electrical impedance tomography with restricted dirichlet-to-neumann map data_ , Inverse Problems 35 (2019), 35005.
  * [28] R.Y. Lay and Q. Li, _Parameter reconstruction for general transport equation_ , SIAM J. Mathematical Analysis 52 (2020), 2734–2758.
  * [29] F. Natterer, _The Mathematics of Computerized Tomography_ , SIAM publications, 2001.
  * [30] R. G. Novikov, _An inversion formula for the attenuated X-ray transformation_ , Ark. Mat. 40 (2002), 145–167.
  * [31] 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 J. Scientific Computing 41 (2019), B929–B952.
  * [32] P. Stefanov and G. Uhlmann, _An inverse source problem in optical molecular imaging_ , Analysis and PDE 1 (2008), 115–126.
  * [33] A. N. Tikhonov, A. V. Goncharsky, V. V. Stepanov, and A. G. Yagola, _Numerical methods for the solution of ill-posed problems_ , Kluwer, London, 1995.
  * [34] M. M. Vajnberg, _Variational method and method of monotone operators in the theory of nonlinear equations_ , Israel Program for Scientific Translations, Jerusalem-London, 1973.