A globally convergent numerical method for a 1-d inverse medium problem with experimental data

In this paper, a reconstruction method for the spatially distributed dielectric constant of a medium from the back scattering wave field in the frequency domain is considered. Our approach is to propose a globally convergent algorithm, which does not require any knowledge of a small neighborhood of the solution of the inverse problem in advance. The Quasi-Reversibility Method (QRM) is used in the algorithm. The convergence of the QRM is proved via a Carleman estimate. The method is tested on both computationally simulated and experimental data.


Introduction
In this paper we develop a globally convergent numerical method for a 1-d inverse medium problem in the frequency domain. The performance of this method is tested on both computationally simulated and experimental data. Propagation of electromagnetic waves is used in experiments. A theorem about the global convergence of our method is proved. Another name for inverse medium problems is Coefficient Inverse Problems (CIPs). We model the process as a 1-d CIP due to some specifics of our data collection procedure, see subsection 6.2. This paper is the first one in which the globally convergent method of [1,2,11,18,19,28,29] is extended to the case of the frequency domain. Indeed, the original version of that method works with the Laplace transform of the time dependent data. Both the 3-d and the 1-d versions of the method of [1,2] were verified on experimental data, see [1,28,29] for the 3-d case and [11,18,19] for the 1-d case. The experimental data here are the same ones as in [11,18,19].
Compared with the previous publications [11,18,19], two additional difficulties occurring here are: (1) we now work with complex valued functions instead of real valued ones and, therefore (2) it is not immediately clear how to deal with the imaginary part of the logarithm of the complex valued solution of the forward problem. On the other hand, we use that logarithm in our numerical procedure. Besides, the previous numerical scheme of [11,18,19] is significantly modified here.
We call a numerical method for a CIP globally convergent if there is a rigorous guarantee that this method reaches at least one point in a sufficiently small neighborhood of the correct solution without any advanced knowledge about this neighborhood. It is well known that the topic of the global convergence for CIPs is a highly challenging one. Thus, similarly with the above cited publications, we prove the global convergence of our method within the framework of a certain approximation. This is why the term "approximate global convergence" was used in above cited references ("global convergence" in short). That approximation is used only on the first iterative step of our algorithm and it is not used in follow up iterations. Besides, this approximation is a quite natural one, since it amounts to taking into account only the first term of a certain asymptotic behavior and truncating the rest. The global convergence property is verified numerically here on both computationally simulated and experimental data.
CIPs are both highly nonlinear and ill-posed. These two factors cause the non-convexity of conventional Tikhonov functionals for these problems. It is well known that typically those functionals have multiple local minima and ravines, see, e.g. numerical examples of multiple local minima in [26] and on Figure 5.3 of [14]. Figure 5.4 of [14] demonstrates an example of a ravine. Hence, there is no guarantee of the convergence of an optimization method to the correct solution, unless the starting point of iterations is located in a sufficiently small neighborhood of that solution. However, it is often unclear in practical scenarios how to reach that neighborhood.
In our inverse algorithm we solve a sequence of linear ordinary differential equations of the second order on the interval x ∈ (0, 1). The peculiarity here is that we have overdetermined boundary conditions for these equations: we have the solution and its first derivative at x = 0 and we have the first derivative of the solution at x = 1. Our attempts to use only boundary conditions at x = 0 did not lead to good reconstruction results. The same observation is in place in Remark 5.1 on page 13 of [18]. Thus, similarly with [11,18,19], we use here the Quasi-Reversibility Method (QRM). The QRM is well suitable to solve PDEs with the overdetermined boundary data. To analyze the convergence rate of the QRM, we use a Carleman estimate.
The QRM was first introduced by Lattes and Lions [20]. However, they have not established convergence rates. These rates were first established in [12,13], where it was shown that Carleman estimates are the key tool for that goal. A survey of applications of Carleman estimates to QRM can be found in [16]. Chapter 6 of [1] describes the use of the QRM for numerical solutions of CIPs. We also mention an active work with the QRM of Bourgeois and Dardè, see, e.g. [4,5,6] for some samples of their publications.
Our experimental data are given in the time domain. To obtain the data in the frequency domain, we apply the Fourier transform. Previously the same data were treated in works [11,18,19] of this group. In that case the Laplace transform was applied. Next, the 1-d version of the globally convergent method of [1,2] was used. Indeed, the original version of this method works with the Laplace transform of the time dependent data. A similar 1-d inverse problem in the frequency domain was solved numerically in [21] by a different method.
The experimental data of this publication were collected by the Forward Looking Radar which was built in the US Army Research Laboratory [23]. The data were collected in the field (as opposed to a laboratory). Thus, clutter was present at the measurement site. This certainly adds some additional difficulties to the imaging problem. The goal of this radar is to detect and identify shallow explosives, such as, e.g. improvised explosive devices and antipersonnel land mines. Those explosives can be located either on the ground surface or a few centimeters below this surface. This radar provides only a single time dependent curve for a target. Therefore, to solve a CIP, we have no choice but to model the 3-d process by a 1-d wave-like PDE. We model targets as inclusions whose dielectric constants are different from the background.
In terms of working with experimental data, the goal in this paper is not to image locations of targets, since this is impossible via solving a CIP with our data, see subsection 6.2 for details. In fact, our goal is to image maximal values of dielectric constants of targets. Our targets are 3-d objects, while we use a 1-d model: since we measure only one time resolved curve for a target. Nevertheless, we show below that our calculated dielectric constants are well within tabulated limits [27].
Of course, an estimate of the dielectric constant is insufficient for the discrimination of explosives from the clutter. On the other hand, the radar community mostly relies on the intensity of the radar image. Therefore, we hope that the additional information about values of dielectric constants of targets might lead in the future to designs of algorithms which would better discriminate between explosives and clutter.
In the 3-d case the globally convergent method of [1,2,28,29] works with the data generated either by a single location of the source or by a single direction of the incident plane wave. The second globally convergent method for this type of measurements was developed and tested numerically in [15]. We also refer to another globally convergent methods for a CIPs, which is based on the multidimensional version of the Gelfand-Levitan method. This version was first developed by Kabanikhin [8] and later by Kabanikhin and Shishlenin [9,10]. Unlike [1,2,28,29], the technique of [8,9,10] works with multiple directions of incident plane waves.
In [11] the performance of the method of [18,19] was compared numerically with the performance of the Krein equation [17] for the same experimental data as ones used here. The Krein equation [17] is close to the Gelfand-Levitan equation [22]. It was shown in [17] that the performance of the Krein equation is inferior to the performance of the method of [18,19] for these experimental data. This is because the solution of the Krein equation is much more sensitive to the choice of the calibration factor than the solution obtained by the technique of [18,19]. On the other hand, it is necessary to apply that factor to those experimental data to make the range of values of the resulting data comparable with the range of values of computationally simulated data.
In section 2 we pose forward and inverse problems and analyze some of their features. In section 3 we study a 1-d version of the QRM. In section 4 we present our globally convergent method. In section 5 we prove the global convergence of our method. In section 6 we present our numerical results for both computationally simulated and experimental data.

Some Properties of Forward and Inverse Problems
It was shown numerically in [3] that the component of the electric field, which is incident upon a medium, dominates two other components. It was also shown in [3] that the propagation of that component is well governed by a wave-like PDE. Furthermore, this finding was confirmed by imaging from experimental data, see Chapter 5 of [1] and [28,29]. Thus, just as in [11,18,19], we use a 1-d wave-like PDE to model the collection of our experimental data of electromagnetic waves propagation.

Formulations of problems
Let c 0 < c 1 be two positive numbers. Let the function c : R → R satisfy the following conditions: Fix the source position x 0 < 0. Consider the generalized Helmholtz equation in the 1-d case, The problem (2.3), (2.4) is the forward problem. Our interest is in the following inverse problem problem: Problem (Coefficient Inverse Problem (CIP)). Fix the source position x 0 < 0. Let [k, k] ⊂ (0, ∞) be an interval of frequencies k. Reconstruct the function β (x) , assuming that the following function g 0 (k) is known (2.6)

Some properties of the solution of the forward problem
In this subsection we establish existence and uniqueness of the forward problem. Even though these results are likely known, we present them here for reader's convenience. In addition, the techniques using in their proofs help us to verify that the function u(x, x 0 , k) never vanishes for x > x 0 , which plays an important role in our algorithm of solving the to CIP. Also, this technique justifies the numerical method for solving the forward problem ( Proof. We prove uniqueness first. Suppose that u 1 and u 2 are two solutions of the problem (2.3), (2.4). Denote U = u 1 − u 2 . Then the function U satisfies Since the function β (x) = 0 outside of the interval (0, 1) , then (2.7) implies that U +k 2 U = 0 for x / ∈ (0, 1) . This, together with (2.8), yields where B 1 and B 2 are some complex numbers depending on x 0 , k. Let R > 1 be an arbitrary number. Multiply both sides of (2.7) by the function U and integrate over the interval (−R, R) using integration by parts. We obtain Hence, the imaginary part of (2.10), −k |B 1 (x 0 , k)| 2 + |B 2 (x 0 , k)| 2 = 0, vanishes, so do B 1 (x 0 , k) and B 2 (x 0 , k). Using (2.9), we obtain U (x, x 0 , k) = 0 for x / ∈ (0, 1) . Due to the classical unique continuation principle, U (x, x 0 , k) = 0, ∀x ∈ R. Thus, u 1 = u 2 .
We now prove existence. Consider the 1-d analog of the Lippman-Schwinger equation with respect to a function P , Fix x 0 < 0 and k > 0. Consider equation (2.11) only for x ∈ (0, 1) . Assume that there exist two functions P 1 , P 2 satisfying (2.11) for x ∈ (0, 1) . Consider their extensions on the set R (0, 1) via the right hand side of (2.11). Then so defined functions P 1 , P 2 satisfy (2.11) for all x ∈ R. Hence, both of them are solutions of the problem (2.3), (2.4). Hence, the above established uniqueness result implies that P 1 ≡ P 2 .
Consider again the integral equation (2.11) for x ∈ (0, 1) and apply the Fredholm alternative. This alternative, combined with the discussion in the previous paragraph, implies that there exists unique solution P ∈ C [0, 1] of equation (2.11). Extending this function for x ∈ R (0, 1) via the right hand side of (2.11), we obtain that there exists unique solution P ∈ C (R) of equation (2.11). Furthermore, the function P − u 0 ∈ C 3 (R) and also this function P is the required solution of the problem Below we consider only such a solution u of the problem (2.3), (2.4) that u − u 0 ∈ C 3 (R): as in Theorem 2.1.

The asymptotic behavior of the function
Furthermore, for any finite interval (a, b) ⊂ R there exists a number γ = γ (a, b, φ) > 0 such that for all x ∈ (a, b) the function u (x, x 0 , k) can be analytically extended with respect k from {k : k > 0} in the half plane C γ = {z ∈ C : Im z < γ}.
Proof. Consider the following Cauchy problem This proof is based on the connection between the solution u of the problem (2.3), (2.4) and the solution u of the problem (2.14), (2.15) via the Fourier transform with respect to t. Here and below the Fourier transform is understood in terms of distributions, see, e.g. the book [33] for the Fourier transform of distributions.
It follows from (2.1), (2.2), (2.16) and (2.21) that the function p ∈ C (R) and has a finite support, p (y) = 0 for y < −x 0 and for y > Recall that φ (x) is the function defined in (2.12). Expression (2.21) in the x−coordinate becomes Let H (z) , z ∈ R be the Heaviside function. It is well known that the solution v (y, t) of the problem (2.19), (2.20) has the following form, see, e.g. Chapter 2 in [24] v (y, t, The backwards substitution y → x transforms the function p (y) in the function φ (x) . Hence, (2.23) implies that the function p (y) in (2.21) is non positive, (2.26) Consider now the operator F of the Fourier transform, In the sense of distributions we have see section 6 in §9 of Chapter 2 of [33]. Consider an arbitrary finite interval (a, b) ⊂ R. Let s, m ≥ 0 be two arbitrary integers such that s + m ≤ 2. Since by (2.22) the function p (y) has a finite support, then (2.26), Lemma 6 of Chapter 10 of the book [32] as well as Remark 3 after that lemma guarantee that functions ∂ s y ∂ m t v (y, t) decay exponentially with respect to t, as long as y ∈ (a, b) . Hence, one can apply the operator F to the functions ∂ s y ∂ m t v (y, t) in the regular sense. The assertion about the analytic extension follows from (2.27) and the exponential decay of the functions ∂ s . The next question is whether the function V satisfies analogs of conditions (2.3), (2.4). Theorem 3.3 of [31] and theorem 6 of Chapter 9 of [32] guarantee that the function V satisfies the following conditions Using (2.25) and the integration by parts, we obtain for k → ∞

Some properties of the solution of the inverse problem
Since the source position x 0 < 0 is fixed, we drop everywhere below the dependence on x 0 in notations of functions. First, we show that, having the function g 0 (k) in (2.6), one can uniquely find the function u x (0, k) . Indeed, for x < 0 conditions (2.3) and (2.4) become for a complex number B (k). Note that by (2.5) and (2.6) Hence, by (2.34) and the definition of the function u, we have Thus, (2.6), (2.35) and (2.36) provide us with an additional data g 1 (k) to solve our CIP, where (2.37) Theorem 2.3 (uniqueness of our CIP). Let φ (x) be the function defined in (2.12). Assume that φ (x) ≤ 0, ∀x ∈ [0, 1] . Then our CIP has at most one solution.
Proof. By Theorem 2.2 one can analytically extend the function g 0 (k) from the interval k ∈ [k, k] in the half plane C γ = {z ∈ C : Im z < γ} , where γ = γ k, k, φ is a positive number. Hence, it follows from (2.36) and (2.37) that we know functions u (0, k) , u x (0, k) for all k ∈ R. Consider the inverse Fourier transform of the function u with respect to k, F −1 (u) . It was established in the proof of Theorem 2.2 that this transform can indeed be applied to the function u and F −1 (u) = u (x, t), where the function u (x, t) is the solution of the problem (2.14), (2.15). Functions are known. Hence, we have obtained the inverse problem for equation (2.14) with initial conditions (2.15) and the data (2.38). It is well known that this inverse problem has at most one solution, see, e.g. Chapter 2 of [24].

A Version of the Quasi-Reversibility Method
As it was pointed out in section 1, we solve ordinary differential equations with over determined boundary conditions using the QRM. In this section, we develop the QRM for an arbitrary linear ordinary differential equation of the second order with over determined boundary conditions. Below for any Hilbert space H its scalar product is ·, · H . For convenience, we use in this section the notation "w" for a generic complex valued function, which is irrelevant to the function w in (2.36). Let functions a (x) and b (x) be in C([0, 1], C) and the function d (x) be in L 2 ([0, 1], C). Let p 0 and p 1 be complex numbers. In this section we construct an approximate solution of the following problem: The QRM for problem (3.1) amounts to the minimization of the following functional Below in this section we will establish existence and uniqueness of the minimizer of J α . We will also show how close that minimizer is to the solution of (3.1) if it exists. We start from the Carleman estimate for the operator d 2 /dx 2 .
3.1 Carleman estimate for the operator d 2 /dx 2 Lemma 3.1 (Carleman estimate). For any complex valued function u ∈ H 2 (0, 1) with u(0) = u (0) = 0 and for any parameter λ > 1 the following Carleman estimate holds where C is a constant independent of u and λ. Proof Moreover, without loss of the generality, we can assume that u is real valued.
Introduce the function v = ue −λx , x ∈ (0, 1). We have A simple calculation yields One the other hand, we have Combine (3.3) and (3.4) and then integrate the result. We obtain Inequality (3.2) follows.

The existence and uniqueness of the minimizer of J α
We first establish the existence and uniqueness for the minimizer of J α .
Let {w n } n≥1 ⊂ W be a minimizing sequence of J α . That means, It is not hard to see that {w n } n≥1 is bounded in H 3 (0, 1). In fact, if {w n } n≥1 has a unbounded subsequence then inf Without the loss of generality, we can assume that {w n } n≥1 weakly converges in H 3 (0, 1) to a function w α and strongly converges to w α in H 2 (0, 1). The function w α belongs to W because W is close and convex. We have The uniqueness of w α is due to the strict convexity of J α . Inequality (3.5) can be verified by the fact that Let r 1 and r 2 be the real and imaginary parts respectively of the complex valued function r. Without confusing, we identify r with the pair of real valued functions (r 1 , r 2 ). Define (3.8) In order to find w α , we find the zero of the Fréchet derivative DJ α of J α (w α ). Since the operators L 1 and L 2 are linear, DJ α is given by The existence of a zero of DJ α follows from the above existence of w α and the uniqueness is, again, deduced from the convexity of J α and W . In our computations, we use the finite difference method to approximate the equation DJ α (w 1 , w 2 ) = 0, together with the condition w 1 + iw 2 ∈ W , as a linear system for w α . The minimizer w α of J α is called the regularized solution of (3.1) [1, 30].

Convergence of regularized solutions
While Theorem 3.1 claims the existence and uniqueness of the regularized solution w α of problem (3.1), we now prove convergence of regularized solutions to the exact solution of this problem, provided that the latter solution exists, see [1,30] for the definition of the regularized solution. It is well known that one of concepts of the Tikhonov regularization theory is the a priori assumption about the existence of an exact solution of an ill-posed problem, i.e. solution with noiseless data [1,30]. Estimate (3.5) is valid for the H 3 (0, 1) −norm and it becomes worse as long as α → 0. However, Theorem 3.2 provides an estimate for the H 2 (0, 1) −norm and the latter estimate is not worsening as α → 0. To prove Theorem 3.2, we use the Carleman estimate of subsection 3.1.
Suppose that there exists the exact solution w * of the problem (3.1) with the exact (i.e. noiseless) data d * ∈ L 2 , p * 0 , p * 1 ∈ C. Let the number δ ∈ (0, 1) be the level of the error in the data, i.e.
Let w α ∈ W be the unique minimizer of the functional J α , which is guaranteed by Theorem 3.1.
Theorem 3.2 (convergence of regularized solutions). Assume that (3.9) holds. Then there exists a constant C 2 > 0 depending only on a L ∞ (0,1) and b L ∞ (0,1) such that the following estimate holds In particular, if α = δ 2 , then the following convergence rate of regularized solutions w α takes place (with a different constant C 2 ) Proof. In this proof, C 2 is a generic constant depending only on a L ∞ (0,1) and b L ∞ (0,1) . Let the function χ ∈ C 2 [0, 1] satisfies the following condition Define the "error" function Since w * is a solution of (3.1), then 1) . (3.14) Denoting v = w α − w * − E, using the test function h = v, and using the Cauchy-Schwarz inequality, we derive in a standard way from (3.9), (3.13) and (3.14) that On the other hand, Lemma 3.1 gives Choosing λ sufficiently large, we obtain  Below w (x, k) is the function defined in (2.36). Since by Theorem 2.2 w (x, k) = 0, ∀x, then we can consider log w (x, k) . Since log z = ln |z|+i arg z, ∀z ∈ C, z = 0, then the natural question is about arg w (x, k) = Im w (x, k) . Hence, we consider the asymptotic behavior at k → ∞ of the function w (x, k). Using (2.5) and Theorem 2.2, we obtain for x > x 0 Obviously arg (1 + O(1/k)) ∈ [−π, π] for sufficiently large k > 0. Hence, we set for sufficiently large k > 0 and for The function log w (x, k) is defined via (4.2) for sufficiently large k. On the other hand, for not large values of k it would be better to work with derivatives of log w (x, k) . Indeed, we would not have problems then with defining arg w (x, k) . Hence, taking k sufficiently large, we define the function φ(x, k) = log w(x, k) as Differentiate (4.3) with respect to k. We have Multiplying both sides of the equation above by exp(−φ(x, k)) gives The function φ(x, k), therefore, defines log w(x, k).
Remark 4.1. It follows from (2.6), (2.36) and (4.1) that g 0 (k) = 1 + O (1/k) as k → ∞. Hence, arg g 0 (k) ∈ [−π, π] for sufficiently large k > 0. Hence, the function log g 0 (k) can be defined similarly with the function v (x, k) in (4.5). Let We call V the "tail function" and this function is unknown. Note that we do not use below the function V . Rather we use only its x−derivatives. Hence, when using these derivatives, we are not concerned with arg w x, k . It easily follows from (2.2), (2.3), (2.5), (2.36) and (2.37) that (4.10) Using (4.4) and (4.9), we obtain Therefore, (4.6)-(4.11) imply that , q x (0, k) = 0. (4.13) We have obtained an integral differential equation (4.12) for the function q with the overdetermined boundary data (4.13). The tail function in (4.12) is also unknown. Hence, to approximate both functions q and V , we need to use not only conditions (4.12), (4.13) but something else as well. Thus, in our iterative procedure, we solve problem (4.12), (4.13), assuming that V is known, and update the function q this way. Then we update the unknown coefficient β (x) . Next, we solve problem (4.9), (4.10) for the function w at k := k and update the tail function via (4.6), (4.7) and (4.8).

Initial approximation V 0 (x) for the tail function
It is important for the above iterative process to properly choose the initial approximation V 0 (x) for the tail function. Since we want to construct a globally convergent method, this choice must not use any advanced knowledge of a small neighborhood of the exact solution c * (x) of our inverse problem.
We now describe how do we choose the initial tail V 0 (x) . It follows from Theorem 2.2 and the definition of V (x) via (4.1)-(4.8) that there exists a function p (x) ∈ C 2 [0, 1] such that Hence, assuming that the number k is sufficiently large, we drop terms O 1/k 2 and O 1/k 3 in (4.14) and set We solve the problem (4.16), (4.17) via the QRM. By the embedding theorem H 2 (0, 1) ⊂ C 1 [0, 1] and f C 1 [0,1] ≤ C f H 2 (0,1) , ∀f ∈ H 2 (0, 1) , where C > 0 is a generic constant. Recall that the function g 1 (k) in (4.17) is linked with the function g 0 (k) as in (2.37). Thus, Theorems 3.1 and 3.2 lead to Theorem 4.1. In this theorem, we use the entire interval k, k rather than just k = k (in (4.18)) for brevity: since we will use this interval below. For k ≥ k, let the exact tail V * (x, k) have the form (4.15). Assume that for k ∈ k, k where δ > 0 is a sufficiently small number, which characterizes the level of the error in the boundary data. Let the function V 0,a (x) ∈ H 3 (0, 1) be the approximate solution of the problem (4.16)-(4.17) obtained via the QRM with α = δ 2 . Then there exists a constants C 3 = C 3 k, c * > 0 depending only on k and c * such that 1. Theorem 4.1 is valid only within the framework of a quite natural approximation (4.15), in which small terms O (1/k 2 ) , O (1/k 3 ) of formulae (4.14) are ignored. We use this approximation only to find the first tail and do not use it in follow up iterations. We believe that the use of this approximation is justified by the fact that the topic of the globally convergent numerical methods for CIPs is a very challenging one.
2. Thus, it follows from Theorem 4.1 that our initial tail function V 0,a x, k provides a good approximation for the exact tail V * x, k already at the start of our iterative process. Hence, setting in (4.11) k = k and recalling (4.8), we conclude that the target coefficient c * (x) is also reconstructed with a good accuracy at the start of our iterative process. The error of the approximation of both V * and c * depends only on the level δ of the error in the boundary data. The latter is exactly what is usually required when solving inverse problems. It is important that when obtaining this approximation for V * , we have not used any advanced knowledge about a small neighborhood of the exact solution c * . In other words, the requirement of the global convergence is in place (see Introduction for this requirement).
3. Even though we obtain good approximations for V * (x, k) and c * (x) from the start, our numerical experience tells us that results improve with iterations in our iterative process described below. A similar observation took place in the earlier above cited works of this group, where the Laplace transform of the time dependent data was used. This is of course due to the approximate nature of (4.15).
4. Even though it is possible to sort of "unite in one" first two conditions (4.18), we are not doing this here for brevity.

5.
In the convergence analysis, we use the form (4.15) for the functions V * and q * only on the first iteration, since this form of functions V ,q is used only on the first iteration of our algorithm.
Below we consider the error parameter η defined as

Numerical method 4.3.1 Equations for q n
Consider a partition of the frequency interval k, k in N subintervals with the step size h, where the number h > 0 is sufficiently small. We assume that the function q (x, k) is piecewise constant with respect to k, q (x, k) = q (x, k n ) for k ∈ [k n , k n−1 ) . For each n = 1, · · · , N and for all x ∈ (0, 1) define q 0 (x) = 0, q n (x) = q(x, k n ), (4.20) Then (4.12) and (4.20)- (4.22) imply that for all n = 1, · · · , N Choose the step size h sufficiently small and ignore terms with h and h 2 . Note that k−k n < h for k ∈ k ∈ [k n , k n−1 ) . Also, we keep in mind that we will iterate with respect to tail functions for each n as well as with respect to n. Thus, we rewrite the last equation as q n,j + 2k 2 n −Q n−1 + V n,j − 2ik n q n,j = −2k n −Q n−1 + V n,j 2 + 2i −Q n−1 + V n,j (4.23) for all x ∈ (0, 1), j = 1, · · · , m for some m > 0. The boundary conditions for q n,j in (4.23) are taken according to those for q in (4.13). Precisely, = ψ 1 n , q n,j (1) = 0.
ii. For s = 1, 2, let v (s) n,j due to an analog of (4.22). iii. Calculate β n,j by (4.26), which will be explained later.
Let w n,j x, k be its solution. Next, using (4.5) and (4.7), set .
In the algorithm, all differential equations are solved via the QRM. Thus, we keep for those "QRM solutions" the same notations for brevity. For simplicity, we assume here that β (x) ≥ 0, although we also work in one case of experimental data with a non-positive function β. Thus, in the algorithm above, we update the function β (x) using (2.1), (2.2), (4.11) and (4.22) as β n,j = min max −v n,j − k 2 n v n,j 2 + 2ik n v n,j , c 0 − 1 , c 1 − 1 . This truncation helps us to get a better accuracy in the reconstructed function β. In fact, it follows from (2.1), (2.2), (4.11) and (4.26) that where β * = c * − 1 is the exact solution of the CIP.

Global Convergence
In this section we prove our main result about the global convergence of the algorithm of the previous section. This method actually has the approximately global convergence property, see the third paragraph of Section 1 and Remarks 4.2. For brevity we assume in this section that j 0 = 1 in the above algorithm. In other words, we assume that we do not perform inner iterations. The case j 0 > 1 can be done similarly. First, we need to introduce some assumptions about the exact solution. Everywhere below the superscript " * " denotes functions which correspond to the exact coefficient c * (x). Denote q * n (x) = q * (x, k n ) . Then q * (x, k) = q * n (x) + O (h) for h → 0 and for k ∈ [k n , k n−1 ) . Set q * 0 (x) ≡ 0.Let the function Q * n−1 (x) be the same as in (4.21), except that functions q j are replaced with q * j . Also, let Then (4.11) and implies that Also, by (4.23) and (4.24) we have for x ∈ (0, 1) Since the number δ characterizes the error in the boundary data and since η > δ is the error parameter introduced in (4.19), then, taking into account (4.24), we set In (5.2) and (5.3) F * n (x) and G * n (x) are error functions, which can be estimated as F * n L 2 (0,1) ≤ M η, G * n L 2 (0,1) ≤ M η, (5.5) where M > 0 is a constant. We also assume that Theorem 5.1. Let conditions of Theorem 4.1 hold. In procedures (i) and (iv) of the algorithm set in the QRM α = η 2 . Assume that the number k > 1 and the number k is so large that in (4.1) |O (1/k)| < 1/2 for k ≥ k for c = c * . Let the function w * x, k ∈ C 2 [0, 1] be the solution of the problem (4.9), (4.10) with the exact coefficient β * (x) = c * (x) − 1 and the exact data g * 0 k , g * 1 k . Let then for n = 1, 2..., N 1 the following estimate holds true Remark 6.1. Thus, this theorem claims that our iteratively found functions β n are located in a sufficiently small neighborhood of the exact solution β * as long as n ∈ [1, N 1 ] . Since this is achieved without any advanced knowledge of a small neighborhood of the exact solution β * , then Theorem 5.1 implies the global convergence of our algorithm, see Introduction. On the other hand, this is achieved within the framework of the approximation of subsection 4.2. Hence, to be more precise, this is the approximate global convergence property, see section 1.1.2 of [1] and section 4 of [18] for the definition of this property. Recall that the number of iterations ( N 1 in our case) can be considered sometimes as a regularization parameter in the theory of ill-posed problems [1,30].
Proof. In addition to (5.10), we will also prove that for n = 1, 2..., N 1 To simplify and shorten the proof, we assume in this proof that we work only with real valued functions. Hence, we replace in two terms of (4.23) "i" with "1" and similarly in two terms of (5.3). The case of complex valued functions is very similar. However, it contains some more purely technical details and is, therefore, more space consuming. We use the mathematical induction method. Denote Using Theorem 4.1, we obtain Hence, by (5.6) and (5.9) Following notations of section 3, denote L n (y) = y + 2k 2 n −Q n−1 + V n − 2k n y , (5.16) L n, * (y) = y + 2k 2 n − Q n−1 * + V * − 2ik n y .
Recall that we find the function w n via solving the problem (4.9), (4.10) with k = k using the QRM. Also, it follows from (2.37) and (4.18) that g 0 k − g * 0 k ≤ η and g 1 k − g * 1 k ≤ M η. Hence, we obtain similarly with (5.18) and (5.19) The function w n = w n − w * is the solution of a QRM problem, which is completely similar with (5.29). Since by (4.26) functions |β n | are uniformly bounded for all n, |β n | ≤ max (|c 0 − 1| , |c 1 − 1|) , then there exists an analog of the constant C 2 of Theorem 3.2, which estimates functions w n for all n as solutions of analogs of problems (5.29). Hence, we can assume that this constant equals M 1 . Using Theorem 3.2, (5.28) and (5.29), we obtain Hence, using (5.9), (5.30) and w 1 = w 1 + w * , we obtain

Numerical results
In all our computation, x 0 = −1 and k ∈ [0.5, 1.5]. We have observed in our computationally simulated data as well as in experimental data that the function |u(x, k)| becomes very small for k > 2. On the other hand, the largest values of |u(x, k)| were observed in some points of the interval k ∈ [0.5, 1.5]. Thus, we assign in all computations k = 1.5, k = 0.5. We have used h = 0.02. Actually in all our computations we go along the interval k ∈ [0.5, 1.5] several times. More precisely, let β (1) (x) be the result obtained in Step 3 of the above globally convergent algorithm. Set c (1) (x) = 1 + β (1) (x). Next, solve the problem the problem (4.9), (4.10) with k := k and β (x) := β (1) (x). Let the function w (1) x, k be its solution. Then we define derivatives of the new tail function V (1) 0 as in (4.25) where w n,j (x, k) is replaced with w (1) x, k . Next, we go to Step 2 and repeat. The process is repeated K = 50 times in our computational program. We choose m 0 such that Our final solution of the inverse problem is c (x) = c m 0 (x) . It is also worth mentioning that in Step 2(b)iii, after calculating β n,j , we replace it by β n,j (x) := 1 l x Ux∩(0,1) β n,j (y)dy where U x is a small neighborhood of x, x ∈ (0, 1) and l x is the length of the interval U x ∩(0, 1). We also use the truncation technique to improve the accuracy of the reconstruction function β n,j (see (4.26)).

Computationally simulated data
In this section, we show the numerical reconstruction of the spatially distributed dielectric constant from computationally simulated data. Let the function c(x) has the form Let the function u(x, k) be the solution of problem (2.3), (2.4). As mentioned in the proof of Theorem 2.1 (see (2.11)), the function u(x, k) satisfies the Lippman-Schwinger equation, This equation can be approximated as a linear system. We have solved that system numerically to computationally simulate the data.  Remark 6.1. In our computer program, we use the linear algebra package, named as Armadillo [25] to solve linear systems. The software is very helpful to speed up the program and to simplify the codes.

Experimental data
The experimental data were collected by the Forward Looking Radar which was built in the US Army Research Laboratory [23]. The device consists of two main parts. The first one, emitter, generates the time resolved electric pulses. The emitter sends out only one component of the electric field. The second part involves 16 detectors. These detectors collect the time resolved backscattering electric signal (voltage) in the time domain. The same component of the electric field is measured as the one which is generated by the emitter: see our comment about this in the beginning of section 2. The step size of time is 0.133 nanosecond. The backscattering data in the time domain are collected when the distance between the radar and the target varies from 20 to 8 meters. The average of these data with respect to both the position of the radar and those 16 detectors is the data on which we have tested our algorithm. To identify horizontal coordinates of the position of the target, Ground Positioning System (GPS) is used. The error in each of horizontal coordinates does not exceed a few centimeters. When the target is under the ground, the GPS provides the distance between the radar and a point on the ground located above the target. As to the depth of a buried target, it is not of a significant interest, since horizontal coordinates are known and it is also known that the depth does not exceed 10 centimeters. We refer to [23] for more details about the data collection process. Publications [11,18,19] contain schematic diagrams of the measurements. Our interest is in computing maximal values of dielectric constants of targets. In one target (plastic cylinder below) we compute the minimal value of its dielectric constant, since its value was less than the dielectric constant of the ground. For each target, the only information the mathematical team (MVK,LHN) had, in addition to just a single experimental curve, was whether it was located in air or below the ground.
We calculate R(x), the relative spatial dielectric constant of the whole structure including the background (air or the ground) and the target. More precisely, where D is the subinterval of the interval (0, 1) , which is occupied by the target. Here c target and c bckgr are values of the function c (x) in target and background respectively. We assume that c bckgr = const. > 0 for each set of experimental data. Hence, c bckgr = 1 if the target is located in air. The ground was dry sand. It is well known that the dielectric constant of the dry sand varies between 3 and 5 [27]. Hence, in the case of buried targets c bckgr ∈ [3,5]. In our mathematical model the time resolved electric signal u(x, x 0 , t) collected by the detectors satisfies the equations (2.14), (2.15) with c(x) being replaced by R(x), where the position x 0 of the source is actually unknown. The latter is one of the difficulties of working with these data. Thus, we set in all our tests x 0 = −1, which is the same as in [11,18,19].
Let R comp (x) be the function R(x) which we compute. Following [18], we define the computed target/background contrast as Since the dielectric constant of air equals 1, then we have R (x) ≥ 1 for targets located in air.
As to the buried targets, we have developed a procedure of the analysis of the original time resolved data, which provides us with the information on which of two cases (6.2) takes place. We refer to Case 1 and Case 2 on page 2944 of [19] for this procedure. In addition, since we had a significant mismatch of magnitudes of experimental and computationally simulated data, we have multiplied, before computations, our experimental data by the calibration number 10 −7 , see [18,19] for details of our choice of this number. There is a significant discrepancy between computationally simulated and experimental data, which was noticed in our earlier publications [11,18,19]. This discrepancy is evident from, comparison of, e.g. Figure 1b with Figure 2b and other similar ones. Therefore, to at least somehow mitigate this discrepancy, we perform a data pre-processing procedure. Besides of the Fourier transform of the time resolved data, we multiply them by a calibration factor, truncate a certain part of the data in the frequency domain and shift the data in the frequency domain, see details below.
The function u(0, x 0 , k), which we have studied in the previous sections, is the Fourier transform of the time resolved data u(0, x 0 , t). The function u(0, x 0 , k) is called "the data in the frequency domain". Observing that |u(0, x 0 , k)| is small when k belongs to a certain interval, we do not analyze u(0, x 0 , k) on that interval. Rather, we only focus on such a frequency interval which contains the major part of the information. Now, to keep the consistency with our study of computationally simulated data, we always force the frequency interval to be [0.5, 1.5]. To do so, we simply shift our data in the frequency domain: compare Figures 2b and 2c, Figures 2f and 2g, Figures 3b and 3c, Figures 3f and 3g and Figures 3j  and 3k.
We consider two cases: targets in air ( Figure 2) and targets buried about a few centimeters under the ground (Figure 3). We had experimental data for total of five (5) targets. The reconstructed dielectric constants of these targets are summarized in Table 1. In this table, computed c comp = R · c bckgr . In tables of dielectric constants of materials, their values are usually given within certain intervals [27]. Now about the intervals of the true c := c true in the 6 th column of Table 1. In the cases when targets were a wood stake and a plastic cylinder, we have taken those intervals from a published table of dielectric constants [27]. The interval of the true c for the case when the target was bush, was taken from [7]. As to the metal targets, it was established in [18] that they can be considered as such targets whose dielectric constants belong to the interval [10,30].   Figure 2: The numerical test to evaluate the "relative" dielectric constant of bush (first row) and wood stake (second row) when they are put in the air. Solid lines on b,c,f,g are real parts and dotted lines are imaginary parts.

Summary
In this paper, we have developed a frequency domain analog of the 1-d globally convergent method of [11,18,19]. We have tested this analog on both computationally simulated and time resolved experimental data. The experimental data are the same as ones used in [11,18,19]. We have modeled the process of electromagnetic waves propagation by the 1-d wave-like PDE. The reason why we have not used a 3-d model, as in, e.g. earlier works of this group on experimental data [1,28,29], is that we had only one time resolved experimental curve for each of our five targets. Our numerical method has the global convergence property. In other words, we have proven a theorem (Theorem 5.1), which claims that we obtain some points in a sufficiently small neighborhood of the exact solution without any advanced knowledge of this neighborhood. Our technique heavily relies on the Quasi Reversibility Method (QRM). The proof of the convergence of the QRM is based on a Carleman estimate. A significant modification of our technique, as compared with [11,18,19], is due to two factors. First, we use the  Fourier transform of time resolved data instead of the Laplace transform in [11,18,19]. Second, when updating tail functions via (4.25), we solve the problem (4.9), (4.10) using the QRM. On the other hand, in [11,18,19] tail functions were updated via solving the "Laplace transform analog" of the problem (2.3), (2.4) as a regular forward problem Since the dielectric constants of targets were not measured in experiments, the maximum what we can do to evaluate our results is to compare them with published data in, e.g. [27]. Results of Table 1 are close to those obtained in [11,18,19]. One can see in Table 1 that our reconstructed values of dielectric constants are well within published limits. We consider the latter as a good result. This is achieved regardless on a significant discrepancy between computationally simulated and experimental data, regardless on a quite approximate nature of our mathematical model and regardless on the presence of clutter at the data collection site. That discrepancy is still very large even after the data pre-processing, as it is evident from comparison of Figures 1b,e with Figures 2c,g and Figures 3c,g,k. Besides, the source position x 0 was unknown but rather prescribed by ourselves as x 0 = −1. Thus, our results indicate a high degree of stability of our method.