\documentstyle[twoside]{article} \input amssym.def \pagestyle{myheadings} \markboth{\hfil Eigenvalues of regular Sturm-Liouville Problems \hfil EJDE--1994/06}% {EJDE--1994/06\hfil H.I. Dwyer \& A. Zettl \hfil} \begin{document} \title{\vspace{-1in}\parbox{\linewidth}{\footnotesize\noindent {\sc Electronic Journal of Differential Equations}, Vol.\ {\bf 1994}(1994), No.\ 06, pp. 1--10. \newline ISSN: 1072-6691. URL: http://ejde.math.swt.edu or http://ejde.math.swt.edu \newline ftp (login: ftp) 147.26.103.110 or 129.120.3.113} \vspace{\bigskipamount} \\ Computing Eigenvalues of Regular Sturm-Liouville Problems \thanks{ {\em 1991 Mathematics Subject Classifications:} 34B24, 34L15, 34L05.\newline\indent {\em Key words and phrases:} Sturm-Liouville problem, numerical computation of eigenvalues.\newline\indent \copyright 1994 Southwest Texas State University and University of North Texas.\newline\indent Submitted April 15, 1994. Published August 23, 1994.} } \date{} \author{ H.I. Dwyer \& A. Zettl} \maketitle % \newcommand{\ddx}{\frac{d}{dx}} % ordinary deriv wrt x % x symbol, used for matrix dimensions, inside math mode \newcommand{\by}{\,\times \,} \newtheorem{thm}{Theorem} \newtheorem{ex}{Example} % \begin{abstract} An algorithm is presented for computing eigenvalues of regular self-adjoint Sturm-Liouville problems with matrix coefficients and arbitrary coupled boundary conditions. \end{abstract} \section*{Introduction} SLEDGE \cite{fupr94}, SLEIGN \cite{bags78},\cite{bail78} (see also \cite{bgkz91}) and the NAG library code \cite{naglib} are all highly effective, state of the art {\em general purpose} codes for computing eigenvalues of regular scalar Sturm-Liouville (SL) problems with {\em separated} boundary conditions. Although there are special purpose codes (e.g. \cite{plum90}) for computing eigenvalues of regular scalar SL problems with periodic boundary conditions there seems to be no general purpose algorithm available for general coupled self-adjoint boundary conditions, even in the scalar case. The purpose of this paper is to present such an algorithm not only for the scalar case but also for SL problems with matrix coefficients and general self-adjoint coupled boundary conditions. Our method is based on a construction which, given an SL problem with a coupled boundary condition, constructs a higher dimensional problem with separated conditions which has exactly the same eigenvalues as the given coupled problem. The general purpose code for SL problems of arbitrary dimension with separated boundary conditions developed by Dwyer in \cite{dwye93} can then be used to compute these eigenvalues. The eigenfunctions of the original problem can also be recovered from the eigenfunctions of the constructed higher dimensional problem. It is worth noting that the eigenvalues and eigenfunctions of the new problem with separated conditions are {\em not} simply approximations for those of the original problem; the correspondence is exact. This may have theoretical implications as well as numerical ones. \section*{Problems and Theorems} A Sturm-Liouville (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} \begin{thm} \label{thm1} Let $p$, $q$, $w$ be complex $m\by m$ matrix functions defined on the interval $[a,b]$, $-\infty < a1$ and probably also not to the case $m=1$ when the boundary conditions are coupled. Our method consists in constructing a new problem of dimension $2m$ (i.e., $m\rightarrow 2m$) which has exactly the same eigenvalues and {\em separated} boundary conditions. On this new problem, the method due to Atkinson, Krall, Leaf, and Zettl \cite{aklz87} as extended and refined by Dwyer \cite{dwye93} can then be applied. \section*{Construction of an Equivalent Separated Problem} We begin with an m-dimensional self-adjoint problem, (\ref{ode}), (\ref{bc1}) satisfying the conditions of Theorem \ref{thm1}. \paragraph{Definition} \begin{enumerate} \item Let $P(x) = \left(\begin{array}{c|c} p(x) & 0 \\ \hline 0 & p((a+b)-x) \end{array}\right) \,\, \mbox{ for all }\,\, x\in [a,b]$ \item Let $Q(x) = \left(\begin{array}{c|c} q(x) & 0 \\ \hline 0 & q((a+b)-x) \end{array}\right) \,\, \mbox{ for all }\,\, x\in [a,b]$ \item Let $W(x) = \left(\begin{array}{c|c} w(x) & 0 \\ \hline 0 & w((a+b)-x) \end{array}\right) \,\, \mbox{ for all }\,\, x\in [a,b]$ \item Let $A_{1}=\left(\begin{array}{c|c} I & -I \\ \hline 0 & 0 \end{array}\right) \,\,\, , \,\,\, A_{2}=\left(\begin{array}{c|c} 0 & 0 \\ \hline I & I \end{array}\right)$ \item Let $B_{1}=\left(\begin{array}{c|c} -D_{2} & C_{2} \end{array}\right) \,\,\, , \,\,\, B_{2}=\left(\begin{array}{c|c} D_{1} & C_{1} \end{array}\right).$ \end{enumerate} Consider the problem on the interval $[c,b]$, where $c=(a+b)/2$, given by \begin{equation} \label{ode2} -(PY')'(x)+Q(x)Y(x)=\lambda W(x)Y(x) \end{equation} subject to: \begin{eqnarray} A_{1}Y(c)+A_{2}(PY')(c) & = & 0 \label{bc2a}\\ B_{1}Y(b)+B_{2}(PY')(b) & = & 0. \label{bc2b} \end{eqnarray} \begin{thm} \label{thm2} Let $m$ be a positive integer, and let Problem~1 be a boundary value problem of dimension~m defined by equations~(\ref{ode}),~(\ref{bc1}) which satisfies the conditions of Theorem~\ref{thm1}. Let Problem~2 be the boundary value problem defined by (6),(7),and (8). Then Problem~2 is a regular self-adjoint Sturm-Liouville problem of dimension 2m on the compact interval $[c,b]$, with separated boundary conditions, and the following statements are true: \begin{enumerate} \item A real value $\lambda$ is an eigenvalue of Problem~1 if, and only if, $\lambda$ is an eigenvalue of Problem~2. \item The multiplicity of $\lambda$ as an eigenvalue of Problem~1 is equal to the multiplicity of $\lambda$ as an eigenvalue of Problem~2. \item Let $\lambda$ be an eigenvalue of Problem~1, and hence also of Problem~2. A function $y$ is an eigenfunction of Problem~1 if, and only if, the function $Y$ defined by \begin{equation} \label{capy} Y(x)=\left(\begin{array}{c} y(x) \\ y((a+b)-x) \end{array}\right) \end{equation} is an eigenfunction of Problem~2. \item Let $\lambda$ be an eigenvalue of multiplicity k for Problem~1, and hence also for Problem~2, and let $ y_{1}, \cdots ,y_{k} $ be linearly independent eigenfunctions of $\lambda$ for Problem~1. Let $Y_{1}, \cdots,Y_{k}$ be defined according to equation (\ref{capy}) using $y_{1},\cdots,y_{k}$ respectively. Then $\{y_{1},\cdots,y_{k} \}$ is an orthonormal set in the Hilbert space $L_{w}^{2}(a,b)$ if, and only if, the set $\{ Y_{1},\cdots ,Y_{k} \}$ is an orthonormal set in the Hilbert space $L_{W}^{2}(c,b)$. \end{enumerate} \end{thm} \paragraph{Proof:} To verify that Problem~2 is a regular, self-adjoint problem, simply verify that all of the constraints are met. We then have: \begin{enumerate} \item Suppose that $\lambda$ is an eigenvalue of Problem~1. Then there is a nontrivial solution $y(x)$ defined on the interval $[a,b]$ such that \begin{enumerate} \item $-(py')'(x)+q(x)y(x)=\lambda w(x)y(x) \,\,\mbox{for all}\,\, x\in [a,b]$ \item $C_{2}y(a)-C_{1}(py')(a)-D_{2}y(b)+D_{1}(py')(b)=0.$ \end{enumerate} Since the interval is $[a,b]$, the differential equation also holds at $(a+b)-x$; for all $x\in [a,b]$ : \[ -(py')'((a+b)-x)+q((a+b)-x)y((a+b)-x)=\lambda w((a+b)-x)y((a+b)-x). \] For $x\in [c,b]$, define \[ Y(x)=\left(\begin{array}{c} y(x) \\ y((a+b)-x) \end{array}\right). \] Clearly $Y$ is nontrivial, and $Y$ and $PY'$ are absolutely continuous, since $y$ is a nontrivial solution. The differential equation of Problem~2 is satisfied by the function~$Y$. At the left endpoint we have \begin{eqnarray*} \lefteqn{A_{1}Y(c)+A_{2}(PY')(c)} \\ & = & \left(\begin{array}{c|c} I & -I \\ \hline 0 & 0 \end{array}\right) \left(\begin{array}{c} y(c) \\ y(c) \end{array}\right) + \left(\begin{array}{c|c} 0 & 0 \\ \hline I & I \end{array}\right) \left(\begin{array}{c} (py')(c) \\ -(py')(c) \end{array}\right) \\ & = & \left(\begin{array}{c} y(c)-y(c) \\ (py')(c)-(py')(c) \end{array}\right) \,\, = \,\, 0 . \end{eqnarray*} At the right endpoint we have \begin{eqnarray*} \lefteqn{B_{1}Y(b)+B_{2}(PY')(b)} \\ & = & \left(\begin{array}{c|c} -D_{2} & C_{2} \end{array}\right) \left(\begin{array}{c} y(b) \\ y(a) \end{array}\right) + \left(\begin{array}{c|c} D_{1} & C_{1} \end{array}\right) \left(\begin{array}{c} (py')(b) \\ -(py')(a) \end{array}\right) \\ & = & -D_{2}y(b)+C_{2}y(a)+D_{1}(py')(b)-C_{1}(py')(a) \,\, = \,\,0 . \end{eqnarray*} Thus the boundary conditions are satisfied, so $Y$ is an eigenfunction for Problem~2 corresponding to the eigenvalue~$\lambda$. The proof of the converse is similar. \item The first part of the proof establishes a correspondence between eigenfunctions of Problem~1 and eigenfunctions of Problem~2. Clearly, a set of eigenfunctions $\{ y_{1},\cdots,y_{k} \}$ is linearly independent if, and only if, the corresponding set of functions $\{ Y_{1},\cdots,Y_{k} \}$ is linearly independent. Thus the number of linearly independent eigenfunctions for Problem~1 is equal to the number for Problem~2; the multiplicities for the eigenvalue $\lambda$ are equal. \item The first part of the proof establishes the relationship between the eigenfunctions. \item Consider two eigenfunctions of Problem~1, $y_{j}$ and $y_{k}$, and the corresponding functions for Problem~2, $Y_{j}$ and $Y_{k}$. The inner product of $y_{j}$ and $y_{k}$ in the Hilbert space $L_{w}^{2}(a,b)$ is \[ _{w} = \int_{a}^{b}\!y_{j}^{*}(t)w(t)y_{k}(t)\,dt. \] The inner product of $Y_{j}$ and $Y_{k}$ in the space $L_{W}^{2}(c,b)$ is $$ _{W} \ = \ \int_{c}^{b}\!Y_{j}^{*}(t)W(t)Y_{k}(t)\,dt $$ \begin{eqnarray*} & = & \int_{c}^{b}\left(\begin{array}{c} y_{j}(t) \\ y_{j}((a+b)-t) \end{array}\right)^{*} \left(\begin{array}{c|c} w(t) & 0 \\ \hline 0 & w((a+b)-t) \end{array}\right) \times \\ && \left(\begin{array}{c} y_{k}(t) \\ y_{k}((a+b)-t) \end{array}\right)\,dt \\ & = & \int_{c}^{b}\!y_{j}^{*}(t)w(t)y_{k}(t)+y_{j}^{*}((a+b)-t)w((a+b)-t)y_{ k}((a+b)-t)\,dt \\ & = & \int_{c}^{b}\!y_{j}^{*}(t)w(t)y_{k}(t)\,dt + \int_{a}^{c}\!y_{j}^{*}(t)w(t)y_{k}(t)\,dt \\ & = & \int_{a}^{b}\!y_{j}^{*}(t)w(t)y_{k}(t)\,dt \,\, = \,\, _{w}. \end{eqnarray*} Since the corresponding inner products are equal, we see that if either set of eigenfunctions is orthonormal, then the other set will be also. $\Box$ \end{enumerate} % % EXAMPLES \section*{Examples} Two examples will be presented. In the first, a comparison is made with results obtained by a technique which is quite different from the one described here. In the second, a problem will be transformed into three distinct problems, sharing the same set of eigenvalues. All three problems produce comparable results, providing evidence that the technique is working as expected. The method used to solve the separated problems produces an interval estimate for the desired eigenvalue. % % PLUM EXAMPLE % In a recent paper \cite{plum90} Plum, using an algorithm based on a homotopy, computes eigenvalues of a scalar periodic problem, Problem~a. This method is in no way similar to the algorithm developed here which transforms Problem~a into a separated problem, Problem~A. \paragraph{Problem a} On the interval $[ 0,\pi]$, \[ -y'' + 100(\cos x)^{2} y = \lambda y \] subject to the periodic boundary conditions \[ y(0) = y(\pi) \,\, , \,\, y'(0)=y'(\pi). \] In the notation we have been using, we have \[ p(x)=1 \,\, , \,\, q(x)=100(\cos x)^{2} \,\, , \,\, w(x)=1 \] with boundary conditions defined by four $2\by 1$ matrices \[ C_{1} = \left(\begin{array}{c} 1 \\ 0 \end{array}\right) \,\, , \,\, C_{2} = \left(\begin{array}{c} 0 \\ 1 \end{array}\right) \,\, , \,\, D_{1} = \left(\begin{array}{c} 1 \\ 0 \end{array}\right) \,\, , \,\, D_{2} = \left(\begin{array}{c} 0 \\ 1 \end{array}\right) .\] \paragraph{Problem A} On the interval $[ \frac{\pi}{2} , \pi]$, \begin{eqnarray*} P(x) & = & \left(\begin{array}{cc} 1 & 0 \\ 0 & 1 \end{array}\right) \,\, , \,\, W(x)\,\, = \,\, \left(\begin{array}{cc} 1 & 0 \\ 0 & 1 \end{array}\right) \\ Q(x) & = & \left(\begin{array}{cc} 100(\cos x)^{2} & 0 \\ 0 & 100(\cos(\pi - x))^{2} \end{array}\right) \end{eqnarray*} with the boundary conditions constants given by \[ A_{1}=\left(\begin{array}{cc} 1 & -1 \\ 0 & 0 \end{array}\right) \,\, , \,\, A_{2}=\left(\begin{array}{cc} 0 & 0 \\ 1 & 1 \end{array}\right) \,\, , \,\, B_{1}=\left(\begin{array}{cc} 0 & 0 \\ -1 & 1 \end{array}\right) \,\, , \,\, B_{2}=\left(\begin{array}{cc} 1 & 1 \\ 0 & 0 \end{array}\right). \] The results given by Plum are quite comparable to those computed by solving Problem~A, as shown in Table~1. The subscripts and superscripts at the end of the numbers in Table~1 and Table~2 below specify the interval in which the given eigenvalue is located; i.e. the superscript gives an upper bound and the subscript gives a lower bound for the eigenvalue. \begin{table}[hbt] \begin{center} \begin{tabular}{|c||r||r|} \hline n & \multicolumn{1}{c||}{Plum} & \multicolumn{1}{c|}{Problem A} \\ \hline 1 & $9.7432204_{0}^{6}$ & $9.7432_{14}^{23}$ \\ \hline 2 & $28.6851_{38}^{41}$ & $28.6851_{31}^{41}$ \\ \hline 3 & $46.47783_{4}^{7}$ & $46.4778_{14}^{23}$ \\ \hline 4 & $62.9864_{86}^{93}$ & $62.9864_{88}^{98}$ \\ \hline 5 & $77.8052_{37}^{44}$ & $77.8050_{61}^{71}$ \\ \hline 6 & $91.80_{07}^{11}$ & $91.801_{071}^{147}$ \\ \hline 7 & $98.975_{4}^{8}$ & $98.975_{372}^{449}$ \\ \hline \end{tabular} \caption{Comparison of Results} \end{center} \end{table} % % COUPLED PROBLEM EXAMPLE % The next example begins with a scalar problem, Problem~b, with a general, coupled boundary condition, with a non-constant, non-periodic coefficient. The second problem, Problem~c, is constructed from the first by \lq\lq folding" the interval in thirds, producing a dimension~3 problem with the same eigenvalues. The third problem, Problem~d, is constructed from the first by a change of variable: $x \mapsto \sqrt{1+x}$. Each problem is solved using the technique described above, involving two problems of dimension~2, and a problem of dimension~6. % % PROBLEM b % \paragraph{Problem b} On the interval $[0,3]$, with coefficients \[ p(x)=2 \,\, , \,\, q(x)=1+x^{2} \,\, , \,\, w(x)=1 \] subject to the boundary conditions \[ \left(\begin{array}{cc} 1 & 2 \\ 0 & 3 \end{array}\right) \left(\begin{array}{c} y \\ py' \end{array}\right) (0) + \left(\begin{array}{cc} 1.5 & 0 \\ 4 & 2 \end{array}\right) \left(\begin{array}{c} y \\ py' \end{array}\right) (3) = 0 . \] % % PROBLEM c % \paragraph{Problem c} On the interval $[0,1]$, with coefficients \begin{eqnarray*} p(x) & = & \left(\begin{array}{ccc} 2 & 0 & 0 \\ 0 & 2 & 0 \\ 0 & 0 & 2 \end{array}\right) \,\, , \,\, w(x) \,\, = \,\, \left(\begin{array}{ccc} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{array}\right) \\ q(x) & = & \left(\begin{array}{ccc} 1+x^{2} & 0 & 0 \\ 0 & 1+(x-2)^{2} & 0 \\ 0 & 0 & 1+(x+2)^{2} \end{array}\right) \end{eqnarray*} subject to the boundary conditions \begin{eqnarray*} \lefteqn{\left(\begin{array}{cccccc} 1 & 0 & 0 & 2 & 0 & 0 \\ 0 & 0 & 0 & 3 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & -1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 1 \end{array}\right) \left(\begin{array}{c} y \\ py' \end{array}\right) (0)}\\ &+& \left(\begin{array}{cccccc} 0 & 0 & 1.5 & 0 & 0 & 0 \\ 0 & 0 & 4 & 0 & 0 & 2 \\ 1 & -1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \end{array}\right) \left(\begin{array}{c} y \\ py' \end{array}\right) (1) = 0 \,.\end{eqnarray*} % % PROBLEM d % \paragraph{Problem d} On the interval $[1,2]$, with coefficients \[ p(x)=\frac{1}{x} \,\, , \,\, q(x)=2x(x^{4}-2x^{2}+2) \,\, , \,\, w(x)=2x \] subject to the boundary conditions \[ \left(\begin{array}{cc} 1 & 2 \\ 0 & 3 \end{array}\right) \left(\begin{array}{c} y \\ py' \end{array}\right) (1) + \left(\begin{array}{cc} 1.5 & 0 \\ 4 & 2 \end{array}\right) \left(\begin{array}{c} y \\ py' \end{array}\right) (2) = 0. \] The corresponding separated problems, Problem~B, Problem~C, and Problem~D were solved; the bisection process was continued in each case until the width of the interval estimate was less than 0.00001. The bisection process began with intervals with integer endpoints in all cases. The results are shown in Table~2. \begin{table}[hbt] \begin{center} \begin{tabular}{|c||r||r||r|} \hline n & \multicolumn{1}{c||}{Problem B} & \multicolumn{1}{c||}{Problem C} & \multicolumn{1}{c|}{Problem D} \\ \hline 1 & $1.6359_{71}^{79}$ & $1.6359_{71}^{79}$ & $1.6359_{71}^{79}$ \\ \hline 2 & $7.9569_{63}^{70}$ & $7.9569_{63}^{70}$ & $7.9569_{63}^{70}$ \\ \hline 3 & $12.8206_{33}^{41}$ & $12.8206_{33}^{41}$ & $12.8206_{33}^{41}$ \\ \hline 4 & $25.2577_{59}^{67}$ & $25.2577_{59}^{67}$ & $25.2577_{59}^{67}$ \\ \hline 5 & $38.4839_{17}^{25}$ & $38.4839_{25}^{32}$ & $38.4839_{17}^{25}$ \\ \hline 6 & $60.1740_{57}^{65}$ & $60.1740_{49}^{57}$ & $60.1740_{19}^{26}$ \\ \hline 7 & $82.2650_{07}^{15}$ & $82.2653_{73}^{81}$ & $82.2652_{44}^{51}$ \\ \hline 8 & $112.7659_{38}^{45}$ & $112.7656_{86}^{94}$ & $112.7644_{65}^{73}$ \\ \hline \end{tabular} \caption{Results} \end{center} \end{table} All of the separated problems were solved using a multi-step numerical integration method. To control run times, the integrator uses a relatively coarse mesh, which we believe is the source of the disagreements between some of these results. % \begin{thebibliography}{9} \bibitem{aklz87} 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} P.B. Bailey, {\em SLEIGN: An Eigenfunction-eigenvalue Code for Sturm-Liouville Problems}, {SAND7-2044}, Sandia Laboratories, Albuquerque, 1978. \bibitem{bags78} P.B. Bailey, M.K. Gordon, L.F. Shampine, {\em Automatic Solution of the Sturm-Liouville Problem}, ACM Trans. Math. Software, {\bf 4}(1978), 193-208. \bibitem{dwye93} H.I. Dwyer, {\em Eigenvalues of Matrix Sturm-Liouville Problems with Separated or Coupled Boundary Conditions}, Doctoral Thesis, Northern Illinois University, 1993. \bibitem{moze94} M. M\"{o}ller, A. Zettl, {\em Semi-boundedness of Ordinary Differential Operators}, Journal of Differential Equations, to appear. \bibitem{bgkz91} P.B. Bailey, B.S. Garbow, H.G. Kaper and A. Zettl, {\em Eigenvalue and eigenfunction computations for Sturm-Liouville problems }, ACM TOMS {\bf 17}(1991), 491-499. \bibitem{plum90} M. Plum, {\em Eigenvalue Inclusions for Second-order Ordinary Differential Operators by a Numerical Homotopy Method}, Journal of Applied Mathematics and Physics (ZAMP), {\bf 41}(March 1990), 205-226. \bibitem{fupr94} S. Pruess and C. Fulton, {\em Mathematical software for Sturm-Liouville problems}, ACM TOMS (1994), to appear. \bibitem{naglib} J. Pryce, {\em D02KEF}, NAG Library Reference Guide. \end{thebibliography} \medskip\noindent {\sc H.I. Dwyer \newline Mathematics Dept., University of Wisconsin, Platteville, WI 53818} \newline E-mail address: dwyer@uwplatt.edu \medskip\noindent {\sc A. Zettl\newline Mathematical Sciences Dept., Northern Illinois University, DeKalb, IL 60115} \newline E-mail address: zettl@math.niu.edu \vspace{3cm} \section*{Addendum} {\bf January 10, 1996}. A remark and a FORTRAN code for the main algorithm of this article have been placed as an appendix. See files \lq\lq appendix.tex" and \lq\lq Dwyer.fort" in the directory EJDE/Volumes/1994/06-Dwyer-Zettl. See also the article \newline H.I. Dwyer, A. Zettl, Eigenvalue Computations for Regular Matrix Sturm--Liouville Problems, Eletrc. J. Diff. Eqns. Vol. 1995(1995), No. 5, pp. 1-13. \end{document}