\documentstyle[twoside]{article}
\input amssym.def
\pagestyle{myheadings}
\markboth{\hfil Eigenvalue Computations \hfil EJDE--1995/05}%
{EJDE--1995/05\hfil H.I. Dwyer \& A. Zettl \hfil}
\begin{document}
\ifx\Box\undefined \newcommand{\Box}{\diamondsuit}\fi
\title{\vspace{-1in}\parbox{\linewidth}{\footnotesize\noindent
{\sc  Electronic Journal of Differential Equations}\newline
Vol. {\bf 1995}(1995), No. 05, pp. 1-13. Published May 2, 1995.\newline
ISSN 1072-6691. URL: http://ejde.math.swt.edu or http://ejde.math.unt.edu
\newline ftp (login: ftp) 147.26.103.110 or 129.120.3.113 }
 \vspace{\bigskipamount} \\
Eigenvalue Computations for Regular Matrix Sturm-Liouville Problems
\thanks{ {\em 1991 Mathematics Subject Classifications:}
34B24, 34L15, 34L05.\newline\indent
{\em Key words and phrases:} Sturm-Liouville problem,
 matrix coefficients, eigenvalues, \newline\indent
 numerical computation of  eigenvalues
\newline\indent
\copyright 1995 Southwest Texas State University  and University of
North Texas.\newline\indent
Submitted: August 25, 1994.} }
\date{}
\author{H.I. Dwyer \& A. Zettl }
\maketitle
%
\begin{abstract}
An algorithm is presented for computing eigenvalues
of regular self-adjoint
Sturm-Liouville~(SL) problems with matrix coefficients and
separated boundary conditions.
\end{abstract}
%
\newcommand{\ddx}{\frac{d}{dx}}   % ordinary deriv wrt x
\newcommand{\ddt}{\frac{d}{dt}}   % ordinary deriv wrt t
\newcommand{\pdx}{\frac{\partial}{\partial x}} % partial deriv wrt x
\newtheorem{thm}{Theorem}
\newtheorem{lemma}[thm]{Lemma}
\newtheorem{defn}{Definition}
\newtheorem{ex}{Example}
% Main body of the paper
%
\section*{Introduction}
In his book {\em Discrete and Continuous Boundary  
Problems}~\cite{atkinson},
F.V.~Atkinson characterizes eigenvalues of matrix Sturm-Liouville  
(SL) problems
in terms of eigenvalues of certain unitary matrices.
Based on this characterization the team of Atkinson, Krall, Leaf,
and Zettl~(AKLZ) developed an algorithm, and constructed a prototype
FORTRAN code for the numerical computation of eigenvalues of SL  
problems.
%
A description of the algorithm, with some test results,
was published in an
Argonne Laboratory Report~\cite{aklz87}.
%
While clearly demonstrating the feasibility
of this algorithm, the tests showed that there were
some difficulties remaining.
%
Dwyer, in his dissertation as a student of Zettl,
refined and significantly extended the algorithm, and corrected
the difficulties in the original code.

In this paper, the completed algorithm for computing the eigenvalues
for a regular Sturm-Liouville problem with separated boundary
conditions is presented.
%
%
%  PROBLEM DEFINITIONS
%
%
\section*{Problem Definitions}
A SL problem consists of the second order linear
ordinary differential equation
\begin{equation} \label{ode}
   -(py')'+qy = \lambda wy \,\,\,\mbox{on} \,\,\,(a,b)
\end{equation}
together with boundary conditions.
For the case when both endpoints $a$, $b$ are regular, these have
the form
\begin{equation} \label{bc1}
 C \left(\begin{array}{c} y(a) \\ (py')(a) \end{array}\right) +
 D \left(\begin{array}{c} y(b) \\ (py')(b) \end{array}\right)
 = \left(\begin{array}{c} 0 \\ 0 \end{array}\right) .
\end{equation}
%
% Theorem
%
\begin{thm} \label{thm1}
  Let $p$, $q$, $w$ be complex $m\,\times\, m$ matrix functions defined
  on the interval $[a,b]$, $-\infty < a<b<\infty$, satisfying the
  following conditions:
  \begin{enumerate}
    \item The matrix $p(t)$ is invertible for almost all $t\in  
(a,b)$.
    \item Each component of $p^{-1}$, $q$, $w$ is in $L^{1}(a,b)$.
    \item The matrices $p(t)$, $q(t)$, $w(t)$ are hermitian for
          almost all $t\in (a,b)$.
    \item The matrices $p(t)$ and $w(t)$ are positive definite
          for almost all $t\in (a,b)$.
  \end{enumerate}
  Assume that the complex constant $2m\,\times\, 2m$ matrices $C$ and $D$  
satisfy:
  \begin{enumerate} \setcounter{enumi}{4}
    \item The $2m\,\times\, 4m$ matrix $(C|D)$ has full rank.
    \item $CJC^{*}=DJD^{*}$, where
          $J=\left(\begin{array}{c|c} 0 & I_{m} \\
                   \hline -I_{m} & 0 \end{array}\right)$ and $I_{m}$  
denotes
          the $m\,\times\, m$ identity matrix.
  \end{enumerate}
  Then the boundary value problem, (\ref{ode}),(\ref{bc1}) is  
self-adjoint;
  its eigenvalues are all real, there are countably many of them
  $\{ \lambda_{n}:n\in N_{0}=\{0,1,2,\cdots\}\}$, and they can be  
ordered
  to satisfy
  \[
     -\infty<\lambda_{0} \leq \lambda_{1} \leq \lambda_{2} \leq  
\cdots \,\, ,
      \,\,\mbox{ with }\,\,
      \lambda_{n} \rightarrow \infty \,\,\mbox{ as }\,\, n  
\rightarrow \infty .
  \]
\end{thm}
\paragraph{Proof.}  See M\"oller and Zettl \cite{moze94}.  $\Box$

If condition (4) is weakened so that $p(t)$ is required to be  
invertible, but is
not required to be positive definite, then a similar result holds,  
except that
the spectrum will not, in general, be bounded below.

We are interested in the special case when the boundary conditions  
are
{\em separated}. This means that the matrices $C$ and $D$ can
each be partitioned
into four $m\,\times\, m$ square matrices, as:
\[  C = \left(\begin{array}{c|c} A_{1} & A_{2}\\
                \hline 0 & 0 \end{array}\right)
    \,\,\, , \,\,\,
    D = \left(\begin{array}{c|c} 0 & 0 \\
                \hline B_{1} & B_{2} \end{array}\right).
\]
In this case, the boundary conditions (2) reduce to
\begin{eqnarray}
    A_{1} y(a) + A_{2} (py'(a)) & = & 0 \label{bca} \\
    B_{1} y(b) + B_{2} (py'(b)) & = & 0 \label{bcb}
\end{eqnarray}
where $A_{1}$, $A_{2}$, $B_{1}$, and $B_{2}$ are complex
constant $m\,\times\, m$ matrices which satisfy
\begin{enumerate}
  \item $rank(A_{1}|A_{2}) = rank(B_{1}|B_{2}) = m$
  \item $A_{1}A_{2}^{*}=A_{2}A_{1}^{*}$
  \item $B_{1}B_{2}^{*}=B_{2}B_{1}^{*}$.
\end{enumerate}
Below we consider the question: How can the eigenvalues $\lambda_{n}$
be computed numerically?
The codes SLEDGE \cite{fupr94}, and SLEIGN  
\cite{bail78},\cite{bags78}
, as well as the codes in the NAG library \cite{naglib},
address this question only for the special case when $m=1$ (scalar  
case).
%
%
%  SHOOTING METHODS FOR SCALAR PROBLEMS
%
%
\section*{Shooting Methods for Scalar Problems}
Before beginning on the matrix problem, we consider the scalar
problem.  The boundary conditions specify the ratio of the values
of the solution $y$ and the quasi-derivative~$py'$ at each of the
endpoints.  We define
\begin{equation}
  z(x)=\frac{y(x)}{(py')(x)}. \label{defz}
\end{equation}
Since the solution is nontrivial, $y$ and $py'$ cannot be  
simultaneously
zero.  However, it is quite likely that $py'$ will be zero at some  
points
along the interval, leading to unbounded values for the ratio~$z$.
To avoid the difficulties of unbounded values, we might use a mapping
which transforms the real line into some bounded set; applying the
mapping to the ratio~$z$ will produce a variable whose initial
value, at endpoint~$a$, is also determined by the boundary condition,
and which avoids the difficulties caused by the division by zero.

If we use the inverse-tangent mapping, the result is the
Pr\"ufer transformation (see~\cite{coddington}).  The ``Pr\"ufer
angle'' is defined by
\[
   \Theta(x) = \arctan(z) = \arctan\left(\frac{y(x)}{(py')(x)}\right) .
\]
This variable satisfies a certain
first order ordinary differential equation.
Using this equation, and the initial value determined by the boundary
condition at the endpoint~$a$, it is possible to compute a numerical
approximation to the resulting value for $\Theta$ at the other
endpoint,~$b$.
Note that the initial value of $\Theta$ depends only on the
boundary condition, but at other points in the interval the value
will depend on the parameter $\lambda$.
The boundary condition at $b$ determines  ``target''
values; if the value of $\Theta$ matches the target, we have found an
eigenvalue.
The success of the search algorithm depends on our ability to do  
three
things:
\begin{enumerate}
  \item Recognize that the value $\lambda$ is an eigenvalue.
  \item Determine whether a trial value of $\lambda$ is ``near'' an
        eigenvalue, and whether it is too high or too low.
  \item Determine {\em which} eigenvalue the trial value is close to.
\end{enumerate}
With appropriate care about which branch of the inverse tangent is  
used,
it can be shown that if $p(t)$ is positive for all $t$ in the  
interval,
then the value $\Theta(b;\lambda)$ is monotone increasing
in $\lambda$, which allows a simple bisection search to be used to
approximate the eigenvalue.  More advanced algorithms
(see~\cite{bgkz91}) use more efficient search
methods, such as Newton's method.

Another possibility is to use a mapping from complex analysis to map  
the real line
onto some bounded curve in the complex plane. The mapping
\[
    \Theta = \frac{1+iz}{1-iz}
\]
takes the real axis onto the unit circle, with $z=\pm \infty$ mapped
to $\Theta = -1$.
This leads somewhat naturally to
\begin{eqnarray*}
  \Theta & = & \frac{1+iz}{1-iz} \\
         & = & \frac{1+i\frac{y}{py'}}{1-i\frac{y}{py'}} \\
         & = & \frac{py'+iy}{py'-iy} \\
         & = & (py'+iy)(py'-iy)^{-1}.
\end{eqnarray*}
Of course, even in the scalar case, and particularly in the matrix  
case,
there are some details which must be addressed; the discussion  
above is
meant only as a motivation for what follows.
%
%
% A SHOOTING METHOD FOR MATRIX PROBLEMS
%
%
\section*{A Shooting Method for Matrix Problems}
Although it is possible to define an inverse-tangent function for
matrix variables, the result is not particularly useful to us here.
The difficulty is that the corresponding sine and cosine
functions do not have
the desirable derivative properties
of their scalar counterparts, introducing difficulties in the
formulation of a first order differential equation.

Instead, we will develop a method based on the second approach
mentioned above.
Let $Y(x)$ be a square matrix function such that
\[ -(pY')'+qY=\lambda wY \]
satisfying the initial conditions at endpoint $a$:
\[ Y(a)=-A_{2}^{*} \,\,\, , \,\,\, (pY')(a)=A_{1}^{*} . \]
The assumptions on the coefficients guarantee that such a solution  
exists, and
is defined on the entire interval $[a,b]$.

We can write an equivalent system of first order equations as  
follows:
let $U=Y$, and $V=pY'$, and we have
\begin{eqnarray*}
  U' & = & p^{-1}V \\
  V' & = & rU \,\, , \,\, \mbox{where} \\
  r(x,\lambda) & = & q(x) - \lambda w(x)
\end{eqnarray*}
subject to:
\begin{eqnarray*}
  U(a) & = & -A_{2}^{*} \\
  V(a) & = & A_{1}^{*}.
\end{eqnarray*}
We are now ready to define the matrix Atkinson \cite{atkinson} calls  
theta.
%
% definition
%
\begin{defn}
  For each $x\in [a,b]$, and for each real $\lambda$, define
  \begin{equation} \label{eqntheta}
     \Theta (x,\lambda)=(V(x)+iU(x))(V(x)-iU(x))^{-1}
  \end{equation}
  for $U$ and $V$ as defined above.
\end{defn}
%
%
The dependence on $\lambda$ arises because $U$ and $V$ are solutions  
to a system
whose coefficients depend on $\lambda$.
%
Having defined the matrix, we must verify that it will exist for a SL  
problem
with separated, self-adjoint boundary conditions.  We will need the  
following
lemmas.
%
% Lemma
%
\begin{lemma} \label{inverts}
  Suppose that $A_{1}$ and $A_{2}$ are complex, square $m\,\times\, m$  
matrices such
  that
  \begin{enumerate}
    \item $rank(A_{1}|A_{2})=m$
    \item $A_{1}A_{2}^{*}=A_{2}A_{1}^{*}$
  \end{enumerate}
  then the matrix $(A_{1}-iA_{2})$ is invertible.
\end{lemma}
  \paragraph{Proof.} We begin by noting that $(A_{1}-iA_{2})$ is invertible if,
  and only if,\newline
  $(A_{1}-iA_{2})(A_{1}-iA_{2})^*$ is invertible.
  \begin{eqnarray*}
     (A_{1}-iA_{2})(A_{1}-iA_{2})^{*}
     & = & (A_{1}-iA_{2})(A_{1}^{*}+iA_{2}^{*}) \\
     & = &  
A_{1}A_{1}^{*}+A_{2}A_{2}^{*}-i(A_{2}A_{1}^{*}-A_{1}A_{2}^{*}) \\
     & = & A_{1}A_{1}^{*}+A_{2}A_{2}^{*}.
  \end{eqnarray*}
  Assume that the sum $A_{1}A_{1}^{*}+A_{2}A_{2}^{*}$ is singular.
  Then there is some nonzero vector $v$ such that
  \begin{eqnarray*}
     (A_{1}A_{1}^{*}+A_{2}A_{2}^{*})v & = & 0 \\
     v^{*}(A_{1}A_{1}^{*}+A_{2}A_{2}^{*})v & = & 0 \\
     v^{*}A_{1}A_{1}^{*}v+v^{*}A_{2}A_{2}^{*}v & = & 0
  \end{eqnarray*}
  which is possible only if $v^{*}A_{1}=0$ and $v^{*}A_{2}=0$.
  However, this would imply that
  \[
      v^{*}(A_{1}|A_{2})=0
  \]
  which is a contradiction since $(A_{1}|A_{2})$ is of full rank, by
  assumption, and $v^{*}\neq 0$.  We conclude that no such $v$  
exists,
  so the sum is nonsingular and $(A_{1}-iA_{2})$ is invertible. $\Box$
%
% Lemma
%
\begin{lemma} (\cite{atkinson}) \label{uvherm}
  For $U$, $V$ as defined above, if $U^{*}V$ is hermitian at any
  $x_{0}\in [a,b]$, then $U^{*}V$ is hermitian for all $x\in [a,b]$.
\end{lemma}
  \paragraph{Proof.} See \cite{atkinson}.   $\Box$
%
% Lemma
%
\begin{lemma} (\cite{atkinson}) \label{uvinverts}
  For $U$, $V$ as defined above, suppose that $U^{*}V$ is hermitian
  for all $x\in [a,b]$, and suppose that $(V-iU)$ is invertible at
  some $x_{0}\in [a,b]$.  Then $(V-iU)$ is invertible for all $x\in  
[a,b]$.
\end{lemma}
  \paragraph{Proof.} See \cite{atkinson}.  $\Box$
%
% Theorem
%
\begin{thm} \label{thetaexists}
  For a regular SL problem with separated, self-adjoint boundary  
conditions,
  the matrix $\Theta(x,\lambda)$ exists for all $x\in [a,b]$, for any  
real
  $\lambda$.
\end{thm}
  \paragraph{Proof.} At the point $x_{0}=a$, we have  
$U(x_{0})^{*}V(x_{0})=-A_{2}A_{1}^{*}$.
  By assumption, this matrix is hermitian.  Applying  
Lemma~\ref{uvherm}, we
  see that $U^{*}V$ is hermitian for all $x\in [a,b]$.  At the point
  $x_{0}=a$, we have $V(x_{0})-iU(x_{0})=(A_{1}-iA_{2})^{*}$, which  
is
  invertible by Lemma~\ref{inverts}.  Finally, by  
Lemma~\ref{uvinverts},
  we see that $(V-iU)$ is invertible for all $x\in [a,b]$, so we  
conclude that
  the matrix $\Theta$ may be computed for every $x\in [a,b]$.  $\Box$

It is easy to verify that the matrix $\Theta(x,\lambda)$ is unitary  
for all
$x\in [a,b]$, for any real value of $\lambda$.

The boundary condition at the endpoint $a$ provides an initial value  
for theta;
in order to compute theta at other points along the interval, we  
develop a
first order differential equation which has theta as a solution.
%
%  theorem
%
\begin{thm} (\cite{atkinson}) \label{thm2}
  For $\Theta(x,\lambda)$ defined above, we have for any real  
$\lambda$
  \[ \pdx \Theta(x,\lambda) = i\Theta(x,\lambda)\Omega(x,\lambda) \]
  where $\Omega(x,\lambda)$ is defined by
  \[ \Omega(x,\lambda)=\frac{1}{2}((\Theta+I)^{*}p^{-1}(x)(\Theta+I)-
     (\Theta-I)^{*}r(x,\lambda)(\Theta-I)). \]
\end{thm}
  \paragraph{Proof.} See \cite{atkinson}.  $\Box$

The first requirement for a method is that we must be able to  
recognize
that a value~$\lambda$ is an eigenvalue.  It is convenient to define  
a
new matrix, constructed by ``scaling'' the matrix $\Theta(x,\lambda)$
by a constant.
%
% definition
%
\begin{defn} \label{defb}
  For $\Theta$ as defined above, let
  \[ B(x,\lambda)=(B_{1}-iB_{2})^{-1}(B_{1}+iB_{2})\Theta(x,\lambda)  
\]
  where $B_{1}$ and $B_{2}$ are the constant matrices used to define  
the
  boundary condition at the endpoint $b$.
\end{defn}
%
%
Note that the factor $(B_{1}-iB_{2})^{-1}$ exists, by  
Lemma~\ref{inverts}.
It is easy to verify that the matrix $B$ is also unitary.
%
% Lemma
%
\begin{lemma} (\cite{aklz87}) \label{nullity}
  For $U$, $V$ as defined above, we have:
  \begin{enumerate}
     \item A real value $\lambda$ is an eigenvalue of the SL~problem
           if, and only if, the matrix \mbox{$B_{1}U(b)+B_{2}V(b)$}
           is singular.
     \item The multiplicity of $\lambda$ as an eigenvalue of the  
SL~problem
           is equal to the dimension of the null space of
           \mbox{$B_{1}U(b)+B_{2}V(b)$}.
  \end{enumerate}
\end{lemma}  
  \paragraph{Proof.} See \cite{aklz87}.  $\Box$

% theorem
%
\begin{thm}
  For $B(x,\lambda)$ as defined above, we have
  \begin{enumerate}
     \item A real value $\lambda$ is an eigenvalue for the SL~problem
        if, and only if, the number 1 is an eigenvalue of the matrix
        $B(b,\lambda)$.
     \item The multiplicity of $\lambda$ as an eigenvalue of the  
SL~problem
        is equal to the multiplicity of the value 1 as an
        eigenvalue of the matrix $B(b,\lambda)$.
  \end{enumerate}
\end{thm}  
  \paragraph{Proof.}  (1) Suppose that 1 is an eigenvalue of $B(b,\lambda)$.
  Then there is a nonzero
  constant vector $c$ such that $B(b,\lambda)c=1c$.  Hence
  \[
    (B_{1}-iB_{2})^{-1}(B_{1}+iB_{2})(V(b)+iU(b))(V(b)-iU(b))^{-1}c=c  . 
  \]
  Define $z=(V(b)-iU(b))^{-1}c$, and note that $z\neq 0$ and
  $c=(V(b)-iU(b))z$.  Hence
  \begin{eqnarray*}
    (B_{1}-iB_{2})^{-1}(B_{1}+iB_{2})(V(b)+iU(b))z & = &  
(V(b)-iU(b))z \\
    (B_{1}+iB_{2})(V(b)+iU(b))z & = & (B_{1}-iB_{2})(V(b)-iU(b))z \\
    2i(B_{1}U(b)+B_{2}V(b))z & = & 0 .
  \end{eqnarray*}
  Since we know $z\neq 0$, we conclude that $B_{1}U(b)+B_{2}V(b)$ is  
singular,
  and the result follows from Lemma~\ref{nullity}.  Conversely, if we  
know that
  $\lambda$ is an eigenvalue of the SL~porblem, we reverse the steps  
to produce
  an eigenvector for $B(b,\lambda)$ corresponding to an eigenvalue of  
1.

  \paragraph{Proof.} (2) The multiplicity of 1 is the number of linearly independent  
choices for the
  vector $c$, and hence the number of linearly independent choices  
for the
  vector $z$, which is the dimension of the null space of
  \mbox{$B_{1}U(b)+B_{2}V(b)$}.  The result is now immediate from
  Lemma~\ref{nullity}. $\Box$

%
The second requirement for a method is that we must be able to  
determine
if a trial value for $\lambda$ is ``close'' to an actual eigenvalue,  
and
whether it is too high or too low.  To do this, we will need the  
following
lemma, which is proved by Atkinson:
%
% lemma
%
\begin{lemma} (\cite{atkinson}) \label{lemma1}
  Suppose that $t$ is a real variable, $\Theta(t)$ is unitary,
  $\Omega(t)$ is hermitian, continuous, and positive definite, and  
suppose
  that $\Theta$ satisfies
  \[ \ddt \Theta=i\Theta\Omega \]
  then as $t$ increases the eigenvalues of $\Theta$ move  
anti-clockwise
  around the complex unit circle.
\end{lemma}
  \paragraph{Proof.} See \cite{atkinson}.  $\Box$

It can be shown (see \cite{atkinson}) that, for any fixed
choice of $x$, the matrix $\Theta(x,\lambda)$ satisfies a (partial)
differential equation in the variable $\lambda$
of the form required by Lemma~\ref{lemma1}.  Since the
matrix $B$ is simply $\Theta$ multiplied by a constant, the same
conclusions hold for $B$; for $x=b$ we have the eigenvalues of  
$B(b,\lambda)$
moving anti-clockwise around the complex unit circle as the
parameter $\lambda$ is increased.  For any trial value of~$\lambda$,  
we
examine the position of the eigenvalues of $B(b,\lambda)$ and we can  
readily
determine if the trial value is ``close'', and whether it is too high  
or
too low.  A simple bisection search procedure will quickly determine  
the
required eigenvalue to any precision we might require.

The third requirement for a method, that we can determine {\em which}
eigenvalue we are close to, is more difficult to satisfy.
We will discuss this issue in the next section.
%
%
%   IMPLEMENTATION DETAILS
%
%
\section*{Implementation Details}
The team of Atkinson, Krall, Leaf, and Zettl (AKLZ) developed a  
prototype
\mbox{FORTRAN} code \cite{aklz87}, based on the algorithm described  
above.
This early code produced excellent results for some trial problems,  
but
had difficulties in certain cases.
A later implementation by Dwyer duplicated the earlier
results, and the difficulties as well.  These difficulties were  
caused by
some details in the way that the algorithm was implemented,
and have since been corrected; the algorithm
itself is perfectly sound.
A corrected FORTRAN code, developed by Dwyer and Zettl, is available  
as
an appendix to the doctoral thesis of Dwyer \cite{dwye93}.

The initial value problem (IVP) for $\Theta$ may be solved using any  
of the
standard techniques.  It is necessary to keep track
of the eigenvalues of the matrix $B(x,\lambda)$ as the variable $x$  
increases
from $a$ to $b$. The numerical method for the IVP solution uses a  
step-size
which is rather small, to increase the accuracy of the computed  
result.
However, since eigenvalue computations are relatively slow, we do not  
wish to
examine the eigenvalues at every step; the code is written to perform  
the
eigenvalue computations only at selected steps.  This modification  
saves
considerable execution time.  The algorithm has been implemented  
using both
the NAG and IMSL linear algebra libraries, with similar results.

The matrix $B(x,\lambda)$ is sampled at equally spaced points along  
the
interval $[a,b]$.  At each point, the eigenvalues of $B$ are computed  
using
a library routine.  We compare the lists of eigenvalues for adjacent  
samples
of $B$ to determine whether the eigenvalues have passed through the  
value~1
on the complex unit circle.  A difficulty arises in this comparison  
if the
two lists for adjacent samples are not returned by the library  
routine in the
same order.  At various points along the interval, the matrix $B$  
will change
enough that the library routine will ``shuffle'' the eigenvalue list;  
the list
will contain eigenvalues which are similar to the previous list, but  
in
a different order.  This re-ordering cannot be avoided by reducing  
the
distance between samples.  It is necessary to correctly match the two  
lists,
before the eigenvalue movements can be tracked.  Since the lists are  
usually
{\em not} shuffled, the matching routine should be designed to  
quickly
recognize this situation.  An algorithm for matching the lists is  
described
in detail in \cite{dwye93}.

The movements of the eigenvalues of $B(x,\lambda)$ are tracked as $x$  
increases
from~$a$ to~$b$.
We keep track of the total number of times that an eigenvalue passes
anti-clockwise through the value~1; we will call this total the {\em  
index}
for the particular choice of~$\lambda$.  As we increase the~$\lambda$  
value, we
will get a ``jump'' in the index value as $\lambda$ passes through an  
eigenvalue
of the SL problem.  The size of the jump indicates the multiplicity  
of
the eigenvalue.
The first (lowest) eigenvalue for the SL problem corresponds to the  
first
jump in the index.  However, the first jump need not be a jump from
index=0 to index=1.
The eigenvalues of $B(x,\lambda)$ do not, in general,
move anti-clockwise as $x$ increases; each
may behave differently, and may move clockwise during some portion
of the interval $[a,b]$.  Thus, it is possible to have a negative  
index value
in some problems.  Other problems may have an index which is never  
zero, since
it is possible that the first eigenvalue may correspond to the jump  
from
index=5 to index=6.

If a lower bound on the eigenvalues is known, then we may solve the  
labeling
difficulty by computing the index corresponding to the lower bound  
and using
this value as an offset.  We may then compute any desired eigenvalue,
$\lambda_{n}$, by searching for the index jump from index=offset+n to
index=offset+n+1.  We continue a bisection search until the  
$\lambda$-interval
$[\lambda_{lo},\lambda_{hi}]$ is sufficiently narrow to meet the user
specified tolerance, where $\lambda_{lo}$ is a value whose index is  
less
than or equal to offset+n, and~$\lambda_{hi}$ is a value whose index  
is
greater than offset+n.
In cases where the spectrum is not bounded below, we use the  
convention
that $\lambda_{0}$ is the first nonnegative eigenvalue. This is  
achieved by
computing the index corresponding to $\lambda=0$, and using that  
value as
offset.

After a bracketing interval has been found using bisection, the  
multiplicity
of the eigenvalue is simply the size of the jump in the index across  
the bracket.
If the user has specified a very tight or
a very loose tolerance on the eigenvalue, the
multiplicity which is computed may be incorrect.
What is actually being computed is the total number of eigenvalues  
which lie
within the bracket.  If the tolerance is larger than the spacing  
between
adjacent eigenvalues, they may both be included within the same  
bracket and
appear as one multiple eigenvalue.  The reported multiplicity will  
thus be
too large.
Conversely, if the tolerance is too tight,
the reported multiplicity may be too low;
due to numerical errors,
an eigenvalue of multiplicity may appear to be a cluster of simple  
eigenvalues,
very close together.  If the tolerance is tighter than the spacing of  
this cluster,
or if the cluster is very near to one of the endpoints of the  
bracket,
some of the simple eigenvalues will be excluded from the bracket,  
resulting in
a multiplicity count which is too low.

Experiments indicate that the best method is to work with a moderate  
tolerance
for initial computations, determine the multiplicity using that  
tolerance, and
then continue the computation to a tighter tolerance.  If the  
multiplicities
do not match, we are alerted that some additional care is required.
The next example was selected to illustrate some of these  
difficulties.
%
%
%      EXAMPLE
%
%
\section*{Example}
We consider the following three-dimensional problem:
%
\begin{ex}
  On the interval $[0,\pi]$, with coefficients
  \[
  P(x)=\left(\begin{array}{ccc} 11 & 6 & 3 \\
                                 6 & 12 & 2 \\
                                 3 & 2 & 1 \end{array}\right)
  \,\, , \,\,
  Q(x)=\left(\begin{array}{ccc} 0 & 0 & 0 \\
                                0 & 0 & 0 \\
                                0 & 0 & 0 \end{array}\right)
  \,\, , \,\,
  W(x)=\left(\begin{array}{ccc} 38 & 24 & 12 \\
                                24 & 18 & 8 \\
                                12 & 8 & 4 \end{array}\right)
  \]
  subject to the condition $Y(0)=Y(\pi)=0$.
\end{ex}
%                                         

It is easy to verify that both $P$ and $W$ are positive definite.
The problem may be
solved without difficulty, so that we know exactly what the  
eigenvalues are;
on this interval they are either integers or integers plus one  
fourth.
The example was constructed so that there are eigenvalues with each  
of
the possible multiplicities.
The first ten eigenvalues are given in Table~1.
%
%
% TABLE 1
%
%
\begin{table}[hbt]
 \begin{center}
  \begin{tabular}{|c|r||c|r|}
    \hline n & $\lambda_{n}$ & n & $\lambda_{n}$ \\
    \hline 0 & 0.25 & 5 & 4.00 \\
    \hline 1 & 1.00 & 6 & 4.00 \\
    \hline 2 & 1.00 & 7 & 6.25 \\
    \hline 3 & 2.25 & 8 & 9.00 \\
    \hline 4 & 4.00 & 9 & 9.00 \\
    \hline
  \end{tabular}
 \caption{Exact Values For Eigenvalues}
 \end{center}
\end{table}
The accuracy which can be achieved
will, of course, depend on
the particular numerical methods used.  A relatively straight-forward
implementation of the algorithm in FORTRAN, using double precision
arithmetic, produced the results shown in Table~2.  The bisection
search produces an interval estimate; the midpoint of the interval is
used as the estimate for the eigenvalue.
Each search began with a random real value.
Each interval was refined until the relative error is less than 0.1\%  
. 

%
%
% TABLE 2
%
%
\begin{table}[hbt]
 \begin{center}
  \begin{tabular}{|c|c|c|c|}
    \hline n & Interval & Midpoint & Mult. \\
    \hline 0 & ( 0.24991 , 0.250968 ) & 0.250479 & 1 \\
    \hline 1 & ( 0.999991 , 1.000968 ) & 1.000479 & 2 \\
    \hline 2 & " & " & \\
    \hline 3 & ( 2.249014 , 2.250968 ) & 2.249991 & 1 \\
    \hline 4 & ( 3.997061 , 4.000968 ) & 3.999014 & 3 \\
    \hline 5 & " & " & \\
    \hline 6 & " & " & \\
    \hline 7 & ( 6.247061 , 6.250968 ) & 6.249014 & 1 \\
    \hline 8 & ( 8.993155 , 9.000968 ) & 8.997061 & 2 \\
    \hline 9 & " & " & \\
    \hline
  \end{tabular}
 \caption{Results, Beginning With Random Real Guess Values}
 \end{center}
\end{table}
If the computations are repeated using an integer as the initial
guess, each eigenvalue will be one of the endpoints
of the bisection interval.
The results of these computations are shown in Table~3.
%
%
% TABLE 3
%
%
\begin{table}[hbt]
 \begin{center}
  \begin{tabular}{|c|c|c|c|}
    \hline n & Interval & Midpoint & Mult. \\
    \hline 0 & ( 0.250000 , 0.250977 ) & 0.250488 & 1 \\
    \hline 1 & ( 0.999023 , 1.000000 ) & 0.999512 & 1 \\
    \hline 2 & ( 1.000000 , 1.000977 ) & 1.000488 & 1 \\
    \hline 3 & ( 2.250000 , 2.251953 ) & 2.250977 & 1 \\
    \hline 4 & ( 3.996094 , 4.000000 ) & 3.998047 & 1 \\
    \hline 5 & ( 4.000000 , 4.003906 ) & 4.001953 & 2 \\
    \hline 6 & " & " & \\
    \hline 7 & ( 6.250000 , 6.253906 ) & 6.251953 & 1 \\
    \hline 8 & ( 8.992188 , 9.000000 ) & 8.996094 & 1 \\
    \hline 9 & ( 9.000000 , 9.007813 ) & 9.003906 & 1 \\
    \hline
  \end{tabular}
 \caption{Results, Beginning With Integer Guess Values}
 \end{center}
\end{table}
Even a small numerical error will split the
multiple eigenvalues across the endpoint, resulting in
errors in the reported multiplicities.
Note that each reported eigenvalue {\em is} within tolerance of
the correct value; it is only the reported multiplicity which is
in error.  Similar errors appear when the tolerance is
specified to be very small.

As a final test, we computed the first and second eigenvalues,
beginning with a random real initial guess, and allowed the
bisection process to continue until an incorrect bisection step
occurred.  Beginning with an interval of length~1, the code
made~32 correct bisections when computing the first eigenvalue,
and~21 correct bisections before the second eigenvalue ``split''
and lost its multiplicity.  These results are shown in Table~4.
%
%
% TABLE 4
%
%
\begin{table}[hbt]
 \begin{center}
  \begin{tabular}{|c|c|c|}
    \hline n & Interval & Midpoint \\
    \hline 0 & ( 0.249999999889 , 0.250000000121 ) &
                 0.250000000005 \\
    \hline 1 & ( 0.999999728102 , 1.000000204939 ) &
                 0.999999966521 \\
    \hline
  \end{tabular}
 \caption{Results of Continued Bisection}
 \end{center}
\end{table}
These results indicate that the values of the eigenvalues may
be computed quite accurately using this algorithm.
However, when computations are continued this far the
reported multiplicity is almost always in error; multiple
eigenvalues are reported as a ``cluster'' of closely spaced simple
eigenvalues.
%
%
%  BIBLIO
%
%
\begin{thebibliography}{99}
\bibitem{atkinson}
  {\sc F.V. Atkinson},
  {\em Discrete and Continuous Boundary Problems},
  Academic Press, New York, 1964.

\bibitem{aklz87}
  {\sc F.V. Atkinson, A.M. Krall, G.K. Leaf, A. Zettl},
  {\em On the Numerical Computation of Eigenvalues of Matrix
    Sturm-Liouville Problems with Matrix Coefficients},
  Argonne National Laboratory Reports, Darien, 1987.
\bibitem{bail78}
  {\sc P.B. Bailey},
  {\em SLEIGN: An Eigenfunction-eigenvalue Code for
    Sturm-Liouville Problems},
  {SAND7-2044}, Sandia Laboratories, Albuquerque, 1978.
\bibitem{bags78}
  {\sc P.B. Bailey, M.K. Gordon, L.F. Shampine},
  {\em Automatic Solution of the Sturm-Liouville Problem},
  ACM Trans. Math. Software, 4, 1978, 193-208.
\bibitem{coddington}
  {\sc E.A. Coddington, N. Levinson},
  {\em Theory of Ordinary Differential Equations},
  Krieger Publishing, Malabar, 1984.
\bibitem{dwye93}
  {\sc H.I. Dwyer},
  {\em Eigenvalues of Matrix Sturm-Liouville Problems
        with Separated or Coupled Boundary Conditions},
  Doctoral Thesis, Northern Illinois University, 1993.
\bibitem{moze94}
  {\sc M. M\"oller, A. Zettl},
  {\em Semi-boundedness of Ordinary Differential Operators},
  Journal of Differential Equations, (to appear).
\bibitem{bgkz91}
  {\sc P.B.Bailey, B.S.Garbow, H.G.Kaper and A.Zettl},
  {\em Eigenvalue and Eigenfunction Computations for Sturm-Liouville  
Problems},
  ACM TOMS 17(1991), 491-499.
\bibitem{fupr94}
  {\sc S.Pruess and C.Fulton},
{\em Mathematical Software for Sturm-Liouville Problems},
  ACM TOMS (1994), (to appear).
\bibitem{naglib}
  {\sc J. Pryce},
  {\em D02KEF},
  NAG Library Reference Guide.
\end{thebibliography}
%
\par\medskip\noindent
{\sc H.I. Dwyer\newline
 Mathematics Department\newline
 University of Wisconsin\newline
 Platteville, Wi. 53818\newline}
 E-mail address:  dwyer@uwplatt.edu
\par\medskip\noindent
{\sc A. Zettl\newline
 Mathematical Sciences Dept.\newline
 Northern Illinois University \newline
 DeKalb, IL. 60115\newline}
 E-mail address:  zettl@math.niu.edu



\end{document}