# A numerical algorithm for solving inverse problems for singular Sturm-Liouville operators.

Abstract

We propose a numerical algorithm for solving inverse problems of spectral analysis for Sturm-Liouville differential operators on the half-line. Moreover some results of numerical experiments are also presented.

AMS subject classification: 65L09, 34A55, 47E05. Keywords: Differential equations, inverse spectral problems, numerical methods.

1. Introduction

In this paper we provide a numerical algorithm for solving inverse problems of spectral analysis for singular Sturm-Liouville differential operators. Inverse spectral problems consist in recovering operators from their spectral characteristics. Such problems often appear in mathematics and various branches of natural sciences and engineering. Inverse problems also play an important role in solving nonlinear evolution equations in mathematical physics. Up to now several theoretical procedures have been created for solving inverse spectral problems (see the monographs [3, 10, 16, 18, 21, 29, 30] and the references therein). Furthermore, a number of works are devoted to numerical aspects for solving inverse spectral problems (see [2, 7, 9,17,20]. In this paper we suggest a new numerical algorithm for solving the inverse problem of recovering a singular Sturm- Liouville operator on the half-line from its spectral function. We also provide some results of numerical experiments resulting from the application of the algorithm.

2. Preliminaries and Theoretical Background

In this section we introduce the main notions and describe briefly necessary facts from the theory of inverse Sturm-Liouville problems on the half-line. For more details on this topic see [16, 18].

We consider a boundary value problem L = L(q(x), h) of the form

-y'' + q(x)y = [lambda]y, x [greater than or equal to] 0, (2.1)

y'(0) - hy(0) = 0, (2.2)

where [lambda] is the spectral parameter, q(x) is a real-valued continuous function, and h is a real number. The function q(x) is called the potential. Denote by [phi](x, [lambda]) the solution of equation (2.1) under the initial conditions [phi](0, [lambda]) = [phi]'(0, [lambda]) = h. For each fixed x [greater than or equal to] 0, the functions [[phi].sup.(v)](x, [lambda]), v = 0, 1, are entire in [lambda] of order 1/2.

Let [sigma]([lambda]), [lambda][member of] R, be the spectral function for L, normalized by the condition [sigma][greater than or equal to](-[infinity]) = 0. The function [sigma]([lambda]) is nondecreasing, continuous from the left, and

[sigma]([lambda]) = 2/[pi][rho] - h + 0(1), [lambda][right arrow] +[infinity],

where [lambda] = [[rho].sup.2]. Moreover, for each function f [member of] [L.sub.2](0,[infinity]) with a finite support, the following Parseval's equality holds:

[[integral].sup.[infinity].sub.0] [[absolute value of f(x).sup.2]] dx = [[integral].sup.[infinity].sub.[-infinity]] [absolute value of F([lambda]).sup.2]] d[sigma]([lambda]), (2.3)

where

F([lambda]) = [[integral].sup.[infinity].sub.0] f(x)[sigma](x, [lambda])dx.

The inverse problem under consideration is formulated as follows.

Inverse Problem 1. Given [sigma]([lambda]), construct q(x) and h.

There exist several methods for solving such inverse spectral problems. Our algorithm is based on the so-called transformation operator method which is due to Marchenko and Levitan [16, 18]. Let us describe briefly the main ideas of the method. For further details see [16, 18].

Take a boundary value problem [??] = L([??](x), [??]) of the same form as L but with different coefficients [??] and [??] (for example, one can take [??](x) [equivalent to] 0, [??] = 0). We agree that if a certain symbol [gamma] denotes an object related to L, then [??] will denote an analogous object related to [??], and [??] = [gamma] - [??]. In particular, if [??](x) [equivalent to] 0, [??] = 0, then

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (2.4)

Below, in our numerical procedure, we will change a "model" boundary value problem [??] in each step such that [??] has a finite support.

Theorem 2.1. The following relation holds:

[psi](x, [lambda]) = [??](x, [lambda]) + [[integral].sup.x.sub.0] A(x, t)[??](t, [lambda])dt, (2.5)

where A(x, t) (which does not depend on [lambda]) is a smooth function, and

q(x) = [??](x) + 2 d/dx A(x, x), h = [??] + A(0, 0). (2.6)

For each fixed x, the function A(x, t) satisfies the following linear integral equation:

A(x, y) + [PHI](x, y) + [[integral].sup.x.sub.0] A(x, t) [PHI](t, y)dt = 0, y < x, (2.7)

where

[PHI](x, y) = [[integral].sup.[infinity].sub.[-infinity]] [??](x, [lambda])[??](y, lambda]) d[??] ([lambda]). (2.8)

The proof of Theorem 2.1 is contained in [16, 18]. The operator

(Af)(x) = f(x) + [[integral].sup.x.sub.0] A(x, t)f(t)dt

is called the transformation operator, and the function A(x, t) is the kernel of the transformation operator.

The solution of Inverse Problem 1 can be found according to the following scheme:

1) Choose a "model" boundary value problem [??].

2) Calculate [PHI](x, y) via (2.8).

3) Find A(x, y) by solving the linear integral equation (2.7).

4) Construct q(x) and h by (2.6).

Now we consider the particular case of a perturbation of only one eigenvalue. Let a "model" boundary value problem [??] be given, and let [??]([lambda]) be the spectral function of [??]. Take numbers [[lambda].sub.0] [member of] R, [[alpha].sub.0] > 0, and denote

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Consider the function

[sigma]([lambda]) = [??]([lambda]) + S([lambda]).

The function [sigma]([lambda]) is the spectral function for a certain boundary value problem L of the form (2.1)-(2.2). Since [??]([lambda]) = S([lambda]), it follows from (2.8) that

[PHI](x, y) = [[alpha].sub.0][??](x, [[lambda].sub.0][??](y, [[lambda].sub.0]). (2.9)

Using (2.7) it is easy to check that here

A(x, y) = - [[alpha].sub.0][??](x, [[lambda].sub.0][??](y, [[lambda].sub.0]) / 1 + [[alpha].sub.0] [[integral].sup.x.sub.0] [[??].sup.2](t, [[lambda].sub.0]) dt. (2.10)

Together with (2.5) this yields

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (2.11)

According to (2.6), the potential q(x) and the coefficient h can be found by the formulae

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (2.12)

We will use the relations (2.9)-(2.12) as the basis of the numerical procedure described in the next section.

3. The Numerical Procedure

In this section we describe a numerical algorithm for solving Inverse Problem 1; moreover we provide results of numerical experiments.

Let the spectral function [sigma]([lambda]) of the boundary value problem L be given. For simplicity, we consider the case when the spectrum of L is bounded from below. Then, without loss of generality, we assume that [sigma]([lambda]) = 0 for [lambda] < 0. Choose a "model" boundary value problem [??] = L([??](x),[??]). Let [??] be the spectral function of [??]. For example, if one takes [??](x) = 0, [??] = 0, then (2.4) is valid.

Choose d > 0, and put [R.sub.k] = [d.sub.k], k [greater than or equal to] 0, and

[[alpha].sub.k] = [??] ([R.sup.2.sub.k]) - [??] ([R.sup.2.sub.k-1]), k [greater than or equal to] 1, (3.1)

where [??] = [sigma] - [??] as before. Denote

[[lambda].sub.k] := [([R.sub.k] - d/2).sup.2], k [greater than or equal to] 1, (3.2)

and consider the step-function

Z([lambda]) = [summation over ([[lambda].sub.k][less than or equal to][lambda] [[alpha].sub.k], [lambda] [greater than or equal to] 0.

We introduce the sequence of the functions [[sigma].sub.k]([lambda]), k [greater than or equal to] 0, by the formulae

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Then

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

For each fixed k [greater than or equal to] 0, the function [[sigma].sub.k]([lambda]) is the spectral function for a certain boundary value problem [L.sub.k] = L(qk(x), [h.sub.k]). In the following algorithm we will construct [L.sub.k] (i.e., [q.sub.k](x), [h.sub.k]) successively for k = 1, 2, 3, ..., using [L.sub.k-1] as a "model" boundary value problem. This allows us to use formulae (2.11)-(2.12) in each step of our calculations. For sufficiently large k (k := M), the pair ([q.sub.M](x), [h.sub.M]) is approximately equal to the desired pair (q(x), h). Let us describe the algorithm in more detail.

Algorithm 1. Let the spectral function [sigma]([lambda]) of the boundary value problem L = L(q(x), h) be given.

1) Choose a "model" boundary value problem [??] = L([??](x), [??]) with the spectral function [??]([lambda]).

2) Take d > 0, [R.sub.k] = dk, k [greater than or equal to] 0, and calculate [[alpha].sub.k], [[lambda].sub.k], k [greater than or equal to] 1, by (3.1)-(3.2).

3) Put [[sigma].sub.0]([lambda]) := [??]([lambda]), [q.sub.0](x) := [??](x), [h.sub.0] := [??], i.e., [L.sub.0] := [??].

4) Perform the following operations successively for k = 1, 2, 3, ...:

(i) Construct [psi]k-1(x) as the solution of the Cauchy problem

-[[psi]''.sub.k-1](x) + [q.sub.k-1](x)[[psi].sub.k-1](x) = [[lambda].sub.k][[psi].sub.k-1](x), [[psi].sub.k-1](O) = 1, [[psi]'.sub.k-1](O, [lambda]) = [h.sub.k-1], (3.3)

and put

[[psi].sub.k-1](x) := 1 + [[alpha].sub.k] [[integral].sup.x.sub.0] [[psi].sup.2.sub.k-1](t)dt.

(ii) Calculate [q.sub.k](x) and [h.sub.k] by the formulae

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (3.4)

We stop our calculations for some k = M, and consider the pair ([q.sub.M](x), [h.sub.M]) as an approximation for the desired pair (q(x), h).

In order to check if ([q.sub.M](x), [h.sub.M]) is sufficiently close to the desired pair, one can use Parseval's equality (2.3) for a sufficiently wide set of the test functions f.

Some remarks on the numerical realization of this algorithm are in order.

1) One can choose [R.sub.k] arbitrary in order to cover the corresponding interval in the [lambda]-axis.

2) In the algorithm we use only the values [sigma]([R.sup.2.sub.k]), k = [??], of the spectral function [sigma]([lambda]) at the points [lambda] = [R.sup.2.sub.k]. Thus, the algorithm can be used without changes in the case when the spectral function is given as a table.

3) In practice we calculate the potential numerically on some interval [0, T] at the points [x.sub.m] = m[delta], 0 [less than or equal to] m [less than or equal to] N with stepsize [delta] = T/N.

4) In each step of the algorithm we solve the Cauchy problem (3.3) numerically. This may be done by standard finite difference methods, with a number of operations of order O(N). Thus, the number of the operations of the whole algorithm is of order O(MN).

5) The derivative [[psi]'.sub.k-1] in (3.4) can be computed analytically with the help of (3.3), namely

[[psi]'.sub.k-1(x) = [h.sub.k-1] + [[integral].sup.x.sub.0] ([q.sub.k-1](t) - [[lambda].sub.k])[[psi].sub.k-1](t)dt.

6) Clearly, [[psi].sub.k-1](x) = [[psi].sub.k-1](x, [[lambda].k]), where [[psi].sub.k-1](x, [lambda]) is the solution of the Cauchy problem

-[[psi]''.sub.k-1](x, [lambda])+[q.sub.k-1](x)[[psi].sub.k-1](x,[lambda]) = [lambda][[psi].sub.k-1](x, [lambda]), [psi].sub.k-1](0, [lambda]) = 1, [[psi]'.sub.k-1](0, [lambda]) = [h.sub.k-1].

For calculating [[psi].sub.k](x, [lambda]), instead of solving the Cauchy problem (3.3), we can also use the following recurrence formula, which follows from (2.11):

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Numerical experiments show that Algorithm 1 works effectively for a wide class of boundary value problems. Below we present two examples for using this procedure.

Example 3.1. Let

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

One can check that [sigma]([lambda]) is the spectral function for a certain boundary value problem L = L(q(x), h) of the form (2.1)-(2.2). The results of the calculations by Algorithm 1 with [??] = 0, [??] = 0, d = 0.6, M = 5, T = p, [pi] = 31, are given in Table 1 and Fig. 1. Parseval's equality (2.3) for the constructed boundary value problem L was checked on the family of the test functions

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Example 3.2. Let

q(x) = 8 / [(1 + 2x).sup.2], h= 0. (3.5)

Then the spectral function can be calculated analytically (see ), namely:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (3.6)

Indeed, the Jost solution e(x, [rho]) for (3.5) has the form [10, p. 144] e(x, [rho]) = 1 - 2/i[rho](1 + 2x)exp(i[rho]x), x [greater than or equal to] 0.

Therefore the characteristic function [DELTA]([rho]) := e'(0, [rho]) is calculated by [DELTA]([rho]) = i[rho] - 2 + 4[(i[rho]).sup.-1].

Since d/d[lambda] [sigma]([lambda]) = V(lambda), [lambda] [greater than or equal to] 0, where

v([lambda)] = [rho] / [pi][[absolute value of [DELTA]([rho]).sup.2] = [[rho].sup.3] / [pi]([[rho].sup.4] - 4 [[rho].sup.2] + 16),

we arrive at (3.6).

We take [sigma]([lambda]) of the form (3.6) as the input data for Algorithm 1. Put [??] = 0, [??] = 0, d = 0.05, M = 200, T = 10, N = 200. The numerical realization of Algorithm 1 gives us a pair ([q.sub.M], [h.sub.M]) such that max [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

[FIGURE 1 OMITTED]

[FIGURE 2 OMITTED]

The approach, suggested in this section, allows one to construct effective numerical algorithms for solving inverse spectral problems for wide classes of differential operators. This approach can be based not only on the transformation operator method, but alternatively also on the method of spectral mappings , which gives us an opportunity to create similar numerical algorithms for arbitrary order differential equations, for systems of differential equations and for other classes of operators. Numerical experiments show that the suggested method is quicker than the methods from [2, 7, 9, 17,20] since it requires less operations. The method suggested in  for solving some inverse scattering problems on the line, requires a similar number of operations, but our approach has a wider area for applications.

In conclusion we note that other aspects of the numerical analysis for inverse spectral problems are reflected in [1, 4-6, 8, 11-15, 19, 22-26, 28, 31].

Acknowledgment

This research was supported in part by DAAD and by the Russian Foundation for Basic Research.

Received February 1, 2007; Accepted February 20, 2007

References

 Alan L. Andrew. Numerov's method for inverse Sturm-Liouville problems, Inverse Problems, 21(1):223-238, 2005.

 David C. Barnes. The inverse eigenvalue problem with finite data, SIAM J. Math. Anal., 22(3):732-753, 1991.

 Richard Beals, Percy Deift, and Carlos Tomei. Direct and inverse scattering on the line, volume 28 of Mathematical Surveys and Monographs, American Mathematical Society, Providence, RI, 1988.

 Liliana Borcea and Vladimir Druskin. Optimal finite difference grids for direct and inverse Sturm-Liouville problems, Inverse Problems, 18(4):979-1001, 2002.

 B. M. Brown,V. S. Samko, I.W. Knowles, and M. Marletta. Inverse spectral problem for the Sturm-Liouville equation, Inverse Problems, 19(1):235-252, 2003.

 B. M. Brown and R. Weikard. A Borg-Levinson theorem for trees, Proc. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci., 461(2062):3231-3243, 2005.

 Khosrow Chadan, David Colton, Lassi Paivarinta, and William Rundell. An introduction to inverse scattering and inverse spectral problems. SIAM Monographs on Mathematical Modeling and Computation. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 1997. With a foreword by Margaret Cheney.

 Y. Chen and V. Rokhlin. On the inverse scattering problem for the Helmholtz equation in one dimension, Inverse Problems, 8(3):365-391, 1992.

 Richard H. Fabiano, Roger Knobel, and Bruce D. Lowe. A finite-difference algorithm for an inverse Sturm-Liouville problem, IMA J. Numer. Anal., 15(1):75-88, 1995.

 G. Freiling and V. Yurko. Inverse Sturm-Liouville problems and their applications, Nova Science Publishers Inc., Huntington, NY, 2001.

 G. M. L. Gladwell. The application of Schur's algorithm to an inverse eigenvalue problem, Inverse Problems, 7(4):557-565, 1991.

 A. A. Gonchar, N. N. Novikova, and G. M. Khenkin. Multipoint Pade approximants in an inverse Sturm-Liouville problem, Mat. Sb., 182(8):1118-1128, 1991.

 Ole H. Hald. The inverse Sturm-Liouville problem and the Rayleigh-Ritz method, Math. Comp., 32(143):687-705, 1978.

 Bui Doan Khanh. A numerical resolution of the Gelfand-Levitan equation, J. Comput. Appl. Math., 72(2):235-244, 1996.

 Roger Knobel and Bruce D. Lowe. An inverse Sturm-Liouville problem for an impedance, Z. Angew. Math. Phys., 44(3):433-450, 1993.

 B. M. Levitan. Inverse Sturm-Liouville problems, VSP, Zeist, 1987. Translated from the Russian by O. Efimov.

 Bruce D. Lowe, Michael Pilant, and William Rundell. The recovery of potentials from finite spectral data, SIAM J. Math. Anal., 23(2):482-504, 1992.

 Vladimir A. Marchenko. Sturm-Liouville operators and applications, volume 22 of Operator Theory: Advances and Applications, Birkhauser Verlag, Basel, 1986. Translated from the Russian by A. Iacob.

 N. N. Novikova. An application of the Pade approximation to the inverse problem of monochromatic vibrosounding, Inverse Problems, 11(1):217-229, 1995.

 J.W. Paine, F. R. de Hoog, and R. S.Anderssen.Onthe correction of finite difference eigenvalue approximations for Sturm-Liouville problems, Computing, 26(2):123-139, 1981.

 Jurgen Poschel and Eugene Trubowitz. Inverse spectral theory, volume 130 of Pure and Applied Mathematics, Academic Press Inc., Boston, MA, 1987.

 Norbert Rohrl.Aleast-squares functional for solving inverse Sturm-Liouville problems, Inverse Problems, 21(6):2009-2017, 2005.

 William Rundell and Paul Sacks. On the determination of potentials without bound state data, J. Comput. Appl. Math., 55(3):325-347, 1994.

 William Rundell and Paul E. Sacks. The reconstruction of Sturm-Liouville operators, Inverse Problems, 8(3):457-482, 1992.

 William Rundell and Paul E. Sacks. Reconstruction techniques for classical inverse Sturm-Liouville problems, Math. Comp., 58(197):161-183, 1992.

 Paul E. Sacks. An iterative method for the inverse Dirichlet problem, Inverse Problems, 4(4):1055-1069, 1988.

 Paul E. Sacks. Reconstruction of steplike potentials, Wave Motion, 18(1):21-30, 1993.

 Catherine Willis. An inverse method for the inverse Dirichlet problem, Inverse Problems, 2:111-130, 1986.

 V. Yurko. Method of spectral mappings in the inverse problem theory, Inverse and Ill-posed Problems Series. VSP, Utrecht, 2002.

 V. A. Yurko. Inverse spectral problems for differential operators and their applications, volume 2 of Analytical Methods and Special Functions. Gordon and Breach Science Publishers, Amsterdam, 2000.

 E. P. Zhidkov and O. V. Kozlova. Continuous analogue of the Newton method in the inverse problem of scattering theory in the presence of eigenfunctions and eigenvalues, Mat. Model., 18(2):120-128, 2006.

G. Freiling Fachbereich Mathematik, Universitat Duisburg-Essen, D-47048 Duisburg, Germany E-mail: gerhard.freiling@uni-due.de

T. Mazur and V. Yurko Department of Mathematics, Saratov University, 410026 Saratov, Russia E-mail: yurkova@info.sgu.ru
```Table 1

x 0 0.1013 0.2027 0.3040 0.4054
q(x) 1.6883 1.9114 2.1345 2.1622 2.0149

x 0.9121 1.0134 1.1148 1.2161 1.3174
q(x) -0.1796 -0.5421 -0.7458 -0.7719 -0.6390

x 1.8242 1.9255 2.0268 2.1282 2.2295
q(x) 0.4171 0.3918 0.2717 0.0881 -0.1163

x 2.7362 2.8376 2.9389 3.0403 3.1416
q(x) -0.1908 0.0051 0.1973 0.3428 0.4083

x 0.5067 0.6081 0.7094 0.8107
q(x) 1.7173 1.3007 0.8072 0.2919

x 1.4188 1.5201 1.6215 1.7228
q(x) -0.3974 -0.1138 0.1466 0.3333

x 2.3309 2.4322 2.5335 2.6349
q(x) -0.2946 -0.4057 -0.4241 -0.3461

x
q(x)
```