Accurate Spectral Collocation Solutions to 2nd-Order Sturm–Liouville Problems

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.

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 self-adjoint problems, while non-self-adjoint problems arise in hydrodynamic and magnetohydrodynamic stability theory (see for instance [1] 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., “discretize-then-solve”. On the contrary, the Chebfun spirit consists in the continuous mode, i.e., “solve-then-discretize” (see [2] 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 [3,4], the NAG’s code SL02F [5,6], SLEIGN and SLEIGN2 [7,8], 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 user-supplied tolerance.

Essentially, the numerical method used in these software packages replaces the coefficients in the equation by a 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 [9] can solve some Schrödinger eigenvalue problems by a constant perturbation method of a higher order. Very recently, this code has been improved (see [10]) but it remains for Schrödinger issues which are outside the scope of this paper.

There is also a class of semi-analytical methods which includes the variational iteration method, the homotopy perturbation method, homotopy analysis method, and Adomian decomposition (see for instance [11]) for solving eigenvalue problems. Their accuracy is far from what spectral collocation methods can provide.

In [12], the authors set up an ambitious method based on the 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 so-called removing technique of independent boundary conditions (see for a review of this technique our monograph [13] 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 the 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, and reliable. All our numerical experiments have been carried out using MATLAB R2020a on an Intel (R) Xeon (R) CPU E5-1650 0 @ 3.20 GHz.

The main purpose of this paper was 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 self-adjoint 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 Sturm–Liouville 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. 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. 2nd-Order Sturm–Liouville Eigenproblems

The 2nd-order SL equation reads

$$\begin{array}{c}{\left(-1\right)}^n{\left(p_n\left(x\right)u^{\left(n\right)}\right)}^{\left(n\right)}+{\left(-1\right)}^{n-1}{\left(p_{n-1}\left(x\right)u^{\left(n-1\right)}\right)}^{\left(n-1\right)}+\dots \hfill \\ \hfill +{\left(p_2\left(x\right)u^{″}\right)}^{″}-{\left(p_1\left(x\right)u^{\prime }\right)}^{\prime }=\lambda u\left(x\right),a<x<b,n\in \mathbb{N},n\ge 2,\end{array}$$

along with separated, (self-adjoint) 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$, $0<k<n-1$, 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 [14]).

The eigenvalues can be ordered in the usual form: ${\lambda }_0\le {\lambda }_1\le {\lambda }_2\le ...,$ such that ${lim}_{k\to \infty }{\lambda }_k=+\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 [15]. 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 [2,16,17,18]. The Chebfun system, in object-oriented 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 [17]. In [16] 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 [2]. Thus, one can get a suggestive image of what they can do working with Chebfun.

Moreover, in ([16] 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 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 [19]. In all spectral collocation methods designed so far, we have used the collocation differentiation matrices from the seminal paper [20]. 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 [21] 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 Orr-Sommerfeld 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 (cut-off parameter) N or the scaling factor (length of integration interval). Only those whose difference or “resolution-dependent drift” is “small” can be believed. In this connection, in the paper [22], the so called absolute (ordinal) drift with respect to the order of approximation has been introduced. We extend this definition in our recent paper [23] 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 [24], the author makes a subtle analysis of spectral methods in solving linear eigenproblems. Among others, he states the so called Boyd’s Eigenvalues Rule-of-Thumb 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 \left(-1,1\right),$$

(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 \infty .$ We comment on the preconditioner introduced in [25]. 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 N-1,$$

applied to ChC as well as to Chebfun discretization ${\mathbf{L}}_{ChC\mathrm{or}Chebfun}^{\left(n\right)}$ of differential operator $L^{\left(n\right)},$ produces matrices $\mathbf{D}{\mathbf{L}}_{ChC\mathrm{or}Chebfun}^{\left(n\right)}$ 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 Eigenproblem-Revisited

Let us consider now the so called Viola’s eigenproblem. It is encountered in porous stability problems (see [1], Chapter 9) and reads

$$\begin{array}{c}\frac{d^2}{d{x}^2}\left[{\left(1-\theta x\right)}^3\frac{d^{2}u}{d{x}^2}\right]=\lambda \left(1-\theta x\right)u,x\in \left(0,1\right),0\le \theta <1,\\ u\left(0\right)=u^{″}\left(0\right)=u\left(1\right)=u^{″}\left(1\right).\end{array}$$

(4)

It is singular as $\theta \to 1^{-}$ in accordance with the definition introduced in Section 2.

By straightforward variational arguments, we have shown in ([13] p. 50) that the 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 Figure 1. Actually, we have obtained

$${\lambda }_1\left(0.98765\right)=8.775218471808549e-01,$$

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 Figure 2 the first four eigenvectors of problem (4) and in Figure 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 uses it and when $\theta :=0.98765$ the degree of approximation grows to more than 80 (see Figure 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.

We have to observe that the problem (4) has been solved by compound matrix method in [26] for $\theta <0.9.$ The author asserts that other methods have to be used in order to resolve the singularity in this problem. We hope that the above analysis sheds some light in this direction.

4.1.2. The Bénard Stability Problem

A simplified form of the Bénard stability problem supplied with self-adjoint boundary conditions reads (see for instance [14,27])

$$\begin{array}{c}{u}^{\left(vi\right)}-\left(2\nu +3\right)u^{\left(iv\right)}+\left({\nu }^2+4\nu +3\right)u^{″}-\left({\nu }^2+2\nu +2\right)u=\lambda u\left(x\right),x\in \left(0,1\right),\\ u\left(0\right)=u\left(1\right)=u^{″}\left(0\right)=u^{″}\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,\dots .$$

All our attempts to solve this problem using Chebfun have failed, so we have resorted to the old $D^2$ strategy.

Thus, we rewrite problem (5) as a homogeneous Dirichlet one attached to a second order differential system, namely

$$\begin{array}{c}{u}^{″}=v\left(x\right),x\in \left(0,1\right),\\ v^{″}=w\left(x\right),\\ w^{″}-\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\mathrm{in}x=0\mathrm{and}x=1.\end{array}$$

(6)

Now we apply to each line the ChC discretization. It leads to the generalized and singular eigenpencil

$$\left(\mathbf{A},\mathbf{B}\right),$$

(7)

where the block matrices are defined by

$$\mathbf{A}:=\left(\begin{array}{ccc}4{\tilde{D}}^{\left(2\right)}& I& Z\\ Z& 4{\tilde{D}}^{\left(2\right)}& I\\ -\left({\nu }^2+2\nu +2\right)I& \left({\nu }^2+4\nu +3\right)I& 4{\tilde{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 ${\tilde{D}}^{\left(2\right)}$ comes from the shift of interval $\left(0,1\right)$ to the canonical Chebyshev interval $\left[-1,1\right]$ and the matrix ${\tilde{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 ${\tilde{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:N-1,2:N-1,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 (N-1).$ 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 [13] pp. 22–23), we see from the inequality

$$Henrici\left(\mathbf{A}\right)=0.300293<Henrici\left({\tilde{D}}^{\left(2\right)}\right)=0.395205,$$

that matrix $\mathbf{A}$ is more normal than ${\tilde{D}}^{\left(2\right)}$.

In our previous paper [28], 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. The first two eigenvalues of Bénard problem (5) for various $\nu$ computed by $D^2$ strategy along with ChC.

$\mathit{\nu }$	${\mathit{\lambda }}_0\left(\mathit{\nu }\right)$	${\mathit{\lambda }}_1\left(\mathit{\nu }\right)$	${\mathit{\lambda }}_0\left(\mathit{\nu }\right)$ according to [14]
$-(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}$

. 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 [14]. 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 Figure 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 [14] as well as [27]. Actually, in [14] 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 Self-Adjoint Eighth-Order Problem

The eigenproblem of the highest order we consider in this paper is the following

$$\begin{array}{c}{u}^{\left(8\right)}\left(x\right)=\lambda u\left(x\right),0<x<1,\hfill \\ \hfill u\left(0\right)=u^{″}\left(0\right)=u^{\left(4\right)}\left(0\right)=u^{\left(6\right)}\left(0\right)=u\left(1\right)=u^{″}\left(1\right)=u^{\left(4\right)}\left(1\right)=u^{\left(6\right)}\left(1\right)=0,\end{array}$$

(8)

with exact eigenvalues ${\lambda }_k={\left(k\pi \right)}^8,k=1,2,\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 Figure 5.

In order to estimate the error with which the eigenvalues were calculated, we display in Figure 6b the relative drift of the first twelve eigenvalues for different approximation orders. As a result that we know 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 [12].

Moreover, this means that we cannot trust more than twelve eigenvalues for this problem.

It is clear that ChC along with the $D^2$ method have the potential to find the first eigenvalues of an 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 Figure 6a. 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 [25,29], as well as in the monograph [30], the following eigenproblem is carefully studied. It consists of the fourth order differential equation

$$u^{\left(iv\right)}+R{u}^{‴}=s{u}^{″},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

$${\left(R^2+4s\right)}^{1/2}\left[cosh\left(R\right)-cosh{\left(R^2+4s\right)}^{1/2}\right]+2ssinh{\left(R^2+4s\right)}^{1/2}=0.$$

(11)

Problems similar to this appear, for example, in linearized stability analysis in fluid dynamics. In [29], the authors noticed spurious eigenvalues when the problem (9) and (10) is solved by ChT method. These spurious eigenvalues appears in the right-half plane suggesting physical instabilities that do not exist.

We have solved the problems (9) and (10) by Chebfun with and without preconditioning. The numerical outcomes are displayed in

Table 2. First two eigenvalues of problems (9) and (10).

i	${\mathit{\lambda }}_{\mathit{i}}$ Chebfun	${\mathit{\lambda }}_{\mathit{i}}$ Preconditioned Chebfun	${\mathit{\lambda }}_{\mathit{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$

. Boldfaced digits in the computed eigenvalues show the extent of agreement with the exact values. Thus, it is very clear that preconditioning Chebfun can considerably improve its accuracy.

4.2.2. A Fourth Order Eigenproblem from Spherical Geometry

In [29], the authors consider the eigenproblem

$$\begin{array}{c}D\left(D-s\right)u=0,0<r_1<r_2,r_1>0,\\ u\left(r_1\right)=u\left(r_2\right)=u^{\prime }\left(r_1\right)=u^{\prime }\left(r_2\right)=0,\end{array}$$

(12)

where the operator D is defined by

$$D\left(u\right):=u^{″}+\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 Figure 7a. It is visible that they satisfy the boundary conditions. Their Chebyshev coefficients are displayed Figure 7b. About the first twenty coefficients of the first four eigenvectors decrease just as abruptly and smoothly.

4.2.3. A Set of Sixth Order Eigenproblems

In [31], the authors consider the following sixth order eigenproblems

$$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^{″}\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 Section 3.4.

Thus, the first four eigenvectors along with their Chebyshev coefficients are depicted in Figure 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 ([25] p. 405).

4.3. Problems with Mixed Boundary Conditions

4.3.1. The Free Lateral Vibration of a Uniform Clamped–Hinged Beam

The fourth order eigenproblem

$$u^{\left(iv\right)}\left(x\right)=\lambda u\left(x\right),0<x<1,u\left(0\right)=u^{\prime }\left(0\right)=u\left(1\right)=u^{″}\left(1\right)=0,$$

(14)

is considered in [32] and is solved by a non conventional spectral collocation method. In this paper, the author shows that the eigenvalues satisfy the transcendental equation

$$tanh\sqrt[4]{\lambda }=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 [12], we report in

Table 3. The first five eigenvalues to problem (14) computed by Chebfun and compared with those provided by Magnus expansion.

j	${\mathit{\lambda }}_{\mathit{j}}$ Chebfun	${\mathit{\lambda }}_{\mathit{j}}$ According to [12]
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$

the latter eigenvalues compared with those provided by Chebfun. A coincidence of at least three decimals can be observed.

The first four eigenvectors to problem (14) computed by Chebfun are displayed in Figure 9. It is perfectly visible that they satisfy the boundary conditions assumed in this problem.

The Chebyshev coefficients of the first four eigenvectors to problem (14) computed by Chebfun are displayed in Figure 10. They decrease smoothly to somewhere around ${10}^{-12}$ which is an argument in favor of the accuracy of numerical results.

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 Euler–Bernouilli theory (see for instance [33]). The equation simply reads

$$u^{\left(iv\right)}={\beta }^{4}u,x\in \left(0,\pi \right),$$

(16)

and is equipped with the following boundary conditions

$$\begin{array}{cc}u\left(0\right)=0,& u^{\prime }\left(0\right)=0,\\ u^{″}\left(\pi \right)=0,& u^{‴}\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 bean 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 problems (16) and (17) are self-adjoint. 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}{min}\frac{{\int }_0^{\pi }{\left(v^{″}\right)}^{2}dx}{{\int }_0^{\pi }{v}^{2}dx},$$

where, roughly, V is a space of continuous functions satisfying the boundary conditions in (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.

The following simple and short Chebfun code solves the problems (16) and (17).

% Cantilevered beam in Euler-Bernouilli 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. The first four eigenvalues of problems (16) and (17) computed by Chebfun compared with numerical solutions to Equation (18).

j	${\mathit{\beta }}_{\mathit{j}}$ by Chebfun	${\mathit{\beta }}_{\mathit{j}}$ Exact Solutions of (18)
1	$\mathbf{5.96}7718563107258e-01$	$\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 first four eigenvalues computed by Chebfun and by Magnus expansion are reported. A satisfactory agreement is observed.

The first four eigenvectors are displayed in Figure 11a 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 rounding-off 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 Figure 12 clearly show 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^{″}+x^{2}u,x\in \left(-\infty ,\infty \right)$$

we will consider 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 \left(-\infty ,\infty \right),$$

and

$$h^3\left(u\right)=-u^{\left(vi\right)}+{\left(3x^2{u}^{″}\right)}^{″}+{\left(\left(8-3x^4\right)u^{\prime }\right)}^{\prime }+\left(x^6-14x^2\right)u,x\in \left(-\infty ,\infty \right).$$

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,\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 limit-point type [15]: the requirement that the eigenfunctions be square integrable suffices as a boundary condition. In [14] 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^{″}=0$ at $x=\pm 100.$ Along with these boundary conditions the eigenproblem becomes self-adjoint.

The first four eigenvectors of the cube of harmonic oscillator computed by SiC are displayed in the upper panels of Figure 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 Figure 14a). It tells us that the first 200 eigenvalues are “good" within an accuracy of approximately ${10}^{-2}.$ It also means the “highest” confirmation of Boyd’s Eigenvalues Rule-of-Thumb (see Section 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 [14], 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 Figure 14b, we display the scalar product of some eigenvectors. They prove that the SiC computed eigenvectors are 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(8-3*(x.^4))*D1)+diag(x.^6-14*(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 Figure 15a). 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 rounding-off plateau.This is the polynomial of the highest degree that Chebfun has used in our numerical experiments. The curves in Figure 15b) 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 Figure 16. It means that the numerical stability is lost for larger X 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 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 second-order 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.

An open problem remains for finding preconditioning methods for the case of hinged boundary conditions or some other types of boundary conditions.

This paper comes shortly after when in another one (see [23]) 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.