Travel time tomography with formally determined incomplete data in 3D

TL;DR

This paper introduces a new globally convergent numerical method for 3D travel time tomography with incomplete data, providing stability estimates and a convergent algorithm.

Contribution

It develops the first globally convergent method with stability estimates for 3D travel time tomography with incomplete data.

Findings

Lipschitz stability estimate derived, ensuring uniqueness.

A weighted convex functional is constructed for the inverse problem.

A gradient projection method is proven to converge globally as data noise diminishes.

Abstract

For the first time, a globally convergent numerical method is developed and Lipschitz stability estimate is obtained for the challenging problem of travel time tomography in 3D for formally determined incomplete data. The semidiscrete case is considered meaning that finite differences are involved with respect to two out of three variables. First, Lipschitz stability estimate is derived, which implies uniqueness. Next, a weighted globally strictly convex Tikhonov-like functional is constructed using a Carleman-like weight function for a Volterra integral operator. The gradient projection method is constructed to minimize this functional. It is proven that this method converges globally to the exact solution if the noise in the data tends to zero.

Figures

Equations

$$\Omega =\left\{\mathbf{x}=\left(x,y,z\right):x,y\in \left(0,1\right),z\in \left(A,A+\sigma \right)\right\}.$$

$$B_A=\left\{\mathbf{x}=\left(x,y,z\right):x,y\in \left(0,1\right),z=A\right\},$$

$$B_{A+\sigma }=\left\{\mathbf{x}=\left(x,y,z\right):x,y\in \left(0,1\right),z=A+\sigma \right\},$$

$$\Gamma =\partial \Omega ╲\left(B_A\cup B_{A+\sigma }\right).$$

$$m\left(\mathbf{x}\right)\ge m_0,\mathbf{x}\in {\mathbb{R}}^3,$$

$$m\left(\mathbf{x}\right)=1,\mathbf{x}\in \left\{z<A\right\},$$

$$m\in C^2\left({\mathbb{R}}^3\right),$$

$$m_z\left(\mathbf{x}\right)\ge 0\text{, }\mathbf{x}\in \overline{\Omega }.$$

$$d\tau =\sqrt{m(\mathbf{x})}\left|d\mathbf{x}\right|,|d\mathbf{x}|=\sqrt{(dx{)}^2+(dy{)}^2+(dz{)}^2}.$$

$$\tau \left(\mathbf{x},{\mathbf{x}}_0\right)=\int {\text{displaylimits}}_{\Gamma \left(\mathbf{x},{\mathbf{x}}_0\right)}\sqrt{m\left(\mathbf{x}\left(s\right)\right)}ds,$$

$$L=\left\{\mathbf{x}=\left(x,y,z\right):x=\alpha \in \left(0,1\right),y=1/2,z=0\right\}.$$

$${\tau }_x^2+{\tau }_y^2+{\tau }_z^2=m\left(\mathbf{x}\right),\mathbf{x}\in {\mathbb{R}}^3,$$

$$\tau \left(\mathbf{x},\alpha \right)=O\left(\left|\mathbf{x}-{\mathbf{x}}_{\alpha }\right|\right),\text{a}\text{s }\left|\mathbf{x}-{\mathbf{x}}_{\alpha }\right|\to 0.$$

$$\sum _{i,j=1}^3\frac{{\partial }^2\ln n(\mathbf{x})}{\partial x_i\partial x_j}{\xi }_i{\xi }_j\ge 0,\forall \xi \in {\mathbb{R}}^3,\forall \mathbf{x}=\left(x_1,x_2,x_3\right)\in \overline{\Omega }.$$

$$\tau \left(\mathbf{x},\alpha \right)=f\left(\mathbf{x},\alpha \right),\forall \mathbf{x}\in \Gamma \cup B_{A+\sigma },\forall \alpha \in \left(0,1\right).$$

$$q\left(\alpha \right)=\sum {\text{displaylimits}}_{n=1}^N{q}_n{\Psi }_n(\alpha ),q_n=\int {\text{displaylimits}}_0^{1}q\left(\alpha \right){\Psi }_n(\alpha )d\alpha .$$

$${\Psi }_n\left(\alpha \right)=\sum {\text{displaylimits}}_{j=1}^N{x}_{nj}{\xi }_j\left(\alpha \right).$$

$$q\left(\alpha \right)=\sum {\text{displaylimits}}_{n=1}^N{q}_n\sum {\text{displaylimits}}_{j=1}^N{x}_{nj}{\xi }_j\left(\alpha \right)=\sum {\text{displaylimits}}_{j=1}^N\left(\sum {\text{displaylimits}}_{n=1}^N{x}_{nj}{q}_n\right){\xi }_j\left(\alpha \right)=\sum {\text{displaylimits}}_{j=1}^N{\tilde{q}}_j{\xi }_j\left(\alpha \right).$$

$${\tau }_z^2\left(\mathbf{x},\alpha \right)\ge \frac{A^2}{A^2+2},\forall \mathbf{x}\in \overline{\Omega },\forall \alpha \in \left[0,1\right].$$

$${\tau }_z\left(\mathbf{x},\alpha \right)>0,\forall \mathbf{x}\in \overline{\Omega }\cup \left\{z\in \left(0,A\right]\right\},\forall \alpha \in \left[0,1\right].$$

$${\tau }_z\left(\mathbf{x},\alpha \right)\ge \frac{A}{\sqrt{A^2+2}},\forall \mathbf{x}\in \overline{\Omega },\forall \alpha \in \left[0,1\right].$$

$$p={\tau }_x(x,y,z,\alpha ),q={\tau }_x(x,y,z,\alpha ),r={\tau }_z(x,y,z,\alpha ).$$

$$\frac{dx}{ds}=\frac{p}{m},\frac{dy}{ds}=\frac{q}{m},\frac{dz}{ds}=\frac{r}{m},$$

$$\frac{dp}{ds}=\frac{m_x}{2m},\frac{dq}{ds}=\frac{m_y}{2m},\frac{dr}{ds}=\frac{m_z}{2m},$$

$$\frac{d\tau }{ds}=\frac{\partial \tau }{\partial x}\frac{dx}{ds}+\frac{\partial \tau }{\partial y}\frac{dy}{ds}+\frac{\partial \tau }{\partial z}\frac{dz}{ds}=p\frac{p}{m}+q\frac{q}{m}+r\frac{r}{m}=1.$$

$$\tau \left(x\left(0\right),y\left(0\right),z\left(0\right),\alpha \right)=0\text{ for }s=0.$$

$$\frac{dx}{dz}=\frac{p}{r},\frac{dy}{dz}=\frac{q}{r},\frac{dp}{dz}=\frac{m_x}{2r},\frac{dq}{dz}=\frac{m_y}{2r},\frac{d{r}^2}{dz}=m_z,\frac{d\tau }{dz}=\frac{m}{r}.$$

$$x{|}_{z=0}=\alpha ,y{|}_{z=0}=1/2,p{|}_{z=0}=p_0,q{|}_{z=0}=q_0,r{|}_{z=0}=r_0,\tau {\mid }_{z=0}=0.$$

$$\tau \left(x,y,z,\alpha \right)=\sqrt{{\left(x-\alpha \right)}^2+{\left(y-1/2\right)}^2+z^2}\text{ for }\left(x,y,z\right)\in \left(0,1\right)\times \left(0,1\right)\times \left(0,A\right).$$

$$r_0=\sqrt{1-p_0^2-q_0^2}\ge 0,$$

Full text

Key words. inverse kinematic problem; Carleman-like estimate for the Volterra operator; Lipschitz stability; convexification; globally convergent numerical method

AMS Classification 35R30

1 Introduction

We construct a globally convergent numerical method for the challenging 3D travel time tomography problem (TTTP) with formally determined incomplete data. The TTTP was first considered by Herglotz [7] in 1905 and then by Wiechert and Zoeppritz [34] in 1907 in the 1D case due to a geophysical application, also, see a detailed description of the 1D case in [24]. However, globally convergent numerical methods for the 3D TTTP with formally determined data were not developed so far. By the definition of, e.g. [1, 2], a numerical method for a nonlinear problem is called globally convergent if there exists a theorem claiming that this method delivers at least one point in a sufficiently small neighborhood of the correct solution without relying on any advanced knowledge of this neighborhood, i.e. a good first guess for the solution is not required.* *

In this paragraph, we indicate those ideas for the TTTP, which are presented here for the first time. More details about the latter statement, including some references, are given in this section below. Since we develop a numerical method, we are allowed to work here with an *approximate mathematical model. *We study the case when the data for the TTTP are both formally determined and incomplete. The TTTP is considered in the semidiscrete form, i.e.  we consider the practically important case of finite differences with respect to two out of three variables, also see Remark 6.1 in section 6. The Lipschitz stability estimate is obtained, which implies uniqueness. Next, the exact solution is constructed via introducing a sequence, which converges to that solution starting from an arbitrary point of a certain bounded set, provided that the noise in the data tends to zero. Since no restrictions are imposed on the size $R$ of that set, then this is a globally convergent numerical method. To construct that method, we introduce a weighted Tikhonov-like functional for the TTTP, which is strictly convex on that bounded set, i.e. we use the so-called “convexification” procedure. However, while in all previous versions of the convexification the so-called Carleman Weight Function (CWF) was applied to differential operators, in the current paper the CWF is applied to a Volterra-like nonlinear integral operator. Our method does not use sophisticated geometrical properties to construct the above sequence. Rather, we straightforwardly minimize the above mentioned globally strictly convex Tikhonov-like functional.

The TTTP is the problem of the recovery of the spatially distributed speed of acoustic waves from first times of arrival of those waves. Another well known term for the TTTP is the “inverse kinematic problem of seismic” [24]. Waves are originated by some point sources located on the boundary of the domain of interest. First times of arrival are recorded on a number of detectors located on that boundary. It is well known that the TTTP is actually a nonlinear problem of the integral geometry, see, e.g. [24]*. The *TTTP has important applications in geophysics [7, 24, 33, 34]. In addition, it was established in [15] that the TTTP arises in the phaseless inverse problem of scattering of electromagnetic waves at high frequencies. The specific TTTP considered here has potential applications in geophysics, checking the bulky baggage in airports, search for possible defects inside the walls, etc..

The “formally determined data” means that the number $m$ of free variables in the data equals to the number $n$ of free variables in the unknown coefficient, $m=n$. If, however, $m>n$, then the data are overdetermined. In our case $n=m=3$. In the 2D case all previously known results for the TTTP work with formally determined data with $m=n=2$. Only complete data for the TTTP were considered in the past. In the 3D case those data were overdetermined with $m=4>n=3$. Complete data for the TTTP are generated by the point source running along the entire surface of the boundary of the domain of interest. Unlike these, our point source runs along an interval of a straight line located outside of the domain of interest. This means that the data are incomplete. The sole purpose of Figure 1 is to illustrate this for a simple case when geodesics are straight lines.

Our approximate mathematical model consists of two items, see Remarks 6.2, 6.3 in the end of section 6 for more details. First, we consider a semidiscrete model. This means that partial derivatives with respect to two out of three variables are written in finite differences. The author believes that this might be valuable for a possible future numerical implementation of the method of this paper. Second, we assume that a certain function generated by the solution of the eikonal equation can be represented as a truncated Fourier series with respect to a special orthonormal basis in the $L_{2}\left(0,1\right)$ space. That basis was first constructed in [16]. Functions of that basis depend only on the position of the source. The number $N\geq 1$ of terms of this series is not allowed to tend to the infinity. We note that such assumptions quite often take place in the theory of ill-posed problems, and corresponding computational results are not affected by these assumptions, see, e.g. [18, 19, 20], where the same basis was used to obtain good quality numerical results.

The TTTP is about the reconstruction of the right hand side of the nonlinear eikonal Partial Differential Equation (PDE) from boundary measurements of the solution of this equation. It is well known that conventional least squares cost functionals for nonlinear inverse problems are non convex. The non convexity, in turn causes the well known phenomenon of multiple local minima and ravines of that functional, see, e.g. [27] for a quite convincing numerical example. Therefore, convergence of any minimization algorithm to the correct solution is in question, unless its starting point is chosen to be in a sufficiently small neighborhood of that solution, which is the so-called small perturbation approach. However, it is unclear in many applications how to actually obtain that sufficiently small neighborhood.

To avoid the phenomenon of local minima, we apply the so-called convexification method. More precisely, we construct a weighted globally strictly convex Tikhonov-like functional for our TTTP. The key element of this functional is the presence in it of the weight function which looks similar with the CWF.

The convexification was first proposed in [11, 12]. However, those initial works lacked some important analytical results, which would allow numerical studies. Such results were recently obtained in [1]. Thus, recent publications [17, 18, 19] contain both the theory and numerical results for the convexification method for some coefficient inverse problems for PDEs, although not for the TTTP. We prove the global convergence to the correct solution of the gradient projection method of the minimization of our functional.

As to the numerical methods for the TTTP in the $n-$D case, $n=2,3$, such a method for the 2D TTTP was published in [28]. Another numerical approach in 3D was published in [35]. Both publications [28, 35] use, at a certain stage, the minimization of a least squares cost functional. Since the convexity of those functionals of [28, 35] is not proven, then the problem of local minima is not addressed there. In both publications [28, 35] complete data are used, and these data are overdetermined ones in the 3D case of [35].

The first global Lipschitz stability and uniqueness result for the TTTP was obtained by Mukhometov in the 2D case [21]. Next, this result was extended in [4, 22, 24] to the $n-$D case, $n\geq 3$. We also refer to the related work [23] for the 2D case. In all these references, the data are complete and the assumption of the regularity of geodesic lines is used. In addition, more recently the question of uniqueness in the 3D case when geodesic lines are not necessarily regular ones was considered in [30]. In the 2D case of [21, 23, 28] the data are formally determined. However, they are overdetermined in the $n-$D case with $n\geq 3$ [4, 22, 24, 30, 35].

As to the Carleman estimates, for the first time they were introduced in the field of Inverse Problems in [5] with the goal of proofs of global uniqueness and stability results for coefficient inverse problems. This idea became quite popular since then with many publications of many authors. To shorten the citation list, we refer here only to, e.g. books [2, 3, 13], the survey [14] and references cited therein. In addition to uniqueness and stability, various modifications of the idea of [5] are currently used for the convexification, see corresponding comments and references above.

All functions considered below are real valued ones. In section 2 we state the TTTP. In section 3 we introduce a special orthonormal basis. In section 4 we estimate from the below a derivative of the solution of the eikonal equation. In section 5 we derive a boundary value problem for a system of coupled nonlinear integral differential equations. In section 6 we rewrite that system in a semidiscrete form. In section 7,  we establish Lipschitz stability and uniqueness. In section 8 we construct the above mentioned globally strictly convex functional and formulate corresponding theorems. These theorems are proven in section 9.

2 Statement of the problem

Below $\mathbf{x}=\left(x,y,z\right)\in\mathbb{R}^{3}.$ Let numbers $A,\sigma=const.>0$. Define the rectangular prism $\Omega\subset\mathbb{R}^{3}$ as

[TABLE]

Denote parts of the boundary $\partial\Omega$ as

[TABLE]

[TABLE]

[TABLE]

Let $n\left(\mathbf{x}\right)$ be the refractive index of the medium at the point $\mathbf{x.}$ Hence, $c\left(\mathbf{x}\right)=1/n\left(\mathbf{x}\right)$ is the sound speed. Denote $m\left(\mathbf{x}\right)=n^{2}\left(\mathbf{x}\right).$ Let the number $m_{0}>0$ be given. We impose the following assumptions on the function $m\left(\mathbf{x}\right):$

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Remark 2.1. We note that the monotonicity condition (8) is not an overly restrictive one. Indeed, it can also be found in the end of section 2 of chapter 3 of the book [24]: see formulas (3.24) and (3.24*′**) there. Also, an analogous monotonicity condition was actually imposed in the 1D case in the originating classical works of Herglotz and Wiechert and Zoeppritz [7, 34], see section 3 of chapter 3 of [24] for a description of their method. We also refer to Figures 5 and 10 in [33] for some geophysical information. *

The function $m\left(\mathbf{x}\right)$ generates the Riemannian metric

[TABLE]

The travel time from the point $\mathbf{x}_{0}\in\mathbb{R}^{3}$ (source) to the point $\mathbf{x}\in\mathbb{R}^{3}$ (receiver) is [24]

[TABLE]

where $\Gamma\left(\mathbf{x},\mathbf{x_{0}}\right)$ is the geodesic line connecting points $\mathbf{x}$ and $\mathbf{x}_{0}$ and $ds$ is the euclidean arc length. We assume that the source $\mathbf{x}_{0}$ runs along an interval $L$ of a straight line located in the plane $\left\{z=0\right\},$

[TABLE]

Hence, $\mathbf{x}_{0}=\mathbf{x}_{\alpha}=\left(\alpha,1/2,0\right),\alpha\in\left(0,1\right).$ Let $\tau\left(\mathbf{x,}\alpha\right)$ be the travel time between points $\mathbf{x}$ and $\mathbf{x}_{\alpha}=\left(\alpha,1/2,0\right).$ Thus, we denote $\tau\left(\mathbf{x,}\alpha\right)=\tau\left(\mathbf{x},\mathbf{x}_{\alpha}\right).$ We denote $\Gamma\left(\mathbf{x},\alpha\right)$ the geodesic line connecting points $\mathbf{x}$ and $\mathbf{x}_{\alpha}.$ It is well known [24] that the function $\tau\left(\mathbf{x},\alpha\right)$ satisfies eikonal equation as the function of $\mathbf{x}$,

[TABLE]

[TABLE]

Everywhere below we assume without further mentioning that the following property holds:

Regularity of Geodesic Lines. *The function *$\tau\left(\mathbf{x},\alpha\right)\in C^{2}\left(\mathbb{R}^{3}\times\left[0,1\right]\right).$ *For any pair of points $\left(\mathbf{x},\mathbf{x}_{\alpha}\right)$ *$\in\overline{\Omega}\times L$ there exists unique geodesic line $\Gamma\left(\mathbf{x},\alpha\right)$ connecting these two points and $\Gamma\left(\mathbf{x},\alpha\right)\cap B_{A}\neq\varnothing$. *In addition, if any geodesic line, which starts at a point *$\mathbf{x}_{\alpha}\in L,$ *intersects *$B_{A},$ *then it intersects it at a single point. Also, it does not go “downwards” in the $z-$direction, but instead intersects $\partial\Omega\diagdown B_{A}$ at another single point, see (1)-(4). In addition, after intersecting $\partial\Omega\diagdown B_{A},$ this line does not “come back” in the domain $\Omega$ but rather goes away from this domain. In other words, this line is not reflected back from any point of its intersection with *$\partial\Omega$.

The following sufficient condition of the regularity of geodesic lines was derived in [25]:

[TABLE]

Travel Time Tomography Problem (TTTP). *Suppose that the function *$m\left(\mathbf{x}\right)$ satisfies conditions (5)-(8). Assume that the following function $f\left(\mathbf{x},a\right)$ is given:

[TABLE]

Determine the function $m\left(\mathbf{x}\right)$ for $\mathbf{x}\in\Omega.$

In other words, by (1)-(4) and (13) the data for the travel time are given for all sources running along the line interval $L$ defined in (11) and for the part $\Gamma\cup B_{A+\sigma}$ of the boundary $\partial\Omega.$ Hence, the data (13) are both formally determined and incomplete.

3 A Special Orthonormal Basis

We now reproduce a special orthonormal basis $\left\{\Psi_{n}\left(\alpha\right)\right\}_{n=0}^{\infty}$ in $L_{2}\left(0,1\right),$ which was constructed in [16]. This basis has the following two properties:

$\Psi_{n}\in C^{1}\left[0,1\right],$ $\forall n=1,2,\dots$ 2. 2.

Let $\left\{,\right\}$ be the scalar product in $L_{2}\left(0,1\right)$. Denote $a_{mn}=\left\{\Psi_{n}^{\prime},\Psi_{m}\right\}.$ Then the matrix $M_{N}=\left(a_{mn}\right)_{m,n=0}^{N}$ should be invertible for any $N=1,2,\dots$

Note that if one would use either the basis of standard orthonormal polynomials or the trigonometric basis, then the first derivative of its first element would be identically zero. Hence, the matrix $M_{N}$ would not be invertible in this case.

We now describe the basis of [16]. Consider the set of functions $\left\{\xi_{n}\left(\alpha\right)\right\}_{n=0}^{\infty}=\{\left(\alpha+a\right)^{n}e^{\alpha}\}_{n=0}^{\infty}$, where $a=const.>0.$ This is a set of linearly independent functions. Besides, this set is complete in the space $L_{2}(0,1)$. After applying the Gram-Schmidt orthonormalization procedure to this set, we obtain the orthonormal basis $\{\Psi_{n}\left(\alpha\right)\}_{n=1}^{\infty}$ in $L_{2}(0,1)$. In fact, the function $\Psi_{n}(\alpha)$ has the form $\Psi_{n}(\alpha)=P_{n}(\alpha+a)e^{\alpha}$, $\forall n\geq 0,$ where $P_{n}$ is a polynomial of the degree $n$. Thus, functions $\Psi_{n}(\alpha)$ are polynomials orthonormal in $L_{2}(0,1)$ with the weight $e^{2\alpha}.$ The matrix $M_{N}$ is invertible since its elements $a_{mn}=\left(\Psi_{n}^{\prime},\Psi_{m}\right)$ are such that $a_{mn}=1$ if $m=n$ and $a_{mn}=0$ if $n<m$.

Consider the function $q\left(\alpha\right)$ in the following form:

[TABLE]

Below we need to impose such a sufficient condition on the vector of coefficients

$q^{N}=\left(q_{0},...,q_{N-1}\right)^{T}$ in the Fourier expansion (14) which would guarantee that the function $q\left(\alpha\right)$ is positive for all $\alpha\in\left[0,1\right].$ Consider vector functions $\xi^{N}\left(\alpha\right)=\left(\xi_{1},...,\xi_{N}\right)^{T}\left(\alpha\right)$ and

$\Psi^{N}(\alpha)=\left(\Psi_{1},...,\Psi_{N-1}\right)^{T}\left(\alpha\right).$ The desired condition is given in Lemma 3.1.

Lemma 3.1. *Let the $N\times N$ matrix $X_{N}$ transforms the vector function $\xi^{N}\left(\alpha\right)$ in the vector function $\Psi^{N}\left(\alpha\right)$ via the Gram-Schmidt orthonormalization procedure, i.e. $X_{N}\cdot\xi^{N}\left(\alpha\right)=\Psi^{N}\left(\alpha\right).$ Let the matrix $X_{N}^{T}$ be the transpose of $X_{N}$. Let $X_{N}^{T}q^{N}=\widetilde{q}^{N}=\left(\widetilde{q}_{0},...,\widetilde{q}_{N-1}\right)^{T}.$ Suppose that all numbers $\widetilde{q}_{0},...,\widetilde{q}_{N-1}>0.$ Then in (14) the function $q\left(\alpha\right)>0$ for all *$\alpha\in\left[0,1\right].$

Proof. It follows from the Gram-Schmidt procedure that elements of the matrix $X_{N}$ are independent on $\alpha.$ Let the raw number $n$ of the matrix $X_{N}$ be $\left(x_{n1},x_{n2},...,x_{nN}\right),n=1,...,N.$ Then

[TABLE]

Hence, by (14)

[TABLE]

Since $\xi_{j}\left(\alpha\right)=\left(\alpha+a\right)^{j}e^{\alpha}>0,$ then (15) implies that $q\left(\alpha\right)>0$ for $\alpha\in\left[0,1\right].$ $\square$

4 Estimate of $\tau_{z}^{2}\left(\mathbf{x},\mathbf{\alpha}\right)$ From the Below

Lemma 4.1. Assume that conditions (5)-(8) hold. Then

[TABLE]

Also,

[TABLE]

Thus,

[TABLE]

Proof. Note that having proven (16) is not enough for our technique: we need to know the sign of the function $\tau_{z}\left(\mathbf{x},\alpha\right),$ i.e. (17), in section 5 (more precisely, in (31)) where we consider the square root of* *$\tau_{z}^{2}\left(\mathbf{x},\alpha\right).$ Denote

[TABLE]

The following equations for geodesic lines can be found on page 66 of [24]:

[TABLE]

[TABLE]

where $s$ is a parameter. Using (19), we obtain* *for $\tau=\tau\left(x\left(s\right),y\left(s\right),z\left(s\right),\alpha\right)$ along a geodesic line

[TABLE]

Set

[TABLE]

Then (22) implies:

$\tau\left(x\left(s\right),y\left(s\right),z\left(s\right),\alpha\right)=s$. Hence, the parameter $s$ coincides with the travel time. In particular, we now rewrite equations (20), (21) in a different form: to have derivatives with respect to $z$ rather than with respect to $s$. Hence, we obtain from (20) and (21):

[TABLE]

[TABLE]

where $\alpha\in\left(0,1\right)$ and $p_{0}$, $q_{0}$ are some given numbers such that $p_{0}^{2}+q_{0}^{2}\leq 1.$ The latter inequality follows from (12) and the fact that by (6) $m\left(x,y,0\right)=1$. Also, by (12) $r_{0}=\pm\sqrt{1-p_{0}^{2}-q_{0}^{2}}.$ To prove that we should take “+” sign in the latter formula, we note that

[TABLE]

Hence, $\tau_{z}=r=z/\tau>0$ for $z\in\left(0,A\right).$ Hence,

[TABLE]

[TABLE]

Suppose that the geodesic line defined by (24) and (25) intersects the part $B_{A}$ of the boundary $\partial\Omega.$ Then the condition of the regularity of geodesic lines implies that there exists a single number $z_{0}=z_{0}\left(p_{0},q_{0},\alpha\right)\in\left(A,A+\sigma\right]$ such that the point $\left(x\left(z_{0}\right),y\left(z_{0}\right),z_{0}\right)\in\partial\Omega\diagdown B_{A}$ and for all numbers $z\in\left(A,z_{0}\right)$ all points $\left(x\left(z\right),y\left(z\right),z\right)$ of that geodesic line belong to $\Omega$. Since by (24) $dr^{2}/dz=m_{z},$ then, using (24) and (28), we obtain

[TABLE]

[TABLE]

which establishes (16). To prove (17), we notice that it follows from (5), (8), (19), the last equation (21) and (23) that

[TABLE]

Estimates (17) and (18) follow immediately from (16), (26) and (29). $\square$

5 A Boundary Value Problem for a System of Nonlinear Coupled Integro

Differential Equations

5.1 A nonlinear integro differential equation

Denote

[TABLE]

By (17)

[TABLE]

Hence, (13) and (31) imply that for all $\alpha\in\left(0,1\right)$

[TABLE]

[TABLE]

[TABLE]

Substituting (30)-(34) in the eikonal equation (12), we obtain the following equation for $\left(x,y,z\right)\in\Omega,\alpha\in\left(0,1\right)$:

[TABLE]

[TABLE]

Differentiating (35) with respect to $\alpha$ and using $\partial_{\alpha}m\left(x,y,z\right)\equiv 0,$ we obtain for $\left(x,y,z\right)\in\Omega,\alpha\in\left(0,1\right)$

[TABLE]

[TABLE]

5.2 Boundary value problem for a system of coupled integro

differential equations

Using (26) and (30), denote for $\left(x,y,\alpha\right)\in\left[0,1\right]^{3}$

[TABLE]

We seek the function $u\left(\mathbf{x},\alpha\right)$ in the form

[TABLE]

[TABLE]

where the function $v\left(\mathbf{x},\alpha\right)$ is unknown for $\mathbf{x}\in\Omega,\alpha\in\left(0,1\right).$ Recall that the part $\Gamma$ of the boundary of the domain $\Omega$ is defined in (4). We need to obtain zero boundary condition at $\Gamma$ for a function associated with the function $v$. To do this, we assume first that there exists a function $g\left(\mathbf{x},\alpha\right)\in H^{1}\left(\Omega\right)$ for every $\alpha\in\left[0,1\right]$ such that

[TABLE]

[TABLE]

Introduce the function $w\left(\mathbf{x},\alpha\right),$

[TABLE]

Then (38)-(42) imply that

[TABLE]

[TABLE]

[TABLE]

We assume that both functions $w\left(\mathbf{x},\alpha\right)$ and $g\left(\mathbf{x},\alpha\right)$ have the form of the truncated Fourier series with respect to the orthonormal basis $\{\Psi_{n}\left(\alpha\right)\},$

[TABLE]

[TABLE]

Here, coefficients $w_{n}\left(\mathbf{x}\right)$ are unknown and coefficients $g_{n}\left(\mathbf{x}\right)$ are known. And similarly for $w_{x},w_{y},w_{x\alpha},w_{y\alpha}$ and the same derivatives of the function $g$. Furthermore, we assume that these functions, being substituted in equation (36), give us zero in its right hand side for $\mathbf{x}\in\Omega,\alpha\in\left(0,1\right).$* *By (44)-(46)

[TABLE]

Denote

[TABLE]

[TABLE]

Remark 5.1. Note that we do not impose an analog of conditions (46), (47) on the function $u(x,y,A,\alpha).$ *The number $N$ of Fourier harmonics in (46), (47) should be chosen in numerical studies. For example, analogs of the series (46) were considered for four different inverse problems in [18, 19, 20, 29]. It was established numerically that for an inverse problem of [18] the optimal $N=8$, for the inverse problem of [19] the optimal $N=3$, for the inverse problem of [20] the optimal $N=15$ and the optimal $N=12$ for the inverse problem of [29]. *

Let

[TABLE]

[TABLE]

Keeping in mind (48) and (49), we define the spaces $C_{N}^{1}\left(\overline{\Omega}\right)$, $C_{N,0}^{1}\left(\overline{\Omega}\right)$ of $N-$D vector functions $W\left(\mathbf{x}\right)$ defined in (49) as

[TABLE]

[TABLE]

Keeping in mind (37)-(52), substitute functions $w$, $g$ and their corresponding first derivatives with respect to $x,y,\alpha$ in equation (36). Next, multiply the resulting equation sequentially by functions $\Psi_{n}\left(\alpha\right),n=1,...,N$ and integrate with respect to $\alpha\in\left(0,1\right)$. Then multiply both sides of obtained system of nonlinear integral differential equations by the matrix $M_{N}^{-1},$ where the matrix $M_{N}$ was introduced in section 3. We obtain

[TABLE]

[TABLE]

where $P$ is the $N-$D vector function,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Thus, we have obtained the desired boundary value problem (53)-(56) for the system of nonlinear coupled integral differential equations. Below we focus on this problem.

5.3 The positivity of the function $\left(u_{0}+w+g\right)\left(\mathbf{x},\alpha\right)$

It follows from (36) and (43) that we need the function $\left(u_{0}+w+g\right)\left(\mathbf{x},\alpha\right)$ to be positive. We discuss this issue in the current section.

Using (16), (30), (37), (43) and (49), define the set of functions $K$ as:

[TABLE]

where $\left(\mathbf{x},\alpha\right)\in\overline{\Omega}\times\left[0,1\right].$ Then by (37) and (43)

[TABLE]

To obtain a sufficient condition guaranteeing (57) in terms of the vector function $W$, we formulate Lemma 5.1.

Lemma 5.1.* Let (46) and (47) hold. Consider the vector function $v^{N}\left(\mathbf{x}\right)=\left(w_{1}+g_{1},w_{2}+g_{2},...,w_{N}+g_{N}\right)^{T}\left(\mathbf{x}\right).$ Let $X_{N}$ be the $N\times N$ matrix of Lemma 3.1. Consider the vector function $\widetilde{v}^{N}\left(\mathbf{x}\right)=X_{N}^{T}\cdot v^{N}\left(\mathbf{x}\right).$ Let $\widetilde{v}^{N}\left(\mathbf{x}\right)=\left(v_{1},...,v_{N}\right)^{T}\left(\mathbf{x}\right).$ Suppose that all functions $v_{n}\left(\mathbf{x}\right)>0$ in $\overline{\Omega}.$ Then the function $w\in K$ and, therefore, (58) holds. Also, the set $K$ is convex*.

Proof. The first part of this lemma follows immediately from Lemma 3.1. We now prove the convexity of the set $K.$ Suppose that two functions $w^{\left(1\right)},w^{\left(2\right)}\in K.$ Let the number $\theta\in\left(0,1\right).$ Then by (57)

[TABLE]

[TABLE]

Summing up these two inequalities, we obtain

[TABLE]

5.4 Applying the multidimensional analog of Taylor formula

We specify in this section how the classical multidimensional analog of Taylor formula [32]* *can be applied to the right hand side of equation (53). Let $R>0$ be an arbitrary number. Denote

[TABLE]

It follows from Lemma 5.1 and (59) that $K\left(R\right)$ is a convex set.

Lemma 5.2. * Let $W^{\left(1\right)}$,$W^{\left(2\right)}\in K\left(R\right)$,let *$G^{\left(1\right)},G^{\left(2\right)}\emph{be}$ the vector functions (50) and $F^{\left(1\right)},F^{\left(2\right)}$ *be the vector functions in (52). Based on (46) and (49), denote $\widetilde{w}\left(\mathbf{x},\alpha\right)=w^{\left(1\right)}\left(\mathbf{x},\alpha\right)-w^{\left(2\right)}\left(\mathbf{x},\alpha\right)$, *$\widetilde{w}_{n}\left(\mathbf{x}\right)=w_{n}^{\left(1\right)}\left(\mathbf{x}\right)-w_{n}^{\left(2\right)}\left(\mathbf{x}\right).$ *Similarly, denote $\widetilde{g}\left(\mathbf{x},\alpha\right)=g^{\left(1\right)}\left(\mathbf{x},\alpha\right)-g^{\left(2\right)}\left(\mathbf{x},\alpha\right),$ where $g^{\left(k\right)}\left(\mathbf{x},\alpha\right)$ corresponds to the vector function *$G^{\left(k\right)}\left(\mathbf{x}\right),k=1,2$ via (47), (50). Also, denote

[TABLE]

* And let functions *$f^{\left(1\right)}=f^{\left(1\right)}\left(x,y,A+\sigma\right)$ and $\widetilde{f}=\widetilde{f}\left(x,y,A+\sigma\right)$ correspond to the vector functions $F^{\left(1\right)}\left(\mathbf{x}\right)$ and $\widetilde{F}\left(\mathbf{x}\right)$ respectively via (51), (52). *Then the following form of the multidimensional Taylor formula is valid: *

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

*where *$\mathbf{x}\in\Omega,i=1,2$ *in $T_{3}$ means that this matrix depends on both vector functions $\left(W^{\left(1\right)},W_{x}^{\left(1\right)},W_{y}^{\left(1\right)},G_{x}^{\left(1\right)},G_{y}^{\left(1\right)},F_{x}^{\left(1\right)},F_{y}^{\left(1\right)}\right)$ and *$\left(W^{\left(2\right)},W_{x}^{\left(2\right)},W_{y}^{\left(2\right)},G_{x}^{\left(2\right)},G_{y}^{\left(2\right)},F_{x}^{\left(2\right)},F_{y}^{\left(2\right)}\right).$ Also, all elements of matrices $T_{j},S_{k},j=0,1,2,3;k=1,2,3,4$ are continuous functions of their variables. Furthermore, the following estimates hold

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where the number

[TABLE]

* depends only on listed parameters. Estimates (61)-(64) mean estimates for each element of corresponding matrices. In terms of the integration with respect to $\alpha,$ matrices *$T_{j},S_{j},j=0,1,2,3$ depend on the integrals of the form

[TABLE]

*where $n,k,m=1,...,N$ and $\mu_{jnkm}$ are continuous functions of their variables. *

Proof. Below $C_{1}>0$ denotes different constants depending only on parameters listed in (64). For $\mathbf{x}\in\Omega$ and $\alpha\in\left(0,1\right)$ consider the functions $\xi\left(\mathbf{x},\alpha\right)$ and $\eta\left(\mathbf{x},\alpha\right)$ defined as

[TABLE]

where functions $w,g$ have the forms (46), (47). Then $\xi\left(w^{\left(2\right)},\mathbf{x},\alpha\right)=\xi\left(w^{\left(1\right)}+\widetilde{w},\mathbf{x},\alpha\right)$ and $\eta\left(w^{\left(2\right)},\mathbf{x},\alpha\right)=\eta\left(w^{\left(1\right)}+\widetilde{w},\mathbf{x},\alpha\right).$

The convexity of the set $K\left(R\right)$ allows us to use the multidimensional analog of Taylor formula in the following form:

[TABLE]

[TABLE]

[TABLE]

where $s_{j}$ are continuous functions of their listed variables. Furthermore,

[TABLE]

Next, substituting (46), (47) and (50) in (66), we obtain

[TABLE]

[TABLE]

where for brevity $\xi_{1}=\xi\left(w^{\left(1\right)},\mathbf{x},\alpha\right).$ Hence, it follows from (65)-(68) that the second line of (56) can be rewritten as

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Similar formulas are obviously valid for the third line of (56). Thus, formulas (56), (65)-(71) imply (60)-(64). $\square$

In Lemma 5.2 we have used second order terms in the Taylor formula. In addition, we will also need the formula which uses only linear terms. The proof of Lemma 5.3 is completely similar to the proof of Lemma 5.2.

Lemma 5.3. Assume that conditions of Lemma 5.2 hold. Then, the following analog of the final increment formula is valid

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $\mathbf{x}\in\Omega$, $i=1,2$, all elements of matrices $Y_{j},j=0,1,2$ are continuous functions of their variables and the following estimates are valid for $t\in\left[z,A+\sigma\right],\left(x,y,t\right),\mathbf{x}\in\overline{\Omega}:$

[TABLE]

[TABLE]

[TABLE]

6 Problem ( 53), ( 54) in the

Semidiscrete Form

We now rewrite equation (53) in the form of finite differences with respect to variables $x,y$ while keeping the continuous derivative with respect to $z$. For brevity we keep the same grid step size $h>0$ in both directions $x,y$. Consider partitions of the intervals $x\in\left(0,1\right),y\in\left(0,1\right)$ in small subintervals of the same length $h$ with $B=1/h$ and corresponding semidiscrete sub-domains of the domains $\Omega$ and $\overline{\Omega},$

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Hence, $\mathbf{x}^{h}\left(z\right)$ is a $z-$dependent vector. Its dimension is $\left(B-2\right)^{2}$ in the case of $\Omega^{h}$ and $\left(B+1\right)^{2}$ in the case of $\Omega_{1}^{h}.$ By (74) only those points $\left(x_{i},y_{j},z\right)\in\Omega^{h}$, which are corresponding interior points of the domain $\Omega.$ On the other hand, in addition to points of $\Omega^{h},$ the semidiscrete domain $\Omega_{1}^{h}$ contains boundary points which belong to the part $\Gamma$ of the boundary $\partial\Omega.$ We assume below that

[TABLE]

Remark 6.1. *Unlike classical forward problems for PDEs, we do not let the grid step size $h$ tend to zero. This is typical for numerical methods for many inverse problems: due to their ill-posed nature, see, e.g. [20, 29]. In other words, the grid step size is often used as the regularization parameter. *

Consider the $N-$D vector function $Q\left(\mathbf{x}\right)=\left(Q_{1},....,Q_{N}\right)^{T}\left(\mathbf{x}\right)$ with $Q_{n}\in C\left(\overline{\Omega}\right).$ We denote $Q^{h}\left(\mathbf{x}^{h}\left(z\right)\right)$ the trace of this vector function on the set $\overline{\Omega}_{1}^{h}$. Thus,

[TABLE]

is the matrix depending on the variable $z\in\left[A,A+\sigma\right]$. Here $Q_{k}^{i,j}\left(z\right)=Q_{k}\left(x_{i},y_{j},z\right),$ where $k=1,...,N.$ Hence, by (46), (47), (51) and (77) the finite difference analogs of functions $w\left(\mathbf{x},\alpha\right)$, $g\left(\mathbf{x},\alpha\right),f\left(x,y,A+\sigma\right)$ are:

[TABLE]

[TABLE]

[TABLE]

Next, by (49), (50), (52) and (78)-(80) the finite difference analogs of matrices $W$, $G$ and $F$ are

[TABLE]

[TABLE]

[TABLE]

For an arbitrary number $z\in\left[A,A+\sigma\right]$ denote

$\Omega_{z}^{h}=\left\{\mathbf{x}^{h}\left(z\right)=\left(x_{i},y_{j},z\right)_{i,j=1}^{B-1},\text{ }z\text{ is fixed}\right\}.$ We introduce semidiscrete functional spaces for matrices like $Q^{h}\left(\mathbf{x}\left(z\right)\right),$

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

We approximate $x,y$ derivatives of the vector function $W\left(\mathbf{x}\right)$ via central finite differences [26]. It is convenient to write this in the equivalent form for the matrix function $W^{h}\left(\mathbf{x}\right)$ as

[TABLE]

[TABLE]

[TABLE]

[TABLE]

See (77) for notations (86) and (87). Using (76) and (84)-(87), we obtain that there exists a constant $C_{2}=C_{2}\left(h_{0},N,\Omega^{h}\right)>0$ depending only on listed parameters such that

[TABLE]

In addition, by embedding theorem $H_{0,N}^{1,h}\left(\Omega_{1}^{h}\right)\subset C_{N,0}^{h}\left(\overline{\Omega}_{1}^{h}\right)$ and there exists a constant $C=C\left(A,\sigma,h_{0},N\right)>0$ such that

[TABLE]

Thus, using (77)-(83), we obtain the following finite difference analog of problem (53), (54):

[TABLE]

[TABLE]

Also, we assume that the vector functions $G^{h}\left(\mathbf{x}^{h}\left(z\right)\right)$ and $F^{h}\left(\mathbf{x}^{h}\left(z\right)\right)$ in (82) and (83) are such that

[TABLE]

As above, let $R>0$ be an arbitrary number. We now introduce the finite difference analogs of sets (57) and (59). Assume that (79) holds. Then

[TABLE]

[TABLE]

It follows from (89), (93) and (94) that $K^{h}\left(R\right)\subset K^{h}\subset C_{N,0}^{h}\left(\overline{\Omega}_{1}^{h}\right).$

Similarly with Lemmata 3.1 and 5.1, Lemma 6.1 provides a sufficient condition imposed on the components of the matrix $w^{h}\left(\mathbf{x}^{h}\left(z\right),\alpha\right),$ which guarantees that

$w^{h}\left(\mathbf{x}^{h}\left(z\right),\alpha\right)\in K^{h}.$

Lemma 6.1.* Let the matrix *$w^{h}=w^{h}\left(\mathbf{x}^{h}\left(z\right),\alpha\right)\in K^{h}.$ Select a triple $\left(i,j,z\right)$ with $i,j=0,...,B,z\in\left[A,A+\sigma\right]$ and *consider the vector *

$v\left(i,j,z\right)=\left(w_{1}+g_{1},w_{2}+g_{2},...,w_{N}+g_{N}\right)^{T}\left(i,j,z\right).$* Let $X_{N}$ be the $N\times N$ matrix of Lemma 3.1. Consider the vector $\widetilde{v}\left(i,j,z\right)=X_{N}^{T}\cdot v\left(i,j,z\right).$ Let *

$\widetilde{v}\left(i,j,z\right)=\left(\widetilde{v}_{1},...,\widetilde{v}_{N}\right)^{T}\left(i,j,z\right).$ Suppose that all numbers $\widetilde{v}_{n}\left(i,j,z\right)>0$ for all $i,j=0,...,B,z\in\left[A,A+\sigma\right].$ Then the function $w^{h}\in K$ and, therefore, by (37) and (43) the following analog of (58) holds for $\mathbf{x}^{h}\left(A+\sigma\right),\mathbf{x}^{h}\left(z\right)\in\overline{\Omega}_{1}^{h}:$

[TABLE]

* Also, the set $K^{h}\left(R\right)$ is convex*.

Proof. The first part of this lemma, the one about the positivity, follows immediately from Lemma 3.1. Consider now two matrices $w^{1,h}\left(\mathbf{x}^{h}\left(z\right),\alpha\right)$,$w^{2,h}\left(\mathbf{x}^{h}\left(z\right),\alpha\right)\in K^{h}\left(R\right).$ Let the number $\theta\in\left[0,1\right].$ Then one can prove completely similarly with the proof of Lemma 5.1 that $\theta w^{1,h}\left(\mathbf{x}^{h}\left(z\right),\alpha\right)+\left(1-\theta\right)w^{2,h}\left(\mathbf{x}^{h}\left(z\right),\alpha\right)\in K^{h}.$ Let $W^{1,h}\left(\mathbf{x}^{h}\left(z\right)\right)$ and $W^{2,h}\left(\mathbf{x}^{h}\left(z\right)\right)$ be two matrices corresponding to matrices $w^{1,h}\left(\mathbf{x}^{h}\left(z\right),\alpha\right)$ and

$w^{2,h}\left(\mathbf{x}^{h}\left(z\right),\alpha\right)$ respectively via (81).* *The triangle inequality and (94) imply that

[TABLE]

[TABLE]

[TABLE]

Lemma 6.2 is a finite difference analog of Lemmata 5.2 and 5.3. The proof is completely similar and is, therefore, omitted.

Lemma 6.2. Assume that (92) holds. Then the direct analogs of Lemmata 5.2 and 5.3, being applied to the right hand side of (90), are true, provided that all functions involved in these lemmata are replaced with their above semidiscrete analogs. *The constant $C_{1}$ in (64) and (73) should be replaced with the constant *$\widetilde{C}_{1}$ depending only on listed parameters, where

[TABLE]

Suppose that we have found such a matrix $W^{h}\left(\mathbf{x}^{h}\left(z\right)\right)\in K^{h}\left(R\right),$ that it solves problem (90), (91). Then, using (38), (42), (43) and (78)-(80), we set:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

. Hence, by (37), (93) and (94)

[TABLE]

Using (30), (31) and (99), we set

[TABLE]

The semidiscrete analogs of formulas (33) and (34) are:

[TABLE]

[TABLE]

Next, using the original eikonal equation (12), we set its semidiscrete form as

[TABLE]

Remark 6.2. *Equation (90) with condition (91) as well conditions (76), (78)-(80) and the assumption that the right hand side of (103) is independent on the parameter $\alpha$ form our approximate mathematical model for the TTTP formulated in section 2. One should expect that in practical computations the left hand side of (103) likely depends on *$\alpha.$ *However, the numerical experience of [18, 19, 20, 29], where the basis $\left\{\Psi_{n}\left(\alpha\right)\right\}_{n=1}^{\infty}$ was successfully used for four different inverse problems, shows that, numerically, one should consider the average value of the left hand side of (103) with respect to *$\alpha.$

Remarks 6.3:

(1) It is well known that the problems like proving the convergence of the numerical methods as ours when in (76), (78)-(80) actual regularization parameters $h_{0}\rightarrow 0$ and $N\rightarrow\infty$ are, generally, extremely challenging ones in the field of inverse problems. The fundamental underlying reason of these challenges is the ill-posed nature of inverse problems. Therefore, we do not analyze this type of convergence here.**

(2) Regardless on item 1, the author believes that the success of numerical studies of [18, 19, 20, 29] indicates that the numerical implementation of the method of this paper will likely lead to good computational results. In this regard, we also mention here the work of Guillement and Novikov [6] for a linear inverse problem as as well as the works of the group of Kabanikhin [8, 9, 10] for nonlinear coefficient inverse problems. In all these four publications truncated Fourier series were used to develop new numerical methods and to test them numerically. Although convergence estimates at $N\rightarrow\infty$ were not provided in [6, 8, 9, 10], numerical results are quite accurate ones.

7 Lipschitz Stability and Uniqueness

Based on (56) as well as on (88) and (89), we consider everywhere below the matrix equation (90), as a system of coupled nonlinear Volterra integral equations whose solution satisfies (91). Denote

[TABLE]

Theorem 7.1 (Lipschitz stability and uniqueness). *Let $R>0$ be an arbitrary number. Consider two matrices $W^{\left(1\right),h}\left(\mathbf{x}^{h}\left(z\right)\right),W^{\left(2\right),h}\left(\mathbf{x}^{h}\left(z\right)\right)\in K^{h}\left(R\right)$. Suppose that these matrices generate two pairs of matrices in (82), (83) $G^{\left(i\right),h}\left(\mathbf{x}^{h}\left(z\right)\right)$ and $F^{\left(i\right),h}\left(\mathbf{x}^{h}\left(A+\sigma\right)\right),i=1,2,$ which satisfy conditions (92).Let *$\Phi^{\left(i\right),h}\left(\mathbf{x}^{h}\left(z\right)\right),i=1,2$ be the corresponding matrices as in (104). Let $m^{\left(1\right),h}\left(\mathbf{x}\left(z\right)\right)$ and $m^{\left(2\right),h}\left(\mathbf{x}\left(z\right)\right)$ *be the corresponding right hand sides of equality (103). Assume that functions *$m^{\left(1\right),h}\left(\mathbf{x}\left(z\right)\right)$ and $m^{\left(2\right),h}\left(\mathbf{x}\left(z\right)\right)$ *are independent on *$\alpha$ (Remarks 6.2). Then there exists a constant

[TABLE]

* depending only on listed parameters such that the following Lipschitz stability estimate holds*

[TABLE]

In particular, if $\Phi^{\left(1\right),h}\left(\mathbf{x}^{h}\left(z\right)\right)\equiv\Phi^{\left(2\right),h}\left(\mathbf{x}^{h}\left(z\right)\right),$ then $m^{\left(1\right),h}\left(\mathbf{x}^{h}\left(z\right)\right)\equiv m^{\left(2\right),h}\left(\mathbf{x}^{h}\left(z\right)\right),$ i.e. uniqueness holds.

Proof. In this proof, $C_{3}>0$ denotes different constants depending on parameters listed in (105). Denote

[TABLE]

It follows from Lemma 6.2, equality (72), estimate (73) of Lemma 5.3, (88), (90), (91), (95)-(99) and (104) that the following inequality with the Volterra integral holds true

[TABLE]

where $\mathbf{x}\left(z\right)\in\overline{\Omega}_{1}^{h}$ and $z\in\left[A,A+\sigma\right].$ Hence, Gronwall’s inequality leads to

[TABLE]

This is the key estimate of this proof. Having this estimate, the target estimate (106) follows immediately from (81)-(83), (88), (95)-(103), (104) and (107). Uniqueness follows from (106). $\square$

8 Weighted Globally Strictly Convex Tikhonov-like Functional

8.1 Estimating an integral

Let $\lambda>0$ be the parameter to be chosen later. We choose the “integral analog” of the CWF as

[TABLE]

Lemma 8.1. *The following estimate holds for all $\lambda>0$ and for every function *$p\in L_{1}\left(A,A+\sigma\right):$

[TABLE]

Proof. Interchanging the integrals, we obtain

[TABLE]

[TABLE]

8.2 The functional

To solve problem (90), (91) numerically, we consider the following minimization problem:

Minimization Problem. *Fix an arbitrary number $R>0$ as well as the gird step size *$h\in\left[h_{0},1\right).$ Let $\gamma>0$ be the regularization parameter. Minimize the functional $J_{\lambda,\gamma}\left(W^{h}\right)$ on the closed set $\overline{K^{h}\left(R\right)},$ where

[TABLE]

[TABLE]

Here we took into account (108). We use the multiplier $e^{-2\lambda A}$ in order to balance two terms in the right hand side of (109). Note that since $R>0$ is an arbitrary number, then we do not impose a smallness condition on the set $\overline{K^{h}\left(R\right)}$ where we search for the solution of problem (90), (91). This is why we are talking below about the global strict convexity and the globally convergent numerical method.

Theorem 8.1 (global strict convexity).* Let $h_{0}$ be the number defined in (76) and let $h\in\left[h_{0},1\right)$. At every point *$W^{h}\in K^{h}\left(2R\right)$ *and for all *$\lambda>0,\gamma\in\left(0,1\right)$ the functional $J_{\lambda,\gamma}\left(W^{h}\right)$ has the Frechét derivative $J_{\lambda,\gamma}^{\prime}\left(W^{h}\right)\in H_{0}^{1,h}\left(\Omega_{1}^{h}\right)$. Furthermore, this derivative is Lipschitz continuous on $K^{h}\left(2R\right)$, i.e. there exists a constant

[TABLE]

* depending only on parameters listed in (110) such that for all $W^{\left(1\right),h},W^{\left(2\right),h}\in K^{h}\left(2R\right)$*

[TABLE]

*where the constant $\overline{C}>0$ depends on the same parameters as ones listed in (110) as well as on *$\lambda.$

In addition, there exists a sufficiently large number $\lambda_{0}>1$ depending on the same parameters as those listed in (110) such that for every $\lambda\geq\lambda_{0}$ and for every $\gamma\in\left(0,1\right)$ the functional $J_{\lambda,\gamma}\left(W^{h}\right)$ is strictly convex on the closed set $\overline{K^{h}\left(R\right)}.$ More precisely, the following estimate holds for all $W^{\left(1\right),h},W^{\left(2\right),h}\in\overline{K^{h}\left(R\right)}$

[TABLE]

[TABLE]

Below $C_{4}$ denotes different constants depending on parameters listed in (110)*. *Theorem 8.2 follows immediately from (111), (112) and Lemma 2.1 of [1]

Theorem 8.2. Suppose that conditions of Theorem 8.1 are in place. Then for every $\lambda\geq\lambda_{0}$ and for every $\gamma\in\left(0,1\right)$ there exists a single minimizer $W_{\min}^{h}\in\overline{K^{h}\left(R\right)}$ of the functional $J_{\lambda,\gamma}\left(W^{h}\right)$ on the set $\overline{K^{h}\left(R\right)}.$ Furthermore,

[TABLE]

Let $P_{\overline{K^{h}\left(R\right)}}:H_{0}^{1,h}\left(\Omega_{1}^{h}\right)\rightarrow\overline{K^{h}\left(R\right)}$ be the projection operator of the space $H_{0}^{1,h}\left(\Omega_{1}^{h}\right)$ on the closed set $\overline{K^{h}\left(R\right)}\subset H_{0}^{1,h}\left(\Omega_{1}^{h}\right).$ Consider an arbitrary point $W^{0,h}\in K^{h}\left(R\right).$ And minimize the functional $J_{\lambda,\gamma}\left(W^{h}\right)$ by the gradient projection method, which starts its iterations at the point $W^{0,h},$

[TABLE]

Here the number $\rho\in\left(0,1\right)$. Theorem 8.3 follows from a combination of Theorems 8.1 and 8.2 with Theorem 2.1 of [1].

**Theorem 8.3. Suppose that conditions of Theorem 8.1 are in place, $\lambda\geq\lambda_{0}$ and $\gamma\in\left(0,1\right).$ Let $W_{\min}^{h}\in\overline{K^{h}\left(R\right)}$ *be the minimizer listed in Theorem 8.2. Then there exists a sufficiently small number $\rho_{0}\in\left(0,1\right)$ depending on the same parameters as ones listed in (110) such that for every $\rho\in\left(0,\rho_{0}\right)$ there exists a number $\eta=\eta\left(\rho\right)\in\left(0,1\right)$ such that the sequence (114) converges to *$W_{\min}^{h},$

[TABLE]

In the regularization theory, the minimizer $W_{\min}^{h}$ of functional (109) is called “regularized solution”, see, e.g. [2, 31]. We now need to show that regularized solutions converge to the exact one when the noise in the data tends to zero. Following the regularization theory, we assume that there exists an exact, i.e. idealized, solution $W^{\ast,h}\in K^{h}\left(R\right)$ of problem (90), (91) with the noiseless data $\Phi^{\ast,h}\left(\mathbf{x}\left(z\right)\right),$

[TABLE]

where $\mathbf{x}^{h}\left(z\right)\in\overline{\Omega}^{h},$ see (104). We assume that there exists the exact, i.e. idealized function $m^{\ast,h}\left(\mathbf{x}^{h}\left(z\right)\right),$ $\mathbf{x}^{h}\left(z\right)\in\overline{\Omega}^{h}$ in (103), which produces the data (116).

Let the number $\delta\in\left(0,1\right)$ be the level of the error in the data $G^{h},F^{h},$ i.e.

[TABLE]

Denote $\widetilde{G}^{h}=G^{\ast,h}-G^{h},\widetilde{F}^{h}=F^{\ast,h}-F^{h}.$ Then (88) and (117) imply that with a constant $C_{2}=C_{2}\left(h_{0},N,\Omega^{h}\right)>0$ depending only on listed parameters the following inequalities hold:

[TABLE]

[TABLE]

Since $\delta\in\left(0,1\right)$ and (117)-(119) hold, then, using the triangle inequality, we replace below the dependence of the constant $C_{4}>0$ on $\left\|\Phi^{h}\right\|_{C_{2N}^{h}\left(\overline{\Omega}_{1}^{h}\right)}$ in (110) with the dependence of $C_{4}$ on $\left\|\Phi^{\ast,h}\right\|_{C_{2N}^{h}\left(\overline{\Omega}_{1}^{h}\right)}.$ Thus, everywhere below $C_{4}>0$ denotes different constants depending on the same parameters as ones listed in (110), except that $\left\|\Phi^{h}\right\|_{C_{2N}^{h}\left(\overline{\Omega}_{1}^{h}\right)}$ is replaced with $\left\|\Phi^{\ast,h}\right\|_{C_{2N}^{h}\left(\overline{\Omega}_{1}^{h}\right)}.$

Consider now the right hand side of equation (90) for the case when $W^{h}$ is replaced with $W^{\ast,h},$ whereas other arguments remain the same. By Lemma 6.2, we can use finite difference analogs of (72) and (73). In addition, we use now (116)-(119). Thus, we obtain

[TABLE]

[TABLE]

[TABLE]

Theorem 8.4. Suppose that conditions of Theorem 8.1 are in place and also that there exists an exact solution $W^{\ast,h}\in K^{h}\left(R\right)$ o*f problem (90), (91) with the noiseless data *$\Phi^{\ast,h}$ as in (116). *Let the $\lambda_{0}>1$ be the number chosen in Theorem 8.1. Fix an arbitrary number *$\lambda=\lambda_{1}\geq\lambda_{0}$ *in the functional *$J_{\lambda,\gamma}\left(W^{h}\right)=J_{\lambda_{1},\gamma}\left(W^{h}\right).$ Just as in the regularization theory, set $\gamma=\gamma\left(\delta\right)=\delta^{2}.$ Let the numbers $\rho\in\left(0,\rho_{0}\right),\rho_{0},\eta=\eta\left(\rho\right)\in\left(0,1\right)$ be the same as in Theorem 8.3. Then the following estimates are valid:

[TABLE]

[TABLE]

[TABLE]

where the functions $m_{n}^{h}\left(\mathbf{x}^{h}\left(z\right)\right)$ are constructed from matrices $W_{n}^{h}$ via the procedure described in section 6 with the final formula (103).

It follows from the above that we need to prove only Theorems 8.1 and 8.4.

9 Proofs of Theorems 8.1 and 8.4

9.1 Proof of Theorem 8.1

Denote $\left(,\right)$ and $\left[,\right]$ the scalar products in the spaces $L_{2,N}^{h}\left(\Omega^{h}\right)$ and $H_{0,N}^{1,h}\left(\Omega_{1}^{h}\right)$ respectively. Let $W^{\left(1\right),h}$ and $W^{\left(2\right),h}$ be two arbitrary points of the set $K^{h}\left(R\right).$ As above, denote $\widetilde{W}^{h}=W^{\left(2\right),h}-W^{\left(1\right),h}.$ Hence, $W^{\left(2\right),h}=W^{\left(1\right),h}+\widetilde{W}^{h}.$ Also, by (94)

[TABLE]

Hence, by (89) and (125)

[TABLE]

It follows from Lemma 6.2 and (60) that

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $i=1,2$ and $T_{j},j=0,1,2,3$ are continuous functions of their variables for $W^{\left(i\right),h}\in\overline{D^{h}\left(R\right)}.$ In addition, by (62) and (126)

[TABLE]

[TABLE]

Denote

[TABLE]

[TABLE]

[TABLE]

Then it follows from (109) and (127)-(131) that

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

In addition, by Lemma 6.2, (60), (61), (89), Lemma 8.1, (126) and (128)-(131) the following estimates hold:

[TABLE]

[TABLE]

Similarly using Cauchy-Schwarz inequality, Lemma 8.1 and assuming that the parameter $\lambda\geq\lambda_{0},$ where $\lambda_{0}$ is a sufficiently large number depending on the same parameters as ones listed in (110), we estimate from the below the sum of the fourth and fifth lines of (132) as:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

It follows from (127), (129) and (130) that the expression in the second line of (132) is linear with respect to $\widetilde{W}^{h},$

[TABLE]

Obviously $\left|Lin\left(\widetilde{W}^{h}\right)\right|\leq C_{4}\left\|\widetilde{W}^{h}\right\|_{H_{0,N}^{1,h}\left(\Omega_{1}^{h}\right)}.$ Hence, $Lin\left(\widetilde{W}^{h}\right):H_{0,N}^{1,h}\left(\Omega_{1}^{h}\right)\rightarrow\mathbb{R}$ is a bounded linear functional. Hence by Riesz theorem there exists a matrix $\Theta\in H_{0,N}^{1,h}\left(\Omega_{1}^{h}\right)$ such that

[TABLE]

Besides, it follows from the above that

[TABLE]

as $\left\|\widetilde{W}^{h}\right\|_{H_{0,N}^{1,h}\left(\Omega_{1}^{h}\right)}\rightarrow 0.$ Hence, $\Theta=J_{\lambda,\gamma}^{\prime}\left(W^{\left(1\right),h}\right)\in H_{0,N}^{1,h}\left(\Omega_{1}^{h}\right)$ is the Frechét derivative of the functional $J_{\lambda,\gamma}$ at the point $W^{\left(1\right),h}\in K^{h}\left(R\right).$ The existence of the Frechét derivative in the larger set $K^{h}\left(2R\right)$ can be proven completely similarly. The Lipschitz continuity property (111) of the Frechét derivative $J_{\lambda,\gamma}^{\prime}$ can be proven similarly with the proof of Theorem 3.1 of [1]. Therefore, we omit this proof for brevity.

Furthermore, using (132)-(136) and recalling that $\widetilde{W}^{h}=W^{\left(2\right),h}-W^{\left(1\right),h},$ we obtain for sufficiently large $\lambda_{0}>1$ and for $\lambda\geq\lambda_{0}$

[TABLE]

[TABLE]

which is the same as the target estimate (112) of this theorem. $\square$

9.2 Proof of Theorem 8.4

Recall that in this theorem we fix and arbitrary number $\lambda=\lambda_{1}\geq\lambda_{0},$ where $\lambda_{0}$ is the number of Theorem 8.1. To indicate the dependence of the functional $J_{\lambda_{1},\gamma}$ on $\Phi^{h},$ we write in this proof $J_{\lambda_{1},\gamma}\left(W^{h},\Phi^{h}\right).$

First, we consider $J_{\lambda_{1},\gamma}\left(W^{\ast,h},\Phi^{\ast,h}\right).$ Since

[TABLE]

then, using (109), we obtain

[TABLE]

Next, by (109), (120), (121), (137) and (138)

[TABLE]

where $Z\left(W^{\ast,h},W^{h},\Phi^{\ast h},\Phi^{h}\right)$ satisfies the following estimate

[TABLE]

Since $\lambda_{1}$ is a fixed arbitrary number such that $\lambda_{1}\geq\lambda_{0}$ and since $\lambda_{0}$ depends on the same parameters as those listed in (110) for $C_{4},$ then recalling that $\gamma=\delta^{2},$ we obtain from (138)-(140)

[TABLE]

Next, (112) implies that

[TABLE]

[TABLE]

Since by (113) $-J_{\lambda_{0},\gamma}^{\prime}\left(W_{\min}^{h},\Phi^{h}\right)\left(W^{\ast,h}-W_{\min}^{h}\right)\leq 0,$ then, using (141) and (142), we obtain

[TABLE]

The first target estimate (122) of this theorem follows from (143). Combining (115) and (122) with the procedure of section 6, which led to (103), we obtain two other target estimates (123) and (124) of Theorem 8.4. $\square$

Acknowledgments. This work was supported by US Army Research Laboratory and US Army Research Office grant W911NF-19-1-0044. The author is grateful to Professor Vladimir G. Romanov for many very fruitful discussions.