\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}