Travel time tomography with formally determined incomplete data in 3D

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.


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 (0, 1) space. That basis was first constructed in [16]. Functions of that basis depend only on the position of the source. The number N ≥ 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.  Figure 1: An illustration for complete and incomplete data in the 2D case, see details in [24]. To simplify, we assume in this figure that the geodesics are straight lines. Thus, we deal in this figure with the data of Radon transform, generated by the function "radon" of MATLAB. a) The true function m (x) to be imaged. b) The complete data of the Radon transform of the function of a). c) The incomplete data of the Radon transform of a) in the case when the source runs along an interval of a straight line, as in this paper below. [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 ≥ 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 ≥ 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 x = (x, y, z) ∈ R 3 . Let numbers A, σ = const. > 0. Define the rectangular prism Ω ⊂ R 3 as Ω = {x = (x, y, z) : x, y ∈ (0, 1) , z ∈ (A, A + σ)} . (1) Denote parts of the boundary ∂Ω as Let n (x) be the refractive index of the medium at the point x. Hence, c (x) = 1/n (x) is the sound speed. Denote m (x) = n 2 (x) . Let the number m 0 > 0 be given. We impose the following assumptions on the function m (x) : 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 (x) generates the Riemannian metric The travel time from the point x 0 ∈ R 3 (source) to the point x ∈ R 3 (receiver) is [24] τ (x, where Γ (x, x 0 ) is the geodesic line connecting points x and x 0 and ds is the euclidean arc length. We assume that the source x 0 runs along an interval L of a straight line located in the plane {z = 0} , Hence, x 0 = x α = (α, 1/2, 0) , α ∈ (0, 1) . Let τ (x,α) be the travel time between points x and x α = (α, 1/2, 0) . Thus, we denote τ (x,α) = τ (x, x α ) . We denote Γ (x, α) the geodesic line connecting points x and x α . It is well known [24] that the function τ (x, α) satisfies eikonal equation as the function of x, Everywhere below we assume without further mentioning that the following property holds: Regularity of Geodesic Lines. The function τ (x, α) ∈ C 2 (R 3 × [0, 1]) . For any pair of points (x, x α ) ∈ Ω × L there exists unique geodesic line Γ (x, α) connecting these two points and Γ (x, α) ∩ B A = ∅. In addition, if any geodesic line, which starts at a point x α ∈ 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 ∂Ω B A at another single point, see (1)-(4). In addition, after intersecting ∂Ω B A , this line does not "come back" in the domain Ω but rather goes away from this domain. In other words, this line is not reflected back from any point of its intersection with ∂Ω.
The following sufficient condition of the regularity of geodesic lines was derived in [25]: Travel Time Tomography Problem (TTTP). Suppose that the function m (x) satisfies conditions (5)- (8). Assume that the following function f (x, a) is given: Determine the function m (x) for x ∈ Ω. 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 Γ∪B A+σ of the boundary ∂Ω. Hence, the data (13) are both formally determined and incomplete.

A Special Orthonormal Basis
We now reproduce a special orthonormal basis {Ψ n (α)} ∞ n=0 in L 2 (0, 1) , which was constructed in [16]. This basis has the following two properties: 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 {ξ n (α)} ∞ n=0 = {(α + a) n e α } ∞ n=0 , 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 {Ψ n (α)} ∞ n=1 in L 2 (0, 1). In fact, the function Ψ n (α) has the form Ψ n (α) = P n (α + a)e α , ∀n ≥ 0, where P n is a polynomial of the degree n. Thus, functions Ψ n (α) are polynomials orthonormal in L 2 (0, 1) with the weight e 2α . The matrix M N is invertible since its elements a mn = (Ψ n , Ψ m ) are such that a mn = 1 if m = n and a mn = 0 if n < m.
Consider the function q (α) in the following form: Below we need to impose such a sufficient condition on the vector of coefficients q N = (q 0 , ..., q N −1 ) T in the Fourier expansion (14) which would guarantee that the func- The desired condition is given in Lemma 3.1.
Proof. It follows from the Gram-Schmidt procedure that elements of the matrix X N are independent on α. Let the raw number n of the matrix X N be (x n1 , x n2 , ..., x nN ) , n = 1, ..., N. Then Hence, by (14) Since ξ j (α) = (α + a) j e α > 0, then (15) Thus, Proof. Note that having proven (16) is not enough for our technique: we need to know the sign of the function τ z (x, α) , i.e. (17), in section 5 (more precisely, in (31)) where we consider the square root of τ 2 z (x, α) . Denote The following equations for geodesic lines can be found on page 66 of [24]: where s is a parameter. Using (19), we obtain for τ = τ (x (s) , y (s) , z (s) , α) along a geodesic line Set Then (22) implies: τ (x (s) , y (s) , z (s) , α) = 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): where α ∈ (0, 1) and p 0 , q 0 are some given numbers such that p 2 0 + q 2 0 ≤ 1. The latter inequality follows from (12) and the fact that by (6) m (x, y, 0) = 1. Also, by (12) To prove that we should take "+" sign in the latter formula, we note that Hence, τ z = r = z/τ > 0 for z ∈ (0, A) . Hence, Suppose that the geodesic line defined by (24) and (25) intersects the part B A of the boundary ∂Ω. Then the condition of the regularity of geodesic lines implies that there exists a single number and for all numbers z ∈ (A, z 0 ) all points (x (z) , y (z) , z) of that geodesic line belong to Ω. Since by (24) dr 2 /dz = m z , then, using (24) and (28), we obtain which establishes (16). To prove (17), we notice that it follows from (5), (8), (19), the last equation (21) and (23) that Estimates (17) and (18) follow immediately from (16), (26) and (29).
Introduce the function w (x, α) , Then (38)-(42) imply that We assume that both functions w (x, α) and g (x, α) have the form of the truncated Fourier series with respect to the orthonormal basis {Ψ n (α)}, Here, coefficients w n (x) are unknown and coefficients g n (x) are known. And similarly for w x , w y , w xα , w yα 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 Denote Remark 5.1. Note that we do not impose an analog of conditions (46), (47) on the function u(x, y, A, α). 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 Keeping in mind (48) and (49), we define the spaces C 1 (37)-(52), substitute functions w, g and their corresponding first derivatives with respect to x, y, α in equation (36). Next, multiply the resulting equation sequentially by functions Ψ n (α) , n = 1, ..., N and integrate with respect to α ∈ (0, 1). Then multiply both sides of obtained system of nonlinear integral differential equations by the matrix M −1 N , where the matrix M N was introduced in section 3. We obtain where P is the N −D vector function, 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 (u 0 + w + g) (x, α) It follows from (36) and (43) that we need the function (u 0 + w + g) (x, α) to be positive. We discuss this issue in the current section.
. Suppose that all functions v n (x) > 0 in Ω. Then the function w ∈ 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 (1) , w (2) ∈ K. Let the number θ ∈ (0, 1) . Then by (57) Summing up these two inequalities, we obtain

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 It follows from Lemma 5.1 and (59) that K (R) is a convex set.
The convexity of the set K (R) allows us to use the multidimensional analog of Taylor formula in the following form: where for brevity ξ 1 = ξ w (1) , x, α . Hence, it follows from (65)-(68) that the second line of (56) can be rewritten as Similar formulas are obviously valid for the third line of (56). Thus, formulas (56), (65)-(71) imply (60)-(64). 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 y , x, y, t G y (x, y, t) dt where x ∈ Ω, 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 ∈ [z, A + σ] , (x, y, t) , x ∈ Ω :

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 ∈ (0, 1) , y ∈ (0, 1) in small subintervals of the same length h with B = 1/h and corresponding semidiscrete sub-domains of the domains Ω and Ω, Hence, x h (z) is a z−dependent vector. Its dimension is (B − 2) 2 in the case of Ω h and (B + 1) 2 in the case of Ω h 1 . By (74) only those points (x i , y j , z) ∈ Ω h , which are corresponding interior points of the domain Ω. On the other hand, in addition to points of Ω h , the semidiscrete domain Ω h 1 contains boundary points which belong to the part Γ of the boundary ∂Ω. We assume below that 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 (x) = (Q 1 , ...., Q N ) T (x) with Q n ∈ C Ω . We denote Q h x h (z) the trace of this vector function on the set Ω h 1 . Thus, is the matrix depending on the variable z ∈ [A, A + σ]. Here Q i,j k (z) = Q k (x i , y j , z) , where k = 1, ..., N. Hence, by (46), (47), (51) and (77) the finite difference analogs of functions w (x, α), g (x, α) , f (x, y, A + σ) are: (79) Next, by (49), (50), (52) and (78)-(80) the finite difference analogs of matrices W , G and F are For an arbitrary number z ∈ [A, , z is fixed . We introduce semidiscrete functional spaces for matrices like Q h (x (z)) , We approximate x, y derivatives of the vector function W (x) via central finite differences [26]. It is convenient to write this in the equivalent form for the matrix function W h (x) as In addition, by embedding theorem Thus, using (77)-(83), we obtain the following finite difference analog of problem (53), (54): Also, we assume that the vector functions G h x h (z) and F h x h (z) in (82) and (83) are such that 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 It follows from (89), (93) and (94) Similarly with Lemmata 3.1 and 5.1, Lemma 6.1 provides a sufficient condition imposed on the components of the matrix w h x h (z) , α , which guarantees that j, z) . Suppose that all numbers v n (i, j, z) > 0 for all i, j = 0, ..., B, z ∈ [A, A + σ] . Then the function w h ∈ K and, therefore, by (37) and (43) the following analog of (58) holds for x h (A + σ) , x h (z) ∈ Ω h 1 : Also, the set K h (R) 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 x h (z) , α ,w 2,h x h (z) , α ∈ K h (R) . Let the number θ ∈ [0, 1] . Then one can prove completely similarly with the proof of Lemma 5.1 that θw 1,h x h (z) , α + (1 − θ) w 2,h x h (z) , α ∈ K h . Let W 1,h x h (z) and W 2,h x h (z) be two matrices corresponding to matrices w 1,h x h (z) , α and w 2,h x h (z) , α respectively via (81). The triangle inequality and (94) imply that 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 C 1 depending only on listed parameters, where x,y,N (Ω) > 0. Suppose that we have found such a matrix W h x h (z) ∈ K h (R) , that it solves problem (90), (91). Then, using (38), (42), (43) and (78)-(80), we set: . Hence, by (37), (93) and (94) Using (30), (31) and (99), we set The semidiscrete analogs of formulas (33) and (34) are: Next, using the original eikonal equation (12), we set its semidiscrete form as 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 α 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 α. However, the numerical experience of [18,19,20,29], where the basis {Ψ n (α)} ∞ n=1 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 α. 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 → 0 and N → ∞ 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 → ∞ were not provided in [6,8,9,10], numerical results are quite accurate ones.
Proof. In this proof, C 3 > 0 denotes different constants depending on parameters listed in (105). Denote Hence, Gronwall's inequality leads to This is the key estimate of this proof. Having this estimate, the target estimate (106) follows immediately from (81)

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 ∈ [h 0 , 1) . Let γ > 0 be the regularization parameter. Minimize the functional J λ,γ W h on the closed set K h (R), where . Here we took into account (108). We use the multiplier e −2λ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 K h (R) 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 ∈ [h 0 , 1). At every point W h ∈ K h (2R) and for all λ > 0, γ ∈ (0, 1) the functional J λ,γ W h has the Frechét derivative J λ,γ W h ∈ H 1,h 0 Ω h 1 . Furthermore, this derivative is Lipschitz continuous on K h (2R), i.e. there exists a constant depending only on parameters listed in (110) such that for all where the constant C > 0 depends on the same parameters as ones listed in (110) as well as on λ.
In addition, there exists a sufficiently large number λ 0 > 1 depending on the same parameters as those listed in (110) such that for every λ ≥ λ 0 and for every γ ∈ (0, 1) the functional J λ,γ W h is strictly convex on the closed set K h (R). More precisely, the following estimate holds for all W (1),h , W (2),h ∈ K h (R) 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 λ ≥ λ 0 and for every γ ∈ (0, 1) there exists a single minimizer W h min ∈ K h (R) of the functional J λ,γ W h on the set K h (R). Furthermore, Consider an arbitrary point W 0,h ∈ K h (R) . And minimize the functional J λ,γ W h by the gradient projection method, which starts its iterations at the point W 0,h , Here the number ρ ∈ (0, 1 Then there exists a sufficiently small number ρ 0 ∈ (0, 1) depending on the same parameters as ones listed in (110) such that for every ρ ∈ (0, ρ 0 ) there exists a number η = η (ρ) ∈ (0, 1) such that the sequence (114) converges to W h min , In the regularization theory, the minimizer W h min 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 * ,h ∈ K h (R) of problem (90), (91) with the noiseless data Φ * ,h (x (z)) , where x h (z) ∈ Ω h , see (104). We assume that there exists the exact, i.e. idealized function (103), which produces the data (116).
Let the number δ ∈ (0, 1) be the level of the error in the data G h , F h , i.e. (88) and (117) imply that with a constant C 2 = C 2 h 0 , N, Ω h > 0 depending only on listed parameters the following inequalities hold: Since δ ∈ (0, 1) and (117)-(119) hold, then, using the triangle inequality, we replace below the dependence of the constant . Thus, everywhere below C 4 > 0 denotes different constants depending on the same parameters as ones listed in (110), except that Φ h Consider now the right hand side of equation (90) for the case when W h is replaced with W * ,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 Theorem 8.4. Suppose that conditions of Theorem 8.1 are in place and also that there exists an exact solution W * ,h ∈ K h (R) of problem (90), (91) with the noiseless data Φ * ,h as in (116). Let the λ 0 > 1 be the number chosen in Theorem 8.1. Fix an arbitrary number λ = λ 1 ≥ λ 0 in the functional J λ,γ W h = J λ 1 ,γ W h . Just as in the regularization theory, set γ = γ (δ) = δ 2 . Let the numbers ρ ∈ (0, ρ 0 ) , ρ 0 , η = η (ρ) ∈ (0, 1) be the same as in Theorem 8.3. Then the following estimates are valid: where the functions m h n x h (z) are constructed from matrices W h n via the procedure described in section 6 with the final formula (103).
It follows from the above that we need to prove only Hence, by (89) and (125) It follows from Lemma 6.2 and (60) that where i = 1, 2 and T j , j = 0, 1, 2, 3 are continuous functions of their variables for W (i),h ∈ D h (R). In addition, by (62)  x h (t) dt, x h (z) ∈ Ω h .
In addition, by Lemma 6.2, (60), (61), (89), Lemma 8.1, (126) and (128)-(131) the following estimates hold: Similarly using Cauchy-Schwarz inequality, Lemma 8.1 and assuming that the parameter λ ≥ λ 0 , where λ 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: .
It follows from (127), (129) and (130) that the expression in the second line of (132) is linear with respect to W h , → 0. Hence, Θ = J λ,γ W (1),h ∈ H 1,h 0,N Ω h 1 is the Frechét derivative of the functional J λ,γ at the point W (1),h ∈ K h (R) . The existence of the Frechét derivative in the larger set K h (2R) can be proven completely similarly. The Lipschitz continuity property (111) of the Frechét derivative J λ,γ can be proven similarly with the proof of Theorem 3.1 of [1]. Therefore, we omit this proof for brevity.

Proof of Theorem 8.4
Recall that in this theorem we fix and arbitrary number λ = λ 1 ≥ λ 0 , where λ 0 is the number of Theorem 8.1. To indicate the dependence of the functional J λ 1 ,γ on Φ h , we write in this proof J λ 1 ,γ W h , Φ h .