A practical method for recovering Sturm-Liouville problems from the Weyl function
aa r X i v : . [ m a t h . C A ] J a n A practical method for recovering Sturm-Liouville problems from theWeyl function
Vladislav V. Kravchenko , Sergii M. Torba Departamento de Matem´aticas, Cinvestav, Unidad Quer´etaro,Libramiento Norponiente
January 25, 2021
Abstract
In the paper we propose a direct method for recovering the Sturm-Liouville potential from theWeyl-Titchmarsh m -function given on a countable set of points. We show that using the Fourier-Legendre series expansion of the transmutation operator integral kernel the problem reduces to aninfinite linear system of equations, which is uniquely solvable if so is the original problem. Thesolution of this linear system allows one to reconstruct the characteristic determinant and henceto obtain the eigenvalues as its zeros and to compute the corresponding norming constants. Asa result, the original inverse problem is transformed to an inverse problem with a given spectraldensity function, for which the direct method of solution from [24] is applied.The proposed method leads to an efficient numerical algorithm for solving a variety of inverseproblems. In particular, the problems in which two spectra or some parts of three or more spectraare given, the problems in which the eigenvalues depend on a variable boundary parameter (includ-ing spectral parameter dependent boundary conditions), problems with a partially known potentialand partial inverse problems on quantum graphs. We consider the inverse problem of recovering the potential of a Sturm-Liouville equation on a finiteinterval together with the boundary conditions from values of the Weyl-Titchmarsh m -function on agiven set of points.A fundamental result of Marchenko [26] states that the Weyl-Titchmarsh function m ( z ) determinesuniquely the potential q ∈ L (0 , π ) as well as the boundary conditions. So some particular inverseproblem is solved, at least in theory, if one finds a way to convert it to a corresponding problem for the m -function. However, several difficulties arise. For some inverse spectral problems the given spectraldata allow one to obtain the values of the m -function at any point z ∈ C . This is the case of theclassical inverse problem with a given spectral density function, formula (2.13) below provides the m -function. While for many problems only the values of the m -function on a countable set of points { z k } ∞ k =0 can be obtained. See Subsection 2.2 and references therein for some examples. Thus one facestwo principle questions. Is it possible to recover uniquely the m -function from the given values m ( z k ), k = 0 , , . . . ? If yes, how can that be done?Since the m -function is meromorphic, the unique recovery is always possible if the sequence { z k } ∞ k =0 contains a subsequence z n k converging to a finite limit and such that the values m ( z n k ) are bounded.This is usually not the case when we start with some inverse problem, the points z k are real and go toinfinity when k → ∞ . Suppose all the points z k are real. A simple sufficient condition for the unique1ecoverability of the m -function from the given values m ( z k ) was given in [9], and a complete answerwas given not so long ago in [17]. Simply speaking, the points z k should be distributed sufficientlydensely, like the points obtained as a union of two spectra.The corresponding results are obtained by using infinite products and analytic continuation. Andas it is stated in the following quote [7, p.106],It is one of the guiding principles of computational analysis that such proofs often lead toelegant mathematics but are useless as far as constructive methods are concerned.The only constructive method for finding the m -function we are aware of is that from [4] (basedon the ideas of [30] and [7, Chapter 3]). This algorithm recovers a certain combination of the integralkernel of the transmutation operator with its derivative by using non-harmonic Fourier series. Howeverit requires additionally a beforehand knowledge of the parameter ω (depending on the average of thepotential q and boundary parameter H ) and numerical realization of this method is expected toconverge slowly (no numerical illustration is provided in [4]).We use Fourier-Legendre series to represent the transmutation integral kernel and its derivative.This idea was originally proposed in [22] and resulted in a Neumann series of Bessel functions (NSBF)representation for the solutions and their derivatives of Sturm-Liouville equations. As a result, aconstructive algorithm for recovering the m -function from its values on a countable set of points isproposed. The problem is reduced to an infinite system of linear equations and the unique solvabilityof this system is proved. We would like to emphasize that no additional information is necessarybeforehand.Moreover, the proposed method leads to an efficient numerical algorithm. Since the NSBF repre-sentation possesses such a remarkable property as a uniform error bound for all real λ >
0, one canobtain any finite set of approximate eigenvalues and corresponding norming constants with a non-deteriorating accuracy. That is, an efficient conversion of the original problem to an inverse problemby a given spectral density function is proposed. The method of solution introduced in [20], [21]and refined in [24] is adapted to recover the potential and the boundary conditions from the spectraldensity function.For classical inverse spectral problems there are many results describing to what extent the po-tential can be recovered from a finite set of spectral data, see, e.g., [14], [28], [33], [31]. For the Weylfunction, on the contrary, we are not aware of any single result. It can be seen that the values of z k , 0 ≤ k ≤ K can be chosen neither from a too small interval nor too sparse. Indeed, taking allthe values z k from one spectrum is insufficient to recover the potential. On the other hand, as wasmentioned in [7, Subsections 3.10.5 and 3.14], when all the values z k are equal to the first eigenvalue λ (for problems with different boundary conditions), the corresponding numerical problem is extremelyill-conditioned. So the proposed algorithm struggles as well when the points z k either belong to asmall interval or are too sparse. Nevertheless, it delivers excellent results for smooth potentials andsufficiently dense points z k , does not require any additional information on the potential or boundaryconditions and can be applied to a variety of inverse problems for which not so many (if any) spe-cialized algorithms are known. So we believe the presented method will be a starting point for thefuture research of the question “to which extent a potential can be recovered from a given finite setof values of the Weyl-Titchmarsh m -function”, as well as of a new class of numerical methods basedon the Weyl-Titchmarsh m -function.The paper is organized as follows. In Section 2 we provide necessary information on the Weylfunction, transmutation operators and Neumann series of Bessel functions representations. In Section3 we show how the inverse problem of recovering the Sturm-Liouville problem from the Weyl functionis reduced to an infinite system of linear algebraic equations, prove its unique solvability and showthe reduction of the original problem to the classical problem of the recovery of a Sturm-Liouville2roblem from a given spectral density function; an algorithm is provided. Finally in Section 4 wepresent several numerical examples illustrating the efficiency and some limitations of the proposedmethod. Let q ∈ L (0 , π ) be real valued. Consider the Sturm-Liouville equation − y ′′ + q ( x ) y = ρ y, x ∈ (0 , π ) , (2.1)where ρ ∈ C , with the boundary conditions y ′ (0) − hy (0) = 0 , y ′ ( π ) + Hy ( π ) = 0 , (2.2)where h and H are arbitrary real constants, and the boundary value problem y ′ (0) − hy (0) = 1 , y ′ ( π ) + Hy ( π ) = 0 . (2.3)If λ = ρ is not an eigenvalue of the spectral problem (2.1), (2.2), then the boundary value problem(2.1), (2.3) possesses a unique solution, which we denote by Φ( ρ, x ). Let M ( λ ) := Φ( ρ, . The functions Φ( ρ, x ) and M ( λ ) are called the Weyl solution and the Weyl function, respectively. See[37, Section 1.2.4] for additional details.By ϕ ( ρ, x ) and S ( ρ, x ) we denote the solutions of (2.1) satisfying the initial conditions ϕ ( ρ,
0) = 1 and ϕ ′ ( ρ,
0) = h (2.4)and S ( ρ,
0) = 0 and S ′ ( ρ,
0) = 1 , (2.5)respectively. Denote additionally for λ = ρ ∆( λ ) := ϕ ′ ( ρ, π ) + Hϕ ( ρ, π ) and ∆ ( λ ) := S ′ ( ρ, π ) + HS ( ρ, π ) . (2.6)Obviously, for all ρ ∈ C the function ϕ ( ρ, x ) fulfills the first boundary condition, ϕ ′ ( ρ, − hϕ ( ρ,
0) = 0, and thus, the spectrum of problem (2.1), (2.2) is the sequence of numbers (cid:8) λ n = ρ n (cid:9) ∞ n =0 such that∆( λ n ) = 0 . It is easy to see thatΦ( ρ, x ) = S ( ρ, x ) + M ( λ ) ϕ ( ρ, x ) (2.7)and M ( λ ) = − ∆ ( λ )∆( λ ) . (2.8)Moreover, the Weyl function is meromorphic with simple poles at λ = λ n , n ≥ emark 2.1. There is another common definition of the Weyl-Titchmarsh m -function [35], [10],[17]. Let v ( ρ, x ) be the solution of (2.1) satisfying v ( ρ, π ) = 1 and v ′ ( ρ, π ) = − H. (2.9) Then the Weyl-Titchmarsh m -function is defined as m ( ρ ) = v ′ ( ρ, v ( ρ, . (2.10) One can verify that m ( ρ ) = h − ∆( ρ )∆ ( ρ ) = h + 1 M ( ρ ) , (2.11) that is, the functions m and M can be easily transformed one into another. Denote additionally α n := Z π ϕ ( ρ n , x ) dx. (2.12)The set { α n } ∞ n =0 is referred to as the sequence of norming constants of problem (2.1), (2.2). Theorem 2.2 (see [37, Theorem 1.2.6]) . The following representation holds M ( λ ) = ∞ X n =0 α n ( λ − λ n ) . (2.13)Let us formulate the inverse Sturm-Liouville problem. Problem 2.3 (Recovery of a Sturm-Liouville problem from its Weyl function) . Given the Weyl func-tion M ( λ ) , find a real valued q ∈ L (0 , π ) , and the constants h , H ∈ R , such that M ( λ ) be the Weylfunction of problem (2.1) , (2.2) . We refer the reader to [37, Theorem 1.2.7] for the result on the unique solvability of Problem 2.3.Having the Weyl function known completely, one can recover the sequence { λ n } ∞ n =0 as the sequenceof its poles, and the sequence { α n } ∞ n =0 can be immediately obtained from the corresponding residuals,the formula (1.2.14) from [37] states that α n = 1Res λ = λ n M ( λ ) = − ( λ n ) · ddλ ∆( λ ) (cid:12)(cid:12)(cid:12)(cid:12) λ = λ n . (2.14)Hence reducing the inverse problem to recovering the Sturm-Liouville problem by its spectral densityfunction. Or, by taking the sequences of poles and zeros of the Weyl function one can reduce Problem2.3 to recovering the Sturm-Liouville problem by two spectra, see [37, Remark 1.2.3].However the situation is more complicated if one looks for a practical method of recovering q , h and H from the Weyl function given only on a countable (even worse, finite) set of points from some(small) interval. The direct reduction to the inverse problems considered, in particular, in [24], is notpossible. Let us formulate the corresponding problem and state some sufficient results for its uniquesolvability. 4 roblem 2.4 (Recovery of a Sturm-Liouville problem from its Weyl function given on a countableset of points) . Given the values M n = M ( z n ) , M n ∈ R ∪ {∞} , of the Weyl function M ( λ ) on a setof points { z n } ∞ n =0 , find a real valued q ∈ L (0 , π ) , and the constants h , H ∈ R , such that the Weylfunction of problem (2.1) , (2.2) takes the values M n at the points z n , n ≥ . It is easy to see that Problem 2.4 is not always uniquely solvable. The following result shows thatchoosing the set of points z n “sufficiently dense” can guarantee the unique solvability. Theorem 2.5 ([9, Theorem 2.1]) . Let a + = max( a, for a ∈ R . Let { z n } ∞ n =0 be a sequence of distinctpositive real numbers satisfying ∞ X n =1 (cid:0) z n − n (cid:1) + n < ∞ . (2.15) Let m and m be Weyl-Titchmarsh m -functions for two Sturm-Liouville problems having the potentials q and q and the constants H and H , correspondingly. Suppose that m ( z n ) = m ( z n ) for all n ∈ N .Then m = m (and hence q = q a.e. on [0 , π ] and H = H ). See also [17, Theorem 1.5] for the necessary and sufficient result.
Below we briefly present some inverse problems, omitting any additional details like the conditions ona spectra ensuring the uniqueness of the solution. They can be consulted in the provided references.
Consider the second set of boundary conditions y (0) = 0 , y ′ ( π ) + Hy ( π ) = 0 . (2.16)Let us denote the spectrum of problem (2.1), (2.16) as (cid:8) ν n = µ n (cid:9) ∞ n =0 .The inverse problem consists in recovering q , h and H from given two spectra, { λ n } ∞ n =0 and { ν n } ∞ n =0 . Note that this problem coincides with the one considered in [24] if one changes x by π − x and h ↔ H . Now the equalities M ( λ n ) = ∞ and M ( ν n ) = 0 , can be easily checked and give the reduction of the two spectra inverse problem to Problem 2.4. There is a class of inverse problems where it is asked to recover the potential from given parts of threespectra [9], four spectra [29], N spectra [16]. All these problems can be regarded as special cases ofthe following inverse problem. Suppose the sequences of numbers { ˜ λ n } ∞ n =0 and { h n } ∞ n =0 are given andsuch that ˜ λ n is an eigenvalue (does not refer to n -th indexed eigenvalue) of the problem (2.1) with theboundary conditions y ′ (0) − h n y (0) = 0 , y ′ ( π ) + Hy ( π ) = 0 . (2.17)The problem is to recover the potential q . We refer the reader to [7, § § m (˜ λ n ) = h n . .2.3 Partially known potential The Hochstadt-Lieberman inverse problem from [15] reads as follows: suppose we are given the spec-trum of the problem (2.1), (2.2) and the potential q is known on the segment [0 , π/ q onthe whole segment [0 , π ]. See also [11], [9] for the case when the potential is known on a differentportion of the segment. If the potential is known on less than half of the segment, more than onespectrum is necessary to recover the potential uniquely. For that reason the general formulation ofthe problem is as follows. Suppose a ∈ (0 , π ) and the potential q is known on [0 , a ]. Moreover, twosequences of numbers are given as in Subsubsection 2.2.2. Recover q on the whole [0 , π ].Let m ( x, ρ ) define the m -function on [ x, π ] with the same boundary condition (2.9), i.e., m ( x, ρ ) = v ′ ( ρ, x ) v ( ρ, x ) . Consider an eigenvalue ˜ λ n = ρ n and the corresponding boundary constant h n . Then m (0 , ˜ λ n ) = m (˜ λ n ) = h n , hence we may consider the following Cauchy problem for equation (2.1) v ( ρ n ,
0) = 1 , v ′ ( ρ n ,
0) = h. Solving it on [0 , a ] (in practice it can be done efficiently for a large set of ˜ λ n using the methoddeveloped in [22]), we obtain values v ( ρ n , a ) and v ′ ( ρ n , a ), hence, m ( a, ˜ λ n ). And by rescaling wereduce the problem of a partially known potential to Problem 2.4. The following spectral parameter dependent boundary condition was considered in [4], [5], [36] f ( λ ) y ′ (0) + f ( λ ) y (0) = 0 , (2.18)where f and f are entire functions, not vanishing simultaneously. The second boundary conditionis the same as in (2.2). The functions f and f are supposed to be known, and an inverse problemconsists in recovering the potential q and the constant H by the given parameter ω (see (2.33)),spectra λ n and corresponding multiplicities m n . See also [34] for a general spectral theory of suchproblems.Supposing for simplicity that all the eigenvalues are simple, i.e., m n = 1 for all n , the problemimmediately reduces to Problem 2.4, since m ( λ n ) = − f ( λ n ) f ( λ n ) . Many inverse problems can be reduced to a problem having boundary condition (2.18). In partic-ular, the Hochstadt-Lieberman problem, the inverse transmission eigenvalue problem, partial inverseproblems for quantum graphs. See [36] and references therein for more details.
The solutions ϕ ( ρ, x ) and S ( ρ, x ) admit the integral representations (see, e.g., [8], [25], [26], [27,Chapter 1], [32], [6]) ϕ ( ρ, x ) = cos ρx + Z x − x K ( x, t ) cos ρt dt, (2.19) S ( ρ, x ) = sin ρxρ + Z x − x K ( x, t ) sin ρtρ dt, (2.20)6here the integral kernel K is a continuous function of both arguments in the domain 0 ≤ | t | ≤ x ≤ π and satisfies the equalities K ( x, x ) = h Z x q ( t ) dt and K ( x, − x ) = h . (2.21)It is of crucial importance that K ( x, t ) is independent of ρ .The following result from [22] will be used. Theorem 2.6 ([22]) . The integral transmutation kernel K ( x, t ) and its derivative K ( x, t ) := ∂∂x K ( x, t ) admit the following Fourier-Legendre series representations K ( x, t ) = ∞ X n =0 β n ( x ) x P n (cid:18) tx (cid:19) , (2.22) K ( x, t ) = ∞ X j =0 γ j ( x ) x P j (cid:18) tx (cid:19) , < | t | ≤ x ≤ π, (2.23) where P k stands for the Legendre polynomial of order k .For every x ∈ (0 , π ] the series converge in the norm of L ( − x, x ) . The first coefficients β ( x ) and γ ( x ) have the form β ( x ) = ϕ (0 , x ) − , γ ( x ) = β ′ ( x ) − h − Z x q ( s ) ds, (2.24) and the rest of the coefficients can be calculated following a simple recurrent integration procedure. Corollary 2.7.
Equality (2.24) shows that the potential q ( x ) can be recovered from β ( x ) . Indeed, q ( x ) = ϕ ′′ (0 , x ) ϕ (0 , x ) = 2 β ′′ ( x )2 β ( x ) + 1 . (2.25) Moreover, the constant h is also recovered directly from β ( x ) , since h = 2 β ′ (0) . The following series representations for the solutions ϕ ( ρ, x ) and S ( ρ, x ) and for their derivatives wereobtained in [22]. Theorem 2.8 ([22]) . The solutions ϕ ( ρ, x ) and S ( ρ, x ) and their derivatives with respect to x admitthe following series representations ϕ ( ρ, x ) = cos ρx + 2 ∞ X n =0 ( − n β n ( x ) j n ( ρx ) , (2.26) S ( ρ, x ) = sin ρxρ + 2 ρ ∞ X n =0 ( − n β n +1 ( x ) j n +1 ( ρx ) , (2.27) ϕ ′ ( ρ, x ) = − ρ sin ρx + (cid:18) h + 12 Z x q ( s ) ds (cid:19) cos ρx + 2 ∞ X n =0 ( − n γ n ( x ) j n ( ρx ) , (2.28) S ′ ( ρ, x ) = cos ρx + 12 ρ (cid:18)Z x q ( s ) ds (cid:19) sin ρx + 2 ρ ∞ X n =0 ( − n γ n +1 ( x ) j n +1 ( ρx ) , (2.29)7 here j k ( z ) stands for the spherical Bessel function of order k (see, e.g., [1]). The coefficients β n ( x ) and γ n ( x ) are those from Theorem 2.6. For every ρ ∈ C all the series converge pointwise. For every x ∈ [0 , π ] the series converge uniformly on any compact set of the complex plane of the variable ρ , andthe remainders of their partial sums admit estimates independent of Re ρ .Moreover, the representations can be differentiated with respect to ρ resulting in ϕ ′ ρ ( ρ, x ) = − x sin ρx + 2 ∞ X n =0 ( − n β n ( x ) (cid:18) nρ j n ( ρx ) − xj n +1 ( ρx ) (cid:19) , (2.30) ϕ ′′ x,ρ ( ρ, x ) = − (cid:18) hx + x Z x q ( s ) ds (cid:19) sin ρx − ρx cos ρx + 2 ∞ X n =0 ( − n γ n ( x ) (cid:18) nρ j n ( ρx ) − xj n +1 ( ρx ) (cid:19) . (2.31)We refer the reader to [22] for the proof and additional details related to this result. The coefficients β n and γ n decay when n → ∞ , and the decay rate is faster for smoother potentials. Namely, thefollowing result is valid. Proposition 2.9 ([23]) . Let q ∈ W p ∞ [0 , π ] for some p ∈ N , i.e., the potential q is p times differentiablewith the last derivative being bounded on [0 , π ] . Then there exist constants c and d , independent of N ,such that | β N ( x ) | ≤ cx p +2 ( N − p +1 / and | γ N ( x ) | ≤ dx p +1 ( N − p − / , N ≥ p + 1 . Let us denote h n := γ n ( π ) + Hβ n ( π ) , n ≥ ω := h + H + 12 Z π q ( t ) dt and ω := ω − h = H + 12 Z π q ( t ) dt. (2.33) Corollary 2.10.
The following representations hold for the characteristic functions ∆ and ∆ ∆( ρ ) = − ρ sin ρπ + ω cos ρπ + 2 ∞ X n =0 ( − n h n j n ( ρπ ) , (2.34)∆ ( ρ ) = cos ρπ + ω sin ρπρ + 2 ρ ∞ X n =0 ( − n h n +1 j n +1 ( ρπ ) , (2.35) and for the derivative with respect to ρddρ ∆( ρ ) = − (1 + πω ) sin ρπ − πρ cos ρπ + 2 ∞ X n =0 ( − n h n (cid:18) nρ j n ( ρπ ) − πj n +1 ( ρπ ) (cid:19) . (2.36) Let G ( x, t ) := K ( x, t ) + K ( x, − t ) = ∞ X n =0 β n ( x ) x P n (cid:18) tx (cid:19) F ( x, t ) = ∞ X n =0 (cid:18) cos ρ n x cos ρ n tα n − cos nx cos ntα n (cid:19) , ≤ t, x < π (2.37)where α n = ( π/ , n > ,π, n = 0 . Then the function G is the unique solution of the Gelfand-Levitan equation, see, e.g., [37, Theorem1.3.1] G ( x, t ) + F ( x, t ) + Z x F ( t, s ) G ( x, s ) ds = 0 , < t < x < π. (2.38)The series (2.37) converges slowly and possesses a jump discontinuity (for ω = 0) at x = t = π .To overcome these difficulties, in [19] another representation for the function F was derived. Namely, F ( x, t ) = ∞ X n =1 (cid:18) cos ρ n x cos ρ n tα n − cos nx cos ntα n + 2 ωπ n (cid:16) x sin nx cos nt + t sin nt cos nx (cid:17)(cid:19) + cos ρ x cos ρ tα − π − ωπ (cid:0) π max { x, t } − x − t (cid:1) . (2.39) ω and ω and of the characteristic functions ∆ and ∆ Let us rewrite (2.8) as M ( ρ )∆( ρ ) + ∆ ( ρ ) = 0and substitute (2.34) and (2.35) into the last expression. We obtain M ( ρ ) − ρ sin ρπ + ω cos ρπ + 2 ∞ X n =0 ( − n h n j n ( ρπ ) ! + cos ρπ + ω sin ρπρ + 2 ρ ∞ X n =0 ( − n h n +1 j n +1 ( ρπ ) = 0 , or ω · M ( ρ ) cos ρπ + ω · sin ρπρ + 2 ∞ X n =0 h n · ( − n M ( ρ ) j n ( ρπ ) + 2 ∞ X n =0 h n +1 · ( − n j n +1 ( ρπ ) ρ = M ( ρ ) ρ sin ρπ − cos ρπ. (3.1)For the case when M ( ρ ) = ∞ , the last equality reduces to ω cos ρπ + 2 ∞ X n =0 h n · ( − n j n ( ρπ ) = ρ sin ρπ. (3.2)9uppose that the Weyl function is known on a countable set of points { z k } ∞ k =0 = { ˜ ρ k } ∞ k =0 , thatis, let the values M k := M ( z k ), M k ∈ R ∪ {∞} be given. Then considering the equality (3.1) for all ρ = ˜ ρ k we obtain an infinite system of linear algebraic equations for the unknown ω , ω and { h n } ∞ n =0 .One can easily see from (2.22) and (2.23) that the functions √ πh n √ n +1 are the Fourier coefficients of thesquare integrable function K ( π, t ) + HK ( π, t ) in the space L ( − π, π ) with respect to the orthonormalbasis { p n ( t ) } ∞ n =0 := (cid:26) √ n + 1 P n ( t/π ) √ π (cid:27) ∞ n =0 . (3.3)Let us introduce new unknowns ξ n = √ πh n √ n + 1 , n ≥ . (3.4)Then { ξ n } ∞ n =0 ∈ ℓ (as a sequence of the Fourier coefficients). Let us rewrite the infinite linear systemof equations as ∞ X n =0 ξ n · ( − n √ n + 2 M k j n (˜ ρ k π ) √ π + ∞ X n =0 ξ n +1 · ( − n √ n + 2 j n +1 (˜ ρ k π ) √ π ˜ ρ k + ω · M k cos ˜ ρ k π + ω · sin ˜ ρ k π ˜ ρ k = M k ˜ ρ n sin ˜ ρ k π − cos ˜ ρ k π, k ∈ N (3.5)(equation ∞ X n =0 ξ n · ( − n √ n + 2 j n (˜ ρ k π ) √ π + ω · cos ˜ ρ k π = ˜ ρ k sin ˜ ρ k π (3.6)should be used if for some k we have M k = ∞ ). Theorem 3.1.
Suppose that the original Problem 2.4 is solvable and the numbers { z k } ∞ k =0 satisfy thecondition (2.15) . Then the infinite system (3.5) possesses a unique ℓ solution { ω, ω , { ξ n } ∞ n =0 } .Proof. Since the original Problem 2.4 is solvable, we can recover the Weyl funcion M and by [37,Theorem 1.2.7] recover the potential q and the constants h and H . From (2.33) we obtain ω and ω , and by Theorem 2.6 we obtain the square integrable with respect to the second variable integralkernel K and its derivative K , as well as the series representations (2.22), (2.23). The coefficients β n and γ n give us the sequence { ξ n } ∞ n =0 ∈ ℓ . One can verify using (2.34) and (2.35) that the obtainednumbers ω , ω and ξ n , n ≥ ℓ solution n ˜ ω, ˜ ω , { ˜ ξ n } ∞ n =0 o of the system (3.5). Let ˜ h n = √ n +1 √ π ˜ ξ n , n ≥
0. Consider the following functions˜ G ( t ) = 2 ∞ X n =0 ˜ h n π P n (cid:18) tπ (cid:19) , (3.7)˜ S ( t ) = 2 ∞ X n =0 ˜ h n +1 π P n +1 (cid:18) tπ (cid:19) , t ∈ [0 , π ] . (3.8)Due to condition { ˜ ξ n } ∞ n =0 ∈ ℓ one has ˜ G ∈ L (0 , π ) and ˜ S ∈ L (0 , π ), hence the following functions˜∆( λ ) := ˜∆( ρ ) = − ρ sin ρπ + ˜ ω cos ρπ + Z π ˜ G ( t ) cos ρt dt, λ ∈ C , (3.9)˜∆ ( λ ) := ˜∆ ( ρ ) = cos ρπ + ˜ ω sin ρπρ + 1 ρ Z π ˜ S ( t ) sin ρt dt, λ ∈ C \ { } (3.10)10re well defined (to specify ρ we use the square root branch with Im √ λ ≥ canbe continuously extended to λ = 0 as well. The resulting functions ˜∆ and ˜∆ are entire functions oforder 1 /
2. Moreover, representations (2.34) and (2.35) hold for the functions ˜∆ and ˜∆ if one changes ω , ω and h n by ˜ ω , ˜ ω and ˜ h n respectively. Let us define˜ M ( λ ) = − ˜∆ ( λ )˜∆( λ ) , λ ∈ C . Then one can verify that˜ M ( z k ) = M k = M ( z k ) , k ≥ . We would like to emphasize here that we do not know if the function ˜ M is the Weyl function corre-sponding to some potential q , so Theorem 2.5 can not be applied directly. Nevertheless, the proof ofthis theorem in [9] requires only that ˜ M is a quotient of two entire functions satisfying some basicasymptotic conditions, which can be easily verified for the functions ˜∆ and ˜∆ . Hence following theproof of Theorem 2.1 from [9] we obtain that ˜ M ≡ M , a contradiction.Now suppose that only a finite set of pairs { z k , M k } Kk =0 := { z k , M ( z k ) } Kk =0 is given. We assumethat the numbers z k are real (the system of equations (3.5) can be formulated for complex numbers z k as well) and ordered: z k < z k +1 , 0 ≤ k ≤ K −
1. Clearly any finite sequence of numbers z k , k ≤ K ,can be extended to an infinite one satisfying the condition (2.15), which is not of a great use sincethe truncated system (3.5) for some particular choices of z k , k ≤ K may not be solved uniquely. Forexample, taking z k = k for all k ≤ K makes it impossible to find the coefficient ω . As a rule ofthumb we can ask that the terms in (2.15) remain small and bounded, even better if they decay as k increases. So we can formulate the following empirical requirement for the numbers z k = ˜ ρ k :˜ ρ k < k z k < k k + 1 , at least starting from some k = K ′ < K . Such requirement is sufficient to deal, for example, with thetwo spectra inverse problem, for which { z k } Kk =0 = { λ k } Kk =0 ∪ { ν k } Kk =0 .Let us consider a truncated system (3.5), [( K − / X n =0 ξ n · ( − n √ n + 2 M k j n (˜ ρ k π ) √ π + [( K − / X n =0 ξ n +1 · ( − n √ n + 2 j n +1 (˜ ρ k π ) √ π ˜ ρ k + ω · M k cos ˜ ρ k π + ω · sin ˜ ρ k π ˜ ρ k = M k ˜ ρ k sin ˜ ρ k π − cos ˜ ρ k π, ≤ k ≤ K. (3.11)Solving this system one obtains approximate values of the parameters ω , ω and the coefficients h , . . . , h K − .The coefficient matrix of the truncated system (3.11) can be badly conditioned. Since we knowthat the coefficients h n decay, see Proposition 2.9, there is no need to look for the same number of theunknowns as the number of points ˜ λ k . One may consider less unknowns and treat the system (3.11)as an overdetermined one. Asking for the condition number of the resulting coefficient matrix to berelatively small one estimates an optimum number M , M ≤ K −
2, of the coefficients h , . . . , h M . Werefer the reader to [24, Subsection 3.2 and Section 5] for additional details.11 .2 Recovery of eigenvalues λ n and norming constants α n Suppose we have the coefficients ω , ω and { h n } ∞ n =0 recovered as described in Subsection 3.1. Thenwe can use the representation (2.34) to calculate the characteristic function ∆( ρ ) for any value of ρ .Recalling that the eigenvalues of the problem (2.1), (2.2) are real and coincide with zeros of ∆( λ ), therecovery of the eigenvalues reduces to finding zeros of the entire function ∆( λ ).Having the spectrum λ k = ρ k , the corresponding norming constants can be found using (2.14),(2.35) and (2.36). Indeed, α k = − ( λ k ) · ddλ ∆( λ ) (cid:12)(cid:12)(cid:12)(cid:12) λ = λ k = − ρ k ∆ ( ρ k ) ddρ ∆( ρ ) (cid:12)(cid:12)(cid:12)(cid:12) ρ = ρ k . Suppose we have only a finite number of coefficients h , . . . , h M . Then approximate eigenvaluesare sought as zeros of the function∆ M ( ρ ) = − ρ sin ρπ + ω cos ρπ + 2 [ M/ X n =0 ( − n h n j n ( ρπ ) , (3.12)and afterwards for each approximate eigenvalue λ k = ρ k the corresponding norming constant can beobtained from α k ≈ (1 + πω ) sin ρ k π − πρ k cos ρ k π + 2 P [ M/ n =0 ( − n h n (cid:16) nρ k j n ( ρ k π ) − πj n +1 ( ρ k π ) (cid:17) ρ k cos ρ k π + 2 ω sin ρ k π + 4 P [( M − / n =0 ( − n h n +1 j n +1 ( ρ k π ) . (3.13) Known the eigenvalues { λ n } ∞ n =0 and the corresponding norming constants { α n } ∞ n =0 , the potential q and the parameters h and H can be recovered by solving the Gelfand-Levitan equation (2.38). In[20] (see also [21]) using representation (2.22) the Gelfand-Levitan equation was reduced to an infinitesystem of linear algebraic equations for the coefficients β n ( x ). In [24] a more accurate modificationof the method was proposed. We present the final result from [24] below.For this subsection and without loss of generality, we suppose that ρ = 0. This always can beachieved by a simple shift of the potential. If originally ρ = 0, then we can consider the new potential e q ( x ) := q ( x ) − ρ . Obviously, the eigenvalues { λ n } ∞ n =0 shift by the same amount, while the numbers h , H and { α n } ∞ n =0 do not change. After recovering e q ( x ) from the shifted eigenvalues, one gets theoriginal potential q ( x ) by adding ρ back.Let us denote e c km ( x ) = − ωx π (cid:18) δ m,k − (2 k − / − δ m,k (2 k − / + δ m,k +1 (2 k + 1 / (cid:19) + ( − k + m ∞ X n =1 (cid:20) j k ( ρ n x ) j m ( ρ n x ) α n − j k ( nx ) j m ( nx ) π (3.14)+ 2 ωπ n (cid:18) xj k ( nx ) j m +1 ( nx ) + xj k +1 ( nx ) j m ( nx ) − k + m ) j k ( nx ) j m ( nx ) n (cid:19)(cid:21) , e C m ( x ) = (cid:18) α − π + 2 ωx π (cid:19) δ m + 2 ωx π δ m + e c m ( x ) , (3.15) e C m ( x ) = 2 ωx π δ m + e c m ( x ) , (3.16) e C km ( x ) = e c km ( x ) , k = 2 , , . . . , m ∈ N , (3.17)12nd e d k ( x ) = − (cid:18) α − π + 4 ωx π − ωxπ (cid:19) δ k − ωx π δ k − ( − k ∞ X n =1 (cid:20) cos ρ n xα n j k ( ρ n x ) − nxπ j k ( nx )+ 2 ωπ n (cid:18) x sin nxj k ( nx ) + x cos nxj k +1 ( nx ) − kn cos nxj k ( nx ) (cid:19)(cid:21) , (3.18)where δ k,m stands for the Kronecker delta and ( k ) m is the Pochhammer symbol.Then the following result is valid. Theorem 3.2 ([24]) . The coefficients β m ( x ) satisfy the system of linear algebraic equations β k ( x )(4 k + 1) x + ∞ X m =0 β m ( x ) e C km ( x ) = e d k ( x ) , k = 0 , , . . . . (3.19) For all x ∈ (0 , π ] and k = 0 , , . . . the series in (3.19) converges. It is of crucial importance the fact that for recovering the potential q as well as the constants h and H it is not necessary to compute many coefficients β m ( x ) (that would be equivalent to approximatethe solution of the Gelfand-Levitan equation) but instead the sole β ( x ) is sufficient for this purpose,see Corollary 2.7.For the numerical solution of the system (3.19) it is natural to truncate the infinite system, i.e., toconsider m, k ≤ N . As was shown in [24], the truncation process possesses some very nice properties.Namely, the truncated system is uniquely solvable for all sufficiently large N , the approximate solutionsconverge to the exact one, the condition numbers of the coefficient matrices are uniformly boundedwith respect to N , and the solution is stable with respect to small errors in the coefficients. Moreover,a reduced number of equations in a truncated system is sufficient for recovering the coefficient β witha high accuracy. However it should be noted that a lot of (approximate) eigenvalues and normingconstants are necessary to obtain values of the series (3.14)–(3.18) accurately enough. As was shownin Subsection 3.2, it is not a problem, thousands of approximate eigendata can be computed. Given a finite set of points { z n } Nn =0 = (cid:8) ˜ ρ n (cid:9) Nn =0 and the corresponding values M n = M ( z n ), 0 ≤ n ≤ N of the Weyl function at these points, the following direct method for recovering the potential q andthe numbers h and H is proposed.1. Solving (overdetermined) system (3.11) find the coefficients ω , ω and h , . . . , h M .2. Compute h = ω − ω .3. Find approximate eigenvalues λ k , k ≤ K as zeros of (3.12). Compute the corresponding normingconstants α k , k ≤ K using (3.13).4. If ρ = 0, perform the shift of the eigenvalues˜ ρ k = q ρ k − ρ , k ≥ ,
13o that 0 becomes the first eigenvalue of the spectral problem. The parameters ω and ω haveto be shifted as well,˜ ω = ω − π ρ , ˜ ω = ω − π ρ . Let us denote the square roots of the shifted eigenvalues by the same expression ρ k , and similarlyfor the shifted parameters ω and ω .5. For a set of points { x l } from (0 , π ], compute the approximate values of the coefficients e C km ( x )and e d k ( x ) for k, m = 0 , . . . , N with the aid of the formulas (3.14)–(3.18) and solve the truncatedsystem (3.19) obtaining thus β ( x ).6. Compute q from (2.25). Take into account that ϕ (0 , x ) is an eigenfunction associated with thefirst eigenvalue λ and hence does not have zeros on [0 , π ] (see, e. g., [2, Theorem 8.4.5]). Thisjustifies the division over ϕ (0 , x ).7. Compute H using (2.33). For this compute the mean of the potential R π q ( t ) dt , and thus, H = ω − Z π q ( t ) dt.
8. Recall that one has to add the original eigenvalue ρ back to the recovered potential to returnto the original potential q . Remark 3.3.
The idea of interval flipping proposed in [24, Subsection 3.6] to improve the recoveryof the potential near x = π can be adapted for the proposed method as well, the formulas (3.39) and(3.40) from [24] are directly applicable with the data obtained on Steps 1–4. We leave the details tothe reader. Remark 3.4.
The proposed method with few modifications can be adapted to the case of Dirichletboundary condition at x = π , corresponding to the value H = ∞ . The Weyl function in such case isgiven by M ∞ ( λ ) = − S ( ρ, π ) ϕ ( ρ, π ) . Following the described steps one can similarly obtain an infinite system of equations for the coefficients β n ( π ) . The coefficients ω and ω are no longer necessary. Squares of the zeros of the function ϕ ( ρ, π ) are the eigenvalues, and the corresponding norming constants can be obtained similarly. Now onehas to solve an inverse problem by the given spectral density function in the case when one boundarycondition is of the Dirichlet type. It can be done with the use of a Gelfand-Levitan equation with adifferent kernel function F , see [25, Chapter 2, § Below, in Section 4 we illustrate the performance of this algorithm with several numerical examples.
The proposed algorithm can be implemented directly, similarly to the algorithm from [24], only severalobservations should be made. On Step 1, to determine the number of unknowns to look for, we usedthe following criteria: we try to recover at least 10 (if less than 30 points z n were given) or at least140 coefficients h n (for more than 30 points z n given); if the condition number of the coefficient matrixof the system (3.11) is less than 1000, we increase the number of recovered coefficients h n until thecondition number surpasses 1000. A significant time saving can be achieved by computing the sphericalBessel functions j m (˜ ρ k π ) using a backwards recursion, see [24], [3], [12] for details.The number K of the eigenvalues to compute on Step 3 can be estimated from the decay of thecoefficients h n . Faster decay means smoother potential, larger value of K . Say, K = 10 for smoothpotentials (coefficients h n decay very fast) and up to K = 10 for non-smooth potentials (coefficients h n decay slowly). For several examples originating from the inverse problems considered in Subsection2.2, one can take an arbitrary value of the parameter h when transforming the Weyl-Titchmarsh m -function into the function M using (2.11). In all these examples we took h = 0. Since the algorithmrecovers two parameters ω and ω which are equal due to the choice h = 0, we used the difference | ω − ω | as some additional indicator of the accuracy of the recovered coefficients. Smaller value | ω − ω | suggests to take more approximate eigenvalues (larger value of K used, as was aforementioned).Similarly to [24], the number of points x l taken on Step 4 should not be large, about one hundreduniformly spaced points works best. We used 8 equations in the truncated system (3.11) in all theexamples. All the computations were performed in Matlab 2017 on an Intel i7-7600U equipped laptopcomputer. We used spline approximation and differentiated the obtained spline on Step 6. In allexamples where eigenvalues or solution of a Cauchy problem were necessary, we applied the methodfrom [22].First we consider several inverse spectral problems which can be reduced to Problem 2.4. In thelast subsection we illustrate that the proposed algorithm can be applied even for complex points z n ,but may fail for some particular choices of the points. Consider an inverse problem of recovering the potential and boundary parameters by two spectra,see Subsubsection 2.2.1. The same problem was considered in [24] with the only difference that theshared boundary condition was at the point 0 (which is not essential since one can apply the changeof variable x ↔ π − x ). Since for this problem the given values of the Weyl function are either 0 or ∞ ,the system (3.1) splits into two independent systems. One (given by (3.2)) coincides with the system(3.8) in [24], but the second system is different from (4.6) from [24]. Moreover, in [24] we recoveredthe parameters ω and ω and obtained larger index eigenvalues from the eigenvalue asymptotics, andused the solutions of the systems only to compute the norming constants. While in this example weapplied the proposed algorithm as is, without separating the system (3.1) or previously obtaining thevalues of ω and ω , and computing all the eigenvalues from the obtained approximate characteristicfunction. Of course, one can not expect obtaining the same accuracy as those of [24].We considered three potentials: smooth q ( x ) = π x exp (cid:0) − xπ (cid:1) (potential from [7], adapted tothe interval [0 , π ]), non-smooth continuous q ( x ) = (cid:12)(cid:12) −| x − | (cid:12)(cid:12) (potential from [18]) and discontinuous q ( x ) = , x ∈ [0 , π ] ∪ [ π , π ) , − xπ + , x ∈ ( π , π ] , xπ − , x ∈ ( π , π ) , , x ∈ [ π , π ) , , x ∈ [ π , π ] . (potential matching the one from [7]). For all potentials we took h = 1 and H = 2.On Figure 1 we present the recovered potential q and its absolute error. 6, 10 or 16 eigenvaluepairs were used for recovery. As one can see, the error almost stabilizes. For 16 eigenvalue pairs the15arameter q was obtained with L (0 , π ) error of 8 . · − , the constants h and H were obtained withthe errors of 4 . · − and 2 . · − .On Figure 2 we present the potential q recovered from 40 and 201 eigenvalue pairs. For 40eigenvalue pairs the potential q was obtained with L (0 , π ) error of 4 . · − , the constants h and H were obtained with the errors of 0 .
017 and 0 . . · − , 0 .
011 and0 . q recovered from 30 and 201 eigenvalue pairs. While for 30eigenvalue pairs the algorithm failed near the interval endpoints, inside the interval (0 . , π − .
5) theresult is acceptable. For 201 eigenvalue pairs the algorithm was able to recover both potential andboundary parameters h and H , L error of the potential was 0 .
24, absolute errors of the boundaryparameters were 0 . . As it is stated in [9], “Two-thirds of the spectra of three spectral problems uniquely determine q .” Inthis subsection we numerically illustrate the solution of such inverse problems.To illustrate the performance of the algorithm we considered two potentials, the same smooth po-tential q and the following potential, possessing a discontinuous second derivative q ( x ) = R x q ( s ) ds ,where q and q were introduced in Subsection 4.1. We took H = 2 and h ∈ { , , } , and eigenvaluesindexed 0 , , , , , , . . . from the first spectrum, 0 , , , , , , , . . . from the second spectrum and1 , , , , , , . . . from the third spectrum.On Figure 4 we present the recovered potential q and its absolute error. The absolute error rapidlydecreases with more eigenvalues used and stabilizes at 24 eigenvalues (8 from each spectra). L errorof the recovered q was 3 . · − , the parameter H had an absolute error of 3 . · − .On Figure 5 we present the recovered potential q and its absolute error. In this case the potential isof finite smoothness, so the algorithm converges slowly. With 120 eigenvalues (40 from each spectrum)we obtained the following errors: L error of the potential was 3 . · − and the error of the parameter H was 1 . · − .In [7, Section 3.14] the following inverse problem was considered. Suppose that the first eigenvalueis known for spectral problems having the same potential but different boundary conditions at x = 0.From countably many values the potential can be recovered uniquely. However, if only finitely manyvalues are given, the situation changes. The corresponding numerical problem is extremely ill-posed,so the method from [30] was not able to solve it, and the modification from [7] was able to recoveronly the general shape of the potential.We tried the proposed algorithm on this problem. For a numerical example we took the potential q from Subsection 4.1, H = 0 and h ∈ { , , , , , , , , , } . On Figure 6, left plot, wepresent the recovered potential. The algorithm failed, it is a case of very small interval for z k . Andthe result does not improve if we take more values of the parameter h . However the situation greatlyimproves if instead of 1 eigenvalue for each spectral problem we take several. On Figure 6, right plot,we present the recovered potential for the inverse problem when the first 3 eigenvalues were takenfrom each spectrum. In this subsection we illustrate the performance of the proposed method applied to solution of inverseproblems with a partially known potential. For comparison we have taken two potentials from [18],a smooth potential q ( x ) = e x − x and a non-smooth continuous potential q ( x ) = | − | x − || . Thepotentials are known on [0 , π/ -10 -8 -6 -4 -2 Recovered q, 6 EV pairsRecovered q, 10 EV pairsRecovered q, 16 EV pairs
Figure 1: On the left: exact (blue line), recovered from 6 eigenvalue pairs (red dots) and recoveredfrom 10 eigenvalue pairs (black dots) potential q from Subsection 4.1. On the right: absolute errorof the recovered potential, 6, 10 or 16 eigenvalue pairs were used. Figure 2: On the left: exact (blue line) and recovered from 40 eigenvalue pairs (black dots); onthe right: exact (blue line) and recovered from 201 eigenvalue pairs (black dots) potential q fromSubsection 4.1. Figure 3: On the left: exact (blue line) and recovered from 30 eigenvalue pairs (black dots); onthe right: exact (blue line) and recovered from 201 eigenvalue pairs (black dots) potential q fromSubsection 4.1. 17 -10 -5 Recovered q, 12 EV usedRecovered q, 18 EV usedRecovered q, 30 EV used
Figure 4: On the left: exact (blue line), recovered from 12 eigenvalues (red dots) and recovered from18 eigenvalues (black dots) potential q from Subsection 4.2. On the right: absolute error of therecovered potential, 12, 18 or 24 eigenvalues were used. -4 -2 Recovered q, 60 EV usedRecovered q, 120 EV used
Figure 5: On the left: exact (blue line), recovered from 60 eigenvalues (red dots) and recovered from120 eigenvalues (black dots) potential q from Subsection 4.2. On the right: absolute error of therecovered potential, 60 or 120 eigenvalues were used. Figure 6: Exact potential q from Subsection 4.2 (blue line) and recovered one (black dots). On theleft: recovered from 10 eigenvalues λ ( h ), h ∈ { , , , , , , , , , } , on the right: recoveredfrom 30 eigenvalues λ ( h ) , λ ( h ) , λ ( h ), h ∈ { , , , , , , , , , } .1818] the Dirichlet boundary conditions are considered, we decided to take mixed boundary conditionswith h = 1 and H = 2 for both potentials.To solve the inverse problem we followed the idea presented in Subsubsection 2.2.3: solve theCauchy problems with initial data v ( ρ n ,
0) = 1, v ′ ( ρ n ,
0) = h on the segment [0 , π/ m -function m ( π/ , ρ n ) = v ′ ( ρ n , π/ v ( ρ n , π/ , transform them to the Weyl function M using (2.11). After a simple rescaling we get Problem 2.4.For the numerical example 40 eigenvalues with 1 missing eigenvalue were considered in [18]. For thepotential q we similarly considered 40 eigenvalues, but with 2 eigenvalues missing: λ and λ . Forthe potential q one eigenvalue λ was missing. We would like to emphasize that we used the proposedmethod “as is”, without recovering missing eigenvalues (which is impossible, since the algorithm hasno information if the points z n correspond to any eigenvalues).On Figures 7 and 8 we present the recovered potentials and their absolute errors. One can appre-ciate an especially excellent accuracy in the case of a smooth potential. z n and some problematic choices The proposed algorithm works even when the numbers z n are complex (however for numerical stabilitythe imaginary parts should be relatively small comparing to the real parts). For illustration weconsidered the potential q from Subsection 4.1, h = 1, H = 2 and the following sets of points z ( k ) n = (cid:18)
15 + n k · i (cid:19) , ≤ n ≤ , k ∈ { , , } . (4.1)The obtained results are presented on Figure 9, left plot. The algorithm was able to recover boththe potential and the boundary conditions, however the accuracy decreases when the points are takenfarther from the real line.It should be mentioned that the proposed algorithm may fail for some particular choices of thepoints z n even for nice potentials. We considered the same problem and the following set of points y (1) n = (cid:18)
14 + n (cid:19) , ≤ n ≤ . The truncated system (3.11) was not ill conditioned, but its solution was nowhere close to the exactone and the situation does not improve when taking more or less unknowns. For example, for 9unknowns ξ , . . . , ξ the parameter h = ω − ω resulted to be 0 . . y (2) n = (cid:16) .
249 + n (cid:17) , ≤ n ≤ . For the latter choice of points, the recovered parameter h resulted to be 0 . { y (1) n } ∞ n =0 and { y (2) n } ∞ n =0 do not satisfy the condition (2.15), so any behavior can be expected for finite subsets. And the set { y (1) n } ∞ n =0 is, in some sense, the worst possible since the points y (1) n are the most distant from anyeigenvalue asymptotics (numbers close either to integers or integers plus one half). The situationchanges completely if we consider the set { } ∪ { y (1) n } ∞ n =0 . This set satisfies the condition (2.15),and the proposed method works for it without any problem, see Figure 9, right plot. The recoveredparameter h was 0 . -12 -10 -8 -6 Figure 7: On the left: exact (blue line) and recovered (black dots) potential q from Subsection 4.3.On the right: absolute error of the recovered potential. Eigenvalues { λ n } n =0 \ { λ , λ } were used.The error of the recovered parameter H is 5 . · − , and the L error of the recovered potential is3 . · − . -4 -3 -2 -1 Figure 8: On the left: exact (blue line) and recovered (black dots) potential q from Subsection 4.3.On the right: absolute error of the recovered potential. Eigenvalues { λ n } n =0 \ { λ } were used. Theerror of the recovered parameter H is 2 . · − , and the L error of the recovered potential is 0 . -12 -10 -8 -6 -4 -2 k=3k=2k=1k=0 -10 -5 y (1)n =(0.25+n/2) y (2)n =(0.249+n/2) {y (1)n } {0} Figure 9: Absolute errors of the recovered potential q from Subsection 4.4. On the left: using thepoints (4.1) with different values of k . On the right: using the points shown on the legend, 0 ≤ n ≤ cknowledgements Research was supported by CONACYT, Mexico via the project 284470. Research of VladislavKravchenko was supported by the Regional mathematical center of the Southern Federal University,Russia.
References [1] M. Abramovitz and I. A. Stegun,
Handbook of mathematical functions , New York: Dover, 1972.[2] F. V. Atkinson,
Discrete and continuous boundary problems , New York: Academic Press, 1964.[3] A. R. Barnett,
The calculation of spherical Bessel and Coulomb functions , in ComputationalAtomic Physics, Electron and Positron Collisions with Atoms and Ions, Berlin: Springer-Verlag,1996, 181–202.[4] N. P. Bondarenko,
Inverse Sturm-Liouville problem with analytical functions in the boundarycondition , Open Mathematics 18 (2020), 512–528.[5] N. P. Bondarenko,
Solvability and stability of the inverse Sturm–Liouville problem with analyticalfunctions in the boundary condition , Math. Meth. Appl. Sci. 43 (2020), 7009–7021.[6] H. Campos, V. V. Kravchenko and S. M. Torba,
Transmutations, L-bases and complete familiesof solutions of the stationary Schr¨odinger equation in the plane , J. Math. Anal. Appl. 389 (2012),no. 2, 1222–1238.[7] Kh. Chadan, D. Colton, L. P¨aiv¨arinta, W. Rundell,
An introduction to inverse scattering andinverse spectral problems . SIAM, Philadelphia, 1997.[8] K. Chadan, P. C. Sabatier,
Inverse problems in quantum scattering theory. Second edition ,Springer-Verlag, New York, 1989.[9] F. Gesztesy, R. del Rio, and B. Simon,
Inverse spectral analysis with partial information onthe potential, III. Updating boundary conditions , Internat. Math. Research Notices 15 (1997),751–758.[10] F. Gesztesy and B. Simon,
A new approach to inverse spectral theory II. General real potentialsand the connection to the spectral measure , Ann. of Math. 152 (2000), 593–643.[11] F. Gesztesy and B. Simon,
Inverse spectral analysis with partial information on the potential, II.The case of discrete spectrum , Trans. Amer. Math. Soc. 352 (2000), 2765–2787.[12] E. Gillman and H. R. Fiebig,
Accurate recursive generation of spherical Bessel and Neumannfunctions for a large range of indices , Comput. Phys. 2 (1988), 62–72.[13] T. N. Harutyunyan and A. M. Tonoyan,
Some properties of the kernel of a transmutation operator ,2021, to appear.[14] H. Hochstadt,
The inverse Sturm–Liouville problem , Comm. Pure Appl. Math., 26 (1973), 715–729.[15] H. Hochstadt and B. Lieberman,
An inverse Sturm-Liouville problem with mixed given data ,SIAM J. Appl. Math. 34 (1978), 676–680.[16] M. Horv´ath,
On the inverse spectral theory of Schr¨odinger and Dirac operators , Trans. Amer.Math. Soc. 353 (2001), 4155–4171.[17] M. Horv´ath,
Inverse spectral problems and closed exponential systems , Ann. of Math. 162 (2005),885–918. 2118] A. Kammanee and C. B¨ockmann,
Determination of partially known Sturm–Liouville potentials ,Appl. Math. Comput. 204 (2008) 928–937.[19] K. V. Khmelnytskaya, V. V. Kravchenko, S. M. Torba,
A representation of the transmutationkernels for the Schr¨odinger operator in terms of eigenfunctions and applications , Appl. Math.Comput. 353 (2019), 274–281.[20] V. V. Kravchenko,
On a method for solving the inverse Sturm–Liouville problem , J. Inverse Ill-posed Probl. 27 (2019), 401–407.[21] V. V. Kravchenko,
Direct and inverse Sturm-Liouville problems: A method of solution ,Birkh¨auser, Cham, 2020.[22] V. V. Kravchenko, L. J. Navarro, S. M. Torba,
Representation of solutions to the one-dimensionalSchr¨odinger equation in terms of Neumann series of Bessel functions , Appl. Math. Comput. 314(2017), 173–192.[23] V. V. Kravchenko and S. M. Torba,
A Neumann series of Bessel functions representation forsolutions of Sturm-Liouville equations , Calcolo 55 (2018), article 11, 23pp.[24] V. V. Kravchenko, S. M. Torba,
A direct method for solving inverse Sturm-Liouville problems ,Inverse Probl. 37 (2021), 015015 (32pp).[25] B. M. Levitan,
Inverse Sturm-Liouville problems , VSP, Zeist, 1987.[26] V. A. Marchenko,
Some questions on one-dimensional linear second order differential operators ,Transactions of Moscow Math. Soc., 1 (1952), 327–420.[27] V. A. Marchenko,
Sturm-Liouville operators and applications: revised edition , AMS Chelsea Pub-lishing, 2011.[28] M. Neher,
Enclosing solutions of an inverse Sturm-Liouville problem with finite data , Computing53 (1994), 379–395.[29] V. Pivovarchik,
An inverse problem by eigenvalues of four spectra , J. Math. Anal. Appl. 396(2012) 715–723.[30] W. Rundell, P. E. Sacks,
Reconstruction techniques for classical inverse Sturm–Liouville problems ,Math. Comput. 58 (1992), 161–183.[31] A. M. Savchuk,
Reconstruction of the potential of the Sturm–Liouville operator from a finite setof eigenvalues and normalizing constants , Math. Notes 99 (2016), 715–728.[32] E. L. Shishkina, S. M. Sitnik,
Transmutations, singular and fractional differential equations withapplications to mathematical physics , Elsevier, Amsterdam, 2020.[33] A. M. Savchuk and A. A. Shkalikov,
Recovering of a potential of the Sturm-Liouville problemfrom finite sets of spectral data , in Spectral Theory and Differential Equations, Amer. Math. Soc.Transl. Ser. 2, 233 (2014), 211–224.[34] A. A. Shkalikov,
Boundary problems for ordinary differential equations with parameter in theboundary conditions , J. Sov. Math. 33 (1986), 1311–1342.[35] B. Simon,
A new approach to inverse spectral theory I. Fundamental formalism , Ann. of Math.150 (1999), 1029–1057.[36] C. F. Yang, N. P. Bondarenko, X. C. Xu,
An inverse problem for the Sturm-Liouville pencil witharbitrary entire functions in the boundary condition , Inverse Problems and Imaging, 14 (2020),153–169.[37] V. A. Yurko,