\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{SL12F -- 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} SL12F 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}, combined with coefficient approximation. \section{ Specification} \bt \hs {\tt SUBROUTINE SL12F(ELAM,EPS,A,B,K,N,USUB1,USUB2,TOL,XMESH,MESH,} \\ \hhs\hqs\hes {\tt \&}\hs\hhs {\tt WORK,IWORK,LWK,ILWK,NXTRAP,IFAIL) } \\ \hs {\tt INTEGER ILWK, LWK(ILWK),IWORK,IFAIL, K, MESH, N, NXTRAP }\\ \hs {\tt {\bf real } ELAM, EPS, A, B, TOL, XMESH(0:MESH), WORK(IWORK)} \\ \hs {\tt EXTERNAL USUB1,USUB2 } \et \section{ Description \label{sec:Descrip}} SL12F 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}$ symmetric and $S_{2,2}$ symmetric positive definite, 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. The eigenvalue indexing used by SL12F is that which ensures that for Sturm-Liouville problems, the eigenvalues will be an infinite sequence \[ \lambda_{0} \leq \lambda_{1} \leq \lambda_{2} \leq \cdots \] With this indexing and the 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 SL12F 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 is performed by first approximating the coefficients $S_{i,j}$ in (\ref{eq:1}) by piecewise constants over a suitable mesh. This enables an eigenvalue approximation to be found. The mesh is then subdivided and Wynn extrapolation applied until the error in the eigenvalue approximation appears to be less than the user's tolerance TOL. \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 SL12F 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 = 3, 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[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. $S_{2,2}$ must be positive definite; if it is not, then SL12F will fail. \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 SL12F. 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 SL12F 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 SL12F 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 SL12F. 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[TOL] -- {\bf real} (input). \begin{description} \item[ On entry:] TOL should be strictly positive. SL12F will perform at most 10 Wynn extrapolations in an attempt to ensure that the error in the eigenvalue approximation returned is less than TOL. \end{description} \item[XMESH] -- {\bf real} array of DIMENSION (0:MESH) (workspace). \item[MESH] -- INTEGER (input/output) \begin{description} \item[On entry:] MESH must specify the first dimension of XMESH as declared in the calling (sub)program. Since XMESH is used to store the mesh used by SL12F, and since the number of mesh intervals required depends on the user's accuracy requirements, it is impossible to specify a value of MESH which will always be sufficient. However for regular problems and modest tolerances (no smaller than $10^{-6}$), MESH=2000 will usually be more than adequate. \item[On exit] MESH specifies the number of mesh intervals actually used. \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*(8+60*N) \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[NXTRAP] -- INTEGER (output). \begin{description} \item[On exit:] NXTRAP specifies the number of extrapolations which were required. This gives some measure of the inherent difficulty of the problem concerned. \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.] An error has occurred while converting the boundary conditions from the user's variables to SL12F's variables. Check the routine USUB2, making sure that the boundary conditions satisfy the rank condition and the conjointness condition. \item[IFAIL = 3.] Too many evaluations of the miss-distanve function have been made without an eigenvalue being located. Try increasing EPS. It is also possible that no eigenvalue exists with the index specified. \item[IFAIL = 4.] The meshing routine is making little progress as the stepsize has become too small. Try increasing the tolerance TOL. \item[IFAIL = 5.] The array XMESH is not large enough. Increase MESH or increase TOL. \item[IFAIL = 6.] The matrix $V+iU$ (see \cite{kn:Marletta}) has become rank deficient. Check the routines USUB1 and USUB2 for coding errors and consider reducing TOL. \item[IFAIL = 7.] The matrix $V+iU$ (see \cite{kn:Marletta}) has become ill-conditioned. Check the routines USUB1 and USUB2 for coding errors and consider reducing TOL. \item[IFAIL = 8.] The eigenvalues of a $\Theta$-matrix (see \cite{kn:Marletta}) cannot be computed accurately, possibly due to ill-conditioning of $V+iU$. Proceed as for IFAIL=7. \item[IFAIL = 9.] The eigenvalues of an H-matrix (see \cite{kn:Marletta}) cannot be computed accurately, possibly due to the presence of large elements. Proceed as for IFAIL=7. \item[IFAIL=10] An H-matrix has become non-Hermitian, possibly due to an accumulation of round-off error. Proceed as for IFAIL=7. \item[IFAIL=11] The coefficient $S_{2,2}$ cannot be diagonalised with sufficient accuracy. Check the routine USUB1 and make sure that $S_{2,2}$ is symmetric. \item[IFAIL=12] The coefficient $S_{2,2}$ is not positive definite. Check the routine USUB1. \end{description} \section{Accuracy} For most problems the eigenvalue error will be less than the user's tolerance TOL. Nevertheless for difficult problems it may be advisable to run SL12F at two different tolerances to obtain an indication of the reliability of the error estimates. \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, Coefficient Approximation. \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 SL12F 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 SL12F EXAMPLE PROGRAM TEXT. C MARK n RELEASE. MARCO MARLETTA COPYRIGHT 1992. C .. PARAMETERS .. INTEGER IWORK,NMAX,MAXMSH PARAMETER (NMAX=10,IWORK=NMAX*(8+60*NMAX),MAXMSH=5000) C .. LOCAL SCALARS .. INTEGER IFAIL,K,MESH,N,NOUT,NXTRAP DOUBLE PRECISION ELAM,EPS,A,B,TOL C .. LOCAL ARRAYS .. INTEGER LWK(NMAX) DOUBLE PRECISION WORK(IWORK),XMESH(0:MAXMSH) C .. EXTERNAL SUBROUTINES .. EXTERNAL USUB1,USUB2 C .. DATA N /4/ DATA NOUT /6/ ELAM = 0.0D0 EPS = 1.0D0 MESH = MAXMSH TOL = 1.D-7 A = 0.1D0 B = 1.0D0 K = 3 IFAIL = 1 CALL SL12F(ELAM,EPS,A,B,K,N,USUB1,USUB2,TOL,XMESH,MESH,WORK, & IWORK,LWK,NMAX,NXTRAP,IFAIL) IF (IFAIL.NE.0) GOTO 100 WRITE(NOUT,9000) K,ELAM 9000 FORMAT(' LAMBDA(',I2,')= ',G18.8) WRITE (NOUT,9030) MESH,NXTRAP 9030 FORMAT(' MESHPOINTS: ',I5,' EXTRAPOLATIONS: ',I3) STOP 100 WRITE(NOUT,9010) IFAIL 9010 FORMAT(' Abnormal exit from SL12F, 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.920732 MESHPOINTS: 171 EXTRAPOLATIONS: 2 \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}