Main Page

The Sturm-Liouville problem

We consider Sturm-Liouville eigenvalue problems on the real line,

$-\frac{d^2u}{dr^2} + Vu = \lambda u,$

where the boundary conditions are that the solutions should be square-integrable. Furthermore we consider potentials $ V $ that are even functions such that $V(r) \to \infty$ as $|r| \to \infty.$

In the "harmonic" case, $ V(r) = r^2, $ the eigenfunctions are the so-called Hermite functions,

$h_n(r) = H_n(r)\exp(-r^2/2), \qquad n = 0, 1,\ldots,$

where $ H_n $ are the Hermite polynomials. The codes in this project explore convergence rates for different potentials $ V $ using the initial sets $h_0, \ldots, h_n$ of Hermite functions for Galerkin bases. In fact, instead of the $ h_n, $ a basis of "normalized" Hermite functions $\phi_n$ satisfying $\int \phi_n^2 = 1$ will be used.

Multiple-precision computations

Presumably the dependence of the convergence rates on the number of basis functions is asymptotic in such a way that it is not numerically observable until $ n$ is beyond a certain minimal threshold. For this reason it is necessary to perform the computations in multiple-precision arithmetic, that is, with floating-point types whose precision can be increased as needed. Two numerical libraries used are:

the MPFR Library, which defines a multiple-precision type in C, and for which there are C++ wrappers which overload the standard arithmetic operations for this type,
the Template Numerical Toolkit, which implements the numerical algorithms from LAPACK for a floating-point type that is passed as a C++ template parameter.

For certain polynomial potentials such as $V(r) = r^\ell$ for $\ell = 2, 4, \ldots,$ the Galerkin discretization of the elliptic operator can be assembled "analytically." For for most potentials, however, the entries of the potential matrix,

$P_{mn} = \int \phi_m \phi_n V~dr,$

must be computed via numerical quadrature. Development of the quadpack++ library was initiated toward this purpose and remains ongoing: this is an implementation of the adaptive routines from QUADPACK for a floating-point type that is passed as a C++ template parameter.

Experiments

The most direct way of examining convergence is to discretize an elliptic operator for two different basis sizes, say $ n $ and $ n+1, $ and to examine the difference $\epsilon_k(n) = \lambda_k(n) - \lambda_k(n+1)$ in the $ k$ -th eigenvalue arising from the respective discretizations. (Assumptions on the potentials ensure that all eigenvalues are non-degenerate, hence such a comparison is possible.) The rate at which $\epsilon_k(n) \to 0$ as $ n $ increases is the quantity of interest.

Since most of the numerical effort is consumed by the computation of the integrals $P_{mn},$ an efficient approach is to assemble just one matrix $ A $ for a reasonably large basis and then to diagonalize each upper left $M\times M$ block, as such any such block is the discretization arising from the first $ M $ basis functions. Programs were written for the potentials

$V_1(r) = |r|^\alpha, \qquad V_2(r) = \exp((1+r^2)^{1/2}).$

Source codes are available by email inquiries.

Internals: the Galerkin discretizations

All functions directly related to the Galerkin discretization are listed under Hermite approximation.

The basis of normalized Hermite functions $\phi_n$ is implemented in Normal_Hermite_Ftns. Assembly of the Galerkin discretization requires computing two kind of integrals,

$K_{mn} = \int \phi_m' \phi_n'~dr, \qquad P_{mn} = \int \phi_m \phi_n V~dr,$

and the respective assembly functions are Hermite_Kinetic and Hermite_Potential. Fortunately there is a closed form expression for the entries $K_{mn}.$

For potentials $ V $ whose growth as $|r| \to \infty$ is dominated by the exponential $\exp(-r^2),$ the integrand of each $P_{mn}$ is non-singular at infinity. This permits the use of "basic" adaptive quadrature on a finite interval after transforming the integral $P_{mn}$ by a variable transform. This is the purpose of the function Potential_QAGI; it defines the variable-transformed integrand.

The integrals $P_{mn}$ have a closed form expression for $V(r) = r^\ell$ for $\ell = 0, 2, \ldots.$ (The matrix is diagonal when $\ell = 2$ since the $ h_n$ are the eigenfunctions.) For experimental purposes the assembly routine Hermite_Moment was written for this special case.

Source codes are available by email inquiries.

References

Gerald Folland, Fourier Analysis and its Applications, Brooks/Cole Publishing (1992).