Recovery of Inhomogeneity from Output Boundary Data

: We consider the Sturm–Liouville equation on a ﬁnite interval with a real-valued integrable potential and propose a method for solving the following general inverse problem. We recover the potential from a given set of the output boundary values of a solution satisfying some known initial conditions for a set of values of the spectral parameter. Special cases of this problem include the recovery of the potential from the Weyl function, the inverse two-spectra Sturm–Liouville problem, as well as the recovery of the potential from the output boundary values of a plane wave that interacted with the potential. The method is based on the special Neumann series of Bessel functions representations for solutions of Sturm–Liouville equations. With their aid, the problem is reduced to the classical inverse Sturm–Liouville problem of recovering the potential from two spectra, which is solved again with the help of the same representations. The overall approach leads to an efﬁcient numerical algorithm for solving the inverse problem. Its numerical efﬁciency is illustrated by several examples.


Introduction
Let q ∈ L 1 (0, L) be real valued. Consider the one-dimensional Schrödinger equation − y + q(x)y = ρ 2 y, x ∈ (0, L), where ρ ∈ C is a spectral parameter, and L > 0 is finite. Let u(ρ, x) be a solution of (1), satisfying some prescribed initial conditions at the origin where a(ρ) and b(ρ) are some known functions that are not identical zeros.
In the present work, we propose a method for solving the following inverse problem.

Problem 1.
Given a(ρ k ), b(ρ k ) and l k := u(ρ k , L) for a number of values of ρ k ∈ C, find q(x).
This problem is of practical importance. For example, consider the following model. A plane wave g(ρ, x) = e −iρx , x < 0, incoming from −∞, interacts with the potential q(x) (i.e., g(ρ, x) solves (1) for x ∈ (0, L)) and is measured at the output, i.e., at the point x = L for a number of frequencies ρ k , k = 1, · · · , m. The potential q(x) is to be recovered from these output boundary data. Thus, we have a(ρ k ) = 1, b(ρ k ) = −iρ k and l k = g(ρ k , L), k = 1, · · · , m, and we need to recover q(x), of course, approximately. Note that g(ρ, x) is a Jost solution for the potential q(x) extended by zero from the segment [0, L] onto the whole axis. Therefore, the problem consists in recovering the potential q(x) from its Jost solution measured at the interface point x = L at several frequencies ρ k .
As another example, consider the problem of the recovery of the potential q(x) from the Weyl function given at a number of points. By Φ(ρ, x) we denote the Weyl solution of (1), which satisfies (1) as well as the boundary conditions Φ (ρ, 0) = 1 and Φ(ρ, L) = 0.
The Problem 1 of recovering q(x) from the data is generally not uniquely solvable even when the set of points ρ k is infinite. For example, when all ρ k are such that ρ 2 k belongs to the Dirichlet-Dirichlet spectrum of (1), the knowledge of M(ρ k ) is not sufficient for recovering q(x). In general, the characterization of such infinite sets of ρ k and of the rest of the data, for which the stated inverse problem is uniquely solvable, is a subtle matter (we refer to [5] for some results in this direction). In the present work, we restrict ourselves to the assumption that for a given infinite set of data of the form (8), the inverse problem is uniquely solvable, and we develop a method for its approximate solution. The overall approach is based on several recent results and ideas.
First of all, we use the Neumann series of Bessel functions (NSBF) representations for solutions of (1) obtained in [10]. They possess important properties, which make them especially convenient for dealing with problems, both direct and inverse, which require considering the spectral parameter ρ admitting a large range of values. In particular, the remainders of the partial sums of the NSBF representations admit estimates independent of (ρ). Roughly speaking, the partial sums approximate the exact solutions equally well for small and for large values of ρ ∈ R. Moreover, the whole information on the potential q(x) is contained already in the first coefficient of the NSBF representation; so, for recovering q(x), it is sufficient to compute the first coefficient alone. Behind these and some other remarkable features of the NSBF representations there stands the fact that they are obtained by expanding into a Fourier-Legendre series the integral kernel of the transmutation (transformation) operator. For the theory of transmutation operators, we refer to [11][12][13].
Second, we develop further the idea proposed in [14,15], which can be formulated as follows. An inverse problem can be efficiently solved by converting the input data into the computed values of the coefficients from the NSBF representations at the endpoint of the interval. This leads to the possibility of computing two different spectra for (1). The knowledge of two spectra and of the NSBF coefficients at the endpoint gives us the possibility of computing the first NSBF coefficient on the whole segment [0, L], which is sufficient for recovering q(x).
The input data of the inverse problem can be of different nature. For example, in [14] they were the spectral data of the spectral problem on a quantum graph, while in [15] the input data were several first eigenvalues from two spectra for (1). In the present work, we show how the endpoint values of a solution of (1), which satisfies some given initial conditions, serve for recovering the potential q(x), following the scheme described above.
In Section 2, we introduce the NSBF representations together with their basic properties. In Section 3, we develop the method for solving Problem 1. In Section 4, we summarize it in the form of an algorithm ready for implementation. In Section 5, we give several illustrations of its numerical performance. Finally, Section 6 contains some concluding remarks.

Preliminaries
By ϕ(ρ, x) and S(ρ, x) we denote the solutions of the equation satisfying the initial conditions The main tool used in the present work is the series representations obtained in [10] for the solutions ϕ(ρ, x) and S(ρ, x) of (9).

Theorem 1 ([10]
). The solutions ϕ(ρ, x) and S(ρ, x) admit the following series representations where j k (z) stands for the spherical Bessel function of order k (see, e.g., [16]). The coefficients g n (x) and s n (x) can be calculated following a simple recurrent integration procedure (see [10]), starting with For every ρ ∈ C, the series converges pointwise. For every x ∈ [0, L], the series converges uniformly on any compact set of the complex plane of the variable ρ, and the remainders of their partial sums admit estimates independent of (ρ).
This last feature of the series representations (the independence of the bounds for the remainders of (ρ)) is a direct consequence of the fact that the representations are obtained by expanding the integral kernels of the transmutation operators (for their theory, we refer to [11][12][13]) into Fourier-Legendre series (see [10]). It is of crucial importance for what follows. In particular, it means that for for all ρ ∈ R, where ε N (x) is a positive function tending to zero when N → ∞. That is, roughly speaking, the approximate solutions ϕ N (ρ, x) and S N (ρ, x) approximate the exact ones equally well for small and for large values of ρ. This is especially convenient when considering direct and inverse spectral problems. Moreover, for a fixed z, the numbers j k (z) rapidly decrease as k → ∞, (see, e.g., [16] [(9.1.62)]). Hence, the convergence rate of the series for any fixed ρ is, in fact, exponential. More detailed estimates for the series remainders depending on the regularity of the potential can be found in [10]. Note that formulas (12) indicate that the potential q(x) can be recovered from the first coefficients of the series (10) or (11). We have and Note that the square roots of the Dirichlet-Dirichlet eigenvalues coincide with the zeros of the function S(ρ, L): S(µ n , L) = 0, n = 1, 2, · · · , while the square roots of the Neumann-Dirichlet eigenvalues coincide with zeros of the function ϕ(ρ, L): ϕ(ρ n , L) = 0, n = 0, 1, · · · .

Solution of Problem 1
From the very beginning, we assume that the set of the numbers ρ k , for which we dispose of the data (8), is finite. Thus, the problem consists in recovering approximately the potential q(x) from the data where in general ρ k and the other data are complex numbers.

Computation of Coefficients g n (L) and s n (L)
In the first step, we use the data (16) (10) and (11) evaluated at the endpoint x = L. For this, let us substitute (10) and (11) into (2). We have for all k = 1, · · · , m. Note that if ρ k = 0, from (10) and (11), we have ϕ(0, x) = 1 + g 0 (x) and S(0, x) = x + xs 0 (x)/3, and (17) takes the form Choosing N ∈ N such that 2(N + 1) ≤ m, we consider the system of linear algebraic equations obtained from (17) by truncating the infinite sums, for k = 1, · · · , m. Note that here and below we do not look to deal with square systems of equations. In computations, a least-squares solution of an overdetermined system (provided by Matlab, which we used for computations in this work) gives very satisfactory results and allows one to make use of all available data, while keeping the number of coefficients relatively small (in practice, N = 4 or 5 may be sufficient).
Solving (18) gives us approximate values of the coefficients {g n (L)} N n=0 and {s n (L)} N n=0 . Here, we come to the problem of converting the knowledge of the coefficients from (10) and (11) evaluated at the endpoint of the interval into the information on q(x) on the whole interval.
The knowledge of the coefficients {g n (L)} N n=0 and {s n (L)} N n=0 allows us to compute the functions ϕ N (ρ, L) and S N (ρ, L) for any value of ρ. Due to the estimates in (13), we know that the accuracy of the approximation of the exact solutions ϕ(ρ, L) and S(ρ, L) by these approximate ones does not deteriorate for large values of ρ ∈ R. So, in fact, we deal with the problem of converting the knowledge of the solutions or more generally of characteristic functions of two Sturm-Liouville problems for (9) at the endpoint into the information on q(x) on the whole interval. This problem naturally arises in the framework of different methods for solving inverse Sturm-Liouville problems. For example, in [17] (see also [18]), this problem was approached via an iterative algorithm based on the use of properties of the corresponding transmutation operator. In the present work, we apply another approach, which extends the one proposed in [14,15], and consists in converting the knowledge of {g n (L)} N n=0 and {s n (L)} N n=0 into an auxiliary inverse problem for (9) consisting in the recovery of q(x) from one spectrum and a sequence of corresponding multiplier constants. We formulate this problem in the next subsection and show how the knowledge of {g n (L)} N n=0 and {s n (L)} N n=0 allows us to reduce the inverse problem to that auxiliary inverse problem.

One Spectrum and Multiplier Constants Inverse Problem
Let α 2 j ∞ j=0 be the eigenvalues of (9) with the Neumann-Dirichlet boundary conditions y (0) = 0 and y(L) = 0.

Problem 2.
Given the Neumann-Dirichlet spectrum α 2 j ∞ j=0 of (9) and the multiplier constants A method for solving the problems of this type was proposed in [19]. Here, we develop the approach proposed in [14,15], where it was applied to slightly different situations. This means that, in particular, we can compute zeros of these two functions. The following result establishes their proximity to zeros of the exact solutions ϕ(ρ, L) and S(ρ, L).

Theorem 2.
For any ε > 0, there exists such N ∈ N that all zeros of the functions ϕ(ρ, L) and S(ρ, L) are approximated by corresponding zeros of the functions ϕ N (ρ, L) and S N (ρ, L), respectively, with errors uniformly bounded by ε, and ϕ N (ρ, L), S N (ρ, L) have no other zeros.
Proof. The proof of this statement is completely analogous to the proof of Proposition 7.1 in [20] and consists in the use of (13), properties of characteristic functions of regular Sturm-Liouville problems, and Rouché's theorem from complex analysis.
Since the zeros of ϕ(ρ, L) and S(ρ, L) coincide with the square roots of the Neumann-Dirichlet and Dirichlet-Dirichlet eigenvalues of (9), the zeros of ϕ N (ρ, L) and S N (ρ, L) for a sufficiently large N give us the approximate values of these singular numbers. We denote the Dirichlet-Dirichlet eigenvalues of (9) as γ 2 j ∞ j=1 . Thus, computing the zeros of ϕ N (ρ, L) and S N (ρ, L), we obtain the approximate numbers α j ∞ j=0 and γ j ∞ j=1 . So, we reduce Problem 1 to a two-spectra inverse problem.
Analogously to solution (11), the solution ψ(ρ, x) admits the series representation where t n (x) are corresponding coefficients, analogous to s n (x) from (11). The multiplier constants can be easily calculated by recalling that ϕ(α j , 0) = 1. Thus, where we take into account that the spherical Bessel functions of odd order are odd functions. The coefficients {t n (0)} N n=0 are computed with the aid of the singular numbers γ j as follows. Since the functions ψ(γ j , x), j = 1, 2, · · · are the Dirichlet-Dirichlet eigenfunctions of (9), we have that ψ(γ j , 0) = 0; hence, This leads to a system of linear algebraic equations for computing the coefficients {t n (0)} N n=0 , which has the form where N D is the number of computed singular numbers γ j , N D ≥ N + 1. Now, having computed {t n (0)} N n=0 , we compute the multiplier constants β j N N j=0 from (22), where N N is the number of the computed singular numbers α j N N j=0 . Next, we use Equation (20) to construct a system of linear algebraic equations for the coefficients g n (x) and t n (x). Indeed, Equation (20) can be written in the form We have as many of these equations as the Neumann-Dirichlet singular numbers α j computed. For computational purposes, we choose some natural number N c -the number of the coefficients g n (x) and t n (x) to be computed. More precisely, we choose a sufficiently dense set of points x m ∈ (0, L), and at every x m , we consider the equations Solving this system of equations, we find g 0 (x m ) and consequently g 0 (x) at a sufficiently dense set of points of the interval (0, L). Finally, with the aid of (14), we compute q(x).

Algorithm
Given the data (8), the algorithm of the solution of Problem 1 consists of the following steps.

2.
Compute a number N N ≥ N of zeros α j N N j=0 of the function ϕ N (ρ, L), which approximate the Neumann-Dirichlet singular numbers, and a number N D ≥ N + 1 of zeros γ j N D j=1 of the function S N (ρ, L), which approximate the Dirichlet-Dirichlet singular numbers of q(x).
In the last step, we used the spline approximation for g 0 (x) and differentiated it twice.

Numerical Examples
The proposed method can be implemented directly using an available numeric computing environment. All the computations were performed in Matlab 2017 on an Intel i7-7600U (Intel, Santa Clara, CA, USA) equipped laptop computer. We start with the simplest example of a constant potential, which is convenient for illustrating the performance of the algorithm in each step by comparing the results with the exact ones.

Example 1.
Consider the potential q 1 (x) = c 2 , x ∈ [0, π], where c > 0 is a constant. Let us consider the problem of recovering q(x) from the Weyl function, that is, from the data of the form (4). For the constant potential, we have We In Table 1, some of the exact and computed (marked by tilde) singular numbers α j are provided together with the absolute error. Similar results were obtained for the Dirichlet-Dirichlet singular numbers, see Table 2.  Thus, the partial sums ϕ 4 (ρ, L) and S 4 (ρ, L) with the corresponding coefficients computed from (18) allow us to compute both spectra with a remarkable accuracy. We emphasize that many more of the higher indices eigenvalues can be computed with a non-deteriorating accuracy, though usually several dozens are sufficient for recovering the potential.
Next, following the algorithm, we computed the coefficients {t n (0)} 4 n=0 from (23) and β j 39 j=0 from (22). In step 5, we chose N c = 10, though the results were similar for an N c chosen between N c = 4 and N c = 20. In Figure 1, the recovered potential is depicted; the maximum absolute error was equal to 0.0012. Choosing more points ρ k and a larger number N leads to more accurate results. For example, with 25 points ρ k distributed uniformly on the same interval [0. 1,10] and N = 9, the maximum absolute error of the recovered potential was 1.09 · 10 −5 . Remarkably, the maximum absolute error of the first 40 Neumann-Dirichlet singular numbers computed was 8.6 · 10 −14 and 5.5 · 10 −14 for the Dirichlet-Dirichlet singular numbers. The numerical solution of the inverse incident plane wave problem of recovering the potential from the data (3) with the same choice of the points ρ k and N led to similarly accurate results: the maximum absolute error of the recovered potential was 1.16 · 10 −5 .
For the same constant potential, we considered the inverse problem with the data corresponding to two spectra (7) with h 1 = −0.5 and h 2 = 0. For the first ten eigenpairs given and N = 9, the potential was recovered with the maximum absolute error 1.34 · 10 −5 , while for five eigenpairs given and N = 4, the maximum absolute error was 0.0035.

Example 2.
Consider the inverse incident plane wave problem for the non-smooth and continuous potential We chose 25 points ρ k distributed uniformly on the segment in the complex plane [−1.5 + 0.1i, 1.5 + 0.1i] and the corresponding data of the type (3). In Figure 2, the result of the solution of the inverse problem is presented. The number of the coefficients was chosen as ten (i.e., N = 9). The maximum absolute error (attained at x = 1) was 0.04.
For the same potential and points ρ k , we generated the data corresponding to a(ρ) = sin ρ and b(ρ) = cos ρ, which led to similar results ( Figure 3).    In Figure 4, the result of the recovery of this potential from the data of type (3) given at 60 points ρ k distributed uniformly on the segment [0.1 + 0.1i, 15 + 0.1i] and N = 24 is shown. The maximum absolute error was approximately 0.07.

Conclusions
A simple method for solving a general inverse problem for the Sturm-Liouville equation, consisting in recovering the potential from the output boundary values of a solution satisfying certain prescribed initial conditions, was developed. The data of the problem were converted into the knowledge of the coefficients of the Neumann series of Bessel functions, representations of two linearly independent solutions evaluated at the endpoint, which allowed us to reduce the problem to a two-spectra inverse Sturm-Liouville problem. This problem was solved by reducing it to a system of linear algebraic equations from which the first coefficient of the Neumann series of Bessel functions representation was found, which led to the recovery of the potential.
The method is simple, direct, accurate, and applicable to a large variety of inverse problems. Its performance was illustrated by numerical examples.

Data Availability Statement:
The data that support the findings of this study are available upon reasonable request.

Conflicts of Interest:
The authors declare no competing interests.