Accurate Spectral Collocation Computation of High Order Eigenvalues for Singular Schr\"odinger Equations

We are concerned with the study of some classical spectral collocation methods as well as with the new software system Chebfun in computing high order eigenpairs of singular and regular Schrodinger eigenproblems. We want to highlight both the qualities as well as the shortcomings of these methods and evaluate them in conjunction with the usual ones. In order to resolve a boundary singularity we use Chebfun with domain truncation. Although it is applicable with spectral collocation, a special technique to introduce boundary conditions as well as a coordinate transform, which maps an unbounded domain to a finite one, are the special ingredients. A challenging set of \hard"benchmark problems, for which usual numerical methods (f. d., f. e. m., shooting etc.) fail, are analyzed. In order to separate \good"and \bad"eigenvalues we estimate the drift of the set of eigenvalues of interest with respect to the order of approximation and/or scaling of domain parameter. It automatically provides us with a measure of the error within which the eigenvalues are computed and a hint on numerical stability. We pay a particular attention to problems with almost multiple eigenvalues as well as to problems with a mixed spectrum.


Introduction
There is clearly an increasing interest to develop accurate and efficient methods of solution to singular Schrödinger eigenproblems. Our main interest here is to compare the capabilities of the new Chebfun package with those of classical spectral methods in solving such problems, showing their capabilities and weaknesses. The latter employ basis functions and/or grid points based on Chebyshev, Laguerre or Hermite polynomials as well as on sinc or Fourier functions. The effort expended by both classes of methods is also of real interest. It can be assessed in terms of the ease of implementation of the methods as well as in terms of computer resources required to achieve a specified accuracy.
For problems on the entire real axis SiC proved to be particularly well suited. Moreover, this method has given excellent results recorded in our contribution [8] and in the works cited there. The so-called generalized pseudospectral (GPS) method, actually the Legendre collocation, is employed in [10] to calculate the bound states of the Hulthén and the Yukawa potentials in quantum mechanics, with special emphasis on higher excited states and stronger couplings. The author uses a two parameter dependent nonlinear transformation in order to map the half-line into the canonical interval [−1, 1] . In contrast, we will use an analogous transformation but which depends on only one parameter. Also, very recently spectral methods based on non-classical orthogonal polynomials have been used in [11] in order to solve some Schrödinger problems connected with Fokker-Planck operator. A particular attention will be paid in this paper to the challenging issue of continuous spectra vs. discrete (numerical) eigenvalues. It is well known that some eigenvalue problems (see for instance the well known text [12]) for differential operators which are naturally posed on the whole real line or the half-line, often lead to some discrete eigenvalues plus a continuous spectrum. Actually, the usual numerical approximation typically involves three processes: 1. reduction to a finite interval; 2. discretization; 3. application of a numerical eigenvalue solver.
Reduction to a finite interval and discretization typically eliminate the continuous spectrum. Even if we do not truncate a priory the domain on which the problem is formulated such an inherent reduction can not be avoided.
It can be argued that, generally speaking, in solving various differential problems, the Chebfun software provides a greater flexibility than the classical spectral methods. This fact is fully true for regular problems.
Unfortunately, in the presence of various singularities, the maximum order of approximation N, of the unknowns can be reached (N ≥ 4000) and then Chebfun issues a message that warns about the possible inaccuracy of the results provided.
We came out of this tangle using modified classical spectral methods. In this way, when we had serious doubts about the accuracy of the solutions given by Chebfun, we managed to establish the correctness of the numerical results.
As a matter of fact, in order to resolve a singularity on the ends of an unbounded integration interval, Chebfun uses only the arbitrary truncation of the domain. Classical spectral methods can also use this method, but it is not recommended. For singular points at finite distances (mainly origin) we will use the so-called removing technique of independent boundary conditions. The boundary conditions at infinity can be enforced using basis functions that satisfy these conditions (Laguerre, Hermite, sinc). An alternative method, which proved to be very accurate, is the Chebyshev collocation (ChC) method in combination with a change of variables (coordinates) which transform the half line into the canonical Chebyshev interval [−1, 1] . Then the removing technique of independent boundary conditions is essential in order to remove the singularities at the end points of integration interval.
A Chebfun code and two MATLAB codes, one for ChC and another for SiC method, are provided in order to exemplify. With minor modifications they could be fairly useful for various numerical experiments.
The structure of this work is as follows. In Section 2 we recall some specific issues for the regular as well as singular Schrödinger eigenproblems. The main comment refers to the notion of mixed spectrum. In Section 3 we review on the Chebfun structure and the classical spectral methods (differentiation matrices, enforcing boundary conditions, etc.). The Section 4 is the central part of the paper. Here we analyze some benchmark problems. In order to separate the "good"from the "bad"eigenvalues we estimate their relative drift with respect to some parameters. The accuracy in computing eigenfunctions is estimated by their departure from orthogonality. We end up with Section 5 where we underline some conclusions and suggest some open problems.

Regular and singular Schrödinger eigenproblems
The Schrödinger equation reads It is a Liouville normal form of a general Sturm-Liouville (SL) equation where λ is proportional with the energy levels of the physical system, q (x) is directly proportional with the potential energy and the "wave function"u may be real or complex such that u u * dx = |u| 2 dx is the probability that the particle under consideration will be "observed"in the interval (x, x + dx) . In problems involving Schrödinger equations, it is customary among chemists and physicists to define the spectrum of this Sturm-Liouville problem as the all eigenvalues λ for which eigenfunctions u exist. The set of isolated points (if any) in this spectrum is called the discrete spectrum; the part (if any) that consists of entire interval is called continuous spectrum. We shall adopt this suggestive terminology here. The equation (1) can be given on a finite, semi-infinite, or infinite interval. Only on a closed and finite interval a ≤ x ≤ b can Schrödinger equation be associated with a regular Sturm-Liouville problem. If the interval of definition is semi-infinite or infinite, or is finite and q (x) vanishes at one or both endpoints, or if q is discontinuous, we can not obtain from (1) a regular Sturm-Liouville problem. In any such case, the Schrödinger equation (1) is called singular. We obtain a singular eigenproblem from a singular Schrödinger equation by imposing suitable homogeneous boundary conditions. They can not always be described by formulae like αu (e) + βu ′ (e) = 0, where e can be the end point a or b. For instance, the condition that u be bounded near a singular end point, which can be finite or ±∞, is a common boundary condition defining a singular eigenproblem. For regular SL problems, it is proved (see for instance [12]) that the spectrum is always discrete, and the eigenfunctions are (trivially) square-integrable. For singular problems the situation is completely different. For instance in their textbook [12] Birkhoff and Rota consider the eigenproblem attached to the free particle equation and show that the spectrum of the free particle is continuous. Some software packages have been designed over time to solve various singular SL problems. The most important would be SLEIGN and SLEIGN2, SLEDGE, SL02F and MATSLISE. The SLDRIVER interactive package supports exploration of a set of SL problems with the four previously mentioned packages. In [13] (see also [14]) the authors designed the software package SLEDGE. They observed that for a class of singular problems their method either fails or converges very slowly. Essentially, the numerical method used in this software package replaces the coefficient function q(x) by step function approximation. Similar behavior has been observed on the NAG code SL02F introduced in [15] and [16] as well as on the packages SLEIGN and SLEIGN2 introduced in [17] and [18]. The MATSLISE code introduced in [19] can solve some Schrödinger eigenvalue problem by a constant perturbation method of a higher order.
The main purpose of this paper is to argue that Chebfun, along with the spectral collocation methods, can be a very feasible alternative to these software packages regarding accuracy, robustness as well as simplicity of implementation. In addition, these methods can compute the "whole"set of eigenvectors and provide some details on the accuracy and numerical stability of the results provided.

Chebfun
The Chebfun system, in object-oriented MATLAB, contains algorithms which amount to spectral collocation methods on Chebyshev grids of automatically determined resolution. Its properties are briefly summarized in [2]. In [1] the authors explain that chebops are the fundamental Chebfun tools for solving ordinary differential (or integral) equations. One may then use them as tools for more complicated computations that may be nonlinear and may involve partial differential equations. This is analogous to the situation in MATLAB itself. The implementation of chebops combines the numerical analysis idea of spectral collocation with the computer science idea of lazy or delayed evaluation of the associated spectral discretization matrices. The grammar of chebops along with a lot of illustrative examples is displayed in the above quoted paper as well as in the text [5]. Thus one can get a suggestive image of what they can do.
In [1] p.12 the authors explain clearly how the Chebfun works, i.e., it solves the eigenproblem for two different orders of approximation, automatically chooses a reference eigenvalue and checks the convergence of the process. At the same time, it warns about the possible failures due to the high non-normality of the analyzed operator (matrix).
Actually, we want to show in this paper that Chebfun along with chebops can do much more, i.e., can accurately solve highly (double) singular Schrödinger eigenproblems.

ChC, LGRC and SiC
In the spectral collocation method the unknown solution to a differential equation is expanded as a global interpolant, such as a trigonometric or polynomial interpolant. In other methods, such as f. e. and f. d., the underlying expansion involves local interpolants such as piecewise polynomials. This means that the accuracy of spectral collocation is superior. For problems with smooth solutions convergence rates are typically of order e −cN or e −c √ N where N is the order of approximation or resolution, i.e., the number of degree of freedom in expansion. In contrast, f. e. or f. d. yield convergence rates that are only algebraic in N , typically of orders N −2 or N −4 . The net superiority of global spectral methods on local methods is discussed in detail in [20]. In all spectral collocation methods designed so far we have used the collocation differentiation matrices from the seminal paper [9]. We preferred this MATLAB differentiation suite for the accuracy, efficiency as well as for the ingenious way of introducing various boundary conditions.
In order to impose (enforce) the boundary conditions we have used two methods that are conceptually different, namely the boundary bordering as well as the basis recombination. A very efficient way to accomplish the boundary bordering is available in [21] and is called removing technique of independent boundary conditions. We have used this technique in the large majority of our papers except [22] where the latter technique has been employed. In the last quoted paper a modified Chebyshev tau method based on basis recombination has been used in order to solve an Orr-Sommerfeld problem with an eigenparameter dependent boundary condition. Even in eigenproblems that contain the spectral parameter in the definition of the boundary conditions we have used the boundary bordering technique (see [23]).
In [24] (see also [25]) we have solved some multiparameter eigenproblems (MEP) which come from separation of variables, in several orthogonal coordinate systems, applied to the Helmholtz, Laplace, or Schrödinger equation. Important cases include Mathieu's system, Lame's system, and a system of spheroidal wave functions. We show that by combining spectral collocation methods, ChC and LGRC, and new efficient numerical methods for solving algebraic MEPs, it is possible to solve such problems both very efficiently and accurately. We improve on several previous results available in the literature, and also present a MATLAB toolbox for solving a wide range of problems.

The drift of eigenvalues
Two techniques are used in order to eliminate the "bad"eigenvalues as well as to estimate the stability (accuracy) of computations. The first one is the drift, with respect to the order of approximation or the scaling factor, of a set of eigenvalues of interest. In a simplified form this concept has been introduced by J. P. Boyd in [26]. The second one is based on the check of the eigenvectors' orthogonality.
In other words, we want to separate the "good"eigenvalues from the "bad"ones, i.e., inaccurate eigenvalues. An obvious way to achieve this goal is to compare the eigenvalues computed for different orders of some parameters such as the approximation order (cut-off parameter) N or the scaling factor. Only those whose difference or "resolution-dependent drift"is "small"can be believed. Actually, in [26] the so called absolute (ordinal) drift with respect to the order of approximation has been introduced.
We extend this definition to the following one. The absolute (ordinal) drift of the jth eigenvalue with respect to the parameter α is defined as where λ (α) j is the jth eigenvalue, after the eigenvalues have been sorted, as computed using a specific value of the parameter. In the most common cases this parameter can be N or c.
The dependence of δ j,absolute,α , j = 1, 2, . . . , N e, where N e is the number of analyzed eigenvalues, on the index (mode) j will be displayed in a log-linear plot. If we divide the right hand side of (2) by λ (α 1 ) j we get the so called relative drift denote by δ j,relative,α .

Numerical benchmark problems and discussions 4.1 A regular Schrödinger eigenproblem
The bounded Coffey-Evans potential reads We attach to equation (1) the homogeneous Dirichlet boundary conditions u (±π/2) = 0 and use β := 30. In spite of being regular, the Coffey-Evans problem is one of the most difficult test problems in the literature because there are very close eigenvalue triplets as β increases.
We have used a short Chebfun code in order to solve this regular eigenproblem. It is available in the next lines. With appropriate changes, this code can be used to analyze any other Schrödinger eigenproblem. The eigenvalues obtained by ChC and Chebfun are extremely close, practically indistinguishable. They also compare very well with those computed in [27] by some coefficient approximation methods of orders 2, 4 and 8 and reported in Table 3 of this paper.
The absolute drift reported in Fig. 1 means that we can compute the first hundred eigenvalues with better accuracy than 10 −10 . Unfortunately no accuracy analysis is reported in [27] (see also [19] and [28]). The eigenvectors are either symmetric (the even ones) or anti-symmetric (the odd ones). The first four of them are depicted in Fig. 2. They look fairly smooth and satisfy the boundary conditions. The Chebyshev coefficients of the first four eigenvectors of Coffey-Evans problem computed by Chebfun are displayed in Fig. 3. These coefficients decrease sharply and smoothly to a rounding-off plateau below 10 −15 . Roughly speaking this means they are computed with the machine precision. For this potential the first eigenvalue λ0 is close to zero (actually we have got λ0 = −6.254 959 429 708 980 · 10 −12 ) and there are very close eigenvalue triplets (λ2, λ3, λ4), (λ6, λ7, λ8), . . . as β increases. The common numeric part of the eigenvalues in the first triplet is 2.31664929 and that of the eigenvalues in the second is a little shorter, i.e. 4.45283.

Chebyshev coefficients
Rounding-off plateau

Two singular Schrödinger eigenproblem on the half line
In this section we study the mixed spectrum of some Schrödinger eigenproblems having a "potential well"dying out at infinity.

Hydrogen atom equation
The first example consists in the equation (1) equipped with the potential along with the boundary conditions The problem (1)-(4)-(5) is clearly singular. We must mention from the beginning that our numerical experiments performed with LGRC and Chebfun together with the truncation of the domain did not produce satisfactory results. In these conditions we have resorted to the mapped ChC method.
In order to implement this method we use the algebraic map which, for each c, transforms the interval [−1, 1] into the half line, and its inverse. The parameter c is free to be tuned for optimum accuracy. The mapping (6) has been introduced in [29] where its practical effects have been discussed. The author observed that the convergence of the Chebyshev expansion is governed by the closeness of the singularities of the function being expanded to the expansion region. The major effect of such mapping is to allow us to move the singularities further away from the expansion region. The mapping may also weaken the effect of the singularities by modifying the strength of the singularity as well as moving it. The value of the scaling parameter c used in our computation was chosen essentially by trial and error along with the drift with respect to this parameter.
In order to write down any second order differential equation in independent variable s, we need the following derivatives: Now it is easy to write the differential equation for u (s) and to attach to this new equation the homogeneous Dirichlet boundary conditions u (±1) = 0. As usual we implement these conditions by deleting the first and the last rows and columns of the collocation matrix attached to the left hand side of the equation (1)- (4). This is the most simplified version of removing technique of the independent boundary conditions introduced in [21]. With (7) the following MATLAB code solve this problem. The first four vectors of the problem are displayed in Fig. 4. It is clear that they satisfy both boundary conditions but in the right neighborhood of origin they have a totally different behavior from the eigenvectors of regular problems, i.e., they vanish out on continuous portions and not in discrete points. However, they clearly approximate square-integrable eigenfunctions and thus confirm some theoretical results proved in [12]. Their Chebyshev coefficients obtained using FCT (fast Chebyshev transformsee [30] for details) are displayed in Fig.5. They decrease sharply and smoothly to some limits, followed by a wide rounding-off plateau. For the first vector (the rightmost one) this limit is around 10 −13 . It increases with the index of the vector. Roughly this means that we cannot hope for a better approximation than something of the order 10 −13 when computing the first eigenvector.
When N := 2048 mapped ChC has found λ0 = −6.250 000 000 166 379· 10 −02 which is a very good approximation of −1/16. The largest negative eigenvalue has been λ66 = −7.865 782 521 027 431 · 10 −07 and the next eigenvalue, the smallest positive has been computed as λ67 = 6.517 056 834 998 433·  In the monograph [12] Sect. 18 it is proved that if the potential q (x) is continuous and has the asymptotic behaviour q (x) = A x + O 1 x 2 as x → ∞, for λ > 0 the spectrum is continuous and the eigenfunctions are not square-integrable and for λ < 0 the spectrum is discrete and the eigenfunctions are square-integrable. The numerical results gathered around this problem plainly confirm this analytical result.
Actually the eigenvector corresponding to the first positive eigenvalue (the sixty-seventh one) is depicted in the left panel of Fig. 7. It is hard to believe that this could approximate a square-integrable eigenfunction. In the right panel of the same figure we displayed the Chebyshev coefficients of the corresponding eigenvector. The oscillations of these coefficients mimic the steep oscillations of the eigenvector.
For the real scaling factor c we have to mention that in case of eigenvalue problems it can be adjusted only on the mathematical basis. Thus, it can be tuned in order to: • improve the decaying rate of the coefficients of spectral expansions.

Potential with a Coulomb type decay
Let's consider now the Schrödinger equation (1) equipped with the potential (8) and supplied with boundary conditions (5). Some eigenvalues computed by ChC method when c := 2 and N := 512 are compared in Table 3 with their counterparts computed by LGRC and perturbation methods. If for the first eigenvalues the coincidence is excellent, the same does not happen for higher indices. That is why in Fig. 8 in a log-linear plot we illustrate the absolute drift in the case of the ChC method. It is observed that, for instance for the tenth mode, the difference between the values of this eigenvalue calculated with two very different orders of approximation N is of the order 10 −10 . This leads us to believe that the eigenvalues computed with ChC are the most accurate. It is also clear from this figure that in the second case, i.e., N1 := 760 and N2 := 1024 the drift oscillates less than in the first case reported in this figure. Table 3: The first five eigenvalues and the tenth one of the Schrödinger problem (1)-(5) when the potential has a Coulomb-type decay (8), computed by three different methods. j λ j computed by LGRC in [8] λ j computed in [28] λ where µ and ν are real parameters. This is a more general (anharmonic) oscillator than the simpler harmonic one. Two behavioral boundary conditions requiring the boundedness of the solutions at large distance, i.e., x → ±∞ are attached to this equation. Actually we impose the conditions We have to observe that these conditions are automatically satisfied in spectral collocation based on Laguerre, Hermite and sinc functions.
SiC with N := 500 and scaling factor h := 0.1 has been used in order to produce the following results. The fist four eigenvectors of Schrödinger eigenproblem (1)-(10) with potential (9) are displayed in the left panel of Fig. 9. Their coefficients are illustrated in the right panel of the same figure. These coefficients symmetrically decrease to approximations of at least 10 −12 which means a reasonable accuracy of the method. We also have computed some high-lying eigenvalues namely λ100, λ150 and λ200, for ν := 1 and µ := 500. They have respectively the numerical values 199.001994801512, 301.001995781805, and 403.001995224433.
In [31] the author solved this problem numerically for general ν and µ. He used a very elegant variational argument of Ritz type, based on Hermite functions, and observed that for large µ and fixed ν the eigenvalues of this problem approximate those of the harmonic oscillator. We have confirmed this observation even for higher index eigenvalues. In the left panel of Fig. 10 we display the relative drift of the first 200 eigenvalues of the problem. It is clear that approximately the first 70 eigenvalues are calculated with an accuracy better than 10 −12 . Eigenvalues with index up to 150 remain at the same accuracy as N ≈ 500. Above this index the accuracy decreases to 10 −3 . Actually, over years, a lot of literature has gathered on this problem. Low order eigenvalues have accurately computed using Runge-Kutta type methods for instance in [32] and [33]. In [34] the author uses Hermite collocation in order to find only the first eigenvalues for some Schrödinger eigenproblem. With respect to the departure from orthogonality, i.e. the distance from zero of the scalar product of two eigenvectors, we display in the right hand side of Fig. 10, in a log linear plot, the absolute values of the of the scalar products of u1 and uj , j = 2, . . . , 200. It is again remarkable that this departure is less than 10 −15 .
It is of some importance to justify our choice for SiC. Unlike all other methods of spectral collocation, where the differentiation matrices are highly non-normal, in SiC these matrices are symmetric or skewsymmetric for even respectively odd values of cut-off parameter N . And this is an important numerical advantage. For instance in the MATLAB code below we use the routine eigs instead of eig which would have been more expensive and slower.
The following very simple MATLAB code has been used. % sort eigenvalues This code is fairly similar with that from Section 4.2 with two differences. First, the Chebyshev differentiation matrices are replaced by sinc differentiation matrices and then the boundary conditions are absent. More exactly, the discrete sinc functions, on which the unknown solution is expanded, are defined by and satisfy the boundary conditions (10). In (11) the nodes x k are equidistant with spacing h and symmetric with respect to the origin.

Concluding remarks and open problems
Regarding the regular problems, Chebfun is unbeatable in terms of accuracy, computation speed, and the information they provide on the accuracy of computational process. It displays the optimal approximation order of unknowns (eigenvectors) and how and to what extent their Chebyshev coefficients decrease. It also specifies the degree to which some boundary conditions are satisfied. As for the singular problems, the situation is not so offering. In this case, Chebfun leaves enough room for the application of usual (classical) spectral methods. For singular eigenproblems on unbounded domains mapped ChC on the half line and SiC, for problems on the real line, perform better than Chebfun even for Chebfun applied on a mapped domain.
However, approaching these problems in parallel, with Chebfun as well as with the classical spectral methods, one can get greater confidence in the accuracy of numerical results. Regarding ChC method, using the socalled absolute or relative drift in terms of some parameters we ensure the numerical stability of the numerical process and can eliminate (numerically) spurious eigenvalues. In other words we get automatically and precisely the accuracy at which a specified set of eigenvalues is computed. The departure of orthogonality of a set of eigenvectors provides a useful hint for the accuracy of the numerical process. We have managed to correctly identify sets of multiple and high indices eigenvalues (triples) and even give some numerical meaning to the notion of continuous spectrum. This second issue obviously remains an open one.
All in all we can say that Chebfun as well as classical spectral methods in various forms, produce more accurate and comprehensive outcomes when solving singular and of course regular eigenproblems than classical methods. However, there remains an open problem, namely that of establishing (possibly automatically) the scaling factor.
The following abbreviations are used in this manuscript: