The main aim of this paper is to show that singular eigenvalue problems, supplied with boundary conditions containing the eigenparameter in their formulation, defined on the finite or unbounded interval, can be accurately solved by spectral collocation.

From the very beginning of our preoccupation with spectral methods, we observed in [8] that the Galerkin method is practically inapplicable to eigenvalue problems containing the eigenparameter in their formulation. This presence hinders the construction of test and trial space functions. Thus, in [8], we designed a modified tau method in order to solve a hydrodynamic stability problem (Orr-Sommerfeld) containing the eigenparameter in the formulation of two boundary conditions.

In this paper, we attempt to argue on along the same line and consider singular eigenvalue problems. The first one comes from elasticity and is defined on a finite interval.

The critical loads for divergence of a clamped-free elastic bar of finite length $l$ and constant flexural rigidity $\alpha$ exposed to an end load $Q$ are determined by the fourth order non-self-adjoint boundary eigenvalue problem

Here, $\lambda:=-\frac{Ql^{2}}{\alpha}$, and $\gamma\varphi_{l}$, with $\gamma\in[0,1]$, is the angle between the end load $Q$ and the $x$-axis and $\varphi_{l}$ denotes the angle between the tangent at the free end of the bar and the $x$-axis (see [20] and [21] and some older papers quoted there). This system is called Petterson-König’s rod, or, if $Q$ is tangential, that is, $\gamma=1$, Pflüger’s rod. The case $\gamma=0$ leads to the well-known Euler-loads.

In the above quoted papers, Tretter provided important analytical results concerning the completeness, minimality, and Riesz basis properties of the corresponding eigenfunctions and associated functions. They are obtained after a suitable linearization of the problem and by means of a detailed asymptotic analysis of the Green’s functions. The function spaces where the above properties hold are described by $\lambda$-independent boundary conditions (see also [22]).

Consequently, we will also consider the boundary $\lambda$-independent eigenvalue problem attached to (1), namely

Particular attention is paid to the way in which boundary conditions are imposed in collocation methods. However, the numerical experiments carried out show that
the critical eigenvalues found by both collocation methods are practically undistinguishable. They are confirmed by the tau method, which is robust with respect to the issue of boundary conditions. The numerical results plainly confirms some laborious analytical studies such as [20-22]. The dependence of the critical values of the eigenparameter on the involved physical parameter is displayed.

On an infinite domain, i.e., on the real positive axis, we first consider the Charney’s baroclinic stability problem. This reads

where $r\in(0,1]$ is the meteorological parameter and $c$ is the eigenparameter. This problem is the so-called point-jet model of the baroclinic and barotropic instability process (see [4] and [5]). This is a second order problem complicated by the singularity in the differential operator, complex eigenvalues, and a semi-infinite domain. Actually, in [5], a more refined version of problem (3) is considered and in [4], Boyd solves a regularized Charney’s problem.

Then, we consider the so-called Fourier eigenvalue problem. It is a SturmLiouville problem which reads

The singularity of this problem comes from the length of its domain. In [15], the authors study the spectrum of this problem and its approximation with eigenvalues from a sequence of regular problems. The right end point is a limit-point.

In the well-known monograph of Gottlieb and Orszag [13] p. 45, the authors express their doubt with respect to the practical value in applications of collocation method based on Laguerre and Hermite polynomials due to the poor resolution of these polynomials. However, Boyd et al. showed in [3] that Laguerre collocation, based on Laguerre functions works fairly well for a singular eigenproblem (the hydrogen atom) on a semi-infinite interval. Also, the numerical results reported in our papers [10-12] advocated that Laguerre and Hermite collocation worked surprisingly well for various non-standard eigenvalue problems, boundary layer type problems or non-linear wave propagation on unbounded domains. On the other hand, in [2], Boyd clearly states that "Because of its simplicity, domain truncation will undoubtedly continue to be used in the real world even if numerical analysts frown."

However, in order to reduce a problem on an infinite interval to one over $[-1,1]$ and to keep the most of the good numerical characteristics of Chebyshev polynomials, in [1], Boyd applies a mapping to these polynomials and defines a new spectral basis, the so-called rational Chebyshev functions on a semi-infinite interval.

In this paper, as well as in our last three quoted above, we have tried to plead that with a comparable computational effort collocation based on Sinc, Laguerre and Hermite functions can provide a competitive alternative for solving boundary value problems on unbounded domains.

The outline of the paper is as follows. In Section 2, we find out the analytical critical values of the eigenparameter in Petterson-König’s rod eigenvalue problem. In Section 3, we review the Chebyshev collocation (ChC) method based on square as well as the rectangular differentiation matrices. These two types of matrices are compared with respect to their normality. We also report the numerical results obtained using both methods in solving the problems (1) and (2). A rough, but fairly useful estimation, of the accuracy of numerical outcomes obtained by ChC based on square differentiation is outlined. In Section 4, we are concerned with solving the singular algebraic eigenvalue problems and establish the order of square ChC method in solving the problem (1) at least with respect to the first eigenvalue. In Section 5, the problem (1) is reconsidered using Chebyshev tau method (ChT). In Section 6, the "point-jet" model of the Charney’s baroclinic stability problem in solved by Laguerre-Gauss-Radau collocation (LGRC). We also comment on the optimal choice of the scaling parameter. In Section 7, we apply LGRC to the Fourier eigenvalue problem and notice a surprising behavior of the first eigenmode. Section 8 contains conclusions.

The general analytical solution to problem (1) has the form

where the real constants $A,B,C$, and $D$ have to be determined by imposing the four boundary conditions from the problem. Accordingly, the general stability criterion for Petterson-König’s rod reads

Thus, for $\gamma=0$, the above equation provides the critical value $\lambda_{cr}=(\pi/2)^{2}$ which confirms the classical results for clamped (fixed)-free column from [25], Sect. 2.1.2. The first eigenvalue $\lambda_{cr}$ of $\lambda$ for intermediate values of $\gamma$ between 0.1 and 0.9 are reported in Table 1. They have been computed out using the MATLAB built-in function fsolve.

We have worked with three methods in order to solve the singular algebraic eigenvalue problems (6) or (7) and (8) (it is clear that the r. h. s. matrices in (6) and (8) are singular). The methods are QZ, Arnoldi, and Jacobi-Davidson. The first two methods are implemented in MATLAB by built-in functions eig and respectively eigs $\left(A^{-1}B\right)$ and the third is available in [16] and [17].

The first three eigenvectors of the problem (1) obtained by solving (6) and of the problem (2) when $\gamma=0.25$ are displayed in Fig. 1. The three pairs of curves are indistinguishable. The first (critical) eigenvalue equals in this case -3.650437 .

In order to quantify the accuracy of ChC solutions, we use the fcgltran routine (see [27]). This routine allows for fast transformation between nodal values at the Chebyshev-Gauss-Lobatto points and spectral coefficients by using the MATLAB built-in functions fft/ifft.

In Fig. 2 in a semilogy plot, we display the values of the spectral coefficients of the first three eigenvectors vs. the order of approximation. This plot suggests an accuracy of order $10^{-11}$.

By classical square differentiation, we also solve the $\lambda$-independent eigenvalue problem (2). We have got extremely close critical values to those obtained solving the the $\lambda$-dependent eigenvalue problem (1). They are displayed in the first two columns of Table 2.

This confirms the analytical results of Tretter [20].
On the other hand, comparing the first and the third columns of Table 2, we conclude that, at least for a few eigenvalues in the lower region of the spectrum, the square ChC and rectangular ChC perform equally well.

It is of some importance to establish the order of ChC method in solving the problem (1) at least with respect to the first eigenvalue. Thus, assuming that the error has the expression

we obtain

However, the Arnoldi iteration proved to the most accurate and stable for large orders N.

In order to evaluate the non-normality of matrices involved in the algebraic eigenvalue problems of this paper, we use the so-called Henrici’s number (see for instance our monograph [10] p. 22). For a square matrix $\mathbf{M}$, it is defined as

where $\mathbf{M}^{*}$ stands for the conjugate transpose of $\mathbf{M}$ and $\varepsilon(\cdot)$ stands for the Frobenius norm of a matrix. We have deduced that

with $H(\mathbf{M})=0$ iff $\mathbf{M}$ is normal, i.e., $\mathbf{M}^{*}\mathbf{M}-\mathbf{MM}^{*}=0$.
It is well known that the non-normality of a matrix is responsible for its high spectral sensitivity, i.e., the sensitivity of its eigenvalues to small perturbations (see for instance Trefethen [19]).

Thus, a short glance at Table 3 shows that, for the problem at hand, the rectangular differentiation does not improve the normality of matrices involved in the eigenproblems but, on the contrary, makes it slightly worse.

In ChT method, the solution is approximated by

and thus its ith order derivatives can be expressed in the following form (see for instance our contribution [10] Sect. 1.1)

For second and fourth order derivatives, the numerical coefficients $a_{k}^{(i)}$ are given by

and

respectively, with $c_{n}$ defined by

In order to incorporate the boundary conditions in the framework of ChT method, the following equalities are useful

Thus, using the relations (11) and (12) along with boundary relations for $k=1,2,3$, we get the following singular eigenvalue problem

where the unknown vector is defined $\mathbf{X}:=\left(a_{0}a_{1}\ldots a_{N}\right)^{T}$ and the structure (sparsity) of matrices $\mathbf{A}$ and $\mathbf{B}$ is provided in Fig. 3.

They look fairly close to upper triangular matrices and are far less populated than the square and rectangular ChC matrices.

Their Henrici’s numbers, i.e., their non-normality, are displayed in Table 3. Actually, the matrices involved in (15) are only slightly worse with respect to normality than those involved in (6) in spite of their sparsity.

We have used the real variant of $JDQZ$ method provided in [16] to find the first six eigenvalues of (15). The imposed target has been 0 and the tolerance is $1e-06$. The results are twice confirmed by using the MATLAB built-in function eig and Arnoldi iterations in the same way as in [9].

It is clear from Fig. 4 that the roundoff plateau settles at a level lower than $10^{-12}$ and thus the magnitude of error in ChT method could be roughly estimated at $10^{-12}$.

Thus, comparing Figs. 2 and 4, we can conclude that ChT method performs better than both square or rectangular ChC at least by one unit of accuracy. Moreover, ChT makes use of less than half the terms in the approximation expansion.

The LGRC has been successfully used for non-linear boundary layer type problems in our recent paper [11]. The Laguerre-Gauss-Radau quadrature nodes and weights are defined for instance in the monograph of Shen et al. [18] p. 243. The set of nodes contains the origin and in ascending order the zeros of the derivative of the generalized Laguerre polynomial of degree $N$. On these nodes, the authors define the differentiation matrix denoted $\mathbf{D}_{N,L}$, and then along with the above weights, they introduce the discrete Laguerre transform and its inverse. The (forward) discrete transform will be fairly useful in order to transfer the information, i.e., values of solution in nodes, from the physical space into the coefficients of expansion in the frequency (phase) space. Thus, we get the LGRC expansion coefficients of the solution. The rough accuracy can be inferred from the order of magnitude of these coefficients.

Actually, the method casts the eigenvalue problem (3) into the following algebraic generalized eigenvalue problem

where

and

The vector $\mathbf{X}$ contains the LGRC nodes and again $\mathbf{I}_{N}$ and $\mathbf{Z}_{N}$ are respectively the identity and the zero matrices. The first row in (16), i.e., in matrices $\mathbf{A}$ and $\mathbf{B}$, stands for the boundary condition in $x_{1}=0$ from the problem (3) and the subsequent rows stand for the differential operators collocated in the nodes $x_{2}$ to $x_{N}$.

The algebraic eigenvalue problem (16) contrasts with the previously considered (6), (8), and (15). The problem (16) is non-singular and the other three are singular. Thus, its solution is straightforward.

An example of the pseudospectrum of Charney’s problem is provided in Fig. 5. We mention that we are interested in the leftmost part of this spectrum. The eigenvalue representing the one unstable mode is visible in this part.

For the scaling factor ranging in the interval [5,10], we have obtained the smoothest curves for the eigenvectors. The reason is that the Laguerre-Gauss-Radau quadrature nodes concentrate on the right of the origin, as the scaling factor increases. Actually, this is the region where a higher resolution is desired.

However, if we consider the phase speed $c$ for the one unstable mode obtained in [1] as the exact one, namely

we can find an optimal value of the scaling parameter. Namely, we compute for several values of scaling parameter $b$ the relative errors for $c$, real(c) and imag(c) with respect to $c_{ex}$ value from (17). Thus, we get Fig. 6. It justifies our choice for the scaling parameter $b$.

The dependence of the real and imaginary parts of the phase speed $c$ on the meteorological parameter $r$ in the Charney’s problem is depicted in Fig. 7 when $N=96$ and $b=5$. This dependence is in perfect accordance with the results in [5] and [1] as well as with the physical requirement that $c$ becomes real and zero when $r=1$.

Also, we have to point out that the numerical algorithm is stable in the sense that the first six digits of the eigenvalue of interest, its real and imaginary parts, remain practically unchanged for values of the cut off parameter $N$ in the interval [64, 128] and the values of scaling parameter in the interval [5, 10].

With the notations from Section 6, the LGRC casts the problem(4) into the algebraic eigenvalue one

where

We shortly report the results obtained using the MATLAB built-in function eig to solve the algebraic eigenvalue problem (18). The first five eigenvalues of (18) are reported in Table 4 and the first eigenvector is depicted in Fig. 8.

In [15], the authors show that $\lambda=-4.0$ is the only eigenvalue of the problem below the essential spectrum. The results reported in Table 4 confirm this analytical result. Moreover, the first eigenvector displayed in Fig. 8 clearly exhibits a boundary layer behavior at the left end. This behavior could be interpreted as a consequence of boundary condition in this point and has not been observed in [15].

Finally, we apply the Laguerre polynomial transform in order to quantify the accuracy of LGRC in solving this problem. Figure 9 shows the behavior of the coefficients of Laguerre expansion for the first three eigenvectors.

As in the case of the fist eigenvalue, the first eigenvector appears to be computed close to machine precision. Unfortunately, for the next eigenvectors this accuracy is dramatically lost.

The implementation of the eigenparameter dependent boundary conditions has proved feasible using the matrix concatenation facilities of MATLAB in both, finite and infinite cases as well as for collocation and for tau methods.

In the specific case of the Petterson-König’s rod eigenvalue problem, we verify numerically that the eigenproblem supplied with $\lambda$ dependent boundary conditions can be fairly well approximated, with respect to the completeness of the eigenfunctions, with a problem supplied with boundary conditions independent of the eigenparameter. Moreover, the presence of $\lambda$ dependent boundary conditions do not affect the accuracy of spectral methods, collocation or tau, based on the Chebyshev polynomials.

However, the ChT method performs somewhat better than ChC based on both rectangular or square differentiation at a "cheaper price", i.e., using fewer than half approximation orders. As the non-normality of matrices involved in ChT can not explain this fact (see Table 3), the only reason for this superiority must be the way in which the boundary conditions are enforced.

Rectangular differentiation neither eases the implementation nor improves the accuracy.

For the singularly perturbed problems on [ $0,\infty$ ), the numerical results provided by LGRC, look reliable and confirm the numerical outcomes obtained by collocation based on rational Chebyshev functions as well as the analytical results.

Acknowledgements The author is grateful to both anonymous referees for their careful reading of the paper and suggestions which resulted in considerable improvement. He is also fairly indebted to Prof. John Boyd for illuminating discussions on the Charney’s problem during our visit to NAG in Oxford with the occasion of Chebyshev Day, 2016. With the same opportunity, Prof. Nick Trefethen provided us hints on rectangular differentiation.

1. 1.

   Boyd, J.P.: Orthogonal rational functions on a semi-infinite interval. J. Comp. Phys. 70, 63-88 (1987)

   Report issue for preceding element
2. 2.

   Boyd, J.P.: Chebyshev domain truncation is inferior to fourier domain truncation for solving problems on an infinite interval. J. Sci. Comput. 3, 109-120 (1988)

   Report issue for preceding element
3. 3.

   Boyd, J.P., Rangan, C., Bucksbaum, P.H.: Pseudospectral methods on a semi-infinite interval with application to the hydrogen atom: a comparison of the mapped Fourier sine method with Laguerre series and rational Chebyshev expansions. J. Comp. Phys. 188, 56-74 (2003)

   Report issue for preceding element
4. 4.

   Boyd, J.P.: Five themes in Chebyshev spectral methods applied to the regularized Charney eigenproblem: extra numerical boundary conditions, a boundary-layer-resolving change of coordinate, parametrzing a curve which is singularat at an endpoint, extending the tau method to log-andpolynomials and finding the roots of a polynomial-and-log approximation. Comput. Math. Appl. 71, 1277-1241 (2016)

   Report issue for preceding element
5. 5.

   Branscome, L.E.: The charney baroclinic stability problem: approximate solutions and modal structures. J. Atmos. Sci., 1393-1409 (1983)

   Report issue for preceding element
6. 6.

   Driscoll, T.A., Hale, N., Trefethen, L.N.: Chebfun Guide. Pafnuty Publications, Oxford (2014)

   Report issue for preceding element
7. 7.

   Driscoll, T.A., Hale, N.: Rectangular spectral collocation. IMA J. Numer. Anal. doi:10.1093/imanum /dru062

   Report issue for preceding element
8. 8.

   Gheorghiu, C.I., Pop, I.S.: A modified Chevyshev-Tau method for a hydrodynamic stability problem. Approximation and Optimization. In: Proceedings of the International Conference on Approximation and Optimization (Romania)-ICAOR Cluj-Napoca, vol. II, pp. 119-126 (1996)

   Report issue for preceding element
9. 9.

   Gheorghiu, C.I., Rommes, J.: Application of the Jacobi-Davidson method to accurate analysis of singular linear hydrodynamic stability problems. Int. J. Numer. Method Fl. (2012)

   Report issue for preceding element
10. 10.

    Gheorghiu, C.I.: Spectral methods for non-standard eigenvalue problems. Fluid and Structural Mechanics and Beyond. Springer Cham Heidelberg New York Dondrecht London (2014)

    Report issue for preceding element
11. 11.

    Gheorghiu, C.I.: Stable spectral collocation solutions to a class of Benjamin Bona Mahony initial value problems. Appl. Math. Comput. 273, 1090-1099 (2016)

    Report issue for preceding element
12. 12.

    Gheorghiu, C.I.: Spectral collocation solutions to systems of boundary layer type. Numer. Algor. doi:10.1007/s11075-015-0083-6

    Report issue for preceding element
13. 13.

    Gottlieb, D., Orszag, S.A.: Numerical analysis of spectral methods: theory and applications. SIAM, Philadelphia (1977)

    Report issue for preceding element
14. 14.

    Goussis, D.A., Pearlstein, A.J.: Removal of infinite eigenvalues in the generalized matrix eigenvalue problems. J. Comput. Phys. 84, 242-246 (1998)

    Report issue for preceding element
15. 15.

    Maohzu, Z., Sun, J., Zettl, A.: The spectrum of singular Sturm-Liouville problems with eigenparameter dependent boundary conditions and its approximation. Results. Math. 63, 1311-1330 (2013)

    Report issue for preceding element
16. 16.

    van Noorden, T., Rommes, J.: Computing a partial generalized real Schur form using the JacobiDavidson method. Numer. Linear Algebra Appl. 14, 197-215 (2007)

    Report issue for preceding element
17. 17.

    Rommes, J.: Arnoldi and Jacobi-Davidson methods for generalized eigenvalue problems $Ax=\lambda Bx$ with $B$ singular. Math. of Comput. 77, 995-1015 (2008)

    Report issue for preceding element
18. 18.

    Shen, J., Tang, T., Wang, L.-L.: Spectral methods. Algorithms, Analysis and Applications Springer Heidelberg Dordrecht London New York (2011)

    Report issue for preceding element
19. 19.

    Trefethen, L.N.: Computation of pseudospectra. Acta Numer. 8, 247-295 (1999)

    Report issue for preceding element
20. 20.

    Tretter, Ch.: Nonselfadjoint spectral problems for linear pencils $N-\lambda P$ of ordinary differential operators with $\lambda$-linear boundary conditions: completeness results. Integr. Equat. Oper. Th. 26, 222-248 (1996)

    Report issue for preceding element
21. 21.

    Tretter, C.h.: A linearization for a class of $\lambda$ - nonlinear boundary eigenvalue problems. J. Math. Analysis Appl. 247, 331-355 (2000)

    Report issue for preceding element
22. 22.

    Tretter, Ch.: Boundary eigenvalue problems for differential equations $N\eta=\lambda P\eta$ with $\lambda$-polynomial boundary conditions. Integr. J. Diff. Eqns. 170, 408-471 (2001)

    Report issue for preceding element
23. 23.

    Xu, K., Hale, N.: Explicit construction of rectangular differential matrices. IMA J. Numer. Anal. doi:10.1093/imanum/drv013

    Report issue for preceding element
24. 24.

    Xu, K., Hale, N.: The Chebyshev points of the first kind. Appl. Numer. Math. 102, 17-30 (2016)

    Report issue for preceding element
25. 25.

    Wang, C.M., Wang, C.Y., Reddy, J.N.: Exact solutions for buckling of structural members. CRC Press (2005)

    Report issue for preceding element
26. 26.

    Weideman, J.A.C., Reddy, S.C.: A MATLAB differentiation matrix suite. ACM T. Math. Software 26, 465-519 (2000)

    Report issue for preceding element
27. 27.

    von Winckel, G.: Fast Chebyshev Transform, MathWorks ^®. File Exchanges (2005)

    Report issue for preceding element