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

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

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

which only partly confirms the above conjecture.

As a validation issue for our computations we have obtained the known value

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.

We have to observe that the problem (4) has been solved by compound matrix method in S0 for 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.

A simplified form of the Bénard stability problem supplied with self-adjoint boundary conditions reads (see for instance GM and LA )

The constant $\nu$ is regarded as a parameter which typically can take the values

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

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

where the block matrices are defined by

and

The factor $4$ in front of ${\tilde{D}}^{\left(2\right)}$ comes from the shift of interval $(0,1)$to the canonical Chebyshev interval $[-1,1]$ 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 Cig14 pp. 22-23) we see from the inequality

that matrix $𝐀$ is more normal than ${\tilde{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 $𝐀$ is non singular and the block matrix $𝐁$ 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.”

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

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.

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

supplied with the clamped boundary conditions

The eigencondition for this problem is

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 right-half plane suggesting physical instabilities that do not exist.

We have solved the problem (9)-(10) by Chebfun with and without preconditioning.

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.

In GTD the authors consider the eigenproblem

where the operator $D$ is defined by

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

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.

In FMB the authors consider the following three eigenproblems of sixth order

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.

The fourth order eigenproblem

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

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.

The first four eigenvectors to problem (14) computed by Chebfun are displayed in Fig. 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 Fig. 10. They decrease smoothly to somewhere around ${10}^{-12}$ which is an argument in favor of the accuracy of numerical results.

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 ZWX ). The equation simply reads

and is equipped with the following boundary conditions

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

Actually the problem (16)-(17) is self-adjoint. Without going into details we will notice that the first eigenvalue of this eigenproblem is the solution of the minimization problem

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.

The following simple and short Chebfun code solves the problem (16)-(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 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 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 Fig. 12 clearly shows that the eigenvectors of this problem computed by Chebfun are orthonormal.

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

we have considered its second and third powers, namely

and

Actually we want to solve the fourth order eigenvalue problem for $h^2\left(u\right),$ namely

and the sixth order eigenvalue problem

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 limit-point 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 self-adjoint.

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 Eigenvalues-Rule-of-Thumb 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 log-linear 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(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 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 rounding-off 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}$.