Abstract
This work is about the use of some classical spectral collocation methods as well as with the new software system Chebfun in order to compute the eigenpairs of some high order Sturm–Liouville eigenproblems. The analysis is divided into two distinct directions. For problems with clamped boundary conditions, we use the preconditioning of the spectral collocation differentiation matrices and for hinged end boundary conditions the equation is transformed into a second order system and then the conventional ChC is applied. A challenging set of “hard” benchmark problems, for which usual numerical methods (FD, FE, shooting, etc.) encounter difficulties or even fail, are analyzed in order to evaluate the qualities and drawbacks of spectral methods. In order to separate ‘‘good” and “bad” (spurious) eigenvalues, we estimate the drift of the set of eigenvalues of interest with respect to the order of approximation N. This drift gives us a very precise indication of the accuracy with which the eigenvalues are computed, i.e., an automatic estimation and error control of the eigenvalue error. Two MATLAB codes models for spectral collocation (ChC and SiC) and another for Chebfun are provided. They outperform the old codes used so far and can be easily modified to solve other problems.
Authors
CalinIoan Gheorghiu
“Tiberiu Popoviciu” Institute of Numerical Analysis, Romanian Academy
Keywords
Sturm–Liouville; clamped; hinged boundary condition; spectral collocation; Chebfun; chebop; eigenpairs; preconditioning; drift; error control
References
see the expanding block below
Cite this paper as:
C.I. Gheorghiu, Accurate spectral collocation solutions to 2ndorder Sturm–Liouville problems, Symmetry, 13 (2021) 3, 385, doi: 10.3390/sym13030385
About this paper
Print ISSN
Not available yet.
Online ISSN
20738994
Google Scholar Profile
soon
Paper (preprint) in HTML form
Accurate Spectral Collocation Solutions to 2nthOrder SturmLiouville Problems.
Abstract
This work is about the use of some classical spectral collocation methods as well as about the use of the new software system Chebfun in order to compute the eigenpairs of some high order SturmLiouville eigenproblems. The analysis is divided into two distinct directions. For problems with clamped boundary conditions we use the preconditioning of the spectral collocation differentiation matrices and for hinged end boundary conditions the equation is transformed into a second order system and then the conventional ChC is applied. A challenging set of “hard”benchmark problems, for which usual numerical methods (FD, FE., shooting etc.) encounter difficulties or even fail, are analysed in order to evaluate the qualities and drawbacks of spectral methods. In order to separate “good”and “bad”(spurious) eigenvalues we estimate the drift of the set of eigenvalues of interest with respect to the order of approximation $N$. This drift gives us a very precise indication of the accuracy with which the eigenvalues are computed, i.e., an automatic estimation and error control of the eigenvalue error. Two MATLAB codes model for spectral collocation (ChC and SiC) and another for Chebfun are provided. They outperform the old codes used so far and can be easily modified to solve other problems.
1 Introduction
Due to the spectacular evolution of advanced programming environments, a special curiosity arose in the numerical analysis of a classical problem, that of accurate solving of high order SL eigenproblems. It seems that quantum mechanics is the richest source of selfadjoint problems, while nonselfadjoint problems arise in hydrodynamic and magnetohydrodynamic stability theory (see for instance S and the vast literature quoted there). The need to compute accurately and efficiently a large set of eigenvalues and eigenfunctions, including those of high index, is now an utmost task.
Our main interest here is to evaluate the capabilities of the new Chebfun package as well as those of conventional spectral methods in meeting these requirements. The latter work in the classical mode, i.e., “discretizethensolve”. On the contrary, the Chebfun spirit consists in the continuous mode, i.e., “solvethendiscretize”(see TBD p.302).
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.
Some FORTRAN software packages have been designed over time to solve various regular and singular SL problems. These seem to be the first attempts to solve numerically (automatically) eigenvalue problems.The most important would be SLEDGE PF ; PFX , the NAG’s code SL02F MP ; PM , SLEIGN and SLEIGN2 BEZ ; BGKZ , and later MATSLISE. The SLDRIVER interactive package supports exploration of a set of SL problems with the first four previously mentioned packages. The SLEDGE, SL02F, SLEIGN2, and NAG’s D02KDF are “automatic”for eigenvalues and not for eigenfunctions. They have built in error estimation and from that they achieve error control. They adjust the accuracy of the discretization so that the delivered eigenvalue has estimated error below a usersupplied tolerance.
Essentially, the numerical method used in this software packages replaces the coefficients in equation by step function approximation. Their most important drawback remains the impossibility to compute the eigenfunctions and a slow convergence in case of some singular eigenproblems.
The MATSLISE code introduced in Ledoux can solve some Schrödinger eigenvalue problem by a constant perturbation method of a higher order. Very recently this code has been improved (see BVD ) but it remains to Schrödinger issues which are outside the scope of this paper.
There is also a class of semianalytical methods which includes the variational iteration method, the homotopy perturbation method, homotopy analysis method and Adomian decomposition (see for instance AS ) for solving eigenvalue problems. Their accuracy is far from what spectral collocation methods can provide.
In PB the authors set up an ambitious method based on Lie group method along with the Magnus expansion in order to solve any order of SL problem with arbitrary boundary conditions.
We believe that spectral collocation methods can contribute to the systematic clarification of some still open issues related to the numeric aspects of SL problems. The most important aspect is how many computed eigenpairs (eigenvalues and eigenfunctions) can we trust when solving a high order SL? This is the outstanding, not completely resolved research issue, we want to address in this paper.
Thus, we will argue that generally Chebfun would provide a greater flexibility in solving various differential problems than the classical spectral methods. This fact is fully true for regular problems. A Chebfun code contains a few lines in which the differential operator is defined along with the boundary conditions and then a subroutine to solve the algebraic eigenproblem. It provides useful information on the optimal order of approximation of eigenvectors and the degree to which the boundary conditions have been satisfied.
Unfortunately, in the presence of various singularities or for problems of higher order than $4$, the maximum order of approximation of the unknowns can be reached ($N\ge 4000$) and then Chebfun issues a message that warns about the possible inaccuracy of the results provided.
Alternative use of conventional spectral collocation methods generally helps to overcome this difficulty.
As a matter of fact, in order to resolve a singularity on one end of the integration interval, Chebfun uses only the truncation of the domain. Classical spectral methods can also use this method, but it is not recommended because much more sophisticated methods are at hand in this case. For singular points at finite distances (mainly origin) we will use the socalled removing technique of independent boundary conditions (see for a review of this technique our monograph Cig14 p. 91). The boundary conditions at infinity can be enforced using basis functions that asymptotically satisfy these conditions (Laguerre, Hermite, sinc).
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. These codes are very easy to implement, efficient, reliable and robust. All our numerical experiments have been carried out using MATLAB R2020a on an Intel (R) Xeon (R) CPU E51650 0 @ 3.20 GHz.
The main purpose of this paper is to argue that Chebfun, along with the spectral collocation methods, can be a very feasible alternative to the above software packages regarding accuracy, robustness as well as simplicity of implementation. In addition, these methods can calculate exactly the “whole”set of eigenvectors approximating eigenfunctions and provide automatic estimation and control of the eigenvalue error. For selfadjoint problems, checking the orthonormality of computed eigenvectors gives us valuable information on the accuracy of the calculation of these vectors.
The structure of this work is as follows. In Section 2 we recall some specific issues for the regular as well as singular SturmLiouville eigenproblems. In Section 3 we review briefly the conventional ChC method as well as Chebfun, the relative drift of a set of eigenvalues and the preconditioning of Chebyshev differentiation matrices. The Section 4 is the core part of our study. By analyzing one set of hinged problems and another one of clamped problems, we want to evaluate the applicability of the two classes of methods as well as their performances in terms of the accuracy of the outcomes they produce. There is also a subsection that contains problems equipped with boundary conditions that are a mixture of these two types, clamped and hinged. We end up with Section 5 devoted to conclusions and open problems.
2 2nthorder SturmLiouville eigenproblems
The 2nthorder SL equation reads

along with separated, (selfadjoint) boundary conditions. We shall assume that all coefficient functions are real valued. The technical conditions for the problem to be non singular are: the interval $(a;b)$ is finite; the coefficient functions ${p}_{k}$, $$, the weight $w$ and $1/{p}_{n}$ are in ${L}^{1}(a,b)$; and the essential infima of ${p}_{n}$ and $w$ are both positive. Under these assumptions, the eigenvalues are bounded below (see for instance GM ).
The eigenvalues can be ordered in the usual form: ${\lambda}_{0}\le {\lambda}_{1}\le {\lambda}_{2}\le \mathrm{\dots},$ such that ${lim}_{k\to \mathrm{\infty}}{\lambda}_{k}=+\mathrm{\infty}.$ In this sequence each eigenvalue has multiplicity at most $n$ (so $k+n>k$ for all $k$). The restriction on the multiplicity arises from the fact that for each $\lambda $ there are at most $n$ linearly independent solutions of the differential equation satisfying either of the endpoint conditions which we shall consider below.
Some of the problems we deal with are also found in the monographic paper Everitt . It contains over 50 challenging examples from mathematical physics and applied mathematics along with a summary of SL theory, differential operators, Hilbert function spaces, classification of interval endpoints and boundary condition functions.
3 Chebfun vs. conventional spectral collocation
3.1 Chebfun
For details on Chebfun we refer to TBD and DBT ; DHT ; DHT1 . The Chebfun system, in objectoriented MATLAB, contains algorithms which amount to spectral collocation methods on Chebyshev grids of automatically determined resolution. This is the main difference compared to conventional spectral methods in which the resolution (order of approximation) is imposed almost arbitrarily. Its properties are briefly summarized in DHT . In DBT the authors explain that chebops are the fundamental Chebfun tools for solving ordinary, partial differential or integral equations.
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 papers as well as in the text TBD . Thus one can get a suggestive image of what they can do working with Chebfun.
Moreover, in DBT 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 nonnormality 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 high order SL problems.
3.2 Spectral collocation methods
Spectral methods have been shown to provide exponential convergence for a large variety of problems, generally with smooth solutions, and are often preferred Cig18 . In all spectral collocation methods designed so far we have used the collocation differentiation matrices from the seminal paper WR . 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 the boundary bordering, which is a simplified variant of the above mentioned removing technique of independent boundary conditions, as well as the basis recombination. We have used the first technique in the large majority of our papers except CigIsp where the latter technique has been employed. In the last quoted paper a modified ChT method based on basis recombination has been used in order to solve an OrrSommerfeld problem with an eigenparameter dependent boundary condition.
Once eigenvectors are calculated in physical space they are transposed into the space of coefficients using FCT. In this way it is possible to estimate the way in which their coefficients decrease.
3.3 The drift of eigenvalues
Two techniques are used in order to eliminate the “bad”eigenvalues as well as to estimate the stability (accuracy) of ChC or Chebfun 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. 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 (cutoff parameter) $N$ or the scaling factor (length of integration interval). Only those whose difference or “resolutiondependent drift”is “small”can be believed. In this connection, in the paper JPB32 the so called absolute (ordinal) drift with respect to the order of approximation has been introduced. We extend this definition in our recent paper Cig21 and will use it without repeating it here.
Whenever the exact eigenvalues of a problem are known, the relative drift is reduced to the relative error.
At this point, the following observation is extremely important. In the highly cited monograph JPB , the author makes a subtle analysis of spectral methods in solving linear eigenproblems. Among others, he states the so called Boyd’s Eigenvalues Rule ofThumb in which he notices that in solving such a problem with a spectral method using $(N+1)$ terms in the truncated spectral series, the lowest $N/2$ eigenvalues are usually accurate to within a few percent, while the larger $N/2$ numerical eigenvalues differ from those of the differential equation by such large amounts as to be useless.
3.4 Preconditioning
To simplify the introduction of a preconditioner we use the differential operator
$${L}^{\left(n\right)}\left(u\right):=\frac{{d}^{n}u}{d{x}^{n}},x\in (1,1),$$  (1) 
subject to clamped boundary conditions
$${u}^{\left(\mu \right)}\left(1\right)=0,0\le \mu \le {l}_{n},$$  (2) 
and
$${u}^{\left(\nu \right)}\left(1\right)=0,0\le \nu \le {r}_{n},$$  (3) 
where $n>1,$ ${l}_{n}$ and ${r}_{n}$ are positive integers such that ${r}_{n}+{l}_{n}+2=n.$
It is well known that for general collocation points the first order differentiation matrix has a condition number of order ${N}^{2}$ and the second order differentiation matrix has a condition number of order ${N}^{4}$ as $N\to \mathrm{\infty}.$ We comment on the preconditioner introduced in HS . These authors show that the preconditioning matrix
$$\mathbf{D}:=diag\left({\left(1+{x}_{k}\right)}^{{l}_{k}+1}{\left(1{x}_{k}\right)}^{{r}_{k}+1}\right),2\le k\le N1,$$ 
applied to ChC as well as to Chebfun discretization ${\mathbf{L}}_{ChC\text{or}Chebfun}^{(n)}$ of differential operator ${L}^{(n)},$ produces matrices ${\mathrm{\mathbf{D}\mathbf{L}}}_{ChC\text{or}Chebfun}^{(n)}$ of an inferior condition number, namely ${N}^{n}.$
4 Numerical experiments
4.1 Hinged ends or simply supported boundary conditions
4.1.1 The Viola’s eigenproblemrevisited
Let us consider now the so called Viola’s eigenproblem. It is encountered in porous stability problems (see S , Ch. 9) and reads
$$  (4) 
It is singular as $\theta \to {1}^{}$ in accordance with the definition introduced in Sect. 2.
By straightforward variational arguments we have shown in Cig14 p. 50, that lowest eigenvalue is positive. In this text we have solved the above problem by ChC using the so called ${D}^{2}$ strategy which involves the change of variables
$$v:={\left(1\theta x\right)}^{3}\frac{{d}^{2}u}{d{x}^{2}}.$$ 
The main deficiency of this strategy is the fact that it produces a lot of numerical spurious eigenvalues (at infinity).
In spite of this, we succeeded in stating the conjecture according to which ${\lambda}_{1}\left(\theta \right),$ the lowest eigenvalue of the problem (4), approaches $1$ as $\theta \to {1}^{}.$
Now, taking advantage of Chebfun we solve directly problem (4). Thus, the dependence of the lowest eigenvalue of the problem (4), computed by Chebfun, on the parameter $\theta $ is depicted in Fig. 1. Actually we have obtained
$${\lambda}_{1}\left(0.98765\right)=8.775218471808549e01,$$ 
which only partly confirms the above conjecture.
As a validation issue for our computations we have obtained the known value
$${\lambda}_{1}\left(0.0\right)={\pi}^{4},$$ 
within an approximation of a thousand.
For the highest computable value of parameter $\theta $ we display in Fig. 2 the first four eigenvectors of problem (4) and in Fig. 3 the Chebyshev coefficients of these eigenvectors.
We have to mention that the singularity in the right end $x=1$ becomes more prominent as the $\theta $ tends to $1$. This is confirmed by increasing the degree of the Chebfun approximation. For instance when $\theta :=0$ only a $25$ degree Chebyshev polynomial it uses and when $\theta :=0.98765$ the degree of approximation grows to more than $80$ (see Fig. 3). It is also worth mentioning that only for $\theta $ growing very close to $1$ a truncation of the domain along with the use of the option splitting have been necessary when Chebfun has been used.
4.1.2 The Bénard stability problem
A simplified form of the Bénard stability problem supplied with selfadjoint boundary conditions reads (see for instance GM and LA )
$$\begin{array}{c}{u}^{\left(vi\right)}\left(2\nu +3\right){u}^{\left(iv\right)}+\left({\nu}^{2}+4\nu +3\right){u}^{\prime \prime}\left({\nu}^{2}+2\nu +2\right)u=\lambda u\left(x\right),x\in (0,1),\\ u\left(0\right)=u\left(1\right)={u}^{\prime \prime}\left(0\right)={u}^{\prime \prime}\left(1\right)={u}^{\left(iv\right)}\left(0\right)={u}^{\left(iv\right)}\left(1\right)=0.\end{array}$$  (5) 
The constant $\nu $ is regarded as a parameter which typically can take the values
$${\nu}_{j}^{\pm}:=\left(1+{j}^{2}{\pi}^{2}\right)\pm {\left(1+{j}^{2}{\pi}^{2}\right)}^{1/2},j=1,2,\mathrm{\dots}.$$ 
All our attempts to solve this problem using Chebfun have failed, so we have resorted to the old ${D}^{2}$ strategy.
$\nu $  ${\lambda}_{0}\left(\nu \right)$  ${\lambda}_{1}\left(\nu \right)$  ${\lambda}_{0}\left(\nu \right)$ according to GM 

$(1+{\pi}^{2})$  $1.000000000102923e+00$  $3.548769279033568e+04$  $1.000005$ 
$(1+4{\pi}^{2})$  $1.000000009534114e+00$  $9.530184561696226e+03$  $1.0001$ 
$(1+{\pi}^{2}){(1+{\pi}^{2})}^{1/2}$  $1.191482998363510e+02$  $2.802486989002433e+04$  $1\times {10}^{7}$ 
$(1+4{\pi}^{2}){(1+4{\pi}^{2})}^{1/2}$  $1.639502291744172e+03$  $1.406538196754713e+04$  $3\times {10}^{5}$ 
Thus we rewrite problem (5) as a homogeneous Dirichlet one attached to a second order differential system, namely
$$\begin{array}{c}{u}^{\prime \prime}=v\left(x\right),x\in (0,1),\\ {v}^{\prime \prime}=w\left(x\right),\\ {w}^{\prime \prime}\left(2\nu +3\right)w+\left({\nu}^{2}+4\nu +3\right)v\left({\nu}^{2}+2\nu +2\right)u=\lambda u\left(x\right),\\ u=v=w=0\text{in}x=0\text{and}x=1.\end{array}$$  (6) 
Now we apply to each line the ChC discretization. It leads to the generalized and singular eigenpencil
$$(\mathbf{A},\mathbf{B}),$$  (7) 
where the block matrices are defined by
$$\mathbf{A}:=\left(\begin{array}{ccc}4{\stackrel{~}{D}}^{\left(2\right)}& I& Z\\ Z& 4{\stackrel{~}{D}}^{\left(2\right)}& I\\ \left({\nu}^{2}+2\nu +2\right)I& \left({\nu}^{2}+4\nu +3\right)I& 4{\stackrel{~}{D}}^{\left(2\right)}\left(2\nu +3\right)I\end{array}\right),$$ 
and
$$\mathbf{B}:=\left(\begin{array}{ccc}Z& Z& Z\\ Z& Z& Z\\ I& Z& Z\end{array}\right).$$ 
The factor $4$ in front of ${\stackrel{~}{D}}^{\left(2\right)}$ comes from the shift of interval $(0,1)$to the canonical Chebyshev interval $[1,1]$ and the matrix ${\stackrel{~}{D}}^{\left(2\right)}$ signifies the second order Chebyshev differentiation matrix with the homogeneous Dirichlet boundary conditions enforced. The matrices $I$ and $Z$ stand respectively for the identity and zeros matrices of the same dimension as ${\stackrel{~}{D}}^{\left(2\right)}.$
The following short MATLAB code has been used to solve (6):
N=256; % order of approximation nu=(1+4*(pi^2))sqrt(1+4*(pi^2)); % parameter \nu [x,D]=chebdif(N,2); D2=D(2:N1,2:N1,2); % differentiation matrices I=eye(size(D2)); Z=zeros(size(D2)); A=[ 4*D2 I Z; Z 4*D2 I; (nu^2+2*nu+2)*I (nu^2+4*nu+3)*I 4*D2(2*nu+3)*I]; B=[Z Z Z; Z Z Z; I Z Z]; % block matrices in pencil k = 8 ; % number of computed eigs E=eigs( @(x)(A\(B*x)), size(A,1),k, ’SM’) % Arnoldi method
When the order of approximation is $N$ both matrices in (7) have order $3\times (N1).$ This tripling of the dimensions of the matrices involved is not a major disadvantage. On the contrary, if we use Henrici’s number as a measure of normality (see for instance our text Cig14 pp. 2223) we see from the inequality
$$ 
that matrix $\mathbf{A}$ is more normal than ${\stackrel{~}{D}}^{\left(2\right)}$.
In our previous paper CigJR we have analyzed various methods to solve singular eigenproblems attached to pencils of the form (7). For the problem at hand we have used the Arnoldi method with the MATLAB sequence eigs(${A}^{1}B$) and with the above code obtained the eigenvalue reported in Table 1. It is very clear that for the values of $\nu $ considered, the block matrix $\mathbf{A}$ is non singular and the block matrix $\mathbf{B}$ is singular and independent of $\nu .$
It is extremely important to point out that for the first two values of the parameter $\nu $ in Table 1 our results are very close to those reported in GM . For the other parameter values this no longer happens. A similar situation occurs even for the second eigenvalue. Then, to decide, over the accuracy of our outcomes we resorted to drift. The relative drift, with respect to the order of approximation $N,$ of the first forty eigenvalues when $\nu :=(1+{\pi}^{2}){(1+{\pi}^{2})}^{1/2}$ is displayed in Fig. 4. It suggests that the first two eigenvalues are computed with an accuracy of at least ${10}^{12}$ and the first forty with an accuracy of at least ${10}^{2}$. This leads us to believe that we have produced much better approximations for these eigenvalues than those reported in GM as well as LA . Actually, in GM the authors use the SLEUTH code and accept that “it is clear that the code is not very accurate on this problem.”
4.1.3 A selfadjoint eighthorder problem
The eigenproblem of the highest order we consider in this paper is the following

(8) 
with exact eigenvalues ${\lambda}_{k}={\left(k\pi \right)}^{8},k=1,2,\mathrm{\dots}.$
All our attempts to solve this problem with Chebfun have failed. The abortion message referred to the extremely small conditioning of the eight order Chebyshev collocation differentiation matrix (of the order ${10}^{40}$).
Instead, the ${D}^{2}$ strategy along with ChC worked well and produced vectors from Fig. 5.
In order to estimate the error with which the eigenvalues were calculated, we display in Fig. 6 b) the relative drift of the first twelve eigenvalues for different approximation orders. Because we know even the exact eigenvalues, we also display the relative errors. It is very clear that the first eigenvalue is computed with better accuracy than ${10}^{12}$. Unfortunately, this means a lower performance by three decimals than that of Magnus expansion reported in Table 10 from PB .
Moreover, this means that we cannot trust more than twelve eigenvalues for this problem.
It is clear that ChC along with ${D}^{2}$ method has the potential to find the first eigenvalues of a SL problem of arbitrary (even) order with good accuracy. In addition to the Magnus method, this strategy calculates its eigenfunctions (eigenvectors) with reasonable accuracy as can be seen in Fig. 6 a). We use FCT to compute the Chebyshev coefficients of the eigenvectors.
4.2 Clamped boundary conditions
4.2.1 A fourth order problem with a third derivative term
With this first example we have to highlight the importance of preconditioning in improving the accuracy of Chebfun. In the papers HS and GTD , as well as in the monograph Fornberg , the following eigenproblem is carefully studied. It consists of the fourth order differential equation
$${u}^{\left(iv\right)}+R{u}^{\prime \prime \prime}=s{u}^{\prime \prime},x\in (1,1),R\in \mathbb{R},$$  (9) 
supplied with the clamped boundary conditions
$$u\left(\pm 1\right)={u}^{\prime}\left(\pm 1\right)=0.$$  (10) 
The eigencondition for this problem is
$${({R}^{2}+4s)}^{1/2}[\mathrm{cosh}\left(R\right)\mathrm{cosh}{({R}^{2}+4s)}^{1/2}]+2s\mathrm{sinh}{({R}^{2}+4s)}^{1/2}=0.$$  (11) 
Problems similar to this appear , for example, in linearized stability analysis in fluid dynamics. In GTD the authors noticed spurious eigenvalues when the problem (9)(10) is solved by ChT method. These spurious eigenvalues appears in the righthalf plane suggesting physical instabilities that do not exist.
$i$  ${\lambda}_{i}$ Chebfun  ${\lambda}_{i}$ Preconditioned Chebfun  ${\lambda}_{i}$ Solution to (11) 

$1$  $\mathbf{9.8}70154876048822e+00$  $\mathbf{9.869604}528925013e+00$  $9.8696044,$ 
$2$  $\mathbf{2.019}216607051227e+01$  $\mathbf{2.0190728}37497370e+01$  $20.1907286$ 
The numerical outcomes are displayed in Table 2. Boldfaced digits in the computed eigenvalues show the extent of agreement with the exact values. Thus is very clear that preconditioning Chebfun we can considerably improve its accuracy.
4.2.2 A fourth order eigenproblem from spherical geometry
In GTD the authors consider the eigenproblem
$$  (12) 
where the operator $D$ is defined by
$$D\left(u\right):={u}^{\prime \prime}+\frac{2}{r}{u}^{\prime}\frac{l\left(l+1\right)}{{r}^{2}},$$ 
and $l$ is a positive integer. They solve this problem by a modified ChT method in order to avoid spurious eigenvalues. We have solved this problem by ChC and obtained for the fist two eigenvalues the numerical values
$$\begin{array}{cc}{s}_{1}=3.947819275687863e+01,& {s}_{2}=8.076297512888706e+01,\end{array}$$ 
which agree up to the fourth decimal with the the true values (determined from the eigencondition).
The first four vectors to problem (12) computed by Chebfun are depicted in Fig. 7 a). It is visible that they satisfy the boundary conditions. Their Chebyshev coefficients are displayed Fig. 7 b). About the first twenty coefficients of the first four eigenvectors decrease just as abruptly and smoothly.
4.2.3 A family of sixth order eigenproblems
In FMB the authors consider the following three eigenproblems of sixth order
$${u}^{\left(vi\right)}\left(x\right)=\lambda {u}^{\left(j\right)}\left(x\right),u\left(\pm 1\right)={u}^{\prime}\left(\pm 1\right)={u}^{\prime \prime}\left(\pm 1\right)=0,j=0,2,4,$$  (13) 
and introduce an extremely simple modification to the ChT method which eliminates the spurious eigenvalues when such high order eigenproblems are solved.
We have tried to solve problem (13) with $j:=4$ by Chebfun but all our attempts failed due to the very small conditioning of matrix involved, i.e., around $O\left({10}^{40}\right).$ The situation became much better with the preconditioner introduced in Subsection 3.4.
Thus, the first four eigenvectors along with their Chebyshev coefficients are depicted in Fig. 8. As it is apparent from the lower panel of this figure the coefficients of degree up to 30 drop sharply to an absolute value below ${10}^{10}$ and then slowly decrease to machine accuracy. This happens at an degree around $120.$
The eigenvalues computed by Chebfun agree up to the first three digits with those provided in HS p.405.
4.3 Problems with mixed boundary conditions
4.3.1 The free lateral vibration of a uniform clampedhinged beam
The fourth order eigenproblem
$$  (14) 
is considered in Mai and is solved by a non conventional spectral collocation method. In this paper the author shows that the eigenvalues satisfy the transcendental equation
$$\mathrm{tanh}\sqrt[4]{\lambda}=\mathrm{tan}\sqrt[4]{\lambda}.$$  (15) 
It is extremely important to observe that neither preconditioning nor ${D}^{2}$ strategy can handle the mixture of boundary conditions in (14). Thus, this problem tests how well the Chebfun can cope with various boundary conditions.
As the eigenvalues computed from (15) are compared with those obtained by Magnus expansion in PB , we report in Table 3 the latter eigenvalues compared with those provided by Chebfun. A coincidence of at least three decimals can be observed.
$j$  ${\lambda}_{j}$ by Chebfun  ${\lambda}_{j}$ According to PB 

$1$  $\mathbf{2.377}373239875730e+02$  $\mathbf{2.377}210675300e+02$ 
$2$  $\mathbf{2.496}524908617440e+03$  $\mathbf{2.496}487437860e+03$ 
$3$  $\mathbf{1.0867}83642364734e+04$  $\mathbf{1.0867}58221697e+04$ 
$4$  $\mathbf{3.17}7977410414838e+04$  $\mathbf{3.17}8009645380e+04$ 
$5$  $\mathbf{7.400}167551416633e+04$  $\mathbf{7.400}084934040e+04$ 
4.3.2 A fourth order eigenproblem with higher order boundary conditions
In order to show again the Chebfun versatility in introducing boundary conditions we consider the following problem called the cantilevered beam in EulerBernouilli theory (see for instance ZWX ). The equation simply reads
$${u}^{\left(iv\right)}={\beta}^{4}u,x\in (0,\pi ),$$  (16) 
and is equipped with the following boundary conditions
$$\begin{array}{cc}u\left(0\right)=0,& {u}^{\prime}\left(0\right)=0,\\ {u}^{\prime \prime}\left(\pi \right)=0,& {u}^{\prime \prime \prime}\left(\pi \right)=0.\end{array}$$  (17) 
The first two boundary conditions state that the beam is clamped in $0$ and the last two state that the beam is free in the right hand end. The eigenvalues $\beta $ satisfy the eigencondition
$$cosh\left(\beta \pi \right)cos\left(\beta \pi \right)+1=0.$$  (18) 
Actually the problem (16)(17) is selfadjoint. Without going into details we will notice that the first eigenvalue of this eigenproblem is the solution of the minimization problem
$${\beta}^{4}=\underset{v\in V}{\mathrm{min}}\frac{{\int}_{0}^{\pi}{\left({v}^{\prime \prime}\right)}^{2}\mathit{d}x}{{\int}_{0}^{\pi}{v}^{2}\mathit{d}x},$$ 
where, roughly, $V$ is a space of continuous functions satisfying the boundary conditions (17). This Ritz formulation, as well as a weak (variational) formulation can be obtained by multiplying the equation with a function $v$ from $V$ and a double integration by parts.
Again, these boundary conditions are not treatable by preconditioning or ${D}^{2}$ strategy.
$j$  ${\beta}_{j}$ by Chebfun  ${\beta}_{j}$ exact solutions of (18) 

$1$  $\mathbf{5.96}7718563107258e01$  $\mathbf{0.596}86$ 
$2$  $\mathbf{1.4941}63617547652e+00$  $\mathbf{1.4941}8$ 
$3$  $\mathbf{2.5002}44462376521e+00$  $\mathbf{2.5002}5$ 
$4$  $\mathbf{3.49999}0154542449e+00$  $\mathbf{3.49999}$ 
The following simple and short Chebfun code solves the problem (16)(17).
% Cantilevered beam in EulerBernouilli theory dom=[0,pi];x=chebfun(’x’,dom); % the domain L = chebop(dom); L.op = @(x,y) diff(y,4); % the operator L.lbc = @(y)[y; diff(y,1)]; % fixed b. c. L.rbc = @(y)[diff(y,2); diff(y,3)];% free b. c. [U,D]=eigs(L,40,’SM’); % first six eigs. % Sorted eigenpairs (eigenvalues and eigenvectors) D=diag(D); [t,o]=sort(D); D=D(o); disp((D.^(1/4))) U=U(:,o);
In Table 4 are reported the first four eigenvalues computed by Chebfun and by Magnus expansion. A satisfactory agreement is observed.
The first four eigenvectors are displayed in Fig. 11 a) and their Chebyshev coefficients are displayed in the same figure panel b). It is clear that Chebfun uses Chebyshev polynomials of slightly lower degree than $20$ and from this level only a sharply decreasing roundingoff plateau follows.
The first four eigenvectors approximating the eigenfunctions are in very good agreement with those exposed in literature. Using the definition of the scalar product of two vectors $u$ and $v$, namely ${u}^{\prime}\ast v$, we can easily check the orthonormality of eigenvectors.
The curves in Fig. 12 clearly shows that the eigenvectors of this problem computed by Chebfun are orthonormal.
4.3.3 The harmonic oscillator and its second and third powers.
We wanted to test our strategy on a problem whose differential equation exhibits stiffness in at least part of the range. Thus, along with the well known harmonic oscillator operator
$$h\left(u\right):={u}^{\prime \prime}+{x}^{2}u,x\in (\mathrm{\infty},\mathrm{\infty})$$ 
we have considered its second and third powers, namely
$${h}^{2}\left(u\right)={u}^{\left(iv\right)}2{\left({x}^{2}{u}^{\prime}\right)}^{\prime}+\left({x}^{4}2\right)u,x\in (\mathrm{\infty},\mathrm{\infty}),$$ 
and
$${h}^{3}\left(u\right)={u}^{\left(vi\right)}+{\left(3{x}^{2}{u}^{\prime \prime}\right)}^{\prime \prime}+{\left(\left(83{x}^{4}\right){u}^{\prime}\right)}^{\prime}+\left({x}^{6}14{x}^{2}\right)u,x\in (\mathrm{\infty},\mathrm{\infty}).$$ 
Actually we want to solve the fourth order eigenvalue problem for ${h}^{2}\left(u\right),$ namely
$${h}^{2}\left(u\right)=\lambda u,$$  (19) 
and the sixth order eigenvalue problem
$${h}^{3}\left(u\right)=\lambda u,$$  (20) 
corresponding to the cube of the harmonic oscillator operator. The eigenvalues of the harmonic oscillator are ${\lambda}_{k}=\left(2k+1\right),$ $k=0,1,2,\mathrm{\dots},$ and those of ${h}^{2}$ and ${h}^{3}$ are the second and the third powers respectively of ${\lambda}_{k}.$ According to the definition for classification of SL problems, given in this Section 2, the eigenproblems (19) and (20) are singular.
We have to observe that no boundary conditions are needed because the problem is of limitpoint type Everitt : the requirement that the eigenfunctions be square integrable suffices as a boundary condition. In GM the problem (20) is solved by a SLEUTH code along with domain truncation. Actually the authors truncate this problem to the interval $(100,100)$, and impose the simplest boundary conditions $u={u}^{\prime}={u}^{\prime \prime}=0$ at $x=\pm 100.$ Along with these boundary conditions the eigenproblem becomes selfadjoint.
The first four eigenvectors of the cube of harmonic oscillator computed by SiC are displayed in the upper panels of Fig. 13. Their sinc coefficients are displayed in the lower panel of the same figure. Roughly speaking it guarantees us an accuracy of at least ${10}^{12}$ in the computation of these eigenvectors.
SiC computes the integers ${\lambda}_{k}$ with an accuracy of at least six digits.
The relative drift, with respect to $N,$ of the first $250$ eigenvalues of the cube of harmonic oscillator computed by SiC are displayed in Fig. 14 a). It tells us that the first $200$ eigenvalues are “good”within an accuracy of approximately ${10}^{2}.$ It means the "highest" confirmation of the EigenvaluesRuleofThumb formulated by Boyd (see Sect. 3.3).
Trying to explain this spectacular phenomenon we cannot forget the fact that the derivation matrices of SiC are symmetric. This leads to normal matrices (operators) whose eigenpairs are properly computable.
If we compare this result with Table 10 from GM , where the best accuracy in computing of the first eigenvalue is ${10}^{2}$, we can speak of a total superiority of SiC method over SLEUTH.
In Fig. 14 b) in a loglinear plot we display the scalar product of some eigenvectors. They prove that the SiC computed eigenvectors are indeed orthonormal. This means that we can trust at least the first $200$ eigenpairs computed by SiC. The following few lines of MATLAB compute the above eigenpairs:
% The sinc differentiation matrices [Weideman & Reddy] N=400;M=6;h=0.1; %Orders of approximation and differentiation and scaling factor [x, D] = sincdif(N, M, h); D1=D(:,:,1);D2=D(:,:,2);D6=D(:,:,6); % The cube of the "harmonic oscillator" operator L=D6+D2*(3*diag(x.^2)*D2)+D1*(diag(83*(x.^4))*D1)+diag(x.^614*(x.^2)); % Finding eigenpairs of L [U,S]=eigs(L,250,0); S=diag(S); [t,o]=sort(S); S=S(o);
Unfortunately, Chebfun along with domain truncation fails in solving the sixth order problem (20) with or without preconditioning. Actually, a warning concerning the very bad conditioning of the matrix is issued.
However, Chebfun behaves fairly well in solving the fourth order problem (19), i.e., computes the corresponding integers with the same accuracy as SiC. We have solved the equation (19) on the truncated interval $(X,X)$ for various $X$ along with the boundary conditions $u(\pm X)={u}^{\prime}(\pm X)=0.$
The Chebyshev coefficients of the first four eigenvectors of eigenproblem (19) computed by Chebfun are displayed in Fig. 15 a). An important aspect must be highlighted, namely the first about $1000$ polynomial coefficients decrease steeply and smoothly to about ${10}^{14}$ after which up to the order of $2500$ follows a wide roundingoff plateau.This is the polynomial of the highest degree that Chebfun has used in our numerical experiments. The curves in Fig. 15 b) show the orthonormality of Chebfun eigenvectors.
The relative drift with respect to the length of integration interval $X$ of the first $250$ eigenvalues to problem (19), when it is solved by Chebfun, is displayed in Fig. 16. It means that the numerical stability is lost for length $X$ larger than $100$ and a set of small eigenvalues can be computed with an accuracy better than ${10}^{9}$.
5 Conclusions and open problems
After analyzing these challenging problems, some firm conclusions can be drawn.
First of all, Chebfun can easily handle any type of boundary condition. This is a significant advantage. Thus, for fourth order eigenproblems, the direct application of Chebfun is versatile in handling various high order boundary conditions and produces reliable outcomes. Furthermore, for problems with clamped boundary conditions both methods, Chebfun as well as ChC, improve their results with at least two, three decimals by preconditioning.
For sixth order eigenproblems, the Chebfun situation is not so encouraging. Its direct application is very uncertain. Matrices whose conditioning order drops to ${10}^{40}$ appear, which most often lead to inaccurate results. For problems of this order or more, subjected to hinged boundary conditions, the reduction to secondorder systems and then the application of the ChC method is the best strategy. In this way we managed to establish a conjecture for the Viola’s problem regarding its lowest eigenvalue.
For fourth order problems on the real line Chebfun along with the truncation of the domain worked fairly well as was the case with the second power harmonic oscillator. As an absolute novelty, we have established in this case the numerical stability with respect to the length of the integration interval. Instead, for the sixth order eigenproblems on the real line the SiC method remains the unique feasible alternative.
In fact, this method is the best in the sense that we can trust the first half of computed eigenpairs. To our knowledge, no software package has reached this performance so far.
As an open problem remains that of finding preconditioning methods for the case of hinged boundary conditions or some other types of boundary conditions.
This paper comes shortly after in another one (see Cig21 ) we have approached, with the same two classes of methods, singular Schröedinger eigenproblems. In this situation we can appreciate that, ChC along with Chebfun, are a better alternative in some respects to other existing methods for a very wide range of eigenproblems. Both compute eigenvectors (approximating eigenfunctions) and by drift estimation demonstrate numerical stability. In addition, the drift with respect to $N$ shows the degree of accuracy up to which a set of eigenvalues is computed. The situation when both types of methods can be applied to the same problem is the ideal one and the one that produces the safest results.
Conflicts of interest: The author declares no conflict of interest.
Sample availability: Samples of the Chebfun codes are available from the author by request.
Abbreviations: The following abbreviations are used in this manuscript:
ChC  Chebyshev collocation method 

ChT  Chebyshev tau method 
${D}^{2}$  strategy to reduce a 2nthorder equation to a second order system 
FCT  fast Chebyshev transform 
FD  finite difference method 
FE  finite element method 
MATSLISE  a MATLAB package for the numerical solution of SL and Schröedinger equations 
SiC  sinc spectral collocation 
SL  SturmLiouville 
SLEDGE  SturmLiouville estimates determined by global errors 
SLEUTH  SturmLiouville Eigenvalues using Theta Matrices 
References
 (1) Straughan, B.: Stability of Wave Motion in Porous Media, Springer Science+Business Media, LLC 2008
 (2) Trefethen, L.N., Birkisson, A., Driscoll, T. A. Exploring ODEs. SIAM, Philadelphia, US, 2018
 (3) Pruess, S., Fulton, C. T. Mathematical Software for SturmLiouville Problem. ACM T. Math. Software 1993, 19, 360–376.
 (4) Pruess, S., Fulton, C. T., Xie, Y. An Asymptotic Numerical Method for a Class of Singular SturmLiouville Problems. ACM SIAM J. Numer. Anal. 1995 32, 1658–1676.
 (5) Marletta, M., Pryce, J.D. LCNO SturmLiouville problems. Computational difficulties and examples. Numer. Math. 1995, 69, 303–320.
 (6) Pryce, J.D., Marletta, M. A new multipurpose software package for Schrödinger and Sturm–Liouville computations. Comput. Phys. Comm. 1991, 62, 42–54.
 (7) Bailey, P. B., Everitt, W. N., Zettl, A. Computing Eigenvalues of Singular SturmLiouville Problems. Results Math. 1991, 20, 391–423.
 (8) Bailey, P.B., Garbow, B., Kaper, H., Zettl, A. Algorithm 700: A FORTRAN software package for SturmLiouville problems. ACM T. Math. Software. 1991, 17, 500–501.
 (9) Ledoux, V., Van Daele, M., Vanden Berghe, G. MATSLISE: A MATLAB Package for the Numerical Solution of SturmLiouville and Schrödinger Equations. ACM T. Math. Software. 2005, 31, 532–554.
 (10) Baeyens, T., Van Daele, M. The fast and accurate computation of eigenvalues and eigenfunctions of timeindependent onedimensional Schrödinger equations. Comput. Phys. Comm. 2021, 258 https://doi.org/10.1016/j.cpc.2020.107568
 (11) Abbasbandy, S., Shirzadi, A. A new application of the homotopy analysis method: Solving the Sturm–Liouville problems. Commun. Nonlinear. Sci. Numer. Simulat. 2011, 16, 112–126.
 (12) Perera, U.; Böckmann, Ch. Solutions of Direct and Inverse EvenOrder SturmLiouville Problems Using Magnus Expansion. Mathematics 2019, 7, 544; . https://doi.org/10.3390/math7060544
 (13) Gheorghiu, C.I. Spectral Methods for NonStandard Eigenvalue Problems. Fluid and Structural Mechanics and Beyond. SpringerVerlag, Cham Heidelberg NewYork Dondrecht London, 2014
 (14) Greenberg, G., Marletta, M. Oscillation Theory and Numerical Solution of Sixth Order SturmLiouville Problems. SIAM J. Numer. Anal. 1998, 35, 20702098.
 (15) Everitt, W.N. A Catalogue of SturmLiouville Differential Equations. In SturmLiouville theory: Past and Present; Amrein, W. O., Hinz, A. M., Hinz, D. B., Eds. Birkhäuser Verlag Basel/Switzerland, 2005; pp.271–331.
 (16) Driscoll, T. A., Bornemann, F., Trefethen, L.N. The CHEBOP System for Automatic Solution of Differential Equations. BIT 2008, 48, 701–723.
 (17) Driscoll, T. A., Hale, N., Trefethen, L.N. Chebfun Guide. Pafnuty Publications, Oxford, UK 2014
 (18) Driscoll, T. A., Hale, N., Trefethen, L.N. Chebfunnumerical computing with functions. http://www.chebfun.org. Accessed 15 November 2019
 (19) Gheorghiu, C.I. Spectral Collocation Solutions to Problems on Unbounded Domains. Casa Cărţii de Ştiinţă Publishing House, ClujNapoca, Romania 2018
 (20) Weideman, J. A. C., Reddy, S. C. A MATLAB Differentiation Matrix Suite. ACM T. Math. Software 2000, 26, 465–519. 26, 465–519
 (21) Gheorghiu, C.I., Pop, I.S. A Modified ChebyshevTau Method for a Hydrodynamic Stability Problem. In Proceedings of the International Conference on Approximation and Optimization; Stancu,D.D., Coman, Gh., Breckner, W. W., Blaga, P., Eds.; Transilvania Press: ClujNapoca, Romania, Volume II, 1996; pp. 119–126.
 (22) Boyd, J. P. Traps and Snares in Eigenvalue Calculations with Application to Pseudospectral Computations of Ocean Tides in a Basin Bounded by Meridians. J. Comput. Phys. 1996, 126, 11–20.
 (23) Gheorghiu, C.I. Accurate Spectral Collocation Computation of High Order Eigenvalues for Singular Schrödinger Equations. Computation 2021, 9, 2. https://doi.org/10.3390/computation9010002
 (24) Boyd, J. P. Chebyshev and Fourier Spectral Methods. Dover Publications, NewYork, US, 2000; pp. 127–158
 (25) Huang, W., Sloan, D.M. The Pseudospectral Method for Solving Differential Eigenvalue Problems. J. Comput. Phys. 1994, 111, 399–409.
 (26) Straughan, B. The Energy Method, Stability, and Nonlinear Convection, SpringerVerlag, NewYork, lnc. 1992; pp. 218–222.
 (27) Lesnic, D., Attili, S. An Efficient Method for Sixthorder SturmLiouville Problems. Int. J. Sci. Technol. 2007, 2, 109–114.
 (28) Gheorghiu, C.I., Rommes, J. Application of the JacobiDavidson method to accurate analysis of singular linear hydrodynamic stability problems. Int. J. Numer. Meth. Fluids 2013, 71, 358–369.
 (29) Gardner, D.R., Trogdon, S.A., Douglass, R.W. A Modified Spectral Tau Method That Eliminates Spurious Eigenvalues. J. Comput. Phys. 1989, 80, 137–167.
 (30) Fornberg. B. A practical guide to pseudospectral methods. Cambridge University Press, Cambridge, UK, 1998; pp.89–90
 (31) McFadden, G.B., Murray, B.T., Boisvert, R.F. Elimination of Spurious Eigenvalues in the Chebyshev Tau Spectral Method. J. Comput. Phys. 1990, 91, 228–239.
 (32) MaiDuy, N. An effective spectral collocation method for the direct solution of highorder ODEs. Commun. Numer. Methods Eng. 2006, 22, 627–642.
 (33) Zhao, S.; Wei, G.W.; Xiang,Y. DSC analysis of freeedged beams by an iteratively matched boundary method. J. Sound. Vib. 2005, 284, 487–493.