Numerical approximations of Sturm–Liouville eigenvalues using Chebyshev polynomial expansions method

Abstract: In this paper, an efficient technique based on the Chebyshev polynomial expansions for computing the eigenvalues of secondand fourth-order Sturm– Liouville boundary value problems is proposed. This technique reduces the given Sturm–Liouville problem to an integral equation. The resulting integral equation is then transformed into the eigenvalue equation by calculating the integrals at the grid points using the Chebyshev expansions. Thus, the required eigenvalues of the given problem are obtained by solving this eigenvalue equation. The excellent performance of this scheme is illustrated through some numerical examples, and comparison with other methods is presented.


Introduction
Sturm-Liouville problems (SLPs) arise throughout engineering mathematics. They arise directly as eigenvalue problems in one space dimension. For example, they describe the vibrational modes of various systems, such as the vibrations of a string or the energy eigenfunctions of a quantum mechanical oscillator, in which case the eigenvalues correspond to the resonant frequencies of vibration or energy levels. They also commonly arise from linear PDEs in several space dimensions when the equations are separable in some coordinate system, such as cylindrical or spherical coordinates.
Mathematicians have studied SLPs for over 200 years. It has a highly developed theory and remains an active area of interest. Therefore, in this work, we consider numerical solutions of both second-and fourth-order SLPs described below.

PUBLIC INTEREST STATEMENT
"Sturm-Liouville problems" are famous boundary value problems that naturally arise when solving certain partial differential equation problems using a "separation of variables" method. Mathematicians have studied these problems for over 200 years. They have a highly developed theory and remain an active area of interest. In this work, a numerical approach is developed to obtain approximate eigenvalues of these problems.
A regular second-order SLP is a second-order linear, homogeneous, ordinary differential equation of the form together with the seperated (each with one boundary point) homogeneous boundary conditions where α i and β i , i = 1, 2, are the real constants. In Equation 1, λ is the unknown eigenvalue pertinent to the differential operator defined by the left-hand side of Equation 1, the interval (a, b) is finite, the functions p, p′, q, and w are in L 1 (a, b) with p(x) > 0 and w(x) > 0 for x ∈ (a, b). This is distinguished from the case when p or w vanishes at some point in the interval [a, b] or when the interval is of infinite length, in which case the problem is called a singular SLP. For a regular second-order SLP, it is found that all eigenvalues are real and non-negative. In addition, there are infinitely many eigenvalues that are discretely distributed and hence do not fill out any interval. More information on the mathematical theory of second-order SLPs may be found in Amrein, Hinz, and Pearson (2005) and Pryce (1993).
A non-singular fourth-order SLP consists of a fourth-order linear ODE of the form together with seperated, self-adjoint boundary conditions specified at both ends of the domain (a, b), two boundary conditions at the end x = a, and another two boundary conditions at the end x = b. Basically, there are three types of boundary conditions and their combinations commonly used with Equation 3 in applications. These boundary conditions are given as y = 0, y � = 0 for clamped end, y = 0, y �� = 0 for hinged (or simply supported) end, and y �� = 0, y �� = 0 for free end. The technical conditions for the problem to be non-singular are: the interval (a, b) is finite; the functions P(x), S(x), Q(x), W(x), and 1∕P(x) are in L 1 (a, b); and the essential infima of P(x) and W(x) are both positive. Under these assumptions, it is well known that the eigenvalues are bounded from below. They can be ordered as 0 ≤ 1 ≤ 2 ≤ ⋯ ≤ k ≤ ⋯, where λ k → ∞ as k → ∞, and where each eigenvalue has multiplicity at most 2. More information on the mathematical theory of fourth-order SLPs may be found in Greenberg (1991) and Greenberg and Marletta (1995).
In most cases, it is not possible to obtain the eigenvalues of the above problems analytically. However, there are various approximate methods as, for example, the weighted residual methods, the variational methods, and finite difference and finite elements methods. Extensive literature reviews for the numerical approximations of second-order and fourth-order SLPs can be found in Yücel (2006) and Yücel and Boubaker (2012), respectively. In order to obtain more efficient numerical results, several ways have been devised in the last years. Recently, El-gamel and Abd El-hady (2013) discussed and compared two useful schemes for second-order Sturm-Liouville eigenvalue problems: differential quadrature method (DQM) and collocation method with sinc functions. Saleh Taher, Malek, and Momeni-Masuleh (2013) proposed a new technique based on Chebyshev differentiation matrix for computing eigenvalues of a general class of regular fourth-order SLPs. Yuan, Ye, Xiao, Kennedy, and Williams (2014) introduced the exact dynamic stiffness vibration method for solving regular second-and fourth-order SLPs. Most recently, Amodio and Settanni (2015) discussed the solution of regular and singular SLPs by means of high-order finite difference schemes. They described a method to define a discrete problem and its numerical solution by means of linear algebra techniques.
In the literature, to the best of the author's knowledge, there is no study on the Chebyshev polynomial expansions (El-gendi's method) (El-gendi, 1969) applied for SLPs. On the other hand, this method is an efficient discretization technique for obtaining accurate numerical solutions of differential, integral, and integro-differential equations. El-gendi's method is based on expanding the unknown function in terms of Chebyshev polynomials. Based on the discrete ortogonality relationships of these polynomials, several methods for solving linear and non-linear ODEs (Fox, 1962) and integral differential equations (Elliott, 1963) were proposed at about the same time. They were found to have considerable advantage over the finite difference methods. Since then, these methods have become standard. They are now part of the larger family of spectral methods (Boyd, 2000). They consist in expanding the unknown function in a series of Chebyshev polynomials, trancating this series, and then substituting the approximation in the actual equation, and finally determining equations for the coefficients. However, El-gendi (1969) described a new method for the numerical solution of differential, integral, and integro-differential equations by computing directly the values of the functions rather than the Chebyshev coefficients. These two approaches are equivalent in the sense that if the function values at some grid points are known, the Chebyshev coefficents for the function can be directly computed. The method has been extended to solve initial value problems in time-dependent quantum field theory and second-order boundary value problems in fluid (Mihaila & Mihaila, 2002). Recently, it has been applied to solve the generalized Kuramoto-Sivashinsky equation (Khater & Temsah, 2008) and convection-diffusion equation (Temsah, 2009).
Our goal in this paper is to extend Chebyshev polynomial expansions, known as El-gendi's method (El-gendi, 1969), to deal with the second-and fourth-order SLPs given above. To achieve this goal, we first consider the numerical approximations of second-order SLP (1) with Drichlet boundary conditions Then, we consider numerical approximations of a simple form of Equation 3 with some specified end conditions to demonstrate directly the use of El-gendi's method on the fourth-order problems. In addition, we choose the fourth-order problems (3) to be squares of the second-order problems (1) for numerical illustrations. This is because if we have a second-order problem with Drichlet boundary conditions given above and the differential operator, then the corresponding fourth-order problem (3) has the differential operator where The boundary conditions, in this case, are The eigenvalues of the fourth-order problem are then simply the squares of the eigenvalues of the second-order problem which means that they are all non-negative (Greenberg & Marletta, 1995).
This paper is organized as follows. In Section 2, we first review the Chebyshev polynomial expansions method (El-gendy's method), then we give some new formulas for the matrix approximation of integrals which are necessary to deal with higher order problems. The method is then applied to second-order and a simple form of the fourth-order SLPs in Section 3. Numerical examples are discussed in Section 4, and some conclusions are drawn in Section 5.

The method of Chebyshev polynomial expansions
In this section, we describe the method of Chebyshev expansion by following closely the procedures outlined in El-gendi (1969). It is assumed that a function y(x), which is continuous and of bounded variation in the interval [-1, 1], can be approximated by where the summation symbol with double prime indicates that the terms with suffixes i = 0 and i = N are to be halved, T i (x) is the Chebyshev polynomials of the first kind of degree i ≤ N, and the coefficients a i are defined by In the equation above, the grid points x j , the so-called Chebyshev points, are given by The approximate formulae (Equation 9) is exact at x = x j given by Equation 11. It is known by the Chebyshev trancation theorem (Boyd, 2000) that the error in approximating y(x) by the sum of its first N terms is bounded by the sum of the absolute values of all the neglected coefficients. A detailed discussion of this theorem and the convergence theory for Chebyshev polynomials can be found in Boyd (2000).
Using the above approximation for the aforementioned function y(x), the integral at the grid points can be approximated as This can be achieved using modern computing tools like Maple and Matlab. For N = 4, the elements of the matrix D, obtained by Maple, are given below. The equation It can be seen that the approximation of the integral (12) based on the Chebyshev expansion (9) can be obtained without computing the Chebyshev coefficients (Equation 10). In other words, the righthand side of Equation 15 gives the approximate values of the integral (12) at the points (13), rather than the Chebyshev coefficients of the integral (10). However, the two approaches are equivalent in the sense that if we know the values of the integral at x k given by Equation 13, its Chebyshev coefficients can be directly computed by using a formula similar to Equation 10. The main advantage of using this approach is that for a certain value of N, the elements of matrix D can be eveluated once and for all.
The following approximations which are necessary for practical applications to be discussed later can be deduced from relation (Equation 15).
where d is an (N + 1) row vector whose elements are the last row of matrix D.
where ẽ is an (N + 1) row vector whose elements are the last row of matrix E.
where s is an (N + 1) row vector whose elements are the last row of matrix DD; D is an (N + 1) × (N + 1) matrix whose elements are given by where h is an (N + 1) row vector whose elements are the last row of matrix H.
Standard Chebyshev polynomials of the first kind T r (x) are valid over the interval [−1, 1]. However, this interval can be normalized using the change of variable x → (x + 1)∕ 2. In this case, the Chebyshev expansion of the integral will be in terms of the shifted Chebyshev polynomials of the first kind defined by which are valid over the interval [0, 1]. All the properties of T * r (x) can be deduced from those of T r (2x−1) (Fox & Parker, 1968). Therefore, when the range is [0, 1], we have the approximation where D = 1 2 D and the right-hand side defines the integral at the grid points The elements of the column y in Equation 25 are also defined at these points. Hence, for practical applications, the operators from Equations 18 to 23 can be similarly defined by replacing the matrix D by D and the points Equation 13 by Equation 26.
As can be seen from the above formulas, the range of x in the computational domain should either be [−1, 1] or [0, 1]. For practical applications, the physical domain is neither [−1, 1] nor [0, 1], but rather [a, b], then for this case, we can perform a coordinate transformation. This can be easily done by remembering that any finite range, a ≤ t ≤ b, can be transformed to the basic range −1 ≤ x ≤ 1 and/ or the range 0 ≤ x ≤ 1 with the change of variables

Applications to SLPs
In this section, we consider the numerical approximations of second-and fourth-order Sturm-Liouville eigenvalues using the matrix representations of the integrals obtained above.

Second-order problem
The general SLP (1) can be easily reduced to the so-called Liouville normal form (or equivalently the Schrödinger equation) We can now substitute the so obtained y′(−1) into Equation 32 to obtain the final integrated form Using the techniques developed in the previous section to calculate integrals, this integral equation can now be transformed into the following eigenvalue equation: where A and B are (N + 1) × (N + 1) matrices whose entries are given by A y =̂B y In the above equation, δ ij is the Kronecker delta, u i = x i + 1 ∕ 2 the elements of an (N + 1) column vector u, s j the elements of (N + 1) row vector s given by Equation 21 with =q, ê ij the elements of (N + 1) × (N + 1) matrix Ê given by Equation 20 with =q, e ij the elements of (N + 1) × (N + 1) matrix E, and ẽ j the elements of (N + 1) row vector ẽ given by Equation 19. It should be noted that proper definition of Equation 36 can be easily obtained when the transformed range is [0, 1]. It is important to note that the first and last columns of B have no effect (since multuplied by zero). Hence, the eigenvalue Equation 35 can be reduced to a (N − 1) × (N − 1) system by eliminating the first and last rows and columns of matrices A and B.
Equation 35 is a generalized eigenvalue problem. When B is non-singular, this equation is mathematically equivalent to B -1 A y =̂y, and when A is non-singular, it is equivalent to A -1 B y = 1 ∕̂ y. Thus, in theory, if one of the matrices A or B is known to be non-singular, the problem could be reduced to a standard eigenvalue problem. However, for this reduction to be satisfactory from the point of view of numerical stability, it is necessary not only that B (or A) should be non-singular, but that it should be well-conditioned with respect to inversion. There are many underlying computational routines that can be used for solving these problems. Once we have the computed ̂ values, they can be used to obtain the eigenvalues of the original problem (28) using

A simple case of the fourth-order problem
Although the numerical results of Equation 3 with some specified end conditions, obtained using the solutions of the second-order problem, will be given in the next section, we consider here a simple case of Equation 3 to demonstrate directly the use of El-gendi's method on the fourth-order problems. Therefore, we put P(x) = W(x) = 1, S(x) = Q(x) = 0, a = 0, and b = 1 in Equation 3. Then, we obtain the following simple equation, the steady-state Euler-Bernoulli equation for the deflection y(x) of a vibrating beam of length 1: As mentioned earlier, there are three types of boundary conditions and their combinations commonly used with this equation in applications. In this work, the formulations are given only for one of them. Other types of boundary conditions can be handled similarly. We assume that the clamped end is at x = 0 and the simply supported end is at x = 1, A coordinate transformation is not necessary for this case since the physical domain is already [0, 1]. By successive integration of Equation 37 from the lower limit 0 to an arbitrary value of x ∈ [0, 1], we obtain (37) where A and B are (N + 1) × (N + 1) matrices whose entries are given by Here, v i = x 3 i − 3x 2 i ∕ 2 are the elements of an (N + 1) column vector v, z i = x 2 i − x 3 i ∕ 4 the elements of an (N + 1) column vector z, h ij the elements of (N + 1) × (N + 1) matrix H given by Equation 22, h j the elements of (N + 1) row vector h given by Equation 23, and ẽ j the elements of (N + 1) row vector ẽ given by Equation 19. For the same reason given in Section 3.1, the eigenvalue equation (44) can be reduced to a (N − 1) × (N − 1) system. A computational routine can now be used to solve Equation 44 for obtaining eigenvalues of the problem.

Numerical results
We present here the results obtained by applying our algorithms to the numerical solution of some second-and fourth-order SLPs in practice. First two problems which have exact solutions are chosen to demonstrate the efficiency and accuracy of the method developed in this work. We will use the relative error ε k which is defined as to measere the performance of the method. In the last problem, we choose the fourth-order problem (3) to be square of the second-order problem (1) as discussed in Section 1. Numerical calculations are carried out using Maple and Matlab. Problem 1. Our first example for the second-order SLP is given by which arises when solving the one-dimensional wave equation u tt −u xx = 0 (the standard model for an oscillating string of length 1 with fixed endpoints) using a "seperation of variables" method. The exact eigenvalues of this problem are k = k 2 2 , k = 1, 2, 3, … . Table 1 lists the relative errors of the obtained results with different number of the Chebyshev points N. It is observed from this table that in order to have good approximations to the first kth eigenvalues, at least 2k Chebyshev points have to be used. It is also observed that the computed values for the lower eigenvalues have a better accuracy than those for the higher eigenvalues. As the number of grid points further increased to above 2k, the accuracy of the obtained results especially for the higher eigenvalues can be further improved as shown in Table 1.
Problem 2. As a second example, we consider the following fourth-order problem which corresponds to the case already discussed in Section 3.2. The exact eigenvalues of problem (48) can be found by solving tanh 1∕ 4 − tan 1∕ 4 = 0.
In Table 2, the relative errors of the first 10 eigenvalues of this problem for several values of the number of Chebyshev points are given. It can be seen that the method developed in this work is very efficient for finding the eigenvalues of the given fourth-order problem.
Problem 3. In this example, we consider the Coffey-Evans equation (Greenberg and Marletta, 1995) which arises in a model of the coupling between dipoles in polarizable liquids such as used in liquid crystal displays (Pryce, 1986). The parameter μ, typically in the range 0-50, measures the strength of coupling.
The lowest eigenvalue λ 0 of this problem is close to zero. Figure 1 shows the computed eigenvalues (λ 1 -λ 19 ) of this problem for several values of the parameter μ, μ = 10, μ = 20, and μ = 30. In Table 3, we compare our results with the available results of Pryce (1986) for λ k where k ≤ 4. It can be observed that our results are in good agreement with those obtained by Pryce (1986). Pryce's paper was concerned with Prüfer shooting methods for finding eigenvalues and corresponding eigenfunctions of the classical SLP (1) with appropriate boundary conditions given by Equation 2. Such methods use a transformation (48) y (4) (x) = y(x), y(0) = y � (0) = 0, y(1) = y �� (1) = 0 (49) −y �� + 2 sin 2 2x − 2 cos 2x y = y, y(− ∕ 2) = y( ∕ 2) = 0 to phase amplitude variables (θ, r) in the (y′, y)-plane (a Prüfer transformation) and solve the resulting equations by a shooting method using an ODE initial value code. Pryce studied the error control of such methods with specific reference to two published codes based on this technique.
The corresponding fourth-order problem for the Equation 49 can be obtained by using Equations 5-9: with the boundary conditions (50) y (4) − 2 2 sin 2 2x − 2 cos 2x y � � + 2 sin 2 2x − 2 cos 2x 2 − 8 2 cos 2 2x − 2 sin 2 2x + cos 2x y = y  As mentioned earlier in Section 1, the eigenvalues of this fourth-order problem are simply the squares of the eigenvalues of the second-order problem (49). The results are shown in Table 4 where the comparisons are made with the available data obtained by Greenberg and Marletta (1995) for μ = 10 and μ = 30. Greenberg and Marletta (1995) developed a shooting method to approximate the eigenvalues and eigenfunctions of the fourth-order problem (3). The problem was first reduced to a Hamiltonian system. They developed an efficient approach to solve this system based on the approximation of the coefficients in the differential equation and a suitable zero-counting algorithm. The zeros which they counted were not zeros of a solution of the fourth-order equation, but nullities (rank deficiencies) of a certain 2 × 2 matrix formed from solutions of the fourth-order equation.

Conclusions
In this paper, we present a method to obtain good approximations to the eigenvalues of second-and fourth-order SLPs defined on finite domains based on a spectral method known as El-gendi's method. This method provides a robust algorithm of computing the integral of a non-singular function defined on a finite domain. Therefore, the method presented in this work converts the given SLP into an integral equation. The resulting integral equation is then transformed into the eigenvalue equation by calculating the integrals at the grid points using the Chebyshev expansions. The required eigenvalues are then obtained by solving this eigenvalue equation. The method is quite general and has some special advantages which were discussed in detail in Boyd (2000). The advantages of El-gendi's method over finite difference methods were also discussed in Boyd (2000) and so will not be repeated here.
Through test examples which have exact solutions, it was found that in order to obtain accurate numerical results for the first kth eigenvalues, at least 2k Chebyshev points have to be used. Comparison with other published works in the literature showed that the method produces highly accurate results for the eigenvalues of the SLPs considered in this work. It may be concluded that the presented method is very powerful and efficient in finding the approximate solutions of SLPs arising in science and engineering.