Reconstruction of potential function in inverse Sturm-Liouville problem via partial data

Article History: Received 26 February 2021 Accepted 09 May 2021 Available 12 May 2021 In this paper, three different uniqueness data are investigated to reconstruct the potential function in the Sturm-Liouville boundary value problem in the normal form. Taking account of Röhrl’s objective function, the steepest descent method is used in the computation of potential functions. To decrease the volume of computation, we propose a theorem to precalculate the minimization parameter that is required in the optimization. Further, we propose a novel time-saving algorithm in which the obligation of using the asymptotics of eigenvalues and eigenfunctions and the appropriateness of selected boundary conditions are also eliminated. As partial data, we take two spectra, the set of the jth elements of the infinite numbers of spectra obtained by changing boundary conditions in the problem, and one spectrum with the set of terminal velocities. In order to show the efficiency of the proposed method, numerical results are given for three test potentials which are smooth, nonsmooth continuous, and noncontinuous, respectively.


Introduction
The inverse Sturm-Liouville (S-L) reconstruction problems consist of the calculation of the potential function from known data. These data which can be referred as spectrum, normalized constants, spectral functions, etc. are based on uniqueness. Although there are many studies on the uniqueness of the inverse S-L problem, it is observed from the literature that only a few studies have been conducted on the reconstruction of the potential function from data that can be obtained experimentally. In this study, we aim to contribute to the literature in this direction.
In 1978, Hald [1] considered the inverse problem of (1) for symmetric potential q with Dirichlet conditions, and by using Rayleigh-Ritz method he reduced the problem consist of Fourier expansion with finite terms to an eigenvalue problem for a matrix. Thus, he showed that the solution for matrix problem converged to the solution for the inverse problem as the dimension of the matrix increased. In 1984, Paine [2] used an algorithm to deal with a similar problem assuming q ∈ C 2 [0, π] in which a perisymmetric tridiagonal matrix obtained from first N eigenvalues. The errors arising from results obtained by using a correction term were analyzed according to the increasing the number N . In 1988, by the use of characteristic values {λ n , ρ n }, Sacks [3] generated the iteration algorithm *Corresponding Author q n+1 (x) = q n (x) + 2g(2x) − 2F (q n )(2x), for inverse S-L problem subject to Dirichlet boundary condition where and w(x, t) was considered to be the solution to related Goursat problem. Then some numerical examples were considered by noticing the convergence of q n to the potential q. In 1992, Lowe et al. [4] investigated the vector q = (q 1 , q 2 , ...) by proposing the potential in the form of for the S-L problem with Dirichlet boundary condition and with also general separable boundary conditions where {φ k (x)} k=2N k=1 = {sin 2πx, cos 2πx, ... , sin 2N πx, cos 2N πx}.
To determine the Fourier coefficients, they took the advantage of Newton's method, proved a convergence theorem of the method, and gave some examples. In the same year, Rundell and Sacks [5] handled two spectra to determine Cauchy data for a transformed hyperbolic partial differential problem, and then they used successive approximation method and Quasi-Newton method respectively. Finally, they considered the same reconstruction problem for different data and controlled effectiveness of their approach on smooth, non-smooth continuous, and noncontinuous test potentials. In 1994, Neher [6] studied the investigation of potential function q when the eigenvalues and symmetric base functions were given for the problem having symmetric potential with Dirichlet condition. His study was based on the determination of a j constants in which the potential function was written in the form of where a = (a j ) ∈ R n . He used Newton's method to find zeros of the function f (a) = (f i (a)) = (λ i (q(x; a)) − v i ) where the eigenvalues were v i ; i = 1, 2, ..., n. In fact, the formulae (2) and (3) used in [4] and [6], respectively, was considered in several earlier papers, notably the paper of O.H. Hald [1]. However, [4] deals with the difficulties arising from the limitation of data available in real applications, while [6] addresses the problem of attempting to enclose q within an interval-valued function. In 1995, Fabiano et al. [7] considered Dirichlet problem with symmetric and general potential, then converted the problem to a matrix equation through a partition of the interval. In this equation, the finite set of eigenvalues and terminal velocities were used to determine the matrix of coefficients. In 2004, Andrew [8] used the same method, which is Modified Newton's method, with Fabiano et al. [7]. The difference between the approaches used in [7] and [8] is that [7] uses a second order finite difference approximation of the differential equation, whereas [8] uses the more accurate Numerov method. The advantages of the latter approach are discussed in [9]. The author generalized the case in 2005 [9], and in 2011 he used a similar method to solve the problem corresponding to a different set of data [10]. In 2003, Brown et al [11] considered a finite number of linear dependence coefficients between some appropriate solutions to the Sturm-Liouville problem and eigenvalues. They used the steepest descent method for the objective functional where {λ n , C n } was the given set of finite data, and the functions u q , v q were the solutions corresponding to λ n for the Sturm-Liouville equation with some special initial conditions. In 2005, Röhrl [12] handled two spectra and used Polak-Ribiere conjugate gradient method to minimize his objective functional given with (4). In his ongoing study [13], he also generalized his work to boundary conditions. In 2007, Rafler and Böckmann [14] modified the Rundell-Sacks method to deal more effectively with potentials having jump discontinuities and gave some numerical examples in L 2 and L ∞ . Additionally, some other methods such as the boundary value method and the finite difference method to find the solutions of the inverse Sturm-Liouville problems can be seen in [15][16][17][18].
This paper is summarized as follows. Section 2.1 discusses the improvement of efficiency of the study given in [12]. In the method described in Sections 2 and 3 of [12], asymptotic formulas were used for the minimization of the objective functional in each iteration so that the original potential was approached as a Fourier series. Here, however, by using an estimate of minimization parameter which is calculated approximately for a random potential, the necessity of using asymptotics is eliminated. So, the original potential is approached as a linear combination of the eigenfunctions which is calculated with the help of an initial value problem. Two important advantages of this statement are expressed below: i) In the reconstruction problems, the errors in calculating the potential increase as the data at hand decreases. Besides, there are two factors that cause errors. The first of these is the error resulting from not being able to solve the direct problem analytically. It is not included in this study. The second is, however, the error that arises from the use of asymptotics in each iteration. As it can be observed from the asymptotics, though this error is small in the calculation of large eigenvalues, one cannot say the same for small eigenvalues. Thus, in the case where the number of experimentally obtained data is limited, one would like to minimize the error which occurs in small eigenvalues. As a result, the elimination of the necessity of the use of asymptotics minimizes the above-mentioned error. In this respect, better results are achieved by taking a small number of data pairs in our study. Moreover, even though increments on the iteration numbers is comprehended as a disadvantage for the case of two spectra, an important gain in terms of time is achieved because the calculation volume is decreased.
ii) The asymptotics of eigenfunctions corresponding to each element of different spectra can be the same (See [12]). Thus, in the calculation, a challenge occurs between the terms included in the gradient. This leads to significant increases in the number of iterations (Please see Figure 1-3 and the first paragraph of page 2013 in [12]). However, the elimination of asymptotics prevents such a challenge. Therefore, a gain in both the number of iterations and the time of calculation is achieved in our paper.
In Section 2.2 and Section 2.3, the process for the two spectra case in the previous section is also applied to the case of two different data sets. The first of these data is McLaughlin-Rundell data which was considered in [19]. When the studies about this problem is surveyed, while almost all of them are about generalizations of the uniqueness problem of [19], this study discusses the numerical solution of the problem. The second one is a spectrum and a set of terminal velocities. Finally, Section 3 is devoted to demonstrating the numerical results obtained for each uniqueness data in Section 2.

Two spectra
We propose to remove the use of asymptotic formulae of eigenvalues and corresponding eigenfunctions while reconstructing the potential. As a result of this suggestion, small eigenvalues and their corresponding eigenfunctions can be used more effectively. Moreover, different eigenvalues do not have to fight against each other when the corresponding asymptotic formulas for eigenfunctions are the same. Indeed, we especially handle the steepest descent method for Röhrl's objective functional.
It is well known that two spectra obtained by changing boundary condition determine the potential uniquely for (1) [20].
} be two spectra of the inverse S-L problem subject to different two boundary conditions for test potential Q. The Röhrl's objective functional is where ω i,j is positive weight constant. Since [12,21], are normalized eigenfunctions which correspond to the eigenvalues λ i,j,q [12]. It is seen from the literature that this type of objective function was firstly considered by Brown et al. [11]. Röhrl handled two spectra instead of Brown's data {λ n , C n } and generalized the problem to determine both the boundary conditions and the potential function [13]. Although it seems that Röhrl's functional is much simple than Brown's, it contains more computation than those in [11] because Röhrl needs to solve the boundary value problem instead of the initial value problem.
Here, it is aimed to minimize functional G(q) step by step so that one can approach potential Q by using iteration functions q n . If I is infinite and the positive weights (iω i,j ) are summable, then the series given by G(q) is convergent (see [12]). Furthermore G(q) = 0 if and only if q = Q from uniqueness [20]. The convergence is obvious if I is finite.
Thus a conjugate gradient algorithm will not get trapped in local minima [12].
By Theorem 1, the minimization process leads to G(q n ) → 0 so that q n → Q. The algorithm to be used is the steepest descent algorithm which is described as follows [22]: Step 0: Choose an initial potential as q 0 and set n = 0, Step 3: Minimize G(q n − h n ∇G(q n )) with respect to h n , Step 4: Set q n+1 := q n − h n ∇G(q n ) and replace n by n + 1, go to step 1.
There are three cases for the asymptotics of eigenvalues according to {α, β} in the eqn. (1) (See [12]). To compute asymptotics of squared normalized eigenfunctions corresponding to eigenvalues with respect to {α, β}, one can consider the study done by Hochstadt [23] which contains asymptotic solution to the normal form Sturm-Liouville differential equation with initial values y(0) = sin α, and α = π 4 , β = − π 4 , then the squared normalized eigenfunctions have the same asymptotics for all i ∈ M and j = 1, 2. These values of {α, β} were considered in [13]. When one chooses the values of {α, β} which have the same asymptotic form for corresponding the squared normalized eigenfunctions, the number of g 2 i,j in the gradient drops to half, and so the coefficients of the same asymptotics g 2 i,j fight each other for each index j. Thus, eliminating the necessity of the use of asymptotics is proposed to overcome such a conflict. Actually, because the asymptotics are not required until Step 3, any solver for eigenvalues and eigenfunctions can be used. Since the eqn. (1) cannot be solved in h n for Step 3, as long as a command that computes the parameter h n , is used, the asymptotics needs to be considered. Therefore, we obtain an estimate for h n by making some omissions: where (a i ) ∈ ℓ 2 and F i,α,β can be seen from asymptotics for eigenvalues (See [12]). Then, for the test potential Q and the iterative potential q n , it can be written as and so we have by neglecting the term A i,j,n = a q n+1 ,j,i − a qn,j,i . Thus by differentiating G(q n+1 ) with respect to the parameter h n and then equating it zero, we obtain Assume ω i,j = 1. Since It could be wondered that how would the calculation be if A i,j,n had not been neglected. For by taking its derivative with respect to the parameter h n , we obtain Since different expressions of g 2 i,j,qn ∇G(q n−1 ) − g 2 i,j,q n+1 ∇G(q n ) for all {i, n} cannot be equal to the same constant − Because we approach the original potential Q by using smooth iterative potentials q n , the small effect of this ratio does not change even if the original potential is not smooth. It should be noted that this neglect has no effect on the eigenvalues or the eigenfunctions. It just affects the number of iteration in the calculation.

McLaughlin-Rundell's data
In this part, we use the data handled by McLaughlin and Rundell in their article [19].
We consider the reconstruction of the potential q when we have a first few eigenvalues for a spectrum by making some modification in their uniqueness theorem. Firstly, we assume that the boundary conditions can be changed sufficiently, i.e., let us consider the problem −y ′′ (x) + (λ + q(x))y(x) = 0, 0 < x < 1, where q ∈ L 2 (0, 1) and β k for k = 1, 2, ... are distinct real numbers. Let λ j (q, β) be the jth eigenvalue of (5) for β instead of β k .
Proof. Since {β k } ∞ k=1 = n {β kn }, the set {λ kn,jn } is an infinite bounded set of real numbers. Therefore it has at least one finite accumulation point. The rest of this proof is similar as in the proof of Theorem 2.
Taking k = k n for each n allows us an extension for Theorem 3. So the eigenvalues can be used more than one for each spectrum. Because the theory is similar to those in the two spectra case (see [12]), we do not give it again here.

One spectrum and set of the terminal velocities
In this part, we use finite subsets of one spectrum {λ n } n≥1 and the set of the terminal velocities κ n = log g ′ n (1,q) g ′ n (0,q) as uniqueness data for Dirichlet problem. It is easy to see that the key is to have derivative of λ n and κ n with respect to potential q n in the steepest descent method. To see the details and to generalize the problem to all separable conditions, one can read the nice book written by Pöschel and Trubowitz [21].
For these data, the new objective functional and its gradient become and where I ⊂ N, respectively. Since y 1,n = y 1 (x, λ n ) and y 2,n = y 2 (x, λ n ) are the solutions to eqn. (6), In other words, Q(x) ∈ L 2 [0, 1] is the global minimum for G(q). Now, let us consider the bilinear form Γ : for the differentiable functions f, g : [0, 1] → R. This transformation is bounded by In particular Γ is continuous on H 1 [12]. Also, it is easy to see that Γ is antisymmetric because of the Wronskian. Some properties for the Wronskian are given as follows [12]: differentiable functions f, g, F, G. ii : For two arbitrary solutions f 1 and f 2 of the eqn. (6) with eigenvalue parameters λ 1 and λ 2 we have Since y 1 (x, λ) and y 2 (x, λ) satisfy the eqn. (6), we obtain and so we have [y 1 , y 2 ] = y 1 (0)y ′ 2 (0) − y ′ 1 (0)y 2 (0) = 1.
Proof. i) A more general proof was given by Röhrl [12].
Theorem 7. If I is finite or iω i is summable, the functional G(q) has no local minima at q with G(q) > 0. In other words, G(q) = 0 ⇔ ∇G(q) = 0.
The computation to find the parameter h n is similar to that in Section 2.1. Since and consequently where Because of g i,qn and therefore On the other hand, for all m, n ∈ N since 1 0 sin(2mπx)(1 − cos(2nπx))dx = 0, for 1 ≤ i ≤ k by using asymptotic formulas. Finally, by substituting (10) and (11) in (9), we obtain By the theory set out in this section, numerical calculations can be performed for three different data sets.

Numerical experiments
In this study numerical computations are obtained using Mathematics 11 with Eigen package [24]. The calculations are conducted on a desktop PC with a processor of Intel(R) Core i5-3470 CPU @3.2 GHz for Windows 7. We consider the initial potential q 0 (x) = 0 and the test potentials Q 1 (x), Q 2 (x), Q 3 (x) as follows: which were used in [11][12][13].

Numerical results for two spectra
To obtain two spectra data we choose α = 0, β 1 = 0 and β 2 = π 2 in the S-L problem (1). The steepest descent method is used to find the potential function for two spectra numerically via the methodology exhibited in Section 2.1. The numerical results that are calculated by our effective algorithm written in Mathematica are demonstrated in figures in which dashed curves show the numerical results while the solid lines show the exact values of potential function. We obtain these results in 381.59 and 3669.09 seconds with 41 and 110 iterations for the first two graphics in Fig. 1, respectively. We have also handled a few different test potential and different boundary conditions. We observed that although the number of iteration or the cost of computation of CPU may increase, there isn't any remarkable difference. However, a very small difference is occurred when we just consider two pairs data. In the numerical calculations of the potentials, the accuracy of G(q n ) ≈ 10 −6 is taken, and Fig. 1 and Fig. 2 are drawn according to this. We also take into account noise for five pairs data by adding random numbers in the interval (−0.01, 0.01) to each eigenvalue, it is demonstrated in Fig. 1c. For uniform noise, we have not seen a remarkable difference. The distribution of G(q n ) as n increases is shown by Fig. 1d.

Numerical results for
McLaughlin-Rundell's data To obtain data, we choose β k = tan π k 2 in the S-L problem (5). Similar to the way that is followed in Section 3.1, to approximate test potential functions with McLaughlin-Rundell's data, we use the methodology given in Section 2.2 with the same algorithm in Section 2.1. The upper bound of the value of G(q n ) is approximately 5 × 10 −5 . The numerical values of q for Q 1 (x) and Q 2 (x) are demonstrated in (a), (b), (c) and (d) of Fig. 3 with the dashed curves while the solid lines show the exact values of test potentials Q 1 and Q 2 . Fig. 3c is plotted to show the effectiveness of our extension on McLaughlin-Rundell's data which is plotted in Fig. 3b. It can be seen also that numerical results in Fig. 3c are very close to test potential and are better than of those in Fig. 3b. It is concluded that the numerical result we obtained is more accurate for smooth test potential Q 1 (x). However, it can be observed from Fig. 3d that for other test potentials we need more than two eigenvalues in calculations to have better approximation. Because we have a few first eigenvalues in real life problems, we can say that it will be more appropriate to prefer extra eigenvalues as auxiliary data with original data.

Numerical results for one spectrum and set of the terminal velocities
In this section, we consider a finite set consisted of the elements from one spectrum and terminal velocities of the S-L problem (6)-(7) to approximate the test potentials Q 1 (x) and Q 2 (x). The upper bound of the value of G(q n ) is approximately 10 −6 in calculations. By using the analogue Mathematica program performed in Section 3.1, we obtain the numerical results which are displayed in Fig. 4. When this data is compared with the first data considered in Section 3.1, it is seen from Fig. 4 that there is no remarkable difference for test potentials Q 1 and Q 2 . However when our computer program runs to obtain required numerical results for Q 3 , the calculation needs some time to get the same results in Section 3.1.

Conclusion
The inverse Sturm-Liouville reconstruction problem is an interesting subject that arises in many physical phenomena and it involves the reconstruction of coefficient functions and boundary conditions by aid of some data which determine these functions uniquely. In this study, we consider inverse S-L problem in the normal form and reconstruct the potential function. Actually, one of the main goal of our study is to eliminate the obligation of using the asymptotics for λ i,qn (α, β), g 2 i,j,qn , κ n (q n ), and ∂κn ∂qn in the steepest descent algorithm by estimating the parameter h n which appears in each iteration when the program runs. The other is to apply our approach to three    problem. Therefore, we think especially this part is worth considering.
Although the reconstruction of the potential by using alternative methods for the data in Section 2.3 appears in the literature, proposed approach in Section 2.1 is firstly considered here. The obtained numerical results are given in Section 3. It can be concluded that the approach for two spectra in Section 3.1 is more effective than other approaches for other data in Section 3.2 and Section 3.3.
The theory in this study can be generalized to other forms of the S-L problems, and it can be used not only for the reconstruction of potential function but also for the recovery of boundary conditions. Recently, there have been studies on fractional Sturm-Liouville problems [25,26]. We think that our approaches will also be useful for fractional Sturm-Liouville problems.