Convexification for the Viscocity Solution for a Coefficient Inverse Problem for the Radiative Transfer Equation

TL;DR

This paper introduces a convexification method for a coefficient inverse problem in radiative transfer, incorporating viscosity solutions and proving stability and convergence, with numerical results showing high efficiency.

Contribution

It is the first to consider viscosity solutions for the boundary value problem in this context and provides stability and convergence analysis using Carleman estimates.

Findings

Proved Lipschitz stability for the boundary value problem.

Established global convergence of the numerical method.

Numerical experiments demonstrate high computational efficiency.

Abstract

A Coefficient Inverse Problem for the radiative transport equation is considered. The globally convergent numerical method, the so-called convexification, is developed. For the first time, the viscosity solution is considered for a boundary value problem for the resulting system of two coupled partial differential equations. A Lipschitz stability estimate is proved for this boundary value problem using a Carleman estimate for the Laplace operator. Next, the global convergence analysis is provided via that Carleman estimate. Results of numerical experiments demonstrate a high computational efficiency of this approach.

Equations

$$0<a<b.$$

$$\Omega =\{x:-B<x_1,...,x_{n-1}<B,a<z<b\},$$

$${\partial }_1\Omega =\left\{x:-B<x_1,...,x_{n-1}<B,z=a\right\},$$

$${\partial }_2\Omega =\left\{x:-B<x_1,...,x_{n-1}<B,z=b\right\},$$

$${\partial }_3\Omega =\left\{x_i=\pm B,z\in \left(a,b\right),i=1,...,n-1\right\}.$$

$${\Phi }_d=\{x_{\alpha }=(\alpha ,0,...,0):\alpha \in [-d,d]\}.$$

$${\Phi }_d\cap \overline{\Omega }=\varnothing .$$

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

$$C_{\sigma }\underset{|x|<\sigma }{\int }\exp \left(\frac{{\left|x\right|}^2}{{\sigma }^2-{\left|x\right|}^2}\right)dx=1.$$

$$f\left(x-x_{\alpha }\right)=0,\forall x\in \overline{\Omega },\forall x_{\alpha }\in {\Phi }_d.$$

$$\overline{B}=\max \left(B,d\right).$$

$$\left\{\begin{array}[]{c}P=\left\{x:\left|x_{1}\right|,...,\left|x_{n-1}\right|<\overline{B},z\in\left(0,b\right)\right\},\\ P_{a}^{+}=P\cap\left\{z>a\right\},\\ P_{a}^{-}=P\diagdown P_{a}^{+}.\end{array}\right.$$

$$\Omega \subseteq P_a^{+}.$$

$$x\in P,\alpha \in \left(-d,d\right).$$

$$\nu (x,\alpha )=\frac{x-x_{\alpha }}{\left|x-x_{\alpha }\right|},x\ne x_{\alpha }.$$

$$\left.\begin{array}[]{c}\nu(x\mathbf{,}\alpha)\cdot\nabla_{x}u(x,\alpha)+a\left(x\right)u(x,\alpha)\\ =\mu_{s}(x)\int\limits_{\Phi_{d}}G(x,\alpha,\beta)u(x,\beta)d\beta+f\left(x-x_{\alpha}\right),\ x\in P,x_{\alpha}\in\Phi_{d},\end{array}\right.$$

$$\displaystyle G(x,\alpha,\beta)\geq 0,\ x\mathbf{\in}\overline{P};\ \alpha,\beta\in\left[-d,d\right],\$$

$$G(x,\alpha ,\beta )\in C^1\left(\overline{P}\times {\left[-d,d\right]}^2\right).$$

$$a\left(x\right)={\mu }_a\left(x\right)+{\mu }_s(x),$$

$${\mu }_a\left(x\right),{\mu }_s(x)\ge 0,{\mu }_a\left(x\right)={\mu }_s(x)=0,x\in P\setminus \Omega ,$$

$${\mu }_a\left(x\right),{\mu }_s(x)\in C\left(\overline{P}\right).$$

$$u(x_{\alpha },\alpha )=0\text{mbox}for{x}_{\alpha }\in {\Phi }_d.$$

$$c\left(x\right)=\exp \left({\int }_{L(x,x_{\alpha })}a(x\left(s\right))ds\right).$$

$$u(x,\alpha )=u_0(x,\alpha )+$$

$$+\frac{1}{c(x,\alpha )}\underset{L(x,x_{\alpha })}{\int }c(x\left(s\right),\alpha ){\mu }_s(x\left(s\right))\left(\underset{{\Phi }_d}{\int }G(x\left(s\right),\alpha ,\beta )u(x\left(s\right),\beta )d\beta \right)ds,$$

$$u_0(x,\alpha )=\frac{1}{c\left(x,\alpha \right)}\underset{L(x,x_{\alpha })}{\int }f\left(x\left(s\right)-x_{\alpha }\right)ds.$$

$$u\left(x,\alpha \right)\ge m>0\text{mbox}for\left(x,\alpha \right)\in {\overline{P}}_a^{+}\times \left[-d,d\right],m=\underset{{\overline{P}}_a^{+}\times \left[-d,d\right]}{\min }{u}_0(x,\alpha ).$$

$$g\left(x,\alpha \right)=u\left(x,\alpha \right),\forall x\in \partial \Omega ,\forall \alpha \in \left(-d,d\right).$$

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

$$\nu (x,\alpha )\cdot {\nabla }_{x}w(x,\alpha )+a\left(x\right)$$

$$=e^{-w\left(x,\alpha \right)}{\mu }_s(x)\underset{{\Phi }_d}{\int }G(x,\alpha ,\beta )e^{w\left(x,\beta \right)}d\beta ,x\in \Omega ,\alpha \in \left(-d,d\right),$$

$$w(x,\alpha ){\mid }_{\partial \Omega }=g_1\left(x,\alpha \right)=\ln g\left(x,\alpha \right).$$

$$\nu (x,\alpha )\cdot {\nabla }_x{w}_{\alpha }(x,\alpha )+{\partial }_{\alpha }\nu (x,\alpha )\cdot {\nabla }_{x}w(x,\alpha )+$$

$$+e^{-w\left(x,\alpha \right)}{w}_{\alpha }\left(x,\alpha \right){\mu }_s(x)\underset{{\Phi }_d}{\int }G(x,\alpha ,\beta )e^{w\left(x,\beta \right)}d\beta -$$

$$-e^{-w\left(x,\alpha \right)}{\mu }_s(x)\underset{{\Phi }_d}{\int }{\partial }_{\alpha }G(x,\alpha ,\beta )e^{w\left(x,\beta \right)}d\beta =0,x\in \Omega ,\alpha \in \left(-d,d\right).$$

$$w_{\alpha }(x,\alpha ){\mid }_{\partial \Omega }=g_2\left(x,\alpha \right)=\frac{g_{\alpha }\left(x,\alpha \right)}{g\left(x,\alpha \right)}.$$

Full text

††: Inverse Problems

, ,

• March 2023

Keywords: Carleman estimate, viscosity solution, radiative transport equation, coefficient inverse problem, Lipschitz stability estimate, convexification, global convergence, numerical studies.

1 Introduction

The radiative transfer equation (RTE) is commonly used in diffusive optics [17]. In the case of light propagation, RTE governs scattering and absorption of photons when they propagate through a diffusive medium, such as, e.g. turbulent atmosphere and biological medium [17]. RTE is also known as the Boltzmann equation [8, 17].

One of the specific applications is in seeing through a turbulent atmosphere. Another attractive application is in optical molecular imaging (OMI) [35] when some optical markers are attached to specific molecules to detect faulty genes. In the single positron emission computed tomography (SPECT) and in photon emission tomography (PET) markers with the high energy are used. Unlike these, markers of OMI emit a relatively low energy near infrared light. Reconstruction of the attenuation coefficient using measurements of the intensity of the emitted light on parts of the human body should lead to the detection and classification of faulty genes. The latter might result in a better diagnostics.

Coefficient Inverse Problems (CIPs) for PDEs are both nonlinear and ill-posed. These two factors cause major challenges in their numerical solutions. The majority of numerical methods for CIPs is based on a very popular procedure of the minimization of least squares cost functionals.In this regard we refer to e.g. [2, 1, 5, 6, 14, 15, 12, 16]. On the other hand, there is no guarantee of the convexity of these functionals, and this might lead to the existence of multiple local minima and ravines, see, e.g. [32].

To avoid the latter, the convexification method was originally proposed in [22, 18] for the case of CIPs for hyperbolic PDEs. The convexification is a numerical development of the idea of [7], which was originally proposed only for the proofs of global uniqueness theorems for multidimensional CIPs and has been explored by many authors since then, see, e.g. [13, 30, 31] as well as [19] for a survey of results as of 2013, [21] for the most recent result, and the book [25, Chapters 2,3]. In [7] the method of Carleman estimates was introduced in the field of CIPs for the first time. Carleman estimates are also the key to the convexification idea.

In the past few years, members of this research group have been working on various applications of the convexification to numerical solutions of various CIPs, see, e.g. [25] for a summary of results as of 2021. In particular, in two most recent publications they have solved by convexification CIPs for two versions of the RTE. [27, 26]. In both these works one obtains first a boundary value problem (BVP) for a nonlinear PDE of the second order, which does not contain the unknown coefficient. Next, the solution of this BVP is represented via a truncated Fourier series with $N>2$ terms with respect to a special orthonormal basis. This basis was originally proposed in [20], also, see [25, section 6.2.3]. As a result, one obtains a new BVP for a system of $N$ coupled nonlinear PDEs of the first order. This BVP is solved via the construction and subsequent minimization of a weighted cost functional. The main element of this functional is the presence of a Carleman Weighted Function (CWF) in it.

The CWF is a function, which is used as a weight in the Carleman estimate for the corresponding PDE operator. The key property of this functional is that it is strictly convex on a convex bounded set of an arbitrary diameter $d>0$ in an appropriate Hilbert space. Since  a smallness condition is not imposed on $d$, then it is natural to call that functional globally strictly convex.

The goal of this work is to numerically solve the above mentioned BVP for that nonlinear PDE via perturbing that PDE operator with the viscosity term $-\varepsilon\Delta,$ where $\varepsilon>0$ is a small number. Then we obtain a BVP for a system of only two coupled nonlinear PDEs of the second order, as opposed to the above case of $N>2$ PDEs . First, we prove the Lipschitz stability estimate for this BVP using a Carleman estimate for the Laplace operator. Next, we construct a globally strictly convex weighted cost functional for the latter BVP. The key element of this functional is the presence of the Carleman Weight Function in it. This function is involved as the weight in the above mentioned Carleman estimate. The proof of the global strict convexity of that functional is the main part of our global convergence analysis. This analysis ends up with the proof of the global convergence of the gradient descent method of the minimization of that functional. Finally, we conduct exhaustive numerical studies, which demonstrate a good performance of our technique. Previously the viscosity solution of the Hamilton-Jacobi equation was numerically computed in [29] by a different version of the convexification method.

Publications [27, 26] represent the first numerical solutions of CIPs for both RTE [27] and its Riemannian version [26]. Previous numerical results were obtained only for the inverse source problems for RTE in [9, 10, 11, 33] and for RRTE in [24]. However, inverse source problems are linear ones, unlike CIPs. There are also some uniqueness and stability results for CIPs for RTE. Since this paper is not concerned with such results, then we refer here only to a limited number of publications on the latter topic [4, 13, 30, 31].

In section 2 we state both forward and inverse problems. To solve the inverse problem, we obtain in in section 3 a boundary value problem for a $2\times 2$ system of nonlinear PDEs with viscosity terms. We prove the Lipschitz stability estimate for this boundary value problem in section 4. In section 5 we construct the convexified Tikhonov-like functional and conduct the convergence analysis for it. Numerical experiments are presented in section 6.

2 Statements of Forward and Inverse Problems

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

[TABLE]

Let $\Omega\subset\mathbb{R}^{n}$ be a rectangular prism with its boundary $\partial\Omega=\partial_{1}\Omega\cup\partial_{2}\Omega\cup\partial_{3}\Omega,$ where

[TABLE]

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

[TABLE]

Hence, $\Phi_{d}$ is a part of the $x_{1}-$axis. By (2.1), (2.2) and (2.6)

[TABLE]

Let $\sigma>0$ be a sufficiently small number. To simplify the presentation, avoid working with singularities. Hence, we consider the following function instead of $\delta\left(x\right)$:

[TABLE]

The constant $C_{\sigma}$ is chosen here such that

[TABLE]

Hence, the function $f\left(x-x_{\alpha}\right)=f\left(x_{1}-\alpha,x_{2},...,x_{n},z\right)\in C^{\infty}\left(\mathbb{R}^{n}\right)$ can be considered as the source function for the source $x_{\alpha}\in\Phi_{d}.$ Let the number $\sigma$ be so small that

[TABLE]

Let

[TABLE]

Introduce the domain $P\subset\mathbb{R}^{n}$ as well as two of its subdomains $P_{a}^{+}$,$P_{a}^{-}\subset P,$

[TABLE]

By (2.2), (2.11) and (2.12)

[TABLE]

Below

[TABLE]

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

[TABLE]

Let $u(x,\alpha)$ be the steady-state radiance at the point $x$ generated by the source function $f\left(x-x_{\alpha}\right).$ The function $u(x,\alpha)$ satisfies the stationary RTE [17]:

[TABLE]

where “$\cdot$” denotes the scalar product in $\mathbb{R}^{n}$, see (2.14). The kernel $G(x,\alpha,\beta)$ of the integral operator in (2.16) satisfies [17, 33]:

[TABLE]

In (2.16),

[TABLE]

where $\mu_{a}\left(x\right)$ and $\mu_{s}(x)$ are the absorption and scattering coefficients respectively. The function $a\left(x\right)$ is the attenuation coefficient. We assume that

[TABLE]

Forward Problem. Let conditions (2.1)-(2.21) hold. Find the function $u(x,\alpha)$ \in C^{1}\Big{(}\overline{P}\times\left[-d,d\right]\Big{)} satisfying equation (2.16) and the initial condition

[TABLE]

Denote

[TABLE]

The following existence and uniqueness theorem for the Forward Problem was proven in [27], and a similar theorem was proven in [26] for the Riemannian analog of RTE:

Theorem 2.1 [27]. Assume that conditions (2.1)-( 2.15) and (2.17)-(2.21) hold. Then *there exists unique solution $u\left(x,\alpha\right)\in C^{1}\left(\overline{P}\times\left[-d,d\right]\right)$ of equation (2.16) with the initial condition (2.22), and the function $u\left(x,\alpha\right)$ has the following form for *$x\in P,\alpha\in\left(-d,d\right):$

[TABLE]

* Furthermore, the following inequality holds:*

[TABLE]

Coefficient Inverse Problem (CIP). Let conditions (2.1 )-(2.21) hold. Let the function $u\left(x,\alpha\right)\in C^{1}\left(\overline{G}\times\left[-d,d\right]\right)$ be the solution of the Forward Problem as in Theorem 2.1. Assume that the attenuation coefficient $a\left(x\right)$ in (2.16) is unknown. Assume that the function $g\left(x,\alpha\right)$ is known,

[TABLE]

*Find the function $a\left(x\right).$ *

3 Boundary Value Problem for a $2\times 2$ System of PDEs With

Viscosity Terms

3.1 Preliminaries

It follows from (2.13), (2.14) and (2.25) that we can introduce a new function $w(x,\alpha)$,

[TABLE]

Substituting (3.1) in (2.16) and (2.26), we obtain

[TABLE]

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

[TABLE]

By (3.3) the boundary condition for $w_{\alpha}(x,\alpha)$ is

[TABLE]

To work with the viscosity solution, we need to figure out the following boundary conditions, see (2.4):

[TABLE]

It follows from (2.6) and (2.15) that $\nu(x\mathbf{,}\alpha)=(\nu_{1}(x\mathbf{,}\alpha)$, $\nu_{2}(x\mathbf{,}\alpha)$, $\cdots$, $\nu_{n-1}(x\mathbf{,}\alpha)$, $\nu_{n}(x\mathbf{,}\alpha))$, where

[TABLE]

Hence,

[TABLE]

By (2.19) and (2.20) $a\left(x\right)=0$ for $x\in\partial_{2}\Omega.$ Hence, (3.2), (3.3) and (3.8)-(3.13) imply:

[TABLE]

3.2 Viscosity solution

Denote

[TABLE]

Based on (3.3), (3.5), (3.7), (3.15)-(3.17), consider the following BVP for the system of viscosity equations with a small parameter $\varepsilon>0$:

[TABLE]

And also

[TABLE]

Therefore, we have obtained the BVP (3.19)-(3.25) with respect to the pair of functions $\left(p_{\varepsilon},q_{\varepsilon}\right)\left(x,\alpha\right)$. By (3.24) and (3.25) this BVP has an overdetermination in the Neumann boundary conditions at $z=b.$ We focus below on the solution of this BVP. Suppose that we have computed a solution $\left(p_{\varepsilon,\mbox{comp}},q_{\varepsilon,\mbox{comp}}\right)\left(x,\alpha\right)$ of BVP (3.19)-(3.25). Then we use (3.2) and (3.18) to compute the target coefficient $a\left(x\right),$

[TABLE]

Remark 3.1. As soon as we got BVP (3.19)-(3.25) for the pair of functions $\left(p_{\varepsilon},q_{\varepsilon}\right),$ we do not require anymore that

[TABLE]

* as in (3.18). In other words, we solve this BVP for a slightly broader class of vector functions $\left(p_{\varepsilon},q_{\varepsilon}\right).$ Nevertheless, it follows from the uniqueness claim of Theorem 4.1 that the solution of this “broader” BVP, if it exists, is still such that (3.27) holds. *

4 Lipschitz Stability Estimate for BVP (3.19)-(3.25)

Introduce the space $H_{2}^{1}\left(\Omega\right)\times L_{2,2}\left(-d,d\right)$ as

[TABLE]

Theorem 4.1 (Lipschitz stability and uniqueness). *Assume that conditions (2.1)-(2.7) and (2.18)-(2.21) hold. Suppose that there exists two pairs of functions *

$\left(p_{\varepsilon,1},q_{\varepsilon,1}\right)\left(x,\alpha\right),\left(p_{\varepsilon,2},q_{\varepsilon,2}\right)\left(x,\alpha\right)\in H^{2}\left(\Omega\right)\times C\left[-d,d\right]$* satisfying equations (3.19) and (3.23) and such that*

[TABLE]

Let

[TABLE]

* Then there exists a constant*

[TABLE]

depending only on listed parameters such that the following Lipschitz stability estimates hold:

[TABLE]

*Let $a_{\varepsilon,1}\left(x\right)$ and $a_{\varepsilon,2}\left(x\right)$ * functions computed via the right hand side of ( 3.26), in which $p_{\varepsilon,\mbox{comp}}(x,\alpha)$ is replaced with $p_{\varepsilon,1}(x,\alpha)$ and $p_{\varepsilon,2}(x,\alpha)$ respectively. Then the following Lipschitz stability estimate holds:

[TABLE]

In particular, suppose that, in addition to (4.2)

[TABLE]

Then (4.4) and (4.5) imply that $\emph{for\ }x\in\Omega,\alpha\in\left(-d,d\right)$

[TABLE]

Below $C_{1}>0$ denotes different constants depending on parameters listed in (4.3). Uniqueness of the BVP (3.19)-(3.25) obviously follows from (4.6).

4.1 Carleman estimate

Prior the proof of Theorem 4.1, we need to prove a Carleman estimate for the operator $\Delta.$  Denote

[TABLE]

Theorem 4.2. Assume that conditions (2.2)-(2.5) hold. There exists a constant $C=C\left(\Omega\right)>0$ and a sufficiently large number $\lambda_{0}=\lambda_{0}\left(\Omega\right)>1$, *both depending only on the domain *$\Omega,$ *such that the following Carleman estimate holds for all  and for all *$\lambda\geq\lambda_{0}:$

[TABLE]

Proof. Everywhere below $C=C\left(\Omega\right)>0$ denotes different constants depending only on the domain $\Omega.$ We assume first that

[TABLE]

Introduce a new function $v\left(x\right),$

[TABLE]

Hence,

[TABLE]

Hence,

[TABLE]

Step 1. Estimate from the below the following term in the second line of (4.12):

[TABLE]

Thus, moving from the function $v$ to the function $u$ via (4.10), we obtain

[TABLE]

Step 2. Estimate from the below the following term in the second line of (4.12):

[TABLE]

Thus, taking

[TABLE]

we obtain

[TABLE]

Summing up (4.13) and (4.15) and taking into account (4.12), we obtain

[TABLE]

**Step 3. **Consider

[TABLE]

Thus,

[TABLE]

Multiply (4.17) by $5\lambda$ and sum up with (4.16). We obtain

[TABLE]

Integrating this over the domain $\Omega$ and using Gauss formula, (2.2)-(2.5), (4.7) and (4.9), we obtain that there exists a sufficiently large number $\lambda_{0}=\lambda_{0}\left(\Omega\right)\geq\lambda_{01}$ and a number $C=C\left(\Omega\right)>0$ such that

[TABLE]

By Cauchy-Schwarz inequality $5\lambda\Delta u\cdot ue^{2\lambda z^{2}}\leq 2.5\left(\Delta u\right)^{2}e^{2\lambda z^{2}}+2.5\lambda^{2}u^{2}e^{2\lambda z^{2}}.$ Substituting this in (4.18) and using density arguments, we obtain the target estimate (4.8). $\ \square$

4.2 Proof of Theorem 4.1

Denote

[TABLE]

Then by (4.2) and (4.7)

[TABLE]

Let the function $f:\mathbb{R}^{k}\rightarrow\mathbb{R}$, $f\in C^{1}\left(\mathbb{R}^{k}\right),k\geq 1.$ It is well known that the following formula is valid

[TABLE]

where “$\circ"$ is the scalar product in $\mathbb{R}^{k}$ and the vector function $\widetilde{f}\left(y_{1},y_{2}\right)$ is such that

[TABLE]

Consider the differences $L_{1}\left(p_{\varepsilon,1},q_{\varepsilon 1}\right)-L_{1}\left(p_{\varepsilon,2},q_{\varepsilon 2}\right)$ and $L_{2}\left(p_{\varepsilon,1},q_{\varepsilon 1}\right)-L_{2}\left(p_{\varepsilon,2},q_{\varepsilon 2}\right).$ Then, using (3.18)-(3.23) and (4.19)-(4.22), we obtain two integral differential inequalities for $x\in\Omega,\alpha\in\left(-d,d\right)$

[TABLE]

Here and below in this proof $C_{2}>0$ denotes different positive constant depending on the same parameters as ones listed in (4.3), except of $\varepsilon.$

Square both sides of each of equations (4.23), (4.24), then multiply by $e^{2\lambda z^{2}}$and then integrate over the domain $\Omega,$ assuming that $\lambda\geq\lambda_{0},$where $\lambda_{0}$ is the parameter of Theorem 4.2. And then sum up two resulting inequalities. Using Cauchy-Schwarz inequality, we obtain

[TABLE]

Recalling (4.20) and applying the Carleman estimate (4.8) to the second line of (4.27), we obtain for all $\lambda\geq\lambda_{0}$

[TABLE]

Choose

[TABLE]

such that $C\varepsilon^{2}\lambda_{1}\geq 2C_{2}$ and then set in (4.28) $\lambda\geq\lambda_{1}.$ We obtain with the constant $C_{1}>0$ as in (4.3)

[TABLE]

Integrate both sides of (4.29) with respect to $\alpha\in\left(-d,d\right).$ Then choose $\lambda_{2}\geq\lambda_{1}$ such that $\lambda_{2}^{3}\geq 2C_{1}.$ Then set $\lambda=\lambda_{2}.$ We obtain

[TABLE]

which is equivalent with (4.4). Estimate (4.5) obviously follows from (4.4).  $\square$

5 Convergence Analysis for the Convexification for BVP (3.19)-(3.25)

Consider an integer $k$,

[TABLE]

where $\left[(n+1)/2\right]$ is the largest integer, which does not exceed $(n+1)/2.$ Then by embedding theorem and (5.1)

[TABLE]

Introduce the space $S,$

[TABLE]

Let $R>0$ be an arbitrary number. Consider the set of 2D vector functions $B\left(g_{1},g_{2},g_{3},g_{4},R\right)$ defined as

[TABLE]

To solve BVP (3.19)-(3.25), we consider

Minimization Problem. Minimize the functional $J_{\lambda,\beta}\left(p_{\varepsilon},q_{\varepsilon}\right)$ on the set $B(g_{1}$, $g_{2},g_{3},g_{4},R),$ where

[TABLE]

In (5.5) $\gamma\in\left(0,1\right)$ is the regularization parameter and the multiplier $e^{-2\lambda b^{2}}$ is introduced to balance two terms in the right hand side of (5.5). Indeed,

[TABLE]

Let $S$ be the space defined in (5.3). We introduce the subspace $S_{0}\subset S$ as:

[TABLE]

5.1 Global strict convexity

Below in section 5 $C_{3}>0$ denotes different constants, all of which depend on the following parameters:

[TABLE]

Recall the definition (4.1)* *of the space $H_{2}^{1}\left(\Omega\right)\times L_{2,2}\left(-d,d\right).$

Theorem 5.1 (global strict convexity). Denote

[TABLE]

* For any $\lambda>0$ functional (5.5) has the Fré chet derivative $J_{\lambda,\gamma}^{\prime}\left(v\right)\in S_{0}$ at every point $v\in\overline{B\left(g_{1},g_{2},g_{3},g_{4},R\right)}.$ This derivative is Lipschitz continuous on $\overline{B\left(g_{1},g_{2},g_{3},g_{4},R\right)},$ *  i.e. there exists a constant $D>0$ such that

[TABLE]

Let $\lambda_{0}\geq 1$ be the number of Theorem 4.2. There exists a sufficiently large number

[TABLE]

such that for every $\lambda\geq\lambda_{3}$ the functional $J_{\lambda,\gamma}\left(v\right)$ is strictly convex on the set $\overline{B\left(g_{1},g_{2},g_{3},g_{4},R\right)},$ i.e. for all $v_{1},v_{2}\in\overline{B\left(g_{1},g_{2},g_{3},g_{4},R\right)}$  the following inequality holds:

[TABLE]

Furthermore, for every $\lambda\geq\lambda_{3}$ there exists unique minimizer $v_{\min,\lambda}$ of the functional $J_{\lambda,\gamma}\left(v\right)$ on the set $\overline{B\left(g_{1},g_{2},g_{3},g_{4},R\right)}$ and the following inequality holds:

[TABLE]

Proof. Let $v_{1}=\left(p_{1},q_{1}\right)$ and $v_{2}=\left(p_{2},q_{2}\right)$ be two arbitrary points of the set $\overline{B\left(g_{1},g_{2},g_{3},g_{4},R\right)}.$ Denote

[TABLE]

By (5.2), (5.3), (5.4), (5.7), (5.8) and (5.14)

[TABLE]

By (3.19) and (5.14)

[TABLE]

Consider now the nonlinear term in the third line of (5.17). By Taylor formula

[TABLE]

where $\xi$ is a point between $p_{1}\left(x,\alpha\right)$ and $p_{1}\left(x,\alpha\right)+h_{1}\left(x,\alpha\right)$ and $\eta$ is a point between $p_{1}\left(x,\beta\right)$ and $p_{1}\left(x,\beta\right)+h_{1}\left(x,\beta\right).$ Hence, using (5.16), (5.18) and (5.19), we obtain

[TABLE]

And also

[TABLE]

Thus, by (5.17), (5.20) and (5.21)

[TABLE]

where $L_{1,\mbox{lin}}\left(p_{1},q_{1,}h\right)$ depends linearly on $h$ and $L_{1,\mbox{nonlin}}\left(p_{1},q_{1,}h\right)$ depends nonlinearly on $h.$ Also, for $x\in\Omega,\alpha\in\left(-d,d\right)$

[TABLE]

Similarly,

[TABLE]

Using (5.5) and (5.22)-(5.29) and recalling (5.9), we obtain

[TABLE]

where $\left[,\right]$ is the scalar product in $S$.

Consider the sum of second and third lines of (5.30),

[TABLE]

It follows from (5.7) and (5.15) that we can consider $I\left(p_{1},q_{1},h\right):S_{0}\rightarrow\mathbb{R}$ as a linear functional of $h\in S_{0}.$ It follows from (5.3), (5.23), (5.24), (5.27) and (5.28) that $I\left(p_{1},q_{1},h\right)$ is a bounded functional. Therefore, by Riesz theorem there exists unique point $\widehat{I}\left(p_{1},q_{1}\right)\in S_{0}$ such that

[TABLE]

It follows from (5.9), (5.14) and (5.23)-(5.32) that

[TABLE]

Therefore $\widehat{I}\left(p_{1},q_{1}\right)=J_{\lambda,\gamma}^{\prime}\left(p_{1},q_{1}\right)$ is the Fréchet derivative of the functional $J_{\lambda,\beta}\left(p,q\right)$ at the point $\left(p_{1},q_{1}\right).$ The Lipschitz continuity property (5.10) is proven similarly with the proof of Theorem 3.1 of [3]. Therefore, we omit this proof here.

Thus, (5.30) can be rewritten as

[TABLE]

Let $RHS$ denotes the right hand side of (5.33). Then using (5.2), (5.4), (5.16), (5.23)-(5.29) and Cauchy-Schwarz inequality, we obtain the following estimate from the below:

[TABLE]

It follows from (4.7), (5.7) and (5.15) that we can apply Carleman estimate (4.8) of Theorem 4.2 to the right hand side of (5.34), and the second line of (4.8) should be zero in this case. Let $\lambda_{0}=\lambda_{0}\left(\Omega\right)\geq 1$ be the number, which was found in Theorem 4.2. We obtain for all $\lambda\geq\lambda_{0}$

[TABLE]

Choose the number $\lambda_{3}\geq\lambda_{0}$ depending on the same parameters as the ones listed in (5.11) and such that $C\varepsilon^{2}\lambda_{3}\geq 2C_{3}.$ Then we obtain

[TABLE]

This, (5.6) and (5.33) imply

[TABLE]

which is equivalent with (5.12).

Given $\lambda\geq\lambda_{3},$ the existence and uniqueness of the minimizer $v_{\min,\lambda}$ of the functional $J_{\lambda,\gamma}\left(v\right)$ on the set $\overline{B\left(g_{1},g_{2},g_{3},g_{4},R\right)}$ as well as inequality (5.13) follow immediately from a combination of either Lemma 2.1 with Theorem 2.1 of [3] or, equivalently, Lemma 5.2.1 and Theorem 5.2.1 of [25]. $\ \ \square$

5.2 The accuracy of the minimizer

In this section we estimate the distance between the minimizer $v_{\min,\lambda}$ which was found in Theorem 5.1, and the exact solution $v^{\ast}$ with the noiseless data of BVP (3.19)-(3.25). In accordance with the regularization theory [34], we assume that there exists a solution $v_{\varepsilon}^{\ast}\left(x,\alpha\right)=\left(p_{\varepsilon}^{\ast},q_{\varepsilon}^{\ast}\right)\left(x,\alpha\right)$ of BVP (3.19)-(3.25) with the noiseless boundary data $g_{1}^{\ast},g_{2}^{\ast},g_{3}^{\ast},g_{4}^{\ast}$ in (3.24), (3.25). By Theorem 4.1 this solution is unique. However, in applications the data (3.24), (3.25) are always given with a noise. Let a small number $\delta\in\left(0,1\right)$ be the level of the noise in the data (3.24), (3.25). We assume that there exist two vector functions $F\left(x,\alpha\right)=\left(F_{1},F_{2}\right)\left(x,\alpha\right)\in S$ and $F^{\ast}\left(x,\alpha\right)=\left(F_{1}^{\ast},F_{2}^{\ast}\right)\left(x,\alpha\right)\in S$ such that

[TABLE]

And we also assume that

[TABLE]

For every vector function $v=\left(p,q\right)\in B\left(g_{1},g_{2},g_{3},g_{4},R\right)$ consider the difference

[TABLE]

Also, denote

[TABLE]

Using (5.4), (5.7) and (5.36)-(5.40), we obtain

[TABLE]

Theorem 5.2. *Let conditions (5.36)-(5.38) and notations (5.39),(5.40) hold. Let the regularization parameter *$\gamma$ be **

[TABLE]

*Let $\lambda_{3}=\lambda_{3}\left(\Omega,d,\varepsilon,R\right)\geq\lambda_{0}\geq 1$ be the number in (5.11). Consider the number *$\lambda_{4},$

[TABLE]

Let $\lambda=\lambda_{4}$ and let $v_{\min,\lambda_{4}}=\left(p_{\min,\lambda_{4}},q_{\min,\lambda_{4}}\right)\in\overline{B\left(g_{1},g_{2},g_{3},g_{4},R\right)}$ be the minimizer of the functional $J_{\lambda_{4},\gamma}\left(v\right)$   on the set $\overline{B\left(g_{1},g_{2},g_{3},g_{4},R\right)}.$  Let $a_{\min,\lambda_{4}}\left(x\right)$ and $a_{\varepsilon}^{\ast}\left(x\right)$ be functions $a\left(x\right),$ which are constructed from $v_{\min,\lambda_{4}}=\left(p_{\min,\lambda_{4}},q_{\min,\lambda_{4}}\right)$ and $v_{\varepsilon}^{\ast}=\left(p_{\varepsilon}^{\ast},q_{\varepsilon}^{\ast}\right)$ respectively via the right hand side of formula (5.11), in which $p_{\varepsilon,\mbox{comp}}$ is replaced with $p_{\min,\lambda}$ and $p_{\varepsilon}^{\ast}$ respectively. Then the following accuracy estimates hold:

[TABLE]

Proof. Let $B_{0}\left(2R\right)$ be the set defined in (5.41). Consider a new functional $I_{\lambda_{4},\gamma}:B_{0}\left(2R\right)\rightarrow\mathbb{R}$ defined as

[TABLE]

Then Theorem 5.2 is applicable to this functional. Let $V_{\min,\lambda_{3}}$ be the minimizer of $I_{\lambda_{4},\gamma}\left(\widetilde{v}\right)$ on the set $\overline{B_{0}\left(2R\right)},$

[TABLE]

Since by (5.41) both vector functions $V_{\min,\lambda_{4}},\widetilde{v}_{\varepsilon}^{\ast}\in S_{0}\left(2R\right),$ then by an obvious analog of (5.12)

[TABLE]

By (5.13)

[TABLE]

Hence, (5.48) implies

[TABLE]

Taking into account dependencies (5.8) and (5.43), we obtain

[TABLE]

By (5.40) and (5.46)

[TABLE]

By (5.5)

[TABLE]

Since $v_{\varepsilon}^{\ast}$ is the exact solution of BVP (3.19)-(3.25), then $L_{1}\left(v^{\ast}\right)=L_{2}\left(v^{\ast}\right)=0.$ Hence, (5.42) and (5.51) imply

[TABLE]

Next, by (5.38), (5.50)-(5.52) and Cauchy-Schwarz inequality

[TABLE]

Hence, using (5.49), we obtain

[TABLE]

Denote

[TABLE]

We have

[TABLE]

Hence,

[TABLE]

Using (5.36), (5.54) and the triangle inequality, we obtain

[TABLE]

Therefore,

[TABLE]

On the other hand, let $v_{\min\lambda_{4}}$ be the minimizer of the functional $J_{\lambda_{3},\gamma}\left(v\right)$ on the set $\overline{B\left(g_{1},g_{2},g_{3},g_{4},R\right)},$ which is found in Theorem 5.1,

[TABLE]

Let $\widetilde{v}_{\min,\lambda_{4}}=v_{\min\lambda_{4}}-F.$ Then $\widetilde{v}_{\min,\lambda_{4}}\in\overline{B_{0}\left(2R\right)}$ and by (5.47)

[TABLE]

However, since by (5.53) and (5.55) $V_{\min,\lambda_{4}}+F\in\overline{B\left(g_{1},g_{2},g_{3},g_{4},R\right)},$ then by (5.56) we should have

[TABLE]

Since the minimizer is unique, then (5.57) and (5.58) imply that $V_{\min,\lambda_{4}}+F=\widehat{V}_{\min,\lambda_{4}}=v_{\min,\lambda_{4}}.$ Thus, (5.54) implies (5.44). Estimate (5.45) obviously follows from (5.11) and (5.44). $\square$

5.3 Global convergence of the gradient descent method

Let $\lambda=\lambda_{4},$ where $\lambda_{4}$ is defined in (5.41). Consider two sets

[TABLE]

We assume now that

[TABLE]

Consider the gradient descent method of the minimization of the functional $J_{\lambda_{4},\gamma}.$ Consider an arbitrary point

[TABLE]

Let $\rho>0$ be a small number. Define the sequence of the gradient descent method as:

[TABLE]

Note that since by Theorem 5.1 $J_{\lambda_{4},\gamma}^{\prime}\left(v_{m-1}\right)\in S_{0}$ for all $m\geq 1,$ then it follows from (5.7) and (5.61) that boundary conditions (3.24), (3.25) are kept the same for all vector functions $v_{m}.$ Theorem 5.3 follows immediately from a combination of Theorem 5.2 with Theorem 6 of [23].

Theorem 5.3. Let $\lambda=\lambda_{4}$ and let conditions of Theorem 5.2, (5.59) and (5.60) hold. Then there exists a sufficiently small number $\rho_{0}\in\left(0,1\right)$ such that for every $\rho\in\left(0,\rho_{0}\right)$ there exists a number $\theta=\theta\left(\rho\right)\in\left(0,1\right)$ such that the sequence $\left\{v_{m}\right\}_{m=0}^{\infty}\subset\widetilde{B}$ and the following convergence estimates for the gradient descent method (5.60), (5.61) hold:

[TABLE]

*where the function $a_{m}\left(x\right)$ is constructed from the vector function $v_{m}=\left(p_{m},q_{m}\right)$ by the right hand side of formula (5.11), in which $p_{\varepsilon,\mbox{comp}}$ is replaced with *$p_{m}.$

Remarks 5.1:

*Since smallness assumptions are not imposed on the number $R>0$ * and since $v_{0}\in\widetilde{B}$ is an arbitrary point, then Theorem 5.3 guarantees the global convergence of the gradient descent method (5.60), (5.61). 2. 2.

Even though the requirement of our theory is that the parameter * $\lambda$ of the Carleman Weight Function $e^{2\lambda z^{2}}$  should be sufficiently large, we have observed in computational experiments of section 6 that $\lambda=5$ is sufficient, which is the same as in two previous publications of this group [27, 26]. Similar observations about reasonable values of $\lambda\in\left[1,3\right]$ were made in other publications about the convexification method [25, 23, 29, 28]. On the other hand, we observe in numerical experiments of section 6 that too large values of *$\lambda=20$ *do not work well, see Figure 3. This is because the Carleman Weight Function changes well too rapidly for *$\lambda=20.$ 3. 3.

Conceptually, the considerations of item 2 are similar with asymptotic theories. Indeed, an asymptotic theory usually claims that if a certain parameter $X$ is sufficiently large, then a certain formula * $Y$ is valid with a good accuracy. However, in a practical computation, which always has a specific ranges of parameters, only numerical experiments can establish reasonable values of *$X,$ * for which $Y$ is valid with a good accuracy. Besides, it is well known that too large values of *$X$ *often do not work well for numerical studies. *

6 Numerical Studies

We have conducted numerical studies in the 2D case. In our numerical testing the domain $\Omega$ and the line $\Phi_{d}$ in (2.2)-(2.6) are:

[TABLE]

We take $\sigma=0.05$ in the function $f_{\sigma}(x)$ in (2.8) and (2.9). As to the kernel $G(x,\alpha,\beta)$ in (2.16), we choose the 2-dimensional Henyey-Greenstein function [17]:

[TABLE]

Here $c_{g}$ represents the ballistic with $c_{g}=0$ and isotropic scattering with $c_{g}=1$ [9, 10, 11], respectively. In this paper, we choose $c_{g}=1/2$.

We have chosen the absorption and scattering coefficients $\mu_{a}(x),\mu_{s}(x)$ in (2.19)-(2.21) as:

[TABLE]

In the numerical tests below, we take $c_{a}=5,10,15,20,30$, and the inclusions with the shape of the letters ‘A’, ‘$\Omega$’ and ‘SZ’.

Remark 6.1. We have intentionally chosen abnormalities with the shapes of letters to demonstrate that our reconstruction technique works well for truly hard cases of non-convex abnormalities containing voids.

By (2.19), (6.3) and (6.6), we have

[TABLE]

Following (6.7), we define the computed inclusion/background contrast as:

[TABLE]

6.1 Data generation

To generate the boundary data (2.26) and then (3.24), (3.25) for our CIP, we have solved the Forward Problem posed in section 2. Using Theorem 2.1, we have solved this problem numerically via the solution of the linear integral equation (2.23) with the condition (2.24).

Consider the partition of the domain $\Omega$ and the line $\Phi_{d}$ in (6.1) with the given mesh sizes $h_{x_{1}},h_{z},h_{\alpha}$:

[TABLE]

Then the set of the discrete points are given as

[TABLE]

We have used the grid step sizes $h_{x_{1}}=h_{z}=h_{\alpha}=1/40.$ To obtain the numerical solution $u(x^{h},\alpha^{h})$ of the Forward Problem, we have solved the corresponding linear algebraic system by the Matlab backslash operator ‘$\backslash$’. This way we have generated the boundary data (2.26). Then, using considerations of subsection 3.1, we have obtained the boundary data (3.24), (3.25).

6.2 Numerical results for the inverse problem

For the inverse problem, we set $h_{x_{1}}=h_{z}=h_{\alpha}=1/20$ to generate the discrete points in (6.9)-(6.13). The discrete form of functional (5.5) is

[TABLE]

where the pair $\left(p^{h},q^{h}\right)(x^{h},\alpha^{h})$ is the pair of functions $\left(p,q\right)\left(x,\alpha\right)$ written on the discrete grid and $L_{1}^{h}\left(p^{h},q^{h}\right)$ and $L_{2}^{h}\left(p^{h},q^{h}\right)$ are operators $L_{1}$ and $L_{2}$ in (3.19) and (3.23), in which differential operators are written in finite differences and integrals are written in discrete forms using the trapezoidal rule.

To numerically solve the Minimization Problem posed in section 5, we have minimized functional (6.14) with respect to the values of discrete functions $p_{\varepsilon}^{h}(x^{h},\alpha^{h})$, $q_{\varepsilon}^{h}(x^{h},\alpha^{h})$ at grid points. The Dirichlet boundary conditions in (3.24) and (3.25) are given as

[TABLE]

By the finite difference method, the Neumann boundary conditions in (3.24) and (3.25) are given as

[TABLE]

We have adopted the Matlab’s built-in optimization toolbox fmincon to minimize the function $J_{\lambda,\gamma}\left(p_{\varepsilon}^{h},q_{\varepsilon}^{h}\right)$ in (6.14) with the boundary conditions (6.15) and (6.16). Here, (6.15) and (6.16) are the constraint conditions used in each iteration of fmincon to ensure that the functions $p_{\varepsilon}^{h}(x^{h},\alpha^{h})$, $q_{\varepsilon}^{h}(x^{h},\alpha^{h})$ at every iteration satisfy the boundary conditions (3.24) and (3.25). The iterations of ** fmincon** were stopped at the iteration number $m$ at which

[TABLE]

see Figure 1.

To solve the minimization problem, we need to provide the starting point $p_{0}^{h},q_{0}^{h}$ for iterations. With the boundary conditions (3.24) and (3.25), for every $\alpha\in[-d,d]$, we have the value of functions $p,q$ on the boundary $\partial\Omega$ via functions $g_{1},g_{2}$. Using the linear interpolations of boundary conditions $g_{1},g_{2}$ with respect to $x_{1}-$direction and $z-$direction, the initial guess $\left(p_{0},q_{0}\right)$ for the pair of functions $\left(p,q\right)$ in the domain $\Omega$ is:

[TABLE]

Then the starting point for the minimization of functional (6.18) is $p_{0}^{h}=p_{0}(x^{h},\alpha^{h})$, $q_{0}^{h}=q_{0}(x^{h},\alpha^{h})$. Even though the first guess does not satisfy required Neumann boundary conditions in (3.24) and (3.25), still all follow up iterations of fmincon satisfy both required boundary conditions: Dirichlet and Neumann, by the constraint conditions (6.15) and (6.16).

We introduce the random noise in the boundary data $g(x,\alpha)$ in (2.26) on the boundary $\partial\Omega$ as follows:

[TABLE]

where $\zeta_{x}$ is the uniformly distributed random variable in the interval $[0,1]$ depending on the point $x\in\partial\Omega$ with $\delta=0.03$ and $\delta=0.05$, which correspond respectively to $3\%$ and $5\%$ noise level. Hence, the random noise is also introduced in boundary conditions $g_{1}\left(x,\alpha\right),g_{2}\left(x,\alpha\right),g_{3}\left(x,\alpha\right)$ and $g_{4}\left(x,\alpha\right)$ in (3.24), (3.25). These functions are defined via $g\left(x,\alpha\right)$ using (3.3), (3.7), (3.15) and (3.16). We now explain how did we differentiate the noisy data for $g\left(x,\alpha\right)$ with respect to $\alpha$ in (3.7) and (3.16). The observation data $g(x,\alpha)$ in (2.26) at the boundary $\partial\Omega$ is generated by the source function $f(x-x_{\alpha})$, whose position $x_{\alpha}$ is determined by the value $\alpha$ in (2.6). Then, for each given $\alpha$, we obtain the corresponding observation data $g(x,\alpha)$ generated by the source function $f(x-x_{\alpha})$ as well as the the sample of the random variable $\zeta_{x}$. Since the samples of the random variable $\zeta_{x}$ for each $\alpha_{k}$ in (6.11) are different, then we use the finite difference method to calculate numerically the derivative of the noisy data $g(x,\alpha)$ with respect to $\alpha$ with the above mentioned grid step size $h_{\alpha}=1/40.$ Results of the differentiation were good enough, and we did not observe instabilities.

Test 1. We consider the coefficient $a(x),$ which corresponds to $\mu_{a}(x)$ in (6.6) with $c_{a}=5$ inside of the letter ‘$A$’. The goal of this test is to find the optimal values of the parameters $\varepsilon$ and $\lambda$ for the minimization problem. Noise (6.19) in the data is not added.

We set $\gamma=0.001$, $\lambda=5$ and perform the numerical tests with different values of $\varepsilon$. The results are displayed in Figure 2. The reconstruction of the header of the letter ‘$A$’ is not good for $\varepsilon\in\left[0.05,0.2\right]$. The reconstruction quality improves when $\varepsilon$ decreases while ranging from 0.2 to 0.01. While $\varepsilon$ varies from 0.01 to 0.001, the difference between the reconstructions is very small. On the other hand, when we choose $\varepsilon=0.0001,$ the reconstruction quality of the top of letter ‘$A$’ becomes worse. In conclusion, although $\varepsilon$ should be small enough, but not too small. Thus, we choose $\varepsilon=0.01$ as an optimal value and use this one in all other tests.

Now we want to select an optimal value of the parameter $\lambda.$ We take $\gamma=0.001,\varepsilon=0.01$ and test values $\lambda=0,2,3,4,5,6,8,20$. Results are presented on Figure 3. The parameter $\lambda$ cannot be neither too small nor too large. The reconstructions are unsatisfactory for $\lambda=0,2,20$. On the other hand, the reconstructions become better when $\lambda$ ranges from 3 to 5, and they are stabilized for $\lambda=5,6,8$. Thus, we choose $\lambda=5$ as the optimal value, see items 2 and 3 of Remarks 5.1 for a relevant discussion.

In summary, we use in the tests below

[TABLE]

Test 2. We take the same values of parameters as listed in (6.20). We consider the coefficient $a(x)$ corresponding to $\mu_{a}(x)$ in (6.6) with $c_{a}=10,15,20,30$ inside of the letter ‘$A$’. Hence, the inclusion/background contrasts in (6.7) are respectively $3:1$, $4:1$, $5:1$ and $7:1$. Noise (6.19) in the data is not added. The results are displayed in Figure 4. The reconstruction quality is good for these four cases, although it slightly deteriorates at $c=20$ and $c=30$. The computed inclusion/background contrasts (6.8) are accurate.

Test 3. We use the same values of parameters as ones in (6.20). We consider the coefficient $a(x)$ corresponding to $\mu_{a}(x)$ in (6.6) with $c_{a}=5$ inside of the letter ‘$\Omega$’. Noise (6.19) in the data is not added. The result is displayed in Figure 5. The reconstruction is quite accurate.

Test 4. We consider the coefficient $a(x)$ corresponding to $\mu_{a}(x)$ in (6.6) with $c_{a}=5$ inside of two letters ‘SZ’, which are two letters in the name of the city (Shenzhen) were the second and the third authors reside. Noise (6.19) in the data is not added. Results are exhibited in Figure 6. The reconstruction is worse than the one for the case of the single letter ‘$\Omega$’ in Figure 5. Nevertheless, the reconstruction is still good and the computed inclusion/background contrasts in (6.8) are both accurate in these two letters.

Test 5. We now consider the noisy data, as in (6.19), with $\delta=0.03$ and $\delta=0.05,$ i.e. with 3% and 5% noise level. We reconstruct the coefficient $a(x)$ with the shape of the letters ‘A’ and ‘$\Omega$’ corresponding to $\mu_{a}(x)$ in (6.6) with $c_{a}=5$ inside of two letters. The results are displayed in Figure 7. In all these four cases, reconstructions of shapes of inclusions and the inclusion/background contrasts in (6.8) are accurate.