Accurate Spectral Collocation Computations of High Order Eigenvalues for Singular Schrödinger Equations-Revisited

In this paper, we continue to solve as accurately as possible singular eigenvalues problems attached to the Schrödinger equation. We use the conventional ChC and SiC as well as Chebfun. In order to quantify the accuracy of our outcomes, we use the drift with respect to some parameters, i.e., the order of approximation N, the length of integration interval X, or a small parameter ε, of a set of eigenvalues of interest. The deficiency of orthogonality of eigenvectors, which approximate eigenfunctions, is also an indication of the accuracy of the computations. The drift of eigenvalues provides an error estimation and, from that, one can achieve an error control. In both situations, conventional spectral collocation or Chebfun, the computing codes are simple and very efficient. An example for each such code is displayed so that it can be used. An extension to a 2D problem is also considered.


Introduction
In the last few years, there has been an increasing interest in investigating the spectrum of the non-selfadjoint Sturm-Liouville (SL) problem of the form − u + q(x)u = λu, −∞ ≤ a < x < b ≤ +∞, λ ∈ C, (1) where the potential q is both complex valued and bounded or real valued and unbounded. In our recent paper [1], we have analyzed the solution of six eigenvalue problems attached to the time independent Schrödinger equation above, from which five have been defined on unbounded domains, the real axis, or the real positive semi-axis. The exception was the Coffey-Evans potential. Some of them were more complicated in the sense that the potential q was additionally singular at the end located at the finite distance, i.e., a. Actually, there are many SL problems in chemistry, physics, astrophysics, engineering, biology, and other fields for which pseudospectral (collocation) methods can be used to an advantage. A lot of such examples are spread and analyzed in the recent monograph [2], where, in addition, a special section is dedicated. The variety of SL problems as well as the difficulties in solving them in closed or numerical form is extensively illustrated in the monograph [3].
Unfortunately, on the date of publication of the above quoted paper, we were not aware of the existence of the reference work [4] in which the authors examine all types of singularity of the real potential q. That is why we will examine here a series of other singular potentials, one of the most interesting being the so-called Weierstrass-Mandelbrot spectrum [5]. We will again rely in this paper on spectral collocation methods based on Chebyshev polynomials and sinc functions as well as on the new computing environment Chebfun.
Thus, for Chebyshev collocation (ChC), we refer among other sources to our contributions [1,6] as well as to the seminal paper [7]. The mapped variant of this method seems to have the origins in the paper [8]. However, we will rely on the explanations given in [1] (see also [9]). It is used for problems formulated on the half real line. For problems on the real axis, we consider that SiC is the most suitable method.
The Chebfun system, in object-oriented MATLAB, contains algorithms which amount to spectral collocation methods on Chebyshev grids of automatically determined resolution. In [10], the authors show that chebops are the fundamental Chebfun tools for solving (ordinary and partial) differential and integral equations, its properties are summarized in [11] and the applicability in solving various problems attached to ODE is exemplified in [12]. Moreover, in [10] (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).
Regardless of the method used, we will use the relative drift to estimate the error with which a specified set of eigenvalues is computed. The relative drift with respect to the order of approximation N or the length of integration interval X was reviewed in the first part of our paper [1]. Now, the drift is extended with respect to a small parameter ε which regularizes a singular problem. This indicator allows us to separate the "good" eigenvalues from the "bad" ones (spurious) and gives us a precise indication of the numerical stability of the used algorithm.
A general theory of the influence of singularities on the accuracy of the computation of pairs of eigenvalue-eigenfunction is far from being rigorously developed. In this context, the main purpose of this work is to analyze, as rigorously as possible, some typical singularities, internal to the field of integration, on its ends or to infinity for unbounded domains.
In other words, we want to offer a realistic test of performance of the Chebfun and some conventional spectral problems in solving such problems.
The paper is organized as follows: In Section 2, we consider the case of complex potential and, in Section 3, we solve some equations with unbounded potential at one endpoint. In the latter section, we devote separate subsections to the classical Dunford and Schwartz eigenproblem, to a problem with logarithmic potential and to a Schrödinger problem involving a quantum-mechanical potential. A brief comparison between ChC and Chebfun is conducted. A Schrödinger equation with an interior singularity is solved by SiC and Chebfun in the first part of Section 4 and then a special topic about the zeros of Bessel's functions. Section 5 is devoted to a regular 2D Schrödinger eigenproblem. The paper ends up with a small section (Section 6) where we summarize our work.
A MATLAB code implementing ChC and another Chebfun code are provided in order to facilitate future numerical experiments. All our numerical experiments have been carried out using MATLAB R2020a on an Intel (R) Xeon (R) CPU E5-1650 0 @ 3.20 GHz.

Complex Schrödinger Potential
For bounded and complex valued potential q and functions u ∈ L 2 [0, ∞), it is known from [13] that it is sufficient to impose only a boundary condition at 0 in order to make (1) into a well-posed eigenvalue problem. The natural choice is Furthermore, as needed, we shall assume appropriate smoothness conditions on the potential q, and we say that the endpoint 0 is a regular one and infinity a singular one for this eigenproblem. Let's consider now that the potential q(x) satisfies the following two conditions: is the Banach space of measurable functions which are bounded outside a set of measure zero; 2. the limit lim x→∞ q(x) =: q ∞ exists and is finite.
In these conditions, in [14], the authors introduce the closed linear operator and prove that: and where q r,min : Let's consider now the eigenproblem (1)-(2) equipped with the pure imaginary potential where α and R are real constants (wave number and Reynolds number, respectively) and We have solved the eigenproblem (1)-(2) by mapped ChC when α := 1, and V(x) := sin(x)exp(−x).
The eigenvalues for various Reynolds numbers are displayed in Figure 1. Their scattering in the complex plane is very similar to that which is only suggested in Figure 1 from [14]. However, a much larger number of eigenvalues is calculated and, more importantly, they all fully meet the conditions (3) and (4). The drift of eigenvalues with respect to order of approximation N is depicted in Figure 2a and the drift with respect to the scaling factor is reported in panel b of the same figure. Both show that no more than approximately 10 eigenvalues can be computed with an accuracy better than 10 −5 . The result in solving the same problem with Chebfun is displayed in Figure 3.

Singular Inverse Square Potential
Let's consider now Equation (1) to which we attach some separated boundary conditions in a and b and the potential q has the form where f and g are analytic functions inside and on a Bernstein ellipse containing the canonical segment [−1, 1], and γ ≥ 0. Problems of this type have applications in engineering and some of them have been solved in [15] by a spectral matrix method.
In Table 1, we compare our results with those reported in [15] (Table 3) and in turn compared with the results obtained by other methods. Unfortunately, the differences are visible. In order to ensure the accuracy of the results obtained by ChC, we report the relative drift with respect to the order of approximation in Figure 4. For the fifteenth eigenvalue, it guarantees an accuracy of the order of 10 −12 . The variation of the magnitude of this eigenvalue in relation to γ is the same, i.e., its magnitude increases with the value of this parameter. We mention that, for the last two values of γ, the point x = −1 is a singular one of nonoscillatory limit-circle type. Relative drift with respect to N The first four eigenvectors of problem (1)-(8) are depicted in Figure 5. It is very clear that they satisfy boundary conditions (8). Their computation is a surplus compared to the work [15] where they are not calculated. Their coefficients decrease smoothly and rapidly to the machine precision as it is apparent from Figure 6 panels (a) and (b). It is extremely important to point out that Figure 6c,d fully confirm the assertion in John Boyd's work [16]. In this paper, it is stated that: the "most useful measure of the rate of convergence of a Chebyshev series is r n := |a n /a n+1 |". Moreover, he states that "r n must decrease to one as n → ∞, if the function is singular at one or both end points ".
We mention at the end of this section that ChC together with a removing technique of independent boundary conditions solved without problems all the singular eigenvalues problems from [15].

The Eigenproblem of Dunford and Schwartz
The generalized hypergeometric equation reads In [17] Example 23, the author attaches to this equation the boundary conditions the first one being of Friedrichs type. We must mention that so far we have not found a suitable technique to introduce this type of boundary conditions when working with ChC. This is why ChC has not even tried for this problem. Anyway, in the proposed benchmark set of singular problems listed in [17], the problem (9)-(10) is considered one of the hardest.
The singularity in origin sinh 2 Thus, the boundary condition in origin comes from the equality [u, v] = 0 where the symbol [ , ] stands for the Wronskian of the functions u and v. The end point b := ∞ is of type limit point non-oscillatory, so we attach the homogeneous Dirichlet boundary condition. Actually, in order to resolve the singularity at ∞, we solve this problem on the truncated interval The best results (the most stable eigenvalues and the smoothest eigenvectors) have been obtained for ε := 1.e − 06 and X := 12. In order to improve the convergence of Chebfun, we use the option 'splitting','on'.
The first ten eigenvalues obtained are reported in Table 2. The first five of them reasonably approximate the exact eigenvalues λ n = −(5 − n) 2 , n = 0, 1, 2, 3, 4. The next five (the second column in Table 2) suggest the continuous spectrum on (0, ∞) in the sense that the eigenvalues are becoming more and more crowded as their index increases. As a precaution not to consume space in excess, we limited ourselves to report only a few eigenvalues from our numerical experiments. The drift with respect to the length X (of the truncated integration interval from (11)) of the discrete eigenvalues of (9)-(10) is depicted in Figure 7. It clearly shows an accuracy of order at least 10 −7 in computation of the first four eigenvalues. The first four eigenvectors of Dunford and Schwartz eigenproblem are displayed in Figure 8. They look plausible viz. the number of their roots.
With its two singularities, this problem has proved to be one of the most difficult and expensive in terms of computing elapsed time. However, Chebfun has provided remarkable numerical results.

A Nasty Eigenproblem with Logarithmic Potential
We will consider in this subsection the eigenproblem (1) supplied with homogeneous Dirichlet boundary conditions and the potential defined by q(x) := ln (x).
The problem is solved in [17] by the SLTSTPAK package on the interval [0, 4], as Example 11, and is considered as a regular problem that looks singular.

ChC Method
The problem contains that nasty q(x), but the classical ChC works extremely fast and accurate. The following simple MATLAB code produces the eigenpairs of the matrix A: We enforce the boundary condition by the removing technique of independent boundary conditions introduced in [18]. Thus, we delete the first and the last rows and columns in the second order Chebyshev differentiation matrix D(:, :, 2) and obtain the matrix Deig.
The nodes x 1 and x N of the Chebyshev-Gauss-Lobatto system are eliminated as slaves and after the problem is solved they are given-back. Actually, we simply add two zeros, one in the first and another in the last position of each eigenvector.

Chebfun Method
In order to avoid the singularity of q(x) in origin, we solve the same problem by Chebfun on the modified (truncated) domain of the form (11), namely [1.e − 08, 4]. The first four eigenvectors are displayed in Figure 9. They look fairly smooth and they do not seem affected at all by the singularity. However, our numerical experiments showed that their coefficients do not decrease smoothly, and this is due to the singularity of q(x). Anyway, no warning has been displayed by Chebfun when the option 'splitting', 'on' has been invoked.

ChC vs. Chebfun
The first and the twenty-fourth eigenvalues of this eigenproblem computed by ChC, Chebfun, and reported in [17] are displayed in Table 3. If, with respect to the first eigenvalue, we observe that they coincide up to the fifth decimal, this does not happen for the value of the twenty-fourth eigenvalue. Thus, in order to separate the "good" eigenvalues from the "bad" ones, i.e., numerically inaccurate eigenvalues, we use again the concept of relative drift, now with respect to the order of approximation N. Actually, the relative drift of eigenvalues computed by ChC and reported in Figure 10 shows a coincidence of at least order 10 −11 for the first twenty-fourth eigenvalue computed by ChC. In other words, it means that the relative drift is resolution independent with respect to N for the outcomes of ChC method, and the method is numerically stable. Inspecting now the Table 3, we can conclude that the eigenvalues computed by ChC and Chebfun are more accurate than those reported in [17], especially for large values of N.   The eigenvectors computed by the two methods look practically identical, so we present in Figure 9 only those obtained through Chebfun. Now, it is clear that the price paid by Chebfun to provide eigenvalues with the same order of accuracy as ChC is the truncation of the domain.

Quantum-Mechanical Potential with Weierstrass Spectrum
We are concerned now in Equation (1) equipped with the Weierstrass-Mandelbrot potential along with boundary conditions In order to avoid singularity in origin, in [5], it is suggested to use the modified potential q(x) := ω 2 − (1/4) /(x + ε) 2 , ω > 0 and ε → 0.
For various values of ω, the outcomes are displayed in Table 4. We think some comments are important at this time. The increase of the number of negative eigenvalues with the increase of parameter ω is in qualitative agreement with the conclusions from the paper [5] Section 5. Unfortunately, we have not been able to prove that the calculated spectrum is indeed a Weierstrass one. This remains an open issue. This problem is listed in the paper [17] but is not actually solved so we do not have a term of comparison. However, our numerical experiments showed that Chebfun is stable with respect to the length X in the range [5,12].

Schrödinger Type Equations with an Interior Singularity
Here, we analyze first the so called Boyd's eigenproblem. It is called singular if q ∈ L[α, β] for any [α, β] ⊂ [a, b] and q / ∈ L[a, b]. For two real numbers a and b, such that −∞ < a < 0 < b < ∞, the potential q(x) is defined by We attach the homogeneous Dirichlet boundary conditions to Equation (1) with q defined in (16).
The corresponding real-valued eigenfunctions {u n , n ∈ N 0 } satisfy the Equation (1) and the Dirichlet boundary conditions (17) and the following properties for n ∈ N 0 : 1.
{u n , n ∈ N 0 } is a complete orthonormal set in the Hilbert space L 2 [a, b]; 7.
degeneracy may occur, but no eigenvalue can have multiplicity greater than two.

Chebfun Results
The Chebfun code has failed to solve this problem in form (1), (16) and (17). However, in order to solve it numerically, we use a regularized version of q(x), namely in the sense of distributions, i.e., lim ε→0 q(x, ε) = −δ 1 x , where, as usual, δ stands for Dirac distribution. With this approximation of the original coefficient q(x), the eigenproblem becomes a regular one. It is temping to hope that the perturbed set of eigenvalues {λ n (ε), n ∈ N 0 } and of eigenvectors {u n (x, ε), n ∈ N 0 } approaches in some sense the corresponding sets of the original problem when ε → 0.
Actually, with Chebfun computations, we have used ε := 1.e − 05 in (18) on the integration interval [− 10,10]. The length of this interval is arbitrary, so, in order to back the results obtained by Chebfun, we use an alternative method in Section 4.2.

Sinc Collocation Results
With the SiC, we have used the order of approximation N in the range 200-500 and the scaling factor h := 0.1. To the parameter, ε in (18) has been assigned the same value, namely 1e − 05.
The first four vectors obtained are displayed in Figure 11. They are continuous in origin but with discontinuous derivatives. They are also mutually orthogonal as it clearly results from the left panel (a) Figure 12.  Thus, the first five items enumerated above are plainly verified. Moreover the spectrum is real, discrete, and converges to ∞ as N → ∞.
In Table 5, we report the first six eigenvalues computed by SiC, Chebfun, and by the SLEIGN package according to [20]. In the latter case, the authors worked with ε := 1.e − 03 on the interval [−10, 10]. However, in the right panel b) of Figure 12, the relative drift of the first eigenvalue is around O(10 −10 ), which means an error of the same order. Unfortunately, the accuracy decreases rather rapidly and thus we can trust only the first five, six eigenvalues. In order to justify the convergence of the eigenvalues calculated with SiC when ε tends to zero, we have estimated the drift of the first 100 eigenvalues with respect to this parameter. This is reported in Figure 13, and it gives us confidence that this convergence does take place. What is interesting and still unexplained is the fact that the largest oscillations appear in the calculation of the first eigenvalues after which the algorithm stabilizes.

Zero-Pairs of Bessel Functions
In [21], the author considers a very interesting generalization to the Bessel's problem.
He tries to answer the question of finding zero-pairs of Bessel functions of order ν, i.e., some numbers a and b for which there exists a Bessel function of order ν which vanishes simultaneously at a and b. We will consider this question from a numerical point of view in searching a deeper understanding of the singularity in origin of Bessel's eigenproblem.
This new problem reads For given τ > 1 (or τ < −1), we have a regular Sturm-Liouville problem with eigenvalues Chebfun confirms this statement and solves this regular problem in no time. It is easy to show that the functions λ n (τ) are analytic for the above values of τ. The author of [21] asks for properties of the analytic continuation of λ n (τ) into the complex plane. For ν ∈ [1/3, 1/2] in the same paper [21], it is shown that λ n (τ), n ∈ N are analytic in the domain C \ [−1, 1] and at τ := ∞, and these functions are also analytic on the segments τ ± i0 for −1 < τ < 1. They can be extended continuously into τ = ±1, respectively.
With the following simple Chebfun code, we have computed some eigenpairs of problem (19) for various τ and ν. They are reported in the first two columns of Table 6. Table 6. The first eight eigenvalues of eigenproblem (19) computed by Chebfun for τ := 0 and τ := 1/4 (first two columns).
In the third column, we report the eigenvalues of the same problem computed by ChC when τ := 0. In both methods, we have assumed ν = 1/3.  It is very important to note that, in solving this eigenproblem both methods, ChC as well as Chebfun, have encountered some difficulties. Thus, for the values of τ close to zero Chebfun displays the following warning: Maximum dimension reached. Solution may not have converged. This warning prompted us to look for a way to validate these results. Consequently, we have solved the problem with the classic ChC. For an unusually large approximation order N := 2028, i.e., half of that used by Chebfun, ChC produced the third column of Table 6.
The drift of the first ten eigenvalues is depicted in Figure 14. This figure suggests a sort of lack of numerical stability and also that we can not hope for a better accuracy than 10 −3 in approximating these pairs of quasi double eigenvalues.   The dependence of the first two sets of eigenvalues λ 0 and λ 1 of problem (19) on τ is depicted in Figure 16. The lower curve corresponds to λ 0 (τ), τ ∈ [0, 1]. This figure has to be seen with a symmetric part on [−1, 0]. In order to smooth these curves we have used the MATLAB pchip code that performs a piecewise Hermite cubic interpolation. In our opinion, this last figure suggests that, in origin, the so-called phenomenon of level crossing of eigenvalues occurs as described in [22] (p. 350) (see also [23]). Due to the scale at which the curves are drawn, the first two eigenvalues λ 0 and λ 1 , from the first column of the Table 6, appear to be superimposed, but actually they are only very nearly, i.e., almost double.
Along with the results displayed in Table 6, we can answer the question formulated in the title of paper [24]. Thusfar, close to the origin multiple eigenvalues are not rare at all for the Bessel problem (19). Moreover, this result is in the spirit of Theorem 6.6 from [24] on degeneracies in the spectra of self-adjoint operators of Schrödinger type.
Other several numerical results are in progress with respect to this challenging problem.

Extension to a Regular 2D Schrödinger Equation
In order to fix ideas, let's consider a 2D Schrödinger eigenproblem. It reads along with homogeneous Dirichlet boundary conditions on the sides of the standard square With respect to the problem (20)- (21), also called the generalized Helmholtz problem, Rayleigh, on physical grounds, and then Hilbert, Weyl, and others proved the so-called Ohm-Rayleigh principle. It states that, on any bounded domain Ω with compact closure, the problem (20)-(21) has a basis of eigenfunctions, and any two eigenfunctions having distinct eigenvalues are orthogonal (see [25]). Moreover, the Laplace operator is separable and in operator notation we can write where A := ∂ 2 /∂x 2 , B := ∂ 2 /∂y 2 and ⊗ stands for the direct or tensor product.
We have solved this problem by ChC. The second order partial derivatives have been computed using the MATLAB function kron as it is suggested in [7]. The main disadvantage of this approach for solving 2D problems is that the order of the involved matrices are large, i.e., it equals N 2 × N 2 (the so-called the curse of dimensionality).
Fortunately, our results plainly confirm these outcomes as it is clear from Table 7. Actually, in order to carry out this comparison, we have worked in the shifted domain [−5.5, 5, 5] × [−5.5, 5.5]. We notice that the discretization involved in [26] leads to matrices of the same large order N 2 × N 2 . We observe that eigenvalues λ 2 , λ 6 and λ 7 are double degenerate. In order to assess the accuracy in computing the eigenenergies, we display in Figure 17 the drift with respect to N of the first 25 eigenvalues.Thus, we can hope for a precision better than 10 −6 for all these eigenvalues. Some eigenvectors are reported in Figure 18. In order to save space and to get clarity, we have used the option off for axis and have interpolated on a finer grid in order to plot.
In fact, we easily calculated the first one hundred eigenpairs using the MATLAB eig function. Figure 18. From upper left to right down, we display the first and the eigenvectors corresponding, respectively, to eigenvalues 2 a , 4, 6 a , 7 a , and 12.

Conclusions and Some Open Problems
In this paper, we have added a series of new singularities to the Schrödinger equation in order to test the qualities of spectral methods in resolving them. A first finding is that both ChC and Chebfun operate accurately in the case of complex potentials. However, as a general remark, we must underline once again the high qualities of both classes of methods used and note that ChC has been successfully extended to a 2D regular problem.
Many other singular SL eigenproblems are waiting to be analyzed, using spectral methods, in order to better understand their possibilities and limitations of application.
A particular open problem that first comes to mind, and for which there are some results in progress, is the simultaneous dependence of the eigenvalues of Bessel's equation on both parameters τ and ν. This could shed more light on the singularity in origin for various problems attached to this classical equation. a second would be a more detailed analysis of the Weierstrass spectrum.
Funding: This research received no external funding.

Acknowledgments:
We would like to thank the editor for the opportunity to publish this paper in the journal Symmetry. We also would like to express our gratitude to both referees for their valuable comments which considerably improved the original form of this paper.

Conflicts of Interest:
The author declares no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: 2D two