A Neumann series of Bessel functions representation for solutions of Sturm-Liouville equations

# A Neumann series of Bessel functions representation for solutions of Sturm-Liouville equations

Vladislav V. Kravchenko and Sergii M. Torba
Departamento de Matemáticas, CINVESTAV del IPN, Unidad Querétaro,
Libramiento Norponiente No. 2000, Fracc. Real de Juriquilla, Querétaro, Qro. C.P. 76230 MEXICO
e-mail: vkravchenko@math.cinvestav.edu.mx, storba@math.cinvestav.edu.mx
Research was supported by CONACYT, Mexico via the projects 166141 and 222478.
###### Abstract

A Neumann series of Bessel functions (NSBF) representation for solutions of Sturm-Liouville equations and for their derivatives is obtained. The representation possesses an attractive feature for applications: for all real values of the spectral parameter the difference between the exact solution and the approximate one (the truncated NSBF) depends on (the truncation parameter) and the coefficients of the equation and does not depend on . A similar result is valid when belongs to a strip . This feature makes the NSBF representation especially useful for applications requiring computation of solutions for large intervals of . Error and decay rate estimates are obtained. An algorithm for solving initial value, boundary value or spectral problems for the Sturm-Liouville equation is developed and illustrated on a test problem.

## 1 Introduction

In the recent work [13] a new representation for solutions of the one-dimensional Schrödinger equation

 −u′′+Q(x)u=ω2u (1.1)

was obtained in terms of so-called (see, e.g., [24], [25] and [2]) Neumann series of Bessel functions (NSBF). The representation possesses an attractive feature for applications: for all the difference between the exact solution and the approximate one (the truncated NSBF) depends on (the truncation parameter) and and does not depend on . A similar result is valid when belongs to a strip . This feature makes the NSBF representation especially useful for applications requiring computation of solutions of (1.1) for large intervals of . For example, as was shown in [13], the NSBF representation allows one to compute hundreds or if necessary even thousands of eigendata with a nondeteriorating for large and remarkable accuracy. In [15] the NSBF representation was extended onto perturbed Bessel equations.

In the present work we derive an NSBF representation for solutions of the Sturm-Liouville equation

 −(p(y)v′)′+q(y)v=ω2r(y)v. (1.2)

The coefficients are assumed to admit the application of the Liouville transformation (see, e.g., [6], [26]). The main result consists in an NSBF representation for solutions of (1.2) and for their derivatives, preserving the same attractive feature described above. Error and decay rate estimates are obtained. An algorithm for solving initial value, boundary value or spectral problems for (1.2) is developed and illustrated on a test problem.

Besides this Introduction the paper contains the following sections. Section 2 presents some well known facts concerning the Liouville transformation together with recent results from [11] showing how the system of formal powers associated with (1.2) is transformed by the Liouville transformation. In Section 3 we recall some relevant results from [13] which are used in Section 4 and Section 5 for obtaining the main results of this work, the NSBF representation for solutions of (1.2) as well as for their derivatives. In Section 6 convenient for computation formulas for the coefficients of the NSBF representations are derived. In Section 7 a computational algorithm based on the NSBF representation is formulated and discussed. Section 8 presents some illustrations of its numerical performance on a test problem admitting an exact solution. In Appendix A we prove error and decay rate estimates of the NSBF representation.

## 2 Preliminaries on the Liouville transformation

### 2.1 Definition of the Liouville transformation

Consider the Sturm-Liouville differential equation

 (p(y)v′)′−q(y)v=−λr(y)v, y∈[A,B], (2.1)

with being a finite interval. Let , and be such that , and for all . Define the mapping by

 l(y):=∫yA{r(s)/p(s)}1/2ds, (2.2)

where . Denote . Then (2.1) is related to the one-dimensional Schrödinger differential equation

 u′′−Q(x)u=−λu,x∈[0,b] (2.3)

by the Liouville transformation of the variables and into and defined as (see, e.g., [6], [26])

 x=l(y),u(x):=ρ(y)v(y) for all y∈[A,B]

and the coefficient is given by the relation

 Q(x)=q(y)r(y)+p(y)4r(y)⎡⎣(p′(y)p(y)+r′(y)r(y))′+34(p′(y)p(y))2+12p′(y)p(y)r′(y)r(y)−14(r′(y)r(y))2⎤⎦=q(y)r(y)−ρ(y)r(y)⎡⎣p(y)(1ρ(y))′⎤⎦′ for all y∈[A,B]. (2.4)

The Liouville transformation can be considered as an operator acting according to the rule

 u(x)=L[v(y)]=ρ(l−1(x))v(l−1(x)).

Let us introduce the following notations for the differential expressions

 B=−d2dx2+Q(x) and C=−1r(y)(ddy(p(y)ddy)−q(y)).

The following proposition summarizes the main properties of the operator .

###### Proposition 2.1.
1. The inverse operator is defined by .

2. The uniform norms of the operators and are and .

3. The operator equality

 BL=LC

is valid on .

### 2.2 Transformation of formal powers

Let be a non-vanishing solution of the equation

 f′′−Q(x)f=0,x∈[0,b], (2.5)

such that

 f(0)=1.

On the existence of a non-vanishing see Remark 2.5 below. Note that in general complex valued solutions are considered even when the coefficient is real valued. Denote .

Consider two sequences of recursive integrals (see [10], [12], [14])

 X(0)≡1,X(n)(x)=n∫x0X(n−1)(s)(f2(s))(−1)nds,n=1,2,…

and

 ˜X(0)≡1,˜X(n)(x)=n∫x0˜X(n−1)(s)(f2(s))(−1)n−1ds,n=1,2,….
###### Definition 2.2.

The families of functions and constructed according to the rules

 φk(x)={f(x)X(k)(x),k\ odd,f(x)˜X(k)(x),k\ evenandψk(x)=⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩˜X(k)(x)f(x),k\ odd,X(k)(x)f(x),k\ even

are called systems of formal powers associated with equation (2.3).

###### Remark 2.3.

The formal powers arise in the spectral parameter power series (SPPS) representation for solutions of (2.3) (see [9], [10], [12], [14]).

Analogously, let us introduce a system of formal powers corresponding to equation (2.1).

Let be a solution of the equation

 (p(y)g′)′−q(y)g=0, y∈[A,B], (2.6)

such that for all (see Remark 2.5). Then the following two families of auxiliary functions are well defined

 ˜Y(0)(y) ≡Y(0)(y)≡1, Y(n)(y) ={n∫yAY(n−1)(s)1g2(s)p(s)ds,n odd,n∫yAY(n−1)(s)g2(s)r(s)ds,n even, ˜Y(n)(y) =⎧⎪⎨⎪⎩n∫yA˜Y(n−1)(s)g2(s)r(s)ds,n odd,n∫yA˜Y(n−1)(s)1g2(s)p(s)ds,n even.

Then similarly to Definition 2.2 we define the formal powers associated to equation (2.1).

###### Definition 2.4.

Let satisfy conditions of Subsection 2.1 and be a non-vanishing solution of (2.6). Then the formal powers associated to equation (2.1) are defined for any as follows

###### Remark 2.5.

The existence of a non-vanishing solution of (2.6) for complex valued and such that , and are continuous on was proved in [14, Remark 5] (see also [4]). Moreover, the only reason for the requirement of the absence of zeros of the functions and is to make sure that the formal powers be well defined. As was shown in [16] this is true even when and/or have zeros, though in that case corresponding formulas are slightly more complicated.

###### Theorem 2.6 ([11]).

Let and be functions satisfying the conditions of Definition 2.4, and hence is a particular solution of (2.5). Then the following relations are valid

 ρ(y)Φn(y)=φn(x)and1ρ(y)Ψn(y)=ψn(x) for all n∈N∪{0},

that is

 φn(x)=L[Φn(y)]andψn(x)=L[Ψn(y)/ρ2(y)].

In order that the equality be fulfilled, must be chosen so that

 g(A)=1ρ(A). (2.7)

We will assume this initial condition to be satisfied.

It is easy to see that for one obtains

 h=√p(A)r(A)(g′(A)g(A)+ρ′(A)ρ(A)). (2.8)

## 3 Solution of the one-dimensional Schrödinger equation

Let , and , denote the solutions of (2.3) satisfying the following initial conditions in the origin

 c(ω,0)=1,c′(ω,0)=h,s(ω,0)=0,s′(ω,0)=ω (3.1)

where is an arbitrary complex number.

###### Theorem 3.1 ([13]).

The solutions and of the equation

 u′′−Q(x)u=−ω2u,x∈(0,b)

 c(ω,x)=cosωx+2∞∑n=0(−1)nβ2n(x)j2n(ωx) (3.2)

and

 s(ω,x)=sinωx+2∞∑n=0(−1)nβ2n+1(x)j2n+1(ωx) (3.3)

where stands for the spherical Bessel function of order , the functions are defined as follows

 βn(x)=2n+12(n∑k=0lk,nφk(x)xk−1), (3.4)

where is the coefficient of in the Legendre polynomial of order . The series in (3.2) and (3.3) converge uniformly with respect to on and converge uniformly with respect to on any compact subset of the complex plane of the variable . Moreover, for the functions

 cN(ω,x)=cosωx+2[N/2]∑n=0(−1)nβ2n(x)j2n(ωx) (3.5)

and

 sN(ω,x)=sinωx+2[(N−1)/2]∑n=0(−1)nβ2n+1(x)j2n+1(ωx) (3.6)

the following uniform estimates hold

 |c(ω,x)−cN(ω,x)|≤√2xεN(x)and|s(ω,x)−sN(ω,x)|≤√2xεN(x) (3.7)

for any , , and

 |c(ω,x)−cN(ω,x)|≤εN(x)√sinh(2Cx)Cand|s(ω,x)−sN(ω,x)|≤εN(x)√sinh(2Cx)C (3.8)

for any , belonging to the strip , . Here is a function satisfying as . These estimates are slight refinements for those presented in [13] using the ideas from [15]. Some estimates for are presented in Appendix A.

Inequalities (3.7) and (3.8) reveal an interesting feature of the representations (3.2) and (3.3): the accuracy of approximation of the exact solutions by the partial sums (3.5) and (3.6) does not depend on meanwhile it belongs to a strip in the complex plane. This feature was tested and confirmed numerically in [13]. An analogous result is valid for the derivatives of the solutions. We formulate it in the following statement.

###### Theorem 3.2 ([13]).

The derivatives of the solutions and with respect to admit the following representations

 c′(ω,x)=−ωsinωx+(h+12∫x0Q(s)ds)cosωx+2∞∑n=0(−1)nγ2n(x)j2n(ωx) (3.9)

and

 s′(ω,x)=ωcosωx+12(∫x0Q(s)ds)sinωx+2∞∑n=0(−1)nγ2n+1(x)j2n+1(ωx) (3.10)

where are defined as follows

 γn(x)=2n+12⎛⎜ ⎜ ⎜ ⎜⎝n∑k=0lk,n(kψk−1(x)+f′(x)f(x)φk(x))xk−n(n+1)2x−12∫x0Q(s)ds−h2(1+(−1)n)⎞⎟ ⎟ ⎟ ⎟⎠. (3.11)

The series in (3.9) and (3.10) converge uniformly with respect to on and converge uniformly with respect to on any compact subset of the complex plane of the variable .

Moreover, for the approximations

 ∘cN(ω,x):=−ωsinωx+(h+12∫x0Q(s)ds)cosωx+2[N/2]∑n=0(−1)nγ2n(x)j2n(ωx)

and

 ∘sN(ω,x):=ωcosωx+12(∫x0Q(s)ds)sinωx+2[(N−1)/2]∑n=0(−1)nγ2n+1(x)j2n+1(ωx)

the following inequalities are valid

 ∣∣c′(ω,x)−∘cN(ω,x)∣∣≤√2xε1,N(x)and∣∣s′(ω,x)−∘sN(ω,x)∣∣≤√2xε1,N(x)

for any , , and

 ∣∣c′(ω,x)−∘cN(ω,x)∣∣≤ε1,N(x)√sinh(2Cx)Cand∣∣s′(ω,x)−∘sN(ω,x)∣∣≤ε1,N(x)√sinh(2Cx)C

for any , belonging to the strip , .

###### Remark 3.3 ([13]).

The functions also can be constructed as solutions of the recurrent equations

 1xnB[xnβn(x)]=2n+12n−3xn−1B[βn−2(x)xn−1],n≥1, (3.12)

with the first functions given by and and initial conditions where .

The functions can be calculated from the equalities

 γ−1(x) =14∫x0Q(s)ds, γ0(x) =β′0(x)−h2−14∫x0Q(s)ds=f′(x)2−h2−14∫x0Q(s)ds, γn(x) =nxβn(x)+β′n(x)+2n+12n−3(γn−2(x)−β′n−2(x)+n−1xβn−2(x)),n≥1. (3.13)

One can use as well.

###### Remark 3.4.

Based on the formulas from Remark 3.3 a stable recurrent integration procedure for computing the functions and was proposed in [13] in the following form.

 ηn(x) =∫x0(tf′(t)+(n−1)f(t))σn−2(t)dt,θn(x)=∫x01f2(t)(ηn(t)−tf(t)σn−2(t))dt, (3.14) σn(x) =2n+12n−3[x2σn−2(x)+cnf(x)θn(x)], (3.15) τn(x) =2n+12n−3[x2τn−2(x)+cn(f′(x)θn(x)+ηn(x)f(x))−(cn−2n+1)xσn−2(x)],n≥1, (3.16)

where if and otherwise.

## 4 Solution of the Sturm-Liouville equation

Let us generalize Theorems 3.1 and 3.2 onto solutions of equation (2.1).

###### Theorem 4.1.

Let the functions , and satisfy the conditions from Subsection 2.1 and be a solution of (2.6) satisfying (2.7) and such that for all . Then two linearly independent solutions and of equation (2.1) for can be written in the form

 v1(ω,y)=cos(ωl(y))ρ(y)+2∞∑n=0(−1)nα2n(y)j2n(ωl(y)) (4.1)

and

 v2(ω,y)=sin(ωl(y))ρ(y)+2∞∑n=0(−1)nα2n+1(y)j2n+1(ωl(y)) (4.2)

with defined by (2.2), the coefficients being defined by the equalities

 αn(y)=2n+12(n∑k=0lk,nΦk(y)lk(y)−1ρ(y)), (4.3)

where are from Definition 2.4. The solutions and satisfy the following initial conditions

 v1(ω,A) =1ρ(A),v2(ω,A)=0, (4.4) v′1(ω,A) =g′(A),v′2(ω,A)=ωρ(A)√r(A)p(A). (4.5)

The series in (4.1) and (4.2) converge uniformly with respect to on and converge uniformly with respect to on any compact subset of the complex plane of the variable .

Moreover, for the functions

 vN1(ω,y)=cos(ωl(y))ρ(y)+2[N/2]∑n=0(−1)nα2n(y)j2n(ωl(y)) (4.6)

and

 vN2(ω,y)=sin(ωl(y))ρ(y)+2[(N−1)/2]∑n=0(−1)nα2n+1(y)j2n+1(ωl(y)) (4.7)

the following estimates hold

 ∣∣v1,2(ω,y)−vN1,2(ω,y)∣∣≤√2l(y)εN(l(y))maxy∈[A,B]1|ρ(y)| (4.8)

for any , , and

 ∣∣v1,2(ω,y)−vN1,2(ω,y)∣∣≤εN(l(y))√sinh(2Cl(y))Cmaxy∈[A,B]1|ρ(y)| (4.9)

for any , belonging to the strip , , where is a function from Theorem 3.1.

###### Proof.

Consider solutions (3.2) and (3.3) of (2.3). The functions

and

are then solutions of (2.1) with . In these representations all the magnitudes except the coefficients can be defined with no reference to equation (2.3). Let us show that the same observation is applicable to the functions

 αn(y):=βn(l(y))ρ(y)=L−1[βn(x)],n=0,1,2,….

Note that

 L−1[xku(x)]=lk(y)u(l(y))ρ(y)=lk(y)L−1[u(x)] (4.10)

for any function and any . Then

 αn(y)=2n+12L−1[n∑k=0lk,nφk(x)xk−1]=2n+12(n∑k=0lk,nL−1[φk(x)]lk(y)−1ρ(y)).

Due to Theorem 2.6 we obtain (4.3) and thus the solutions and of (2.1) have the form (4.1), (4.2).

Estimates (4.8) and (4.9) follow from (3.7) and (3.8) by taking into account Proposition 2.1.

Equalities (4.4) and (4.5) are obtained directly from (3.1). Indeed, due to statement 1. from Proposition 2.1, and that gives us (4.4). From the definition of the Liouville transformation we have the equality for and being related by , and thus,

 vy=1ρ(√rpux−ρyv), (4.11)

from where (4.5) are derived. ∎

###### Remark 4.2.

Inequalities (4.8) and (4.9) show that the representations (4.1) and (4.2) preserve the important property of the representations (3.2) and (3.3): the accuracy of approximation of the exact solutions by the partial sums (4.6) and (4.7) does not depend on meanwhile it belongs to a strip in the complex plane.

###### Remark 4.3.

Similarly to the coefficients (see Remark 3.3) the functions can be constructed as solutions of the recurrent equations

 1ln(y)C[ln(y)αn(y)]=2n+12n−3ln−1(y)C[αn−2(y)ln−1(y)] (4.12)

obtained by applying to (3.12) and taking into account (4.10) and the equality . The first functions of the sequence are given by

 α−1=L−1[β−1]=12ρ% andα0=L−1[β0]=12(g−1ρ). (4.13)

The initial conditions satisfied by have the form

 Σn(A)=Σ′n(A)=0. (4.14)

This sequence of equations for leads to a stable recursive integration procedure which is proposed in Section 6.

## 5 Representation of the derivatives of the solutions v1 and v2

Due to (4.11) for and with the aid of (3.9), (3.10) and (4.1), (4.2) we obtain

 v′1(ω,y)=√r(y)p(y)(1ρ(y)(G1(y)cos(ωl(y))−ωsin(ωl(y)))+2∞∑n=0(−1)nμ2n(y)j2n(ωl(y)))−ρ′(y)ρ(y)v1(ω,y) (5.1)

and

 v′2(ω,y)=√r(y)p(y)(1ρ(y)(G2(y)sin(ωl(y))+ωcos(ωl(y)))+2∞∑n=0(−1)nμ2n+1(y)j2n+1(ωl(y)))−ρ′(y)ρ(y)v2(ω,y) (5.2)

where

 μn(y):=γn(l(y))ρ(y)=L−1[γn(x)],n=0,1,2,…,

and the functions and have the form (cf. [11, Remark 4.7])

 G1(y):=h+12∫l(y)0Q(s)ds=h+12∫yA1(pr)1/4(q(pr)1/4−[p{(pr)−1/4}′]′)(s)ds=h+ρρ′2r∣∣∣yA+12∫yA[qρ2+(ρ′)2r](s)ds (5.3)

and

 G2(y):=12∫l(y)0Q(s)ds=G1(y)−h. (5.4)

Let us obtain a formula for calculating which would not involve magnitudes related to equation (2.3) but only those related to (2.1). Consider (3.11). Direct calculation gives us the relation

 f′(x)f(x)=√p(y)r(y)(g′(y)g(y)+ρ′(y)ρ(y)).

Hence, due to Theorem 2.6,

 kψk−1(x)+f′(x)φk(x)/f(x)=kΨk−1(y)ρ(y)+ρ(y)√p(y)r(y)(g′(y)g(y)+ρ′(y)ρ(y))Φk(y).

Thus,

 μn(y)=2n+12ρ(y)(n∑k=0lk,nlk(y)(kΨk−1(y)ρ(y)+ρ(y)√p(y)r(y)(g′(y)g(y)+ρ′(y)ρ(y))Φk(y))−n(n+1)2l(y)−G2(y)−h2(1+(−1)n)). (5.5)

This is a direct formula for calculating in terms of the functions and . A recurrent formula for can be obtained from (3.13). We have

 μ−1(y) =G2(y)2ρ(y), (5.6) μ0(y) =√p(y)r(y)(α′0(y)+ρ′(y)ρ(y)α0(y))−G1(y)2ρ(y)=√p(y)r(y)(g(y)ρ(y))′2ρ(y)−G1(y)2ρ(y), (5.7)

and

 μn(y)=nαn(y)l(y)+√p(y)r(y)(α′n(y)+ρ′(y)ρ(y)αn(y))+2n+12n−3(μn−2(y)−√p(y)r(y)(α′n−2(y)+ρ′(y)ρ(y)αn−2(y))+n−1l(y)αn−2(y)),n=1,2,…. (5.8)

One may also use

 μ1(y)=α1(y)l(y)+√p(y)r(y)(α′1(y)+ρ′(y)ρ(y)α1(y))−3G2(y)2ρ(y). (5.9)

## 6 Computation formulas for the coefficients

Besides the direct formulas (