Approximate recovery of the Sturm–Liouville problem on a half-line from the Weyl function

2.1 Definitions and inverse problem setting

Let q ⁢ ( x ) be a real-valued function, Lebesgue-measurable and such that

We consider the Schrödinger equation

where λ = ρ 2 , ρ ∈ ℂ + ¯ is a spectral parameter and ℂ + := { ρ ∈ ℂ : Im ⁡ ρ > 0 } . The solution y ⁢ ( x ) is sought to satisfy the boundary condition

where h is some fixed real number. Under the conditions imposed, the Sturm–Liouville problem (2.2)–(2.3) possesses at most a finite set of simple eigenvalues, which are the values ρ 2 , which admit nontrivial solutions of (2.2)–(2.3) belonging to ℒ 2 ⁢ ( 0 , ∞ ) . The eigenvalues, if they exist, are all negative (see, e.g., [23, 24, 26, 31]). By φ ⁢ ( ρ , x ) and S ⁢ ( ρ , x ) we denote the solutions of (2.2) satisfying the initial conditions

and

For all ρ ∈ ℂ + ¯ , equation (2.2) possesses a unique solution e ⁢ ( ρ , x ) (the Jost solution) such that for ν = 0 , 1 the following asymptotic relations hold:

uniformly in ℂ + ¯ (see, e.g., [9, 26, 31]). Thus, for every ρ ∈ ℂ + fixed, the function e ⁢ ( ρ , x ) belongs to ℒ 2 ⁢ ( 0 , ∞ ) . When Im ⁡ ρ = 0 , it is only bounded, and e ⁢ ( - ρ , x ) = e ⁢ ( ρ , x ) ¯ . For ρ = i ⁢ τ , τ > 0 , the Jost solution e ⁢ ( ρ , x ) is real valued. Moreover, for large values of | ρ | the following asymptotic relation holds:

uniformly with respect to x ≥ 0 (see [31, p. 129]), where

For all x ∈ [ 0 , ∞ ) fixed, both e ⁢ ( ρ , x ) and e ′ ⁢ ( ρ , x ) are analytic in ℂ + and continuous in ℂ + ¯ (see [31, p. 206]).

Denote

and

The function Φ ⁢ ( ρ , x ) (called the Weyl solution) satisfies equation (2.2), as well as the conditions

Note that Φ ⁢ ( ρ , x ) is defined uniquely by (2.2) and (2.9).

Denote

We call M ⁢ ( ρ ) the Weyl function. Note that for ρ ∈ ℝ , M ⁢ ( - ρ ) = M ⁢ ( ρ ) ¯ .

The singularities of M ⁢ ( ρ ) , ρ ∈ ℂ + , if they exist, coincide with the square roots of the eigenvalues of (2.2)–(2.3), while zeros of M ⁢ ( ρ ) , ρ ∈ ℂ + , if they exist, coincide with the square roots of the eigenvalues of (2.2) with the Dirichlet boundary condition

In both cases the eigenvalues are strictly negative numbers.

We have an obvious relation, which is a consequence of (2.4), (2.5) and (2.8),

We are interested in the following inverse problem.

A fairly complete theory of this problem can be found in [31]. The main result of the present work is a practical method for the reconstruction of q ⁢ ( x ) and h from the knowledge of M ⁢ ( ρ ) for ρ ≥ 0 and of its zeros. The problem is of interest in many applications. Depending on a concrete physical model, the Weyl function, for example, may represent a ratio of the components of an electromagnetic field inciding perpendicularly onto a layered medium (see [26, p. 3]) or be defined by a reflection coefficient in a scattering problem (see below). Moreover, M ⁢ ( ρ ) is closely related to the spectral function and to the scattering data in such a way that it can be computed from each of these data (see the work of A.G. Ramm who considers the I-function, which is 1 / M ⁢ ( ρ ) for h = 0 [26, 27] and the book of V. A. Yurko [31]). For solving (IP1), V. A. Yurko developed a method of spectral mappings which allows one to reduce the problem to the solution of an integral equation involving integration along an infinite contour in a complex plane. To our best knowledge the method has not been realized numerically.

2.2 Relation of (IP1) to inverse scattering problem

Note that M ⁢ ( ρ ) is closely related to the reflection coefficient. Consider equation (2.2) on a whole axis ( - ∞ , ∞ ) with a potential q ⁢ ( x ) , which is identical zero for x < 0 , and satisfies (2.1). Then for a plane wave e i ⁢ ρ ⁢ x incoming from - ∞ , the solution of the scattering problem has the form

where R ⁢ ( ρ ) and T ⁢ ( ρ ) are the reflection and transmission coefficients, respectively. Additionally, the continuity of u ⁢ ( ρ , x ) and u ′ ⁢ ( ρ , x ) at the origin leads to the relations

and

Thus,

and hence

Thus, if h is known, the problem of the reconstruction of q ⁢ ( x ) from the Weyl function is equivalent to the reconstruction from the reflection coefficient. The uniqueness result for this problem was obtained in [2].

In particular, when h = 0 , we have

4.1 Reconstruction of e ⁢ ( ρ , 0 ) and Δ ⁢ ( ρ )

Given M ⁢ ( ρ ) for ρ ≥ 0 and the location of its zeros (h is unknown), in a first step, we are interested in recovering e ⁢ ( ρ , 0 ) and Δ ⁢ ( ρ ) . Denote

and

Here ρ k = i ⁢ τ k , τ k > 0 , denote zeros of M ⁢ ( ρ ) , if they exist.

Thus, e ⁢ ( ρ , 0 ) is recovered from the knowledge of M ⁢ ( ρ ) for ρ ≥ 0 and of the location of its zeros. In its turn, Δ ⁢ ( ρ ) is obtained from

4.2 Reconstruction of q ⁢ ( x ) and h

For recovering q ⁢ ( x ) and h from e ⁢ ( ρ , 0 ) and Δ ⁢ ( ρ ) , we make use of identity (2.12). It can be written as

Substituting series representations (3.1), (3.2) and (3.5), we obtain the equality

Here e ⁢ ( ρ , 0 ) and Δ ⁢ ( ρ ) are known, computed in the first step.

Now, let us fix x ≥ 0 and consider a sufficiently large number K of the points ρ k ∈ ℂ + ¯ . To compute approximate values of a number of the coefficients { g n ⁢ ( x ) , s n ⁢ ( x ) } n = 0 N 1 , { a n ⁢ ( x ) } n = 0 N 2 , from (4.10) we obtain a finite system of linear algebraic equations

Note that the number of equations K can be taken arbitrarily large, while the number of the unknowns ( 2 ⁢ ( N 1 + 1 ) + N 2 + 1 ), due to estimates (3.3), (3.7), (3.8) can be kept relatively small, especially for a convenient choice of the points ρ k belonging to a strip 0 ≤ Im ⁡ ρ ≤ C , with C > 0 being not too large. Solving this system for a sufficiently dense set of points x j ∈ ( 0 , b ] on an arbitrary interval of interest ( 0 , b ] ⊂ ( 0 , ∞ ) , we find g 0 ⁢ ( x ) in [ 0 , b ] ( g 0 ⁢ ( 0 ) = 0 ). Now, due to the first equality from (3.4), we have that φ ⁢ ( 0 , x ) = g 0 ⁢ ( x ) + 1 , and thus,

Example 1.

Consider the potential

From [3, Example 6.1] we have

and hence for h = 0 we obtain that

and

It is easy to see that no eigenvalue exists, and M ⁢ ( 0 ) ≠ 0 .

In the first step we recovered e ⁢ ( ρ , 0 ) and Δ ⁢ ( ρ ) from M ⁢ ( ρ ) numerically, by using (4.1). For this, the integral from (4.1) was computed by the modified 6 point Newton–Cottes integration formula (see [22] for details) on a large interval ( - T , T ) with T = 3000 and for K = 501 points ρ k with Im ⁡ ρ k = 1 and Re ⁡ ρ k distributed logarithmically equally spaced on the segment [ 0.001 ,  1001 ] . That is, Re ⁡ ρ k = 10 α k with α k , k = 1 , … , K , being distributed uniformly on [ log ⁡ ( 0.001 ) , log ⁡ ( 1001 ) ] ≈ [ - 3 ,  3 ] . Thus, e ⁢ ( ρ , 0 ) was computed with the aid of (4.1), and then Δ ⁢ ( ρ ) was computed by (4.8).

In the second step, system (4.11) was considered at 251 points x ∈ [ 0 , π ] and for the same set of the points ρ k . Note that, since all the unknown { g n ⁢ ( x ) , s n ⁢ ( x ) } n = 0 N 1 , { a n ⁢ ( x ) } n = 0 N 2 are real, for each ρ k , equation (4.11) gives us two equations: with real and imaginary parts of matrix coefficients. Finally, numbers N 1 and N 2 were chosen automatically. For every x j ∈ [ 0 , π ] we compute the rank of each of the three thirds of the matrix of the system (matrix consisting of the first N 1 + 1 columns, matrix consisting of the columns N 1 + 2 , … , 2 ⁢ N 1 + 2 , and matrix consisting of the columns 2 ⁢ N 1 + 3 , … , 2 ⁢ N 1 + N 2 + 3 ), which is the number of singular values of the matrix that are larger than a default tolerance (which was chosen as 10 - 8 ). Then the number of the columns remained in each third of the matrix is taken equal to the corresponding rank, this means that for different points x j , a different number of the unknown NSBF coefficients is computed, so that, in fact, system (4.11) is considered in the form