We study the distribution of eigenspectra for operators of the form with self-adjoint boundary conditions on both bounded and unbounded interval domains. With integrable potentials , we explore computational methods for calculating spectral density functions involving cases of discrete and continuous spectra where discrete eigenvalue distributions approach a continuous limit as the domain becomes unbounded. We develop methods from classic texts in ODE analysis and spectral theory in a concrete, visually oriented way as a supplement to introductory literature on spectral analysis. As a main result of this study, we develop a routine for computing eigenvalues as an alternative to , resulting in fast approximations to implement in our demonstrations of spectral distribution.

Introduction

We follow methods of the texts by Coddington and Levinson [1] and by Titchmarsh [2] (both publicly available online via archive.org) in our study of the operator and the associated problem

(1)

where on the interval with real parameter and boundary condition

(2)

for fixed , where . For continuous (the set of absolutely integrable functions on ), we study the spectral function associated with (1) and (2) using two main methods: First, following [1], we approximate by step functions associated with related eigenvalue problems on finite intervals for some sufficiently large positive ; then, we apply asymptotic solution estimates along with an explicit formula for spectral density [2]. For some motivation and clarification of terms, we recall a major application: For certain solutions of (1) and (2) and for any (the set of square-integrable functions on ), a corresponding solution to (1) may take the form

where

(in a sense described in Theorem 3.1 of Chapter 9 [1]); here, is said to be a spectral transform of . By way of such spectral transforms, the differential operator may be represented alternatively in the integral form

where induces a measure by which (roughly, the set of square-integrable functions when integrated against ) and by which Parseval’s equality holds. Typical examples are the complete set of orthogonal eigenfunctions for and the corresponding Fourier sine transform in the limiting case (cf. Chapter 9, Section 1 [1]).

For a fixed, large finite interval , we consider the problem (1), (2) along with the boundary condition

(3)

(), which together admit an eigensystem with correspondence

where the eigenvalues satisfying and where the eigenfunctions form a complete basis for . Since the associated spectral function is a step function with jumps at the various , we first estimate these by way of a related equation arising from Prüfer (phase-space) variables and compute the corresponding jumps .

Then, we use interpolation to approximate the continuous spectral function using data from a case of large at points and using

(4)

imposing the condition for all .

We compare our results with those of a well-known formula [2] appropriate to our case on , which we outline as follows: For fixed , let be the solution to (1) with boundary values

for which the asymptotic formula

(5)

holds as . Then we have

(6)

from Section 3.5 [2].

Finally, in the last section, we apply the above techniques to extend our study to operators on large domains and on , where spectral matrices take the place of spectral functions as a matrix analog of spectral transforms on these types of intervals (cf. equation (5.5) [1]). The techniques are described in detail below, but it is of particular interest that our computations uncover an interesting pattern in a discrete-spectrum case, as we are forced to reformulate our approach according to certain eigen-subspaces involved: our desired spectral approximations are resolved by way of an averaging procedure in forming Riemann sums.

Various sections of Chapters 7–9 [1] (see also [3] and related articles) present useful introductory discussion applied to material presented in this article; yet, with our focus on equations (1)–(6), one may proceed given basic understanding of Riemann–Stieltjes integration along with knowledge of ordinary differential equations and linear algebra, commensurate with (say) the use of and .

An Eigenvalue Estimator

We compute eigenvalues by first computing solutions on to the following, arising from Prüfer variables (equation 2.4, Chapter 8 [1]):

(7)

Here, , where is a nontrivial solution to (1), (2) and (3) and satisfies

(8)

for positive integers . We interpolate to approximate such solutions as an efficient means to invert (8) in the variable . And we use the following function on (7) throughout this article.

Consider an example with , , and potential for parameter with , , in the case , .

We create an interpolation approximation for eigenvalues .

It is instructive to graphically demonstrate the theory behind this method. Here, we consider the eigenvalues as those values of where the graph of intersects the various lines as we use to find (or ), our maximum index , depending on .

We choose these boundary conditions so that we may compare our results with those of applied to the corresponding problem (1) and (2) using .

We now compare and contrast the methods in this case. The percent differences of the corresponding eigenvalues are all less than 0.2%, even within our limits of accuracy.

In contrast, our interpolation method allows some direct control of which eigenvalues are to be computed, whereas (in the default setting) outputs a list up to 39 values, starting from the first. Moreover, our method admits nonhomogeneous boundary conditions, where admits only homogeneous conditions, Dirichlet or Neumann.

Spectral Density: Discrete Approximation

We proceed to build our approximate spectral density function for the problem (1) and (2) on with the same potential as above. We compute eigenvalues likewise but now on a larger interval for and with nonhomogeneous boundary conditions, say given by , (albeit does not depend on ).

We compute eigenvalues via our interpolation method and compute a minimum (or ) as well as a maximum index so as to admit only positive eigenvalues; is supported on and negative eigenvalues result in dubious approximations by .

We now compute the values .

Fitting Method

We now apply the method of [2] as outlined in equation (6). We use to include data from an interval near the endpoint that includes at least one half-period of the period of the fitting functions and .

The function may return non-numerical results among the first few, in which case we recommend that either or be readjusted or that be set large enough to disregard such results.

We now compare our results of the discrete and continuous (asymptotic fit) spectral density approximations.

We compare the results by plotting percent differences, all being less than 0.1%.

Check with Exact Calculation

We chose as above because, in part, the solutions can be computed in terms of well-known (modified Bessel) functions. Replacing by , for , the solutions are linear combinations of

(9)

From asymptotic estimates (cf. equation 9.6.7 [4]), we see that the former is dominant and the latter is recessive as when . Then, from Chapter 9 [1], equation 2.13 and Theorem 3.1, we obtain the density function by computing

(10)

where is a solution as above and is a solution with boundary values , . (Here, is commonly known as the Titchmarsh–Weyl -function.) In the following code, we produce the density function in exact form by replacing functions from (9), the dominant by 1 and the recessive by 0, to compute the inside limit and thereafter simply allowing to be real.

We likewise compare the exact formula for the continuous spectrum with the discrete results, noting that the exact graph appears to essentially be the same as that obtained by our asymptotic fitting method (not generally expecting the fits to be accurate for small !).

Extension to Unbounded Domains: A Proof of Concept

For the operator we now extend our study to large domains in the discrete-spectrum case and to the domain in the continuous-limit case. We choose an odd function potential of the form for positive constants , . We focus on the spectral density associated with specific boundary values at and an associated pair of solutions to (1): namely, we consider expansions in the pair and such that

(11)

We apply the above computational methods to the analytical constructs from Chapter 5 [1] in both the discrete and continuous cases. First, for the discrete case, we compute spectral matrices associated with self-adjoint boundary-value problems and the pair as in (11): We estimate eigenvalues for an alternative two-point boundary-value problem on for (moderately) large to compute the familiar jumps of the various components . These components induce measures that appear in the following form of Parseval’s equality for square-integrable functions on (taken in a certain limiting sense):

(real-valued case). Second, we compute the various densities as limits as by the formulas

(12)

where and are certain limits of -functions, related to equation (10), but for our ODE problem on domains and , respectively. The densities are computed by procedures more elaborate than (6), as discussed later. Then, we compare results of the discrete case like in (4), approximating

(13)

Discrete Case

After choosing (self-adjoint) boundary conditions (of which the limits happen to be independent)

(14)

on an interval , we estimate eigenvalues and compute coefficients , from the linear combinations

for the associated orthonormal (complete) set of eigenfunctions ; , whereby

(real-valued case). Here, the functions result by normalizing eigenfunctions satisfying (14) so that we obtain

We are ready to demonstrate. Let us choose , and , (arbitrary). Much of the procedure follows as above, with minor modification, as we include to obtain the values and (the next result may take around three minutes on a laptop).

We now approximate the density functions by plotting where

(15)

(for certain ) as we compute the difference quotients at the various jumps, over even and odd indices separately, and assign the corresponding sums to the midpoints of corresponding intervals .

We give the plots below, in comparison with those of the continuous spectra, and give a heuristic argument in the Appendix as to why this approach works.

Continuous Case

First, we apply the asymptotic fitting method using the solutions and . Here, we have to compute full complex-valued formulas for the corresponding -functions (cf. Section 5.7 [2]) where a slight modification of the derivation of , via a change of variables and a complex conjugation, results in (See Appendix).

We now compare the result of the discrete and asymptotic fitting methods for the elements .

Appendix

We have deferred some discussion on our use of , comparison of eigenvalue computations, discrete eigenspace decomposition and Weyl -functions to this section.

First, we have used to suppress messages warning that some solutions may not be found. From Chapter 8 [1], we expect unique solutions since the functions are strictly increasing. We have also used to suppress various messages from and other related functions regarding small values of to be expected with short-range potentials and large domains.

Second, our formulation of and the midpoints as in (15) arises from a decomposition of the eigenspace by even and odd indices. We motivate this decomposition by an example plot of the values , where the dichotomous behavior is quite pronounced, certainly for large .

We are thus inspired to compute the quotients over even and odd indices separately. Then, we consider, say, a relevant expression from Parseval’s equality: for appropriate Fourier coefficients , , associated with respective solutions , we write

We suppose that and for the corresponding transforms in the limit . Of course, a rigorous argument is beyond the scope of this article.

Finally, we elaborate on the calculations of the -functions and : Given the asymptotic expressions

as (resp.), we follow Section 5.7 of [2], making changes as needed, with a modification via complex conjugation (, say) for to arrive at

Acknowledgments

The author would like to thank the members of MAST for helpful and motivating discussions concerning preliminary results of this work in particular and Mathematica computing in general.

References

[1]	E. A. Coddington and N. Levinson, Theory of Ordinary Differential Equations, New York: McGraw-Hill, 1955. archive.org/details/theoryofordinary00codd.
[2]	E. C. Titchmarsh, Eigenfunction Expansions Associated with Second-Order Differential Equations, 2nd ed., London: Oxford University Press, 1962. archive.org/details/eigenfunctionexp0000titc.
[3]	E. W. Weisstein. “Operator Spectrum” from MathWorld–A Wolfram Web Resource. mathworld.wolfram.com/OperatorSpectrum.html.
[4]	M. Abramowitz and I. A. Stegun, eds., Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, New York: Wiley, 1972.
C. Winfield, “From Discrete to Continuous Spectra,” The Mathematica Journal, 2019. https://doi.org/10.3888/tmj.21-3.

About the Author

C. Winfield holds an MS in physics and a PhD in mathematics and is a member of the Madison Area Science and Technology amateur science organization, based in Madison, WI.

Christopher J. Winfield
Madison Area Science and Technology
3783 US Hwy. 45
Conover, WI 54519
cjwinfield2005@yahoo.com