\documentstyle{article}
\setlength{\headheight}{0.0in}
\setlength{\textheight}{8.75in}
\setlength{\textwidth}{6.0in}
\setlength{\oddsidemargin}{0.15in}
\newcommand{\hs}{\hspace{1.0in}}
\newcommand{\hhs}{\hspace{0.5in}}
\newcommand{\hqs}{\hspace{0.25in}}
\newcommand{\hes}{\hspace{0.125in}}
\newcommand{\bt}{\begin{tabbing}}
\newcommand{\et}{\end{tabbing}}
\newcommand{\be}{\begin{equation}}
\newcommand{\ee}{\end{equation}}
\newcommand{\ds}{\displaystyle}
\newcommand{\ver}{\verbatim}
\newtheorem{lemma}{Lemma}
\begin{document}
\title{SL11F -- MARCO FORTRAN Library Routine Document}
\author{}
\date{}
\maketitle
{\scriptsize NOTE: before using this routine, please read the
appropriate implementation document to check the interpretation of {\bf bold
italicised } terms and other implementation-dependent details. The routine
name may be precision dependent.}
\section{ Purpose}
SL11F finds a specified eigenvalue of a differential equation eigenvalue
problem posed in the form of an eigen-parameter dependent Hamiltonian
system with suitable separated boundary conditions. The method used is
a shooting method based on matrix oscillation theory, using the algorithm
described in \cite{kn:Marletta}.
\section{ Specification}
\bt
\hs {\tt SUBROUTINE SL11F(ELAM,EPS,A,B,K,N,TOL,USUB1,USUB2,KNOBS,} \\
\hhs\hqs\hes {\tt \&}\hs\hhs {\tt WORK,IWORK,LWK,ILWK,IFAIL) } \\
\hs {\tt INTEGER ILWK, LWK(ILWK),IWORK,IFAIL, K, KNOBS, N }\\
\hs {\tt {\bf real } ELAM, EPS, A, B, TOL, WORK(IWORK)} \\
\hs {\tt EXTERNAL USUB1,USUB2 }
\et
\section{ Description \label{sec:Descrip}}
SL11F finds a specified eigenvalue $\lambda_{k} $ of a Hamiltonian eigenvalue problem of
the form
\be
\left( \begin{array}{c} u' \\ v' \end{array} \right) =
\left(\begin{array}{ll} S_{1,1}(x,\lambda) & S_{1,2}(x,\lambda) \\
S_{1,2}^{T}(x,\lambda) & S_{2,2}(x,\lambda) \end{array} \right)
\left( \begin{array}{cc} u \\ v \end{array} \right), \;\;\; x\in(a,b) \label{eq:1} \ee
\be A_{1}^{T}u(a) + A_{2}^{T}v(a) = 0, \label{eq:2} \ee
\be B_{1}^{T}u(b) + B_{2}^{T}v(b) = 0, \label{eq:3} \ee
where the $S_{i,j}$ are real $n\times n$ matrix functions of $x$ and $\lambda$, with
$S_{1,1}$ and $S_{2,2}$ symmetric, and $A_{1}$, $A_{2}$, $B_{1}$ and $B_{2}$ are
$n\times n$ matrices satisfying the conjointness conditions
\be A_{1}^{T}A_{2}-A_{2}^{T}A_{1} = 0, \;\;\; B_{1}^{T}B_{2}-B_{2}^{T}B_{1} = 0, \ee
and the rank conditions
\be \mbox{rank}(A_{1}|A_{2}) = n, \;\;\; \mbox{rank}(B_{1}|B_{2}) = n. \label{eq:rank} \ee
Here $(A_{1}|A_{2})$ denotes the $n\times 2n$ matrix whose first $n$ columns are the
columns of $A_{1}$ and whose $n+1$st to $2n$th columns are the columns of $A_{2}$.
There is no a-priori reason for the eigenvalue problem which we have just posed to have
any solutions. For a number of well-known cases (see, e.g., Atkinson \cite{kn:Atkinson})
the existence of eigenvalues is established. These cases include Hamiltonian
systems derived from $2m$-th order linear formally self-adjoint differential
operators.
The eigenvalue indexing used by SL11F is that which ensures that for Sturm-Liouville
problems and for problems derived from $2m$-th order linear self-adjoint operators,
the eigenvalues will be an infinite sequence
\[ \lambda_{0} \leq \lambda_{1} \leq \lambda_{2} \leq \cdots \]
For other problems, possibly with nonlinear dependence of (\ref{eq:1})
on the eigenparameter $\lambda$, it cannot be guaranteed that an eigenvalue will exist
with any particular index. It is possible that the first few eigenvalues will be missing,
or that there will be only finitely many eigenvalues; it is even possible for the eigenvalues
to form a sequence which extends both to $-\infty$ and to $+\infty$; such a sequence requires
an indexing $(\lambda_{k})_{k=-\infty}^{+\infty}$.
The shooting method used by SL11F integrates a differential equation to find a function
from which the position of the eigenvalues may be deduced by a rootfinding process.
The integration of the differential equation requires a tolerance TOL, which is supplied
by the user. The accuracy of the eigenvalue approximation obtained depends on TOL,
although the error in the eigenvalue approximation may well exceed TOL. It is therefore
recommended that SL11F be called with two different values of TOL in order to assess the
accuracy of the eigenvalue approximations obtained.
\section{References}
See the bibliography at the end of this document.
\section{ Parameters}
\begin{description}
\item[ELAM] -- {\bf real} scalar (input/output).
\begin{description}
\item[On entry:] ELAM must
specify an initial approximation to the eigenvalue $ \lambda_{k} $ which is
sought. This need not be at all accurate since the routine will always find an
approximation to $\lambda_{k} $ even if this is not the eigenvalue closest to
the user's initial approximation.
\item[On exit:] ELAM will contain the approximation to $\lambda_{k}$.
\end{description}
\item[EPS] -- {\bf real} scalar (input).
\begin{description}
\item[On entry:] EPS must be set to a {\bf strictly positive} order-of-maginitude
estimate of the error in the initial estimate ELAM. Over-estimates are definitely
preferable to under-estimates.
\end{description}
\item[A] -- {\bf real} (input).
\begin{description}
\item[On entry:] A must specify the left-hand endpoint $a$ of the interval $(a,b)$
on which the problem is posed.
\end{description}
\item[B] -- {\bf real} (input).
\begin{description}
\item[On entry:] B must specify the right-hand endpoint $b$ of the interval $(a,b)$
on which the problem is posed. B $>$ A is required.
\end{description}
\item[K] -- INTEGER (input).
\begin{description}
\item[On entry:] K must specify the index $k$ of the eigenvalue sought. To check for a
sequence of eigenvalues decreasing to $-\infty$, call SL11F with a large negative value of
K (e.g. -100) and a large value of EPS (e.g. EPS=10); if failure occurs with IFAIL = 2,
then such a sequence is unlikely.
\end{description}
\item[N] -- INTEGER (input).
\begin{description}
\item[On entry:] N must specify the dimension of the matrices $S_{i,j}$ in (\ref{eq:1}).
\end{description}
\item[TOL] -- {\bf real} (input).
\begin{description}
\item[ On entry:] TOL should be strictly positive. SL11F will use TOL as an error
control parameter during the shooting process. Reducing TOL should reduce the
error in the eigenvalue approximation.
\end{description}
\item[USUB1] -- SUBROUTINE, supplied by the user. USUB1 must supply the
values of the coefficient matrices $S_{i,j}$ in (\ref{eq:1}).
The specification of USUB1 is as follows.
\bt
\hs {\tt SUBROUTINE USUB1(X,ELAM,S11,S12,S21,S22,N)} \\
\hs {\bf real} {\tt X, ELAM} \\
\hs {\tt INTEGER N} \\
\hs {\bf real} {\tt S11(N,N), S12(N,N), S21(N,N), S22(N,N) }
\et
\begin{description}
\item[X] -- {\bf real} (input).
\begin{description}
\item[On entry:] X contains the value of a point in $(a,b)$. X must not be
changed by USUB1.
\end{description}
\item[ELAM] -- {\bf real} (input).
\begin{description}
\item[On entry:] ELAM contains the current value of $\lambda$. ELAM must
not be changed by USUB1.
\end{description}
\item[S11] -- {\bf real} ARRAY of DIMENSION (N,N) (output).
\begin{description}
\item[On exit:] S11 must contain the value of $S_{1,1}$ at the point X for the
current value of $\lambda$ specified by ELAM.
\end{description}
\item[S12] -- {\bf real} ARRAY of DIMENSION (N,N) (output).
\begin{description}
\item[On exit:] S12 must contain the value of $S_{1,2}$ at the point X for the
current value of $\lambda$ specified by ELAM.
\end{description}
\item[S21] -- {\bf real} ARRAY of DIMENSION (N,N) (output).
\begin{description}
\item[On exit:] S21 must contain the value of $S_{1,2}^{T}$ at the point X for the
current value of $\lambda$ specified by ELAM.
\end{description}
\item[S22] -- {\bf real} ARRAY of DIMENSION (N,N) (output).
\begin{description}
\item[On exit:] S22 must contain the value of $S_{2,2}$ at the point X for the
current value of $\lambda$ specified by ELAM.
\end{description}
\item[N] -- INTEGER (input).
\begin{description}
\item[On entry:] N specifies the value of the variable N in the user's call to
SL11F. N must not be changed by USUB1.
\end{description}
\end{description}
\item[USUB2] -- SUBROUTINE, supplied by the user. This subroutine supplies the
boundary conditions for the problem as specified by equations
(\ref{eq:2},\ref{eq:3}). The specification of USUB2 is as follows.
\bt
\hs {\tt SUBROUTINE USUB2(IEND,ISING,X,U,V,NDIM,N,ELAM) } \\
\hs {\tt INTEGER IEND, NDIM, N } \\
\hs {\tt LOGICAL ISING } \\
\hs {\bf real }{\tt X, U(NDIM,N), V(NDIM,N), ELAM }
\et
\begin{description}
\item[IEND] -- INTEGER (input).
\begin{description}
\item[On entry:] IEND specifies the end of the interval
$ (a,b) $ at which the boundary conditions are to be imposed. If IEND
= 0 then boundary conditions are required at the end $ x = a$; if IEND = 1
then boundary conditions are required at $ x = b$. IEND must not be
changed by the routine UUSB1.
\end{description}
\item[ISING] -- LOGICAL (output).
\begin{description}
\item[On exit:] ISING must have the value .FALSE. This parameter is
included to allow for a future upgrade of SL11F to handle singular
problems.
\end{description}
\item[X] -- {\bf real} scalar (input).
\begin{description}
\item[On entry:] X specifies the actual endpoint at which the
boundary condition is required (i.e. the value A or B). Like
ISING, this parameter is included to allow for a future upgrade
of SL11F to handle singular problems.
\end{description}
\item[U] -- {\bf real} ARRAY of DIMENSION (NDIM,N) (output).
\begin{description}
\item[On exit:] U must specify the matrix $A_{2}$ if IEND=0
or $B_{2}$ if IEND=1.
\end{description}
\item[V] -- {\bf real} ARRAY of DIMENSION (NDIM,N) (output).
\begin{description}
\item[On exit:] V must specify the matrix $-A_{1}$ if IEND=0
or $-B_{1}$ if IEND=1.
\end{description}
\item[NDIM] -- INTEGER (input).
\begin{description}
\item[On entry:] NDIM specifies the first dimensions of the
arrays U and V. The value of NDIM must not be changed by USUB2.
\end{description}
\item[N] -- INTEGER (input).
\begin{description}
\item[On entry:] N specifies the value of the variable $N$ in the user's call to
SL11F. N must not be changed by USUB2.
\end{description}
\item[ELAM] -- {\bf real} scalar (input).
\begin{description}
\item[On entry:] ELAM specifies the current value of $\lambda$, allowing
for $\lambda$-dependent boundary conditions. The value of ELAM must not
be changed by USUB2.
\end{description}
\end{description}
\item[KNOBS] -- INTEGER (input).
\begin{description}
\item[On entry:] KNOBS allows the user to increase the maximum number of
miss-distance function evaluations used to find the eigenvalue. The
maximum number of miss-distance function evaluations is max(60,KNOBS).
In general, it is recommended that KNOBS be assigned the value 0.
\end{description}
\item[WORK] -- {\bf real} array of DIMENSION IWORK (workspace).
\item[IWORK] -- INTEGER (input).
\begin{description}
\item[On entry:] IWORK must specify the dimension of WORK as declared in the
calling (sub)program; IWORK $\geq $ N*(16+73*N)+65.
\end{description}
\item[LWK] -- INTEGER array of DIMENSION ILWK (workspace).
\item[ILWK] -- INTEGER (input).
\begin{description}
\item[ON entry:] ILWK must specify the dimension of LWK as declared in the
calling (sub)program; ILWK $\geq $ N.
\end{description}
\item[IFAIL] -- INTEGER (input/output).
\begin{description}
\item[On entry:] IFAIL must be set to -1, 0 or 1. For users not familiar with
this parameter (described in Chapter P01 of the NAG Library documentation) the
recommended value is 0.
\item[On exit:] IFAIL = 0 unless the routine detects an error (see Section \ref{section:6}).
\end{description}
\end{description}
\section{ Error Indicators and Warnings \label{section:6}}
Errors detected by the routine:
\begin{description}
\item[IFAIL = 1.]
A parameter error on input. This error may be trigerred by all of the following:
N $<$ 2, TOL $\leq$ 0, EPS $\leq$ 0 and B $\leq$ A, IWORK or LWK less than
the minimum values specified above.
If IFAIL was -1 or 0 on entry then an error message will give more details.
\item[IFAIL = 2.]
Too many attempts have been made to find the required eigenvalue. Try increasing
the parameter KNOBS (to a value greater than 60) or increase TOL.
\item[IFAIL = 3.]
An error has occurred in F02BJF when converting the boundary conditions supplied
by the user into the variables used by SL11F. This may occur if the matrices
$(A_{1}|A_{2})$ or $(B_{1}|B_{2})$ are almost rank deficient, or if they have
some extremely large entries. Check the routine USUB2 for coding errors.
\item[IFAIL = 4.]
The boundary conditions most recently supplied do not appear to satisfy the
rank condition (\ref{eq:rank}). Check the routine USUB2 for coding errors.
\item[IFAIL = 5.]
An $H$-matrix (see \cite{kn:Marletta}) has occurred whose exponential cannot
be accurately computed. Check the routine USUB1 for coding errors and consider
reducing the value of TOL.
\item[IFAIL = 6.]
A failure has occurred on a call to F02AJF; the eigenvalues of the central
miss-distance unitary matrix $\Theta_{R}^{*}\Theta_{L}$ (see \cite{kn:Marletta})
cannot be computed accurately. Check the routine USUB1 for coding errors and
consider reducing the value of TOL.
\item[IFAIL = 7.]
An $H$-matrix (see \cite{kn:Marletta}) has occurred which cannot be
accurately diagonalised. Check the routine USUB1 for coding errors and consider
reducing the value of TOL.
\item[IFAIL = 8.]
Too many unsuccessful attempts have been made to find a suitable initial
stepsize to integrate the shooting equations. Check the routine USUB1 for
coding errors and consider changing the value of TOL (increasing TOL is
more likely to yield success).
\item[IFAIL = 9.]
The integration of the shooting equations has ground to a halt because
a suitable stepsize cannot be found for the next step. Check the routine
USUB1 for coding errors and consider changing the value of TOL (increasing
TOL is more likely to yield success).
\end{description}
\section{Accuracy}
The error in the eigenvalue approximation obtained is found empirically to be
a modest multiple of TOL for most problems. However the user is always advised
to run SL11F at two different tolerances to obtain a more reliable indication
of the accuracy which is being achieved.
\section{Run times}
The run times are proportional to $N^{3}$ and increase with decreasing TOL.
\section{ Keywords}
Eigenvalue Problems, Hamiltonian Equations, Sturm-Liouville Problems,
Matrix Oscillation Theory, Shooting Methods.
\section{ Example}
We consider the problem
\[ -y'' + Q(x)y = \lambda y, \;\;\; x\in (0.1,1], \;\;\; y(0.1) = y(1) = 0, \]
where $Q(x)$ is the following $4\times 4$ matrix:
\be Q_{i,j}(x) = \frac{1}{\max(i,j)}\cos x + \frac{\delta_{i,j}}{x^{i}}. \label{eq:exeq} \ee
This may be written in the form of a Hamiltonian system with $S_{1,1}=\lambda I - Q$,
$S_{1,2}=S_{2,1}=0$, and $S_{2,2}=I$, where $I$ denotes the identity matrix.
The symbol $\delta_{i,j}$ in (\ref{eq:exeq}) is the usual Kronecker $\delta$:
\[ \delta_{i,j} = \left\{ \begin{array}{ll} 0 & \mbox{for $i\neq j$} \\
1 & \mbox{for $i=j$}
\end{array}{ll} \right. \]
We ask SL11F to compute an approximation to $\lambda_{3}$.
\subsection{ Program Text}
\noindent {\scriptsize WARNING: This {\bf double precision} program may
require amendment for certain implementations. The results produced may not be
the same. If in doubt, please seek further advice. }
\begin{verbatim}
C SL11F EXAMPLE PROGRAM TEXT.
C MARK n RELEASE. MARCO MARLETTA COPYRIGHT 1992.
C .. PARAMETERS ..
INTEGER IWORK,NMAX
PARAMETER (NMAX=10,IWORK=NMAX*(16+73*NMAX)+65)
C .. LOCAL SCALARS ..
INTEGER IFAIL,K,KNOBS,N,NOUT
DOUBLE PRECISION ELAM,EPS,A,B,TOL
C .. LOCAL ARRAYS ..
INTEGER LWK(NMAX)
DOUBLE PRECISION WORK(IWORK)
C .. EXTERNAL SUBROUTINES ..
EXTERNAL USUB1,USUB2
C ..
DATA N /4/
DATA NOUT /6/
ELAM = 0.0D0
EPS = 1.0D0
TOL = 1.D-7
A = 0.1D0
B = 1.0D0
K = 3
IFAIL = 1
CALL SL11F(ELAM,EPS,A,B,K,N,TOL,USUB1,USUB2,KNOBS,WORK,
& IWORK,LWK,NMAX,IFAIL)
IF (IFAIL.NE.0) GOTO 100
WRITE(NOUT,9000) K,ELAM
9000 FORMAT(' LAMBDA(',I2,')= ',G18.8)
STOP
100 WRITE(NOUT,9010) IFAIL
9010 FORMAT(' Abnormal exit from SL11F, IFAIL = ',I2)
WRITE(NOUT,9020) ELAM
9020 FORMAT(' ELAM = ',G18.8)
STOP
END
SUBROUTINE USUB1(X,ELAM,S11,S12,S21,S22,N)
INTEGER I,J,N
DOUBLE PRECISION X,ELAM,S11(N,N),S12(N,N),S21(N,N),S22(N,N)
DO 10 J=1,N
DO 20 I=1,N
S12(I,J) = 0.0D0
S21(I,J) = 0.0D0
S22(I,J) = 0.0D0
S11(I,J) = -COS(X)/DBLE(MAX(I,J))
20 CONTINUE
S22(J,J) = 1.0D0
S11(J,J) = S11(J,J) + ELAM - X**(-J)
10 CONTINUE
RETURN
END
SUBROUTINE USUB2(IEND,ISING,X,U,V,NDIM,N,ELAM)
INTEGER I,J,IEND,NDIM,N
LOGICAL ISING
DOUBLE PRECISION ELAM,X,U(NDIM,N),V(NDIM,N)
ISING = .FALSE.
DO 10 J=1,N
DO 20 I=1,N
U(I,J) = 0.0D0
V(I,J) = 0.0D0
20 CONTINUE
V(J,J) = 1.0D0
10 CONTINUE
C REMARK: WE COULD ACTUALLY ASSIGN V TO BE ANY NON-SINGULAR MATRIX.
RETURN
END
\end{verbatim}
\subsection{Results}
\begin{verbatim}
LAMBDA( 3) = 26.920694
\end{verbatim}
\begin{thebibliography}{99}
\bibitem{kn:Atkinson} F.V. Atkinson, {\em Discrete and Continuous Boundary
Value Problems}, Academic Press, New York (1964).
\bibitem{kn:Marletta} M. Marletta, {\em Numerical Solution of Eigenvalue
Problems for Hamiltonian Systems}, University of Leicester Mathematics
Department Technical Report 1992/17.
\end{thebibliography}
\end{document}