Eigenvalue Computation of Regular 4th Order Sturm-Liouville Problems

Ahmad Alalyani

Al Baha University, Al Baha, KSA.

DOI: 10.4236/am.2019.109056                

Abstract

In this paper we present and test a numerical method for computing eigenvalues of 4th order Sturm-Liouville (SL) differential operators on finite intervals with regular boundary conditions. This method is a 4th order shooting method based on Magnus expansions (MG4) which use MG4 shooting as the integrator. This method is similar to the SLEUTH (Sturm-Liouville Eigenvalues Using Theta Matrices) method of Greenberg and Marletta which uses the 2nd order Pruess method (also known as the MG2 shooting method) for the integrator. This method often achieves near machine precision accuracies, and some comparisons of its performance against the well-known SLEUTH software package are presented.

Keywords

4th Order Sturm-Liouville Problem, Magnus Methods (MG4), Regular Boundary Conditions

1. Introduction

In this paper1 we consider the self-adjoint differential operators which arise from the 4th order differential equation

$L\left(y\right)=y^{\left(4\right)}\left(x\right)-{\left(s\left(x\right)y^{\prime }\left(x\right)\right)}^{\prime }+q\left(x\right)y\left(x\right)=\lambda y\left(x\right)\mathrm{,}a\le x\le b$ (1)

when separated, self-adjoint boundary conditions are imposed at each of the two regular endpoints $x=a$ and $x=b$.

We make the assumptions:

*2010 Mathematics Subject Classication 34L15; 34L16; 65L15; 70J50.

¹The content of this paper is related to my Ph.D. dissertation [1] .

1) $q\left(x\right)$ is continuous on $\left[a\mathrm{,}b\right]$.

2) $s\left(x\right)$ and $s^{\prime }\left(x\right)$ are continuous on $\left[a\mathrm{,}b\right]$.

3) $-\infty <a<b<\infty$.

Under these assumptions both endpoints a and b are regular endpoints. The most general separated, self-adjoint boundary conditions which can be imposed at $x=a$ and $x=b$ are

$\left(\begin{array}{c}{l}_1\left(y\right)\\ l_2\left(y\right)\end{array}\right)=A_1\left(\begin{array}{c}y\left(a\right)\\ y^{″}\left(a\right)\end{array}\right)+A_2\left(\begin{array}{c}-y^{‴}\left(a\right)+s\left(a\right)y^{\prime }\left(a\right)\\ -y^{\prime }\left(a\right)\end{array}\right)=\left(\begin{array}{c}0\\ 0\end{array}\right)$ (2)

and

$\left(\begin{array}{c}{r}_1\left(y\right)\\ r_2\left(y\right)\end{array}\right)=B_1\left(\begin{array}{c}y\left(b\right)\\ y^{″}\left(b\right)\end{array}\right)+B_2\left(\begin{array}{c}-y^{‴}\left(b\right)+s\left(b\right)y^{\prime }\left(b\right)\\ -y^{\prime }\left(b\right)\end{array}\right)=\left(\begin{array}{c}0\\ 0\end{array}\right)$ (3)

where $A_1\mathrm{,}{A}_2$ and $B_1\mathrm{,}{B}_2$ are any choice of real, $2\times 2$ matrices satisfying the properties

$A_1{A}_2^{\text{T}}-A_2{A}_1^{\text{T}}=0$ (4)

$A_1{A}_1^{\text{T}}+A_2{A}_2^{\text{T}}=I_2.$ (5)

and

$B_1{B}_2^{\text{T}}-B_2{B}_1^{\text{T}}=0$ (6)

$B_1{B}_1^{\text{T}}+B_2{B}_2^{\text{T}}=I_2.$ (7)

The above boundary conditions can be shown to be equivalent to the general forms of boundary conditions used by Everitt [2] (in his PhD dissertation on 4th order Sturm-Liouville problems), Fulton [3] and many others.

The domain of the maximal operator $L_1$ associated with the Equation (1) on the closed interval $\left[a\mathrm{,}b\right]$ is

$D\left(L_1\right)=\left\{f\in L_2\left(a,b\right):f,f^{\prime },f^{″},f^{‴}\in A{C}_{loc}\left(a,b\right),Lf\in L_2\left(a,b\right)\right\},$ (8)

where $A{C}_{loc}$ is the space of functions which are absolutely continuous on compact subsets of $\left(a\mathrm{,}b\right)$. The self-adjoint operators associated with Equation (1) are then obtained by restricting $D\left(L_1\right)$ by two boundary conditions at the left endpoint and two boundary conditions at the right endpoint as in (2) and (3), namely

$D\left(L_{A_1,A_2,B_1,B_2}\right)=\left\{f\in D\left(L_1\right):l_1\left(f\right)=0,l_2\left(f\right)=0,r_1\left(f\right)=0,r_2\left(f\right)=0\right\}$ (9)

The Green’s formula for the 4th order equation is

${\displaystyle {\int }_a^b}\left(zLy-yLz\right)\text{d}x={{\left[y\mathrm{,}z\right]}_{\left(x\right)}|}_a^b\mathrm{,}$ (10)

where the bilinear concomitant is defined as

${\left[y,z\right]}_{\left(x\right)}:=\left|\begin{array}{cc}{y}^{‴}\left(x\right)& z^{‴}\left(x\right)\\ y\left(x\right)& z\left(x\right)\end{array}\right|-\left|\begin{array}{cc}{y}^{″}\left(x\right)& z^{″}\left(x\right)\\ y^{\prime }\left(x\right)& z^{\prime }\left(x\right)\end{array}\right|+s\left(x\right)\left|\begin{array}{cc}y\left(x\right)& z\left(x\right)\\ y^{\prime }\left(x\right)& z^{\prime }\left(x\right)\end{array}\right|.$ (11)

Using this definition, and the boundary conditions (2) and (3), it can be shown that the operators $L_{A_1,A_2,B_1,B_2}$ on $L_2\left(a\mathrm{,}b\right)$ are symmetric; that is, for all $f\mathrm{,}g\in D\left(L_1\right)$ :

$\left(f\mathrm{,}Lg\right)-\left(Lg\mathrm{,}f\right)={\displaystyle {\int }_a^b}f\cdot \left(Lg\right)-g\cdot \left(Lf\right)\text{d}x=0.$ (12)

This paper is devoted to the 4th order shooting method based on Magnus Expansions (MG4) for computation of eigenvalues of 4th order SL problems of the type (1), (2), (3) having regular endpoints.

For 2nd order SL problems, on both regular and singular intervals, there are several well developed software packages for eigenvalue and eigenfunction computations: SLEIGN [4] , SLEIGN 2 [5] , SLEDGE [6] [7] , SLO2FM [8] [9] [10] , MATSLISE [11] [12] [13] . The best source of information where their capabilities and general performance is discussed, is the book, Numerical Solution of Sturm-Liouville Problems, of John Pryce [8] .

For the 4th order SL equation, the only reliable software package for eigenvalue and eigenfunction computations is ACM Algorithm 775: SUBROUTINE SLEUTH, produced by L. Greenberg and M. Marletta in 1997 [14] . This package restricts attention to problems on finite intervals with regular endpoints; at this writing there is still no readily available software for singular endpoints. The SLEUTH code handles eigenvalues and eigenfunctions only for SL problems with regular endpoints, and it is capable of handling a wide variety of possible boundary conditions. The Greenberg-Marletta algorithm is based on an integration scheme using piecewise trigonometric hyperbolic splines (the Pruess method), also known as the MG2 shooting method. This was also the underlying integration scheme in the SLEDGE package. The SLEUTH code is based on using formulas of Greenberg [14] [15] [16] for the number of eigenvalues less than $\lambda$.

Since eigenvalue/eigenfunction calculations for the 4th order equation have been tackled by many other methods we give here a brief overview of some of the existing competitive methods. The most prominent approaches to date, and those which continue to receive much attention, are as follows:

1) Extended Sampling Method (ESM) which relies on the classical Whittaker-Shannon-Kotelnikov sampling theorem [17] (also used for 2nd order SL problems).

2) Fliess Series Method [18] [19] [20] which represents the solution of an IVP for the 4th order Equation (1) in terms of iterated integrals involving the coefficient functions, $s\left(x\right)$ and $q\left(x\right)$.

3) Chebyshev Method [21] - [26] which approximates solutions of the 4th order Equation (1) by Chebyshev polynomials.

4) Boubaker Polynomials Expansion Scheme (BPES) [27] which approximates solutions of of the 4th order Equation (1) using Boubaker polynomials and utilizes a Differential Quadrature Method (DQM).

5) Spectral Parameter Power Series (SPPS) Method [28] [29] [30] which expands the solutions of the SL equation (both second and 4th order) in a convergent Taylor expansion in the eigenparameter $\lambda$. This has recently proven to have much advantages both for theoretical problems and numerical computations for SL equations.

Selection of Test Problems

To investigate the performance of the method, we make the following selection of test problems. These problems are the square of a 2nd order SL problem.

1) The square of the 2nd order Bessel equation

$Ly=y^{\left(4\right)}-{\left(s\left(x\right)y^{\prime }\right)}^{\prime }+q\left(x\right)y=\lambda y,x\in \left[1,5\right]$ (13)

where

$s\left(x\right)=\frac{-2}{4x^2},$ (14)

and

$q\left(x\right)=\frac{1}{16x^4}+\frac{3}{2x^4}.$ (15)

2) The square of the 2nd order Modified Harmonic Oscillator equation

$Ly=y^{\left(4\right)}-{\left(s\left(x\right)y^{\prime }\right)}^{\prime }+q\left(x\right)y=\lambda y,x\in \left[1,5\right]$ (16)

where

$s\left(x\right)=2\left(x^2+x^4\right),$ (17)

and

$q\left(x\right)={\left(x^2+x^4\right)}^2-\left(2+12x^2\right)\mathrm{.}$ (18)

3) The square of the 2nd order equation

$Ly=y^{\left(4\right)}-{\left(s\left(x\right)y^{\prime }\right)}^{\prime }+q\left(x\right)y=\lambda y,x\in \left[0,\pi \right]$ (19)

where

$s\left(x\right)=2\left(\cos \left(x\right)+2\cos \left(2x\right)+3\cos \left(3x\right)\right),$ (20)

and

$q\left(x\right)={\left(\cos \left(x\right)+2\cos \left(2x\right)+3\cos \left(3x\right)\right)}^2-\left(-\cos \left(x\right)-8\cos \left(2x\right)-27\cos \left(3x\right)\right).$ (21)

4) The square of the 2nd order Coffey-Evan equation with $\beta =10$

$Ly=y^{\left(4\right)}-{\left(s\left(x\right)y^{\prime }\right)}^{\prime }+q\left(x\right)y=\lambda y,x\in \left[\frac{-\pi }{2},\frac{\pi }{2}\right]$ (22)

where

$s\left(x\right)=2\left({\beta }^2{\sin }^2\left(2x\right)-2\beta \cos \left(2x\right)\right),$ (23)

and

$q\left(x\right)={\left({\beta }^2{\sin }^2\left(2x\right)-2\beta \cos \left(2x\right)\right)}^2-\left(8{\beta }^2\cos \left(4x\right)+8\beta \cos \left(2x\right)\right).$ (24)

5) The square of the 2nd order Legendre equation

$Ly=y^{\left(4\right)}-{\left(s\left(x\right)y^{\prime }\right)}^{\prime }+q\left(x\right)y=\lambda y,x\in \left[0,\frac{\pi }{4}\right]$ (25)

where

$s\left(x\right)=\frac{1}{2}{\sec }^2\left(x\right),$ (26)

and

$q\left(x\right)=\left(\frac{{\sec }^4\left(x\right)}{16}\right)-\left({\sec }^2\left(x\right){\tan }^2\left(x\right)+\frac{{\sec }^4\left(x\right)}{2}\right).$ (27)

Problem 5, the Legendre squared equation, arises from changes of variables to the non-LNF form discussed in [31] .

2. The MG4 Shooting Method Associated with the 4th Order Sturm-Liouville Equation

In this section we describe an implementation of the MG4 shooting technique for the 4th order SL Equation (1) on regular intervals with $s\left(x\right)\mathrm{,}q\left(x\right)$ continuous. The Equation (1) can be converted to the 1st order system (Atkinson [32] , pp. 323-324) $Y^{\prime }\left(x\mathrm{,}\lambda \right)=A\left(x\right)\cdot Y\left(x\mathrm{,}\lambda \right)$, where

$Y^{\prime }\left(x\mathrm{,}\lambda \right)=\left(\begin{array}{cccc}0& 0& 0& -1\\ 0& 0& -1& -s\left(x\right)\\ q\left(x\right)-\lambda & 0& 0& 0\\ 0& -1& 0& 0\end{array}\right)Y\left(x\mathrm{,}\lambda \right)=A\left(x\right)\cdot Y\left(x\mathrm{,}\lambda \right)\mathrm{,}$ (1)

and

$Y\left(x\mathrm{,}\lambda \right)=\left(\begin{array}{c}y\left(x\mathrm{,}\lambda \right)\\ y^{″}\left(x\mathrm{,}\lambda \right)\\ -y^{‴}\left(x\mathrm{,}\lambda \right)+s\left(x\right)y^{\prime }\left(x\mathrm{,}\lambda \right)\\ -y^{\prime }\left(x\mathrm{,}\lambda \right)\end{array}\right)\mathrm{.}$ (2)

Remark 2.1 Currently the most reliable software package for eigenvalues and eigenfunctions of the 4th order Sturm-Liouville equation with regular endpoints is ACM Algorithm 775: SUBROUTINE SLEUTH, produced by L. Greenberg and M. Marletta in 1997 [14] , which is available from NETLIB at ORNL. The SLEUTH (Sturm-Liouville Eigenvalues Using Theta Matrices) code employs an MG2 approximation for the solution $Y\left(x\right)={\text{e}}^{\sigma \left(x\right)}{Y}_0$, (see [33] , p. 283), on each mesh interval of a Hamiltonian system similar to the system (1) (the order of the derivatives in the above $Y\left(x\mathrm{,}\lambda \right)$ vector being slightly different). The SLEUTH code is based on using formulas of Greenberg [15] [16] [34] [35] for the number of eigenvalues less than $\lambda$, so the eigenvalue algorithm is quite different than the method we are proposing here for eigenvalue computation for (1) using (1).

For the IVP, $Y^{\prime }\left(x\mathrm{,}\lambda \right)=A\left(x\right)\cdot Y\left(x\mathrm{,}\lambda \right)$, $Y\left(0\right)=I$, where A is a constant matrix is also basic to Magnus methods for (1). We introduce the following lemma and the definitions of the Lie-group and Lie-algebra (see [36] , Prop. 7.2.3 and Def. 7.2.4).

Definition 2.1 $SL\left(\mathrm{4,}ℝ\right)$ is the Lie-group and defined as:

$SL\left(\mathrm{4,}ℝ\right)\mathrm{<mo>\; :\; </mo>}=\left\{A\mathrm{<mo>\; |\; </mo>}A\text{is}4\times 4\text{matrixwithrealentriesand}\text{det}\left(A\right)=1\right\}\mathrm{.}$

Definition 2.2 $sl\left(\mathrm{4,}ℝ\right)$ is the Lie-algebra and defined as:

$sl\left(\mathrm{4,}ℝ\right)\mathrm{<mo>\; :\; </mo>}=\left\{A\mathrm{<mo>\; |\; </mo>}A\text{is}4\times 4\text{matrixwithrealentriesand}\text{tr}\left(A\right)=0\right\}\mathrm{.}$

Lemma 2.1 If $X\in sl\left(\mathrm{4,}ℝ\right)$, then $\exp \left(X\right)\in SL\left(4,ℝ\right)$, i.e.

$\det \left(\exp \left(X\right)\right)=1.$ (3)

Remark 2.2 For any constant matrix $X\in sl\left(\mathrm{4,}ℝ\right)$, it follows from this lemma that the solution

$Y\left(t\mathrm{,}\lambda \right)={\text{e}}^{Xt}$ (4)

of the IVP,

$Y^{\prime }\left(t,\lambda \right)=X\left(t\right)\cdot Y\left(t,\lambda \right),$ (5)

$Y\left(0\right)=I,$

lies in the Lie Group $SL\left(\mathrm{4,}ℝ\right)$.

The Magnus methods originate (see [37] ) with the observation that an analytical solution of (1) with initial condition $Y\left(0\right)=Y_0$ can be written as, (see [33] , p. 283),

$Y\left(x\right)={\text{e}}^{\sigma \left(x\right)}{Y}_0$ (6)

where

$\begin{array}{c}\sigma \left(x\right)={\displaystyle {\int }_0^x}A\left(k\right)\text{d}k+\frac{1}{2}{\displaystyle {\int }_0^x}\left[A\left(k\right)\mathrm{,}{\displaystyle {\int }_0^k}A\left(\xi \right)\text{d}\xi \right]\text{d}k\\ +\frac{1}{4}{\displaystyle {\int }_0^x}\left[A\left(k\right)\mathrm{,}{\displaystyle {\int }_0^k}\left[A\left(\xi \right)\mathrm{,}{\displaystyle {\int }_0^{\xi }}A\left(\eta \right)\text{d}\eta \right]\text{d}\xi \right]\text{d}k\\ +\frac{1}{12}{\displaystyle {\int }_0^x}\left[\left[A\left(k\right)\mathrm{,}{\displaystyle {\int }_0^k}A\left(\xi \right)\text{d}\xi \right]\mathrm{,}{\displaystyle {\int }_0^k}A\left(\eta \right)\text{d}\eta \right]\text{d}k+\cdots \end{array}$ (7)

and where the square brackets denote the matrix commutator and are defined as:

$\left[A,B\right]:=A\cdot B-B\cdot A$ (8)

The MG4 method is a well known 4th order method obtained by truncation of the above Magnus series, together with evaluation of the A matrix in (1) at two gaussian points $A_1$ and $A_2$ :

For the Hamiltonian system (1), we put

$A_1:=A\left(x_n+h\left(\frac{\sqrt{3}-1}{2\sqrt{3}}\right)\right)=A\left(x^{\prime }\right),$ (9)

$A_2:=A\left(x_n+h\left(\frac{\sqrt{3}+1}{2\sqrt{3}}\right)\right)=A\left(x^{″}\right),$ (10)

(meaning that $q\left(x\right)$ and $s\left(x\right)$ in (1) are to be evaluated at $x=x^{\prime }$ and $x=x^{″}$ in the nth mesh interval $\left[x_n\mathrm{,}{x}_{n+1}\right]$ ). Then the MG4 method of Iserles and Norsett ( [38] , p. 1012), and ( [33] , p. 288) for (1) takes the form (for a fixed step size h)

$Y_{n+1}=\exp \left(h\tilde{A}\right)\cdot Y_n,$ (11)

where the transfer matrix for passing from $x_n$ to $x_{n+1}$ is $M:=\exp \left(h\tilde{A}\right)$ with

$\begin{array}{l}\tilde{A}:=\frac{1}{2}\left(A_1+A_2\right)+\frac{1}{4\sqrt{3}}\left[A_2,A_1\right]h\\ =\frac{1}{2}\left(A\left(x^{\prime }\right)+A\left(x^{″}\right)\right)+\frac{1}{4\sqrt{3}}h\left(A\left(x^{″}\right)A\left(x^{\prime }\right)-A\left(x^{\prime }\right)A\left(x^{″}\right)\right)\\ =\left(\begin{array}{cccc}0& 0& 0& -1\\ \frac{-\sqrt{3}h\left(q_1-q_2\right)}{12}& \frac{-\sqrt{3}h\left(s_1-s_2\right)}{12}& -1& \frac{-s_1-s_2}{2}\\ \frac{q_1+q_2}{2}-\lambda & 0& 0& \frac{\sqrt{3}h\left(q_1-q_2\right)}{12}\\ 0& -1& 0& \frac{\sqrt{3}h\left(s_1-s_2\right)}{12}\end{array}\right),\end{array}$ (12)

where the square bracket $\left[A_2\mathrm{,}{A}_1\right]$ denotes the matrix commutator and is defined as:

$\left[A_2,A_1\right]=A_2{A}_1-A_1{A}_2,$ (13)

and

$q_1=q\left(x^{\prime }\right),q_2=q\left(x^{″}\right),s_1=s\left(x^{\prime }\right),s_2=s\left(x^{″}\right).$

The four eigenvalues of $\tilde{A}$ are

${\lambda }_1=-{\left[A-\frac{\sqrt{B}}{2}\right]}^{\frac{1}{2}},$ (14)

${\lambda }_2=-{\left[A+\frac{\sqrt{B}}{2}\right]}^{\frac{1}{2}},$ (15)

${\lambda }_3=-{\lambda }_1,$ (16)

${\lambda }_4=-{\lambda }_2,$ (17)

where

$A=\frac{s_1+s_2}{4}+\frac{h^2\left(s_1^2+s_2^2\right)}{96}-\frac{h^2{s}_1{s}_2}{48},$ (18)

and

$\begin{array}{c}B=4\lambda +\frac{h^4\left(s_1^4+s_2^4\right)}{2304}-\frac{h^4\left(s_1^3{s}_2+s_1{s}_2^3\right)}{576}+\frac{h^4{s}_1^2{s}_2^2}{384}+\frac{h^2\left(s_1^3+s_2^3-s_1^2{s}_2-s_1{s}_2^2\right)}{48}\\ +\frac{s_1^2+s_2^2}{4}+\frac{s_1{s}_2}{2}-2\left(q_1+q_2\right).\end{array}$ (19)

Eigenvalues of $\tilde{A}$ matrix:

Let us define

$D:=A-\frac{\sqrt{B}}{2},$ (20)

and

$E:=A+\frac{\sqrt{B}}{2}.$ (21)

Then the following four cases of eigenvalues of $\tilde{A}$ arise, involving both complex and real eigenvalues:

Case 1: $B>0$ and ( $D>0,E>0$ ).

${\lambda }_1=-{\left|D\right|}^{\frac{1}{2}},$ (22)

${\lambda }_2=-{\left|E\right|}^{\frac{1}{2}},$ (23)

${\lambda }_3={\left|D\right|}^{\frac{1}{2}},$ (24)

${\lambda }_4={\left|E\right|}^{\frac{1}{2}}.$ (25)

Case 2: $B>0$ and ( $D<0,E<0$ ).

${\lambda }_1=-i{\left|D\right|}^{\frac{1}{2}},$ (26)

${\lambda }_2=-i{\left|E\right|}^{\frac{1}{2}},$ (27)

${\lambda }_3=i{\left|D\right|}^{\frac{1}{2}},$ (28)

${\lambda }_4=i{\left|E\right|}^{\frac{1}{2}}.$ (29)

Case 3: $B>0$ and ( $D<0,E>0$, where $D<E$ ).

${\lambda }_1=-i{\left|D\right|}^{\frac{1}{2}},$ (30)

${\lambda }_2=-{\left|E\right|}^{\frac{1}{2}},$ (31)

${\lambda }_3=i{\left|D\right|}^{\frac{1}{2}},$ (32)

${\lambda }_4={\left|E\right|}^{\frac{1}{2}}.$ (33)

Case 4: $B<0$.

${\lambda }_1=-{\left[A-\frac{i\sqrt{\left|B\right|}}{2}\right]}^{\frac{1}{2}},$ (34)

${\lambda }_2=-{\left[A+\frac{i\sqrt{\left|B\right|}}{2}\right]}^{\frac{1}{2}},$ (35)

${\lambda }_3=-{\lambda }_1,$ (36)

${\lambda }_4=-{\lambda }_2.$ (37)

Remark 2.3 It follows from (12) that

$trace\left(\left(x-x_n\right)\tilde{A}\right)=0,$ (38)

so that $\left(x-x_n\right)\tilde{A}\in sl\left(\mathrm{4,}ℝ\right)$. Also we observe that on diagonalization we have (in all the above 4 cases),

$\begin{array}{c}\det \left(\exp \left[\left(x-x_n\right)\tilde{A}\right]\right)=\det \left(P\right)\cdot \det \left(\exp \left[\left(x-x_n\right)D\right]\right)\cdot \det \left(P^{-1}\right)\\ =\exp \left(\left(x-x_n\right){\displaystyle \underset{j=1}{\overset{4}{\sum }}}{\lambda }_j\right)=1,\end{array}$ (39)

where $\left\{{\lambda }_j\right\},j=1,2,3,4$, are the eigenvalues of $\tilde{A}$. Hence on each mesh interval,

$Y\left(x,\lambda \right)=\exp \left[\left(x-x_n\right)\tilde{A}\right]\cdot Y_n\left(x,\lambda \right)$ (40)

is a solution of

$Y^{\prime }\left(x,\lambda \right)=\tilde{A}\left(x\right)\cdot Y\left(x,\lambda \right)$ (41)

which remains in the Lie Group, $SL\left(\mathrm{4,}ℝ\right)$.

We consider the SL problem for Equation (1) with the following choices of Dirichlet boundary conditions at the left and right endpoints (compare (2) and (3)).

$L\left(y\right)=y^{\left(4\right)}\left(x\right)-{\left(s\left(x\right)y^{\prime }\left(x\right)\right)}^{\prime }+q\left(x\right)y\left(x\right)=\lambda y\left(x\right),a\le x\le b,$ (42)

$A_1\left(\begin{array}{c}y\left(a\right)\\ y^{″}\left(a\right)\end{array}\right)+A_2\left(\begin{array}{c}-y^{‴}\left(a\right)+s\left(a\right)y^{\prime }\left(a\right)\\ -y^{\prime }\left(a\right)\end{array}\right)=\left(\begin{array}{c}y\left(a\right)\\ y^{″}\left(a\right)\end{array}\right)=\left(\begin{array}{c}{l}_1\left(y\right)\\ l_2\left(y\right)\end{array}\right)=\left(\begin{array}{c}0\\ 0\end{array}\right),$ (43)

$B_1\left(\begin{array}{c}y\left(b\right)\\ y^{″}\left(b\right)\end{array}\right)+B_2\left(\begin{array}{c}-y^{‴}\left(b\right)+s\left(b\right)y^{\prime }\left(b\right)\\ -y^{\prime }\left(b\right)\end{array}\right)=\left(\begin{array}{c}y\left(b\right)\\ y^{″}\left(b\right)\end{array}\right)=\left(\begin{array}{c}{r}_1\left(y\right)\\ r_2\left(y\right)\end{array}\right)=\left(\begin{array}{c}0\\ 0\end{array}\right).$ (44)

where

$A_1=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right),$ (45)

$A_2=\left(\begin{array}{cc}0& 0\\ 0& 0\end{array}\right)$ (46)

and

$B_1=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right),$ (47)

$B_2=\left(\begin{array}{cc}0& 0\\ 0& 0\end{array}\right).$ (48)

The left boundary conditions are implemented by fixing initial conditions for two solutions $Y_1\left(x\mathrm{,}\lambda \right)$ and $Y_2\left(x\mathrm{,}\lambda \right)$ of (1) at $x=a$, namely we define solutions $Y_1\left(x\mathrm{,}\lambda \right)$ and $Y_2\left(x\mathrm{,}\lambda \right)$ at $x=a$ by requiring

$\begin{array}{c}\left(Y_1\left(a,\lambda \right),Y_2\left(a,\lambda \right)\right)=\left(\begin{array}{cc}{y}_1\left(a,\lambda \right)& y_2\left(a,\lambda \right)\\ {y^{″}}_1\left(a,\lambda \right)& {y^{″}}_2\left(a,\lambda \right)\\ -{y^{‴}}_1\left(a,\lambda \right)+s\left(a\right){y^{\prime }}_1\left(a,\lambda \right)& -{y^{‴}}_2\left(a,\lambda \right)+s\left(a\right){y^{\prime }}_2\left(a,\lambda \right)\\ -{y^{\prime }}_1\left(a,\lambda \right)& -{y^{\prime }}_2\left(a,\lambda \right)\end{array}\right)\\ =\left(\begin{array}{cc}0& 0\\ 0& 0\\ 1& 0\\ 0& 1\end{array}\right)=\left(\begin{array}{c}0\\ I_2\end{array}\right).\end{array}$ (49)

Then the corresponding solutions $\left\{y_1\left(x\mathrm{,}\lambda \right)\mathrm{,}{y}_2\left(x\mathrm{,}\lambda \right)\right\}$ of (1) automatically satisfy the boundary conditions (43) at $x=a$. Using the solutions $\left\{Y_1\left(x\mathrm{,}\lambda \right)\mathrm{,}{Y}_2\left(x\mathrm{,}\lambda \right)\right\}$ of (1) defined by (49) we define the $2\times 2$ matrices

$U\left(x,\lambda \right):=\left(\begin{array}{cc}{y}_1\left(x\mathrm{,}\lambda \right)& y_2\left(x\mathrm{,}\lambda \right)\\ {y^{″}}_1\left(x\mathrm{,}\lambda \right)& {y^{″}}_2\left(x\mathrm{,}\lambda \right)\end{array}\right),$ (50)

and

$V\left(x\mathrm{,}\lambda \right)\mathrm{<mo>\; :\; </mo>\; <mo>\; =\; </mo>}\left(\begin{array}{cc}-{y^{‴}}_1\left(x\mathrm{,}\lambda \right)+s\left(x\right){y^{\prime }}_1\left(x\mathrm{,}\lambda \right)& -{y^{‴}}_2\left(x\mathrm{,}\lambda \right)+s\left(x\right){y^{\prime }}_2\left(x\mathrm{,}\lambda \right)\\ -{y^{\prime }}_1\left(x\mathrm{,}\lambda \right)& -{y^{\prime }}_2\left(x\mathrm{,}\lambda \right)\end{array}\right)\mathrm{.}$ (51)

The two-dimensional subspace,

$\Im =\left\{y\left(x,\lambda \right):y\left(x,\lambda \right)=A{y}_1\left(x,\lambda \right)+B{y}_2\left(x,\lambda \right)\text{for}\text{some}\text{choice}\text{of}A\text{and}B\right\},$ (52)

where $y_1\left(x\mathrm{,}\lambda \right)$ and $y_2\left(x\mathrm{,}\lambda \right)$ are solutions of (1) corresponding to the $Y_1\left(x\mathrm{,}\lambda \right)$ and $Y_2\left(x\mathrm{,}\lambda \right)$ solutions of the IVP, (1) and (49), of the four-dimensional solution space of (1) satisfies the boundary conditions (43) at $x=a$.

We have the following theorem giving necessary sufficient conditions (N.S.Cs) for the eigenvalues of the SL problem for Equation (1) which has boundary conditions at $x=a$ (43) and boundary conditions at $x=b$ (44).

Theorem 2.1 1) A N.S.C. for $\lambda \in \left(-\infty \mathrm{,}\infty \right)$ to be an eigenvalue of the SL problem for Equation (1) with boundary conditions at $x=a$ (43) and boundary conditions at $x=b$ (44), and having multiplicity one, is:

$\left|\begin{array}{cc}{r}_1\left(y_1\right)& r_1\left(y_2\right)\\ r_2\left(y_1\right)& r_2\left(y_2\right)\end{array}\right|=\left|\begin{array}{cc}{y}_1\left(b,\lambda \right)& y_2\left(b,\lambda \right)\\ {y^{″}}_1\left(b,\lambda \right)& {y^{″}}_2\left(b,\lambda \right)\end{array}\right|=0$ (53)

and

$rank\left(\begin{array}{cc}{r}_1\left(y_1\right)& r_1\left(y_2\right)\\ r_2\left(y_1\right)& r_2\left(y_2\right)\end{array}\right)=rank\left(\begin{array}{cc}{y}_1\left(b,\lambda \right)& y_2\left(b,\lambda \right)\\ {y^{″}}_1\left(b,\lambda \right)& {y^{″}}_2\left(b,\lambda \right)\end{array}\right)=1.$ (54)

2) A N.S.C. for $\lambda \in \left(-\infty \mathrm{,}\infty \right)$ to be an eigenvalue of the SL problem for Equation (1) with boundary conditions at $x=a$ (43) and boundary conditions at $x=b$ (44), and having multiplicity one, is:

$\left|\begin{array}{cc}{r}_1\left(y_1\right)& r_1\left(y_2\right)\\ r_2\left(y_1\right)& r_2\left(y_2\right)\end{array}\right|=\left|\begin{array}{cc}{y}_1\left(b,\lambda \right)& y_2\left(b,\lambda \right)\\ {y^{″}}_1\left(b,\lambda \right)& {y^{″}}_2\left(b,\lambda \right)\end{array}\right|=0$ (55)

and

$rank\left(\begin{array}{cc}{r}_1\left(y_1\right)& r_1\left(y_2\right)\\ r_2\left(y_1\right)& r_2\left(y_2\right)\end{array}\right)=rank\left(\begin{array}{cc}{y}_1\left(b,\lambda \right)& y_2\left(b,\lambda \right)\\ {y^{″}}_1\left(b,\lambda \right)& {y^{″}}_2\left(b,\lambda \right)\end{array}\right)=0.$ (56)

Proof Let $\left\{y_1\left(x\mathrm{,}\lambda \right)\mathrm{,}{y}_2\left(x\mathrm{,}\lambda \right)\right\}$ be the unique solutions of the 4th order Equation (1) which are defined by the initial conditions at $x=a$ and $\left\{Y_1\left(x\mathrm{,}\lambda \right)\mathrm{,}{Y}_2\left(x\mathrm{,}\lambda \right)\right\}$ be the corresponding solutions of the Hamiltonian system (1):

$\left(\begin{array}{cc}{y}_1\left(a,\lambda \right)& y_2\left(a,\lambda \right)\\ {y^{″}}_1\left(a,\lambda \right)& {y^{″}}_2\left(a,\lambda \right)\\ -{y^{‴}}_1\left(a,\lambda \right)+s\left(a\right){y^{\prime }}_1\left(a,\lambda \right)& -{y^{‴}}_2\left(a,\lambda \right)+s\left(a\right){y^{\prime }}_2\left(a,\lambda \right)\\ -{y^{\prime }}_1\left(a,\lambda \right)& -{y^{\prime }}_2\left(a,\lambda \right)\end{array}\right)=\left(\begin{array}{cc}0& 0\\ 0& 0\\ 1& 0\\ 0& 1\end{array}\right)$ (57)

By fixing (57) the two-dimensional space $\Im$ is fixed by the $4\times 2$ matrix (57). The constants in the (57) matrix were chosen to ensure that the boundary conditions at $x=a$ (43) was satisfied, so we know that

$\left(\begin{array}{c}{l}_1\left(y_1\right)\\ l_2\left(y_1\right)\end{array}\right)=\left(\begin{array}{c}{y}_1\left(a,\lambda \right)\\ {y^{″}}_1\left(a,\lambda \right)\end{array}\right)=\left(\begin{array}{c}0\\ 0\end{array}\right),$ (58)

$\left(\begin{array}{c}{l}_1\left(y_2\right)\\ l_2\left(y_2\right)\end{array}\right)=\left(\begin{array}{c}{y}_2\left(a,\lambda \right)\\ {y^{″}}_2\left(a,\lambda \right)\end{array}\right)=\left(\begin{array}{c}0\\ 0\end{array}\right),$ (59)

that is, that both $y_1\left(x\mathrm{,}\lambda \right)$ and $y_2\left(x\mathrm{,}\lambda \right)$ satisfy the boundary conditions at $x=a$ (43). Also, of course, the space $\Im$ of solutions spanned by these two solutions $\left\{y_1\left(x\mathrm{,}\lambda \right)\mathrm{,}{y}_2\left(x\mathrm{,}\lambda \right)\right\}$ of (1) satisfies the boundary conditions at $x=a$ (43), that is

$l_1\left(A{y}_1\left(x,\lambda \right)+B{y}_2\left(x,\lambda \right)\right)=0,$ (60)

$l_2\left(A{y}_1\left(x\mathrm{,}\lambda \right)+B{y}_2\left(x\mathrm{,}\lambda \right)\right)=0$ (61)

for all $A\mathrm{,}B\in ℂ$. It remains only to apply the boundary conditions at $x=b$ (44), i.e. to require that

$r_1\left(y\right)=y\left(b,\lambda \right)=0,$ (62)

and

$r_2\left(y\right)=y^{″}\left(b,\lambda \right)=0$ (63)

for boundary conditions at $x=b$ (44). Hence, we find

$r_1\left(A{y}_1\left(x,\lambda \right)+B{y}_2\left(x,\lambda \right)\right)=A{r}_1\left(y_1\left(x,\lambda \right)\right)+B{r}_1\left(y_2\left(x,\lambda \right)\right)=0$ (64)

and

$r_2\left(A{y}_1\left(x,\lambda \right)+B{y}_2\left(x,\lambda \right)\right)=A{r}_2\left(y_1\left(x,\lambda \right)\right)+B{r}_2\left(y_2\left(x,\lambda \right)\right)=0$ (65)

as the requirement for solutions in $\Im$ to also satisfy the boundary conditions at $x=b$ (44). But the Equations (64) and (65) can have a nonzero solution for $A\mathrm{,}B$ if and only if

$\left|\begin{array}{cc}{r}_1\left(y_1\right)& r_1\left(y_2\right)\\ r_2\left(y_1\right)& r_2\left(y_2\right)\end{array}\right|=\left|\begin{array}{cc}{y}_1\left(b,\lambda \right)& y_2\left(b,\lambda \right)\\ {y^{″}}_1\left(b,\lambda \right)& {y^{″}}_2\left(b,\lambda \right)\end{array}\right|=0$ (66)

Hence (66) is a N.S.C. condition for $\lambda$ to be an eigenvalue of the SL problem for Equation (1). The multiplicity of the eigenvalue is defined as the number of linearly independent solutions of (1) which satisfy both boundary conditions at $x=a$ and both boundary conditions at $x=b$. Since the dimension of $\Im$ is two, this can be at most two. For multiplicity one, we must have

$rank\left(\begin{array}{cc}{r}_1\left(y_1\right)& r_1\left(y_2\right)\\ r_2\left(y_1\right)& r_2\left(y_2\right)\end{array}\right)=rank\left(\begin{array}{cc}{y}_1\left(b,\lambda \right)& y_2\left(b,\lambda \right)\\ {y^{″}}_1\left(b,\lambda \right)& {y^{″}}_2\left(b,\lambda \right)\end{array}\right)=1,$ (67)

and in this case the solution of the two Equations (64) and (65) will be

$\left(\begin{array}{c}A\\ B\end{array}\right)=\left(\begin{array}{c}\frac{-r_1\left(y_2\right)}{r_1\left(y_1\right)}\\ 1\end{array}\right),\text{if}{r}_1\left(y_1\right)\ne 0$ (68)

and

$\left(\begin{array}{c}A\\ B\end{array}\right)=\left(\begin{array}{c}1\\ 0\end{array}\right),\text{if}{r}_1\left(y_1\right)=0.$ (69)

In this case the eigenfunction, which is unique up to a constant multiple, is

$y\left(x,\lambda \right)=B\left(\frac{-r_1\left(y_2\right)}{r_1\left(y_1\right)}\cdot y_1\left(x,\lambda \right)+y_2\left(x,\lambda \right)\right),\text{if}{r}_1\left(y_1\right)\ne 0$ (70)

or

$y\left(x,\lambda \right)=A{y}_1\left(x,\lambda \right),\text{if}{r}_1\left(y_1\right)=0.$ (71)

For multiplicity two, we need to require

$rank\left(\begin{array}{cc}{r}_1\left(y_1\right)& r_1\left(y_2\right)\\ r_2\left(y_1\right)& r_2\left(y_2\right)\end{array}\right)=rank\left(\begin{array}{cc}{y}_1\left(b,\lambda \right)& y_2\left(b,\lambda \right)\\ {y^{″}}_1\left(b,\lambda \right)& {y^{″}}_2\left(b,\lambda \right)\end{array}\right)=0,$ (72)

and in this case $y_1\left(x\mathrm{,}\lambda \right)$ and $y_2\left(x\mathrm{,}\lambda \right)$ would be linearly independent eigenfunctions of (1). $■$

More generally, the general case of boundary conditions (3) would give

$f\left(\lambda \right):=\det \left(B_{1}U\left(b,\lambda \right)+B_{2}V\left(b,\lambda \right)\right)=0$ (73)

as a N.S.C. for $\lambda$ to be an eigenvalue of (1) with boundary conditions (44) and (3). More generally, the boundary conditions (2) could be handled by changing the initial conditions (57) appropriately.

3. Description of the MG4 Algorithm

We obtained the $4\times 4$ transfer matrix

$M=\exp \left(h\tilde{A}\right)$ (1)

by doing the following steps:

1) Calculate the eigenvalues and the eigenvectors of $\tilde{A}$.

2) Diagonalize

$\tilde{A}=P\cdot D\cdot P^{-1},$ (2)

where P denotes the matrix of eigenvectors of $\tilde{A}$, and D denotes the diagonal matrix of the four eigenvalues of $\tilde{A}$ (12).

3) Put

$M=P\left(\exp \left(hD\right)\right)P^{-1}.$ (3)

This gives a $4\times 4$ transfer matrix with 4 cases of eigenvalues of $\tilde{A}$. (The matrix elements could be reduced to expressions involving sinh, cosh, sin and cos functions). Our MG4 code implements the above transfer matrix $M=\exp \left(h\tilde{A}\right)$ by doing the matrix multiplication (3) numerically. Here we describe the implementation of the MG4 method for computing the eigenvalues of the SL problem for Equation (1) with the choices of Dirichlet boundary conditions (43) and (44) at the left and right endpoints. To impose the boundary conditions (44) on $\Im$, we integrate the IVP, (1) and (49), from $x=a$ to $x=b$ using the MG4 method on the $4\times 2$ matrix $Y\left(x\mathrm{,}\lambda \right)$. Then when $Y\left(b\mathrm{,}\lambda \right)$ has been computed, the boundary conditions (44),

$y\left(b,\lambda \right)=A{y}_1\left(b,\lambda \right)+B{y}_2\left(b,\lambda \right)=0$ (4)

$y^{″}\left(b,\lambda \right)=A{y^{″}}_1\left(b,\lambda \right)+B{y^{″}}_2\left(b,\lambda \right)=0,$ (5)

will be satisfied for some choices of real constants A and B, not both zero, if and only if

$f\left(\lambda \right):=\det \left(U\left(b,\lambda \right)\right)=\left|\begin{array}{cc}{y}_1\left(b,\lambda \right)& y_2\left(b,\lambda \right)\\ {y^{″}}_1\left(b,\lambda \right)& {y^{″}}_2\left(b,\lambda \right)\end{array}\right|=0.$ (6)

The computation is performed using an initial uniform mesh, applying bisection method with initial upper and lower bounds for a given eigenvalue ${\lambda }_n$, and then doubling the number of mesh points by bisecting the mesh to generate a Richardson h⁴-extrapolation table over successively bisected meshes. Then the extrapolated eigenvalue is selected when the eigenvalue extrapolation error satisfies a tolerance test.

3.1. Description of h⁴-Richardson’s Extrapolation

If $q\left(x\right)\mathrm{,}s\left(x\right)\in C^{\infty }\left[a\mathrm{,}b\right]$, we can assume that if MG4 method is applied we will have for each choice of h:

$\hat{\lambda }\left(h\right)\approx {\lambda }_{Exact}+{\tau }_1{h}^4+{\tau }_2{\left(h^4\right)}^2+\cdots +{\tau }_m{\left(h^4\right)}^m$ (7)

for some choices of real constants ${\tau }_1\mathrm{,}{\tau }_2\mathrm{,}\cdots \mathrm{,}{\tau }_m$. For each $m=1,2,3,\cdots$

$\hat{\lambda }\left(h\right)-{\lambda }_{Exact}={\tau }_1{h}^4+{\tau }_2{\left(h^4\right)}^2+\cdots +{\tau }_m{\left(h^4\right)}^m+O\left({\left(h^4\right)}^{m+1}\right)=O\left(h^4\right).$ (8)

Putting

$x_1=h_0^4$ (9)

$x_2=h_1^4={\left(\frac{h_0}{2}\right)}^4$ (10)

$⋮$

$x_i=h_{i-1}^4={\left(\frac{h_0}{2^{i-1}}\right)}^4,$ (11)

in Neville’s algorithm ( [39] , 2.1.2.5b, p. 42),

$T_{ik}\left(x\right)=T_{i,k-1}\left(x\right)+\frac{T_{i,k-1}\left(x\right)-T_{i-1,k-1}\left(x\right)}{\frac{x-x_{i-k}}{x-x_i}-1},1\le k\le i,i=0,1,2,\cdots ,n,$ (12)

we find

$T_{I-1,J-1}\left(0\right)=T_{I-1,J-2}\left(0\right)+\frac{T_{I-1,J-2}\left(0\right)-T_{I-2,J-2}\left(0\right)}{{16}^{J-1}-1}$ (13)

where we have taken

$k=J-1,i=I-1,x=0,x_i=h_{I-1}^4,x_{i-k}=h_{I-J}^4.$

Applying Neville’s algorithm generates the h⁴-Richardson’s extrapolation table for the eigenvalue computation. Defining

${\hat{\lambda }}_{IJ}:=T_{I-1,J-1}\left(0\right)$ (14)

we have from (13) that

${\hat{\lambda }}_{I,J}={\hat{\lambda }}_{I,J-1}+\frac{{\hat{\lambda }}_{I,J-1}-{\hat{\lambda }}_{I-1,J-1}}{{16}^{J-1}-1},$ (15)

where

${\hat{\lambda }}_{11}=\text{computed value for}h=h_0,$

${\hat{\lambda }}_{21}=\text{computed value for}h=\frac{h_0}{2},$

$⋮$

${\hat{\lambda }}_{n1}=\text{computed value for}h=\frac{h_0}{2^n}.$

here the second term in (15) is the extrapolation error. The first column of the extrapolation table, that is, the eigenvalues ${\hat{\lambda }}_{11}\mathrm{,}{\hat{\lambda }}_{21}\mathrm{,}\cdots$, are computable quantities. The columns two, three, four, five, $\cdots$, are generated from column one using (15).

3.2. Computing Large Eigenvalues by Using the Invariant-Imbedding Variables

In a manner similar to Greenberg and Marletta in their SLEUTH code (see [14] , Section 3.2, pp. 461-462), we applied the change to “invariant-imbedding” variables, $V{U}^A,U{V}^A,\det U$ and $\det V$, in our code in order to provide a good stable integration scheme. We generated the “invariant-imbedding” variables by using the $2\times 2$ matrices $U\left(x\mathrm{,}\lambda \right)$ and $V\left(x\mathrm{,}\lambda \right)$ which we defined in (50) and (51) for eigenvalue computation of the SL problem.

4. Numerical Results and Discussion

In this section we give some numerical outputs for each of the 5 test problems in Section 1.1, and compare with the comparable SLEUTH outputs. The 5 test problems are squares of 2nd order SL equations. For such problems the choice of Dirichlet boundary conditions for the 2nd order problem, generates, by squaring, a 4th order SL problem whose eigenvalues are the squares of the 2nd order SL problem, Greenberg and Marletta ( [14] pp. 478-481). For example, we consider in section 5 the following 2nd order SL problems

$-y^{″}+\left(\frac{-1}{4x^2}\right)y=\lambda y$ (1)

$y\left(1\right)=0,y\left(5\right)=0$

$-y^{″}+\left(x^2+x^4\right)y=\lambda y$ (2)

$y\left(1\right)=0,y\left(5\right)=0$

$-y^{″}+\left(\cos \left(x\right)+2\cos \left(2x\right)+3\cos \left(3x\right)\right)y=\lambda y$ (3)

$y\left(0\right)=0,y\left(\pi \right)=0$

$-y^{″}+\left({\beta }^2{\sin }^2\left(2x\right)-2\beta \cos \left(2x\right)\right)y=\lambda y$ (4)

$y\left(\frac{-\pi }{2}\right)=0,y\left(\frac{\pi }{2}\right)=0$

$-y^{″}+\frac{1}{4}{\sec }^2\left(x\right)y=\lambda y$ (5)

$y\left(0\right)=0,y\left(\frac{\pi }{4}\right)=0$

Squaring the self-adjoint operator corresponding to (1)-(5) gives the 4th order self-adjoint operator corresponding to the 4th order problems in Tables 1-5, respectively. Accordingly, the eigenvalues of the problems in Tables 1-5 are the squares of the eigenvalues of the 2nd order SL problems (1)-(5), respectively. Tables 1-5 give outputs of MG4 and SLEUTH codes on the test problems 1, 2, 3, 4 and 5 respectively, with the choices of Dirichlet boundary conditions. In these tables, we list the SLEUTH and MG4 outputs to 17 digits. The number of these digits which are correct is always a key issue in assessing the performance of a numerical algorithm. Since the exact eigenvalues of these 4th order SL problems are the squares of the exact eigenvalues of the 2nd order SL problems (1)-(5), we computed the eigenvalues of the 2nd order SL problems (1)-(5) and computed their squares to provide a benchmark against which we can compare MG4 and SLEUTH algorithm outputs. The purpose for the following tables is to make comparisons at reasonably high accuracy; so we ran the MG4 and SLEUTH codes with the tolerance parameter TOL = 10⁻¹². Outputs of the SLEDGE code (Sturm-Liouville Estimates Determined by Global Error Control) of Pruess and Fulton [6] for the problems (1)-(5) were obtained for the eigenvalue indices listed in Tables 1-5, and then their squares were computed to generate the second columns of Tables 1-5. Since SLEDGE is known to compute eigenvalues quite reliably to the user-requested accuracy, we used the squares of the SLEDGE eigenvalues as the benchmark for MG4 and SLEUTH codes in Tables 1-5. The error characterization $E_{sledge}$ in columns 4 and 6 of Tables 1-5 was computed as

$E_{sledge}=(\begin{array}{l}\left|{\lambda }_n^{\left(SLEDGE\right)}-{\lambda }_n\right|\mathrm{,}\text{if}{\lambda }_n<1\\ \left|\frac{{\lambda }_n^{\left(SLEDGE\right)}-{\lambda }_n}{{\lambda }_n^{\left(SLEDGE\right)}}\right|\mathrm{,}\text{if}{\lambda }_n\ge \mathrm{1,}\end{array}$ (6)

where ${\lambda }_n^{\left(SLEDGE\right)}$ are the SLEDGE-squared eigenvalues and ${\lambda }_n$ are eigenvalues obtained by the SLEUTH and MG4 codes. This represents the absolute or relative eigenvalue errors of each code relative to the benchmark values obtained from SLEDGE. The obtained absolute or relative eigenvalue errors of each code are used to measure the performance of each method. In Table 1, we see that the error characterization $E_{sledge}$ for the eigenvalues ${\lambda }_0$, ${\lambda }_{20}$ and ${\lambda }_{100}$ obtained by the SLEUTH are slightly larger than the error characterization $E_{sledge}$ for the eigenvalues ${\lambda }_0$, ${\lambda }_{20}$ and ${\lambda }_{100}$ achieved by the MG4 method.

Table 1. Eigenvalues of Problem 1.

Table 2. Eigenvalues of Problem 2.

Table 3. Eigenvalues of Problem 3.

Table 4. Eigenvalues of Problem 4.

Table 5. Eigenvalues of Problem 5.

In Table 2, we see that the error characterization $E_{sledge}$ for the eigenvalues ${\lambda }_0$, ${\lambda }_{50}$ and ${\lambda }_{100}$ obtained by the SLEUTH are slightly larger than the error characterization $E_{sledge}$ for the eigenvalues ${\lambda }_0$, ${\lambda }_{50}$ and ${\lambda }_{100}$ achieved by the MG4 method. In Table 3, we see that the error characterization $E_{sledge}$ for the eigenvalues ${\lambda }_0$ and ${\lambda }_{100}$ achieved by the MG4 method are quite comparable to the error characterization $E_{sledge}$ for the eigenvalues ${\lambda }_0$ and ${\lambda }_{100}$ obtained by the SLEUTH. But, the error characterization $E_{sledge}$ for the eigenvalue ${\lambda }_{50}$ obtained by the SLEUTH is slightly larger than the error characterization $E_{sledge}$ for the eigenvalue ${\lambda }_{50}$ achieved by the MG4 method. In Table 4, we see that the error characterization $E_{sledge}$ for the eigenvalue ${\lambda }_{100}$ achieved by the MG4 method is quite comparable to the error characterization $E_{sledge}$ for the eigenvalue ${\lambda }_{100}$ obtained by the SLEUTH. But, the error characterization $E_{sledge}$ for the eigenvalues ${\lambda }_2$ and ${\lambda }_{50}$ obtained by the SLEUTH are slightly larger than the error characterization $E_{sledge}$ for the eigenvalues ${\lambda }_2$ and ${\lambda }_{50}$ achieved by the MG4 method. In Table 5, we see that the error characterization $E_{sledge}$ for the eigenvalues ${\lambda }_0$, ${\lambda }_8$, ${\lambda }_{30}$ and ${\lambda }_{100}$ obtained by the SLEUTH are slightly larger than the error characterization $E_{sledge}$ for the eigenvalues ${\lambda }_0$, ${\lambda }_8$, ${\lambda }_{30}$ and ${\lambda }_{100}$ achieved by the MG4 method.

Remark 5.1 The machine precision, obtained from the FORTRAN routine EPSLON, on the desktop computer (with Pentium 4 processors) used to obtain the following outputs of the SLEUTH and SLEDGE codes was MACHEPS = 2.22D-16.

5. Conclusion

In this paper we have presented the MG4 algorithm of Iserles et al. [33] [40] , for eigenvalue computation of regular 4th order Sturm-Liouville problems. Applying the change to “Invariant-Imbedding” variables in a manner similar to Greenberg and Marletta in their SLEUTH code [14] provides good stabilization for the MG4 algorithm, and this resulted in its capability for achieving high accuracy, very often on the order of machine precision. The MG4 algorithm appears to be competitive to the SLEUTH code.

Acknowledgements

The author would like to thank Al Baha University for their financial and moral support. The author would also like to express his sincere gratitude to the following people:

1) Professor Charles Fulton and Dr. Steven Pruess for supplying a FORTRAN 90 version of the SLEDGE code.

2) Professor M. Marletta for supplying a FORTRAN 77 version of the SLEUTH code based on use of the BLAS underlying the LAPACK software, and for suggesting a modification to it which allowed us to run SLEUTH with high accuracy requests.

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Alalyani, A. (2019) Numerical Methods for Eigenvalue Computation of Fourth-Order Self-Adjoint Ordinary Differential Operators. PhD Dissertation, Florida Institute of Technology, Melbourne.
[2]	Everitt, W.N. (1957) The Sturm-Liouville Problem for Fourth Order Differential Equations. The Quarterly Journal of Mathematics, 8, 146-160. https://doi.org/10.1093/qmath/8.1.146
[3]	Fulton, C. (1989) The Bessel-Squared Equation in the Lim-2, Lim-3, Im-4 Cases. The Quarterly Journal of Mathematics, 40, 423-456. https://doi.org/10.1093/qmath/40.4.423
[4]	Bailey, P., Gordon, M. and Shampine, L. (1978) Automatic Solution of the Sturm-Liouville Problem. ACM Transactions on Mathematical Software, 4, 193-208. https://doi.org/10.1145/355791.355792
[5]	Bailey, P.B., Everitt, W.N. and Zettl, A. (2001) The SLEIGN2 Sturm-Liouville Code. ACM Transactions on Mathematical Software, 21, 143-192. https://doi.org/10.1145/383738.383739
[6]	Pruess, S. and Fulton, C. (1993) Mathematical Software for Sturm-Liouville Problems. ACM Transactions on Mathematical Software, 19, 360-376. https://doi.org/10.1145/155743.155791
[7]	Fulton, C. and Pruess, S. (1998) The Computation of Spectral Density Functions for Singular Sturm-Liouville Problems Involving Simple Continuous Spectra. ACM Transactions on Mathematical Software, 34, 107-129. https://doi.org/10.1145/285861.285867
[8]	Pryce, J.D. (1993) Numerical Solution of Sturm-Liouville Problems. Oxford University Press, Oxford.
[9]	Pryce, J.D. and Marletta, M. (1992) Automatic Solution of Sturm-Liouville Problems Using the Pruess Method. Journal of Computational and Applied Mathematics, 39, 57-78. https://doi.org/10.1016/0377-0427(92)90222-J
[10]	Pryce, J. (1986) Error Control of Phase-Function Shooting Methods for Sturm-Liouville Problems. IMA Journal of Numerical Analysis, 6, 103-123. https://doi.org/10.1093/imanum/6.1.103
[11]	Ledoux, V., Van Daele, M. and Vanden Berghe, G. (2004) CP Methods of Higher Order for Sturm-Liouville and Schrödinger Equations. Computer Physics Communications, 162, 153-165. https://doi.org/10.1016/j.cpc.2004.07.001
[12]	Ledoux, V., Van Daele, M. and Vanden Berghe, G. (2005) MATSLISE: A Matlab Package for the Numerical Solution of Sturm-Liouville and Schrödinger Equations. ACM Transactions on Mathematical Software, 31, 532-554. https://doi.org/10.1145/1114268.1114273
[13]	Ixaru, L. (2000) CP Methods for the Schrödinger Equation. Journal of Computational and Applied Mathematics, 125, 347-357. https://doi.org/10.1016/S0377-0427(00)00478-7
[14]	Greenberg, L. and Marletta, M. (1997) Algorithm 775: The Code SLEUTH for Solving Fourth-Order Sturm-Liouville Problems. ACM Transactions on Mathematical Software, 23, 453-493. https://doi.org/10.1145/279232.279231
[15]	Greenberg, L. and Marletta, M. (1995) Oscillation Theory and Numerical Solution of Fourth-Order Sturm-Liouville Problems. IMA Journal of Numerical Analysis, 15, 319-356. https://doi.org/10.1093/imanum/15.3.319
[16]	Greenberg, L. (1991) An Oscillation Method for Fourth-Order, Self-Adjoint, Two-Point Boundary Value Problems with Nonlinear Eigenvalues. SIAM Journal on Mathematical Analysis, 22, 1021-1042. https://doi.org/10.1137/0522067
[17]	Chanane, B. (2010) Accurate Solutions of Fourth Order Sturm-Liouville Problems. Journal of Computational and Applied Mathematics, 234, 3064-3071. https://doi.org/10.1016/j.cam.2010.04.023
[18]	Chanane, B. (1998) Eigenvalues of Fourth Order Sturm-Liouville Problems Using Fliess Series. Journal of Computational and Applied Mathematics, 96, 91-97. https://doi.org/10.1016/S0377-0427(98)00086-7
[19]	Chanane, B. (2002) Fliess Series Approach to the Computation of the Eigenvalues of Fourth-Order Sturm-Liouville Problems. Applied Mathematics Letters, 15, 459-563. https://doi.org/10.1016/S0893-9659(01)00159-8
[20]	Chanane, B. (1998) Eigenvalues of Sturm-Liouville Problems Using Fliess Series. Applicable Analysis, 69, 233-238. https://doi.org/10.1080/00036819808840659
[21]	Mihaila, B. and Mihaila, I. (2002) Numerical Approximations Using Chebyshev Polynomial Expansions: El-gendi’s Method Revisited. Journal of Physics A: Mathematical and General, 35, 731-746. https://doi.org/10.1088/0305-4470/35/3/317
[22]	El-Gamel, M. and Sameeh, M. (2012) An Efficient Technique for Finding the Eigenvalues of Fourth-Order Sturm-Liouville Problems. Applied Mathematics, 3, 920-925. https://doi.org/10.4236/am.2012.38137
[23]	El-Gamel, M. and Abd El-hady, M. (2013) Two Very Accurate and Efficient Methods for Computing Eigenvalues of Sturm-Liouville Problems. Applied Mathematical Modeling, 37, 5039-5046. https://doi.org/10.1016/j.apm.2012.10.019
[24]	Fox, L. (1962) Chebyshev Methods for Ordinary Differential Equations. The Computer Journal, 4, 318-331. https://doi.org/10.1093/comjnl/4.4.318
[25]	Saleh Taher, A.H., Malek, A. and Momeni-Masuleh, S.H. (2013) Chebyshev Differential Matrices for Efficient Computation of the Eigenvalues of Fourth-Order Sturm-Liouville Problems. Applied Mathematical Modeling, 37, 4634-4642. https://doi.org/10.1016/j.apm.2012.09.062
[26]	Ycel, U. (2015) Numerical Approximations of Sturm-Liouville Eigenvalues Using Chebyshev Polynomial Expansions Method. Cogent Mathematics, 2, Article ID: 1045223. https://doi.org/10.1080/23311835.2015.1045223
[27]	Ycel, U. and Boubaker, K. (2012) Differential Quadrature Method (DQM) and Boubaker Polynomials Expansion Scheme (BPES) for Efficient Computation of the Eigenvalues of Fourth-Order Sturm-Liouville Problems. Applied Mathematical Modelling, 36, 158-167. https://doi.org/10.1016/j.apm.2011.05.030
[28]	Khmelnytskaya, K., Kravchenko, V. and Baldenebro-Obeso, J. (2012) Spectral Parameter Power Series for Fourth-Order Sturm-Liouville Problems. Applied Mathematics and Computation, 219, 3610-3524. https://doi.org/10.1016/j.amc.2012.09.055
[29]	Kravchenko, V. and Porter, R. (2010) Spectral Parameter Power Series for Sturm-Liouville Problems. Mathematical Methods in the Applied Sciences, 33, 459-468. https://doi.org/10.1002/mma.1205
[30]	Kravchenko, V. and Porter, R. (2015) Eigenvalue Problems, Spectral Parameter Power Series, and Modern Applications. Mathematical Methods in the Applied Sciences, 38, 1945-1969. https://doi.org/10.1002/mma.3213
[31]	Fulton, C. and Krall, A. (1983) Self-Adjoint 4th Order Boundary Value Problem in the Lim-4 Case. Proceedings of Symposium on Ordinary Differential Equations and Operators, Dundee, Lecture Notes in Mathematics 1032, 240-256. https://doi.org/10.1007/BFb0076800
[32]	Atkinson, F. (1964) Discrete and Continuous Boundary Problems. Academic Press, New York. https://doi.org/10.1063/1.3051875
[33]	Iserles, A., Marthinsen, A. and Nrsett, S.P. (1999) On the Implementation of the Method of Magnus Series for Linear Differential Equations. BIT Numerical Mathematics, 39, 281-304. https://doi.org/10.1023/A:1022393913721
[34]	Greenberg, L. and Marletta, M. (1997) The Validity of Richardson Extrapolation in Coefficient Approximation Methods for Fourth Order and Higher Order Self-Adjoint ODE Eigenvalue Problems. Technical Report, Department of Mathematics and Computer Science, University of Leicester, Leicester.
[35]	Greenberg, L. (1991) A Prüfer Method for Calculating Eigenvalues of Self-Adjoint Systems of Ordinary Differential Equations: Parts 1 and 2. Technical Report TR91-24, University of Maryland at College Park, College Park.
[36]	Gallier, J. (2005) Manifolds, Lie Groups, Lie Algebras, Riemannian Manifolds, with Applications to Computer Vision and Robotics [Lecture Notes]. Department of Computer and Information Science, University of Pennsylvania, Philadelphia.
[37]	Magnus, W. (1954) On the Exponential Solution of Differential Equations for a Linear Operator. Communications on Pure and Applied Mathematics, 7, 649-673. https://doi.org/10.1002/cpa.3160070404
[38]	Iserles, A. and Nrsett, S.P. (1999) On the Solution of Linear Differential Equations in Lie Groups. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 357, 983-1019. https://doi.org/10.1098/rsta.1999.0362
[39]	Stoer, J. and Bulirsch, R. (1993) Introduction to Numerical Analysis. Second Edition, Springer-Verlag, New York. https://doi.org/10.1007/978-1-4757-2272-7
[40]	Iserles, A., Munthe-Kaas, H.Z., Norsett, S.P. and Zanna, A. (2005) Lie-Group Methods. Acta Numerica, 9, 1-148. https://doi.org/10.1017/S0962492900002154