In this paper1 we consider the self-adjoint differential operators which arise from the 4th order differential equation
when separated, self-adjoint boundary conditions are imposed at each of the two regular endpoints and .
We make the assumptions:
*2010 Mathematics Subject Classication 34L15; 34L16; 65L15; 70J50.
1The content of this paper is related to my Ph.D. dissertation  .
1) is continuous on .
2) and are continuous on .
Under these assumptions both endpoints a and b are regular endpoints. The most general separated, self-adjoint boundary conditions which can be imposed at and are
where and are any choice of real, matrices satisfying the properties
The above boundary conditions can be shown to be equivalent to the general forms of boundary conditions used by Everitt  (in his PhD dissertation on 4th order Sturm-Liouville problems), Fulton  and many others.
The domain of the maximal operator associated with the Equation (1) on the closed interval is
where is the space of functions which are absolutely continuous on compact subsets of . The self-adjoint operators associated with Equation (1) are then obtained by restricting by two boundary conditions at the left endpoint and two boundary conditions at the right endpoint as in (2) and (3), namely
The Green’s formula for the 4th order equation is
where the bilinear concomitant is defined as
Using this definition, and the boundary conditions (2) and (3), it can be shown that the operators on are symmetric; that is, for all :
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  , SLEIGN 2  , SLEDGE   , SLO2FM    , MATSLISE    . 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  .
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  . 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    for the number of eigenvalues less than .
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  (also used for 2nd order SL problems).
2) Fliess Series Method    which represents the solution of an IVP for the 4th order Equation (1) in terms of iterated integrals involving the coefficient functions, and .
3) Chebyshev Method  -  which approximates solutions of the 4th order Equation (1) by Chebyshev polynomials.
4) Boubaker Polynomials Expansion Scheme (BPES)  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    which expands the solutions of the SL equation (both second and 4th order) in a convergent Taylor expansion in the eigenparameter . 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
2) The square of the 2nd order Modified Harmonic Oscillator equation
3) The square of the 2nd order equation
4) The square of the 2nd order Coffey-Evan equation with
5) The square of the 2nd order Legendre equation
Problem 5, the Legendre squared equation, arises from changes of variables to the non-LNF form discussed in  .
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 continuous. The Equation (1) can be converted to the 1st order system (Atkinson  , pp. 323-324) , where
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  , which is available from NETLIB at ORNL. The SLEUTH (Sturm-Liouville Eigenvalues Using Theta Matrices) code employs an MG2 approximation for the solution , (see  , p. 283), on each mesh interval of a Hamiltonian system similar to the system (1) (the order of the derivatives in the above vector being slightly different). The SLEUTH code is based on using formulas of Greenberg     for the number of eigenvalues less than , so the eigenvalue algorithm is quite different than the method we are proposing here for eigenvalue computation for (1) using (1).
For the IVP, , , 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  , Prop. 7.2.3 and Def. 7.2.4).
Definition 2.1 is the Lie-group and defined as:
Definition 2.2 is the Lie-algebra and defined as:
Lemma 2.1 If , then , i.e.
Remark 2.2 For any constant matrix , it follows from this lemma that the solution
of the IVP,
lies in the Lie Group .
The Magnus methods originate (see  ) with the observation that an analytical solution of (1) with initial condition can be written as, (see  , p. 283),
and where the square brackets denote the matrix commutator and are defined as:
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 and :
For the Hamiltonian system (1), we put
(meaning that and in (1) are to be evaluated at and in the nth mesh interval ). Then the MG4 method of Iserles and Norsett (  , p. 1012), and (  , p. 288) for (1) takes the form (for a fixed step size h)
where the transfer matrix for passing from to is with
where the square bracket denotes the matrix commutator and is defined as:
The four eigenvalues of are
Eigenvalues of matrix:
Let us define
Then the following four cases of eigenvalues of arise, involving both complex and real eigenvalues:
Case 1: and ( ).
Case 2: and ( ).
Case 3: and ( , where ).
Case 4: .
Remark 2.3 It follows from (12) that
so that . Also we observe that on diagonalization we have (in all the above 4 cases),
where , are the eigenvalues of . Hence on each mesh interval,
is a solution of
which remains in the Lie Group, .
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)).
The left boundary conditions are implemented by fixing initial conditions for two solutions and of (1) at , namely we define solutions and at by requiring
Then the corresponding solutions of (1) automatically satisfy the boundary conditions (43) at . Using the solutions of (1) defined by (49) we define the matrices
The two-dimensional subspace,
where and are solutions of (1) corresponding to the and solutions of the IVP, (1) and (49), of the four-dimensional solution space of (1) satisfies the boundary conditions (43) at .
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 (43) and boundary conditions at (44).
Theorem 2.1 1) A N.S.C. for to be an eigenvalue of the SL problem for Equation (1) with boundary conditions at (43) and boundary conditions at (44), and having multiplicity one, is:
2) A N.S.C. for to be an eigenvalue of the SL problem for Equation (1) with boundary conditions at (43) and boundary conditions at (44), and having multiplicity one, is:
Proof Let be the unique solutions of the 4th order Equation (1) which are defined by the initial conditions at and be the corresponding solutions of the Hamiltonian system (1):
By fixing (57) the two-dimensional space is fixed by the matrix (57). The constants in the (57) matrix were chosen to ensure that the boundary conditions at (43) was satisfied, so we know that
that is, that both and satisfy the boundary conditions at (43). Also, of course, the space of solutions spanned by these two solutions of (1) satisfies the boundary conditions at (43), that is
for all . It remains only to apply the boundary conditions at (44), i.e. to require that
for boundary conditions at (44). Hence, we find
as the requirement for solutions in to also satisfy the boundary conditions at (44). But the Equations (64) and (65) can have a nonzero solution for if and only if
Hence (66) is a N.S.C. condition for 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 and both boundary conditions at . Since the dimension of is two, this can be at most two. For multiplicity one, we must have
and in this case the solution of the two Equations (64) and (65) will be
In this case the eigenfunction, which is unique up to a constant multiple, is
For multiplicity two, we need to require
and in this case and would be linearly independent eigenfunctions of (1).
More generally, the general case of boundary conditions (3) would give
as a N.S.C. for 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 transfer matrix
by doing the following steps:
1) Calculate the eigenvalues and the eigenvectors of .
where P denotes the matrix of eigenvectors of , and D denotes the diagonal matrix of the four eigenvalues of (12).
This gives a transfer matrix with 4 cases of eigenvalues of . (The matrix elements could be reduced to expressions involving sinh, cosh, sin and cos functions). Our MG4 code implements the above transfer matrix 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 , we integrate the IVP, (1) and (49), from to using the MG4 method on the matrix . Then when has been computed, the boundary conditions (44),
will be satisfied for some choices of real constants A and B, not both zero, if and only if
The computation is performed using an initial uniform mesh, applying bisection method with initial upper and lower bounds for a given eigenvalue , and then doubling the number of mesh points by bisecting the mesh to generate a Richardson h4-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 h4-Richardson’s Extrapolation
If , we can assume that if MG4 method is applied we will have for each choice of h:
for some choices of real constants . For each
in Neville’s algorithm (  , 22.214.171.124b, p. 42),
where we have taken
Applying Neville’s algorithm generates the h4-Richardson’s extrapolation table for the eigenvalue computation. Defining
we have from (13) that
here the second term in (15) is the extrapolation error. The first column of the extrapolation table, that is, the eigenvalues , are computable quantities. The columns two, three, four, five, , 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  , Section 3.2, pp. 461-462), we applied the change to “invariant-imbedding” variables, and , in our code in order to provide a good stable integration scheme. We generated the “invariant-imbedding” variables by using the matrices and 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 (  pp. 478-481). For example, we consider in section 5 the following 2nd order SL problems
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−12. Outputs of the SLEDGE code (Sturm-Liouville Estimates Determined by Global Error Control) of Pruess and Fulton  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 in columns 4 and 6 of Tables 1-5 was computed as
where are the SLEDGE-squared eigenvalues and 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 for the eigenvalues , and obtained by the SLEUTH are slightly larger than the error characterization for the eigenvalues , and 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 for the eigenvalues , and obtained by the SLEUTH are slightly larger than the error characterization for the eigenvalues , and achieved by the MG4 method. In Table 3, we see that the error characterization for the eigenvalues and achieved by the MG4 method are quite comparable to the error characterization for the eigenvalues and obtained by the SLEUTH. But, the error characterization for the eigenvalue obtained by the SLEUTH is slightly larger than the error characterization for the eigenvalue achieved by the MG4 method. In Table 4, we see that the error characterization for the eigenvalue achieved by the MG4 method is quite comparable to the error characterization for the eigenvalue obtained by the SLEUTH. But, the error characterization for the eigenvalues and obtained by the SLEUTH are slightly larger than the error characterization for the eigenvalues and achieved by the MG4 method. In Table 5, we see that the error characterization for the eigenvalues , , and obtained by the SLEUTH are slightly larger than the error characterization for the eigenvalues , , and 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.
In this paper we have presented the MG4 algorithm of Iserles et al.   , 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  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.
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.