NAG Library Routine Document
D02KAF
1 Purpose
D02KAF finds a specified eigenvalue of a regular secondorder Sturm–Liouville system defined on a finite range, using a Pruefer transformation and a shooting method.
2 Specification
INTEGER 
K, IFAIL 
REAL (KIND=nag_wp) 
XL, XR, BCOND(3,2), TOL, ELAM, DELAM 
EXTERNAL 
COEFFN, MONIT 

3 Description
D02KAF finds a specified eigenvalue
$\stackrel{~}{\lambda}$ of a Sturm–Liouville system defined by a selfadjoint differential equation of the secondorder
together with boundary conditions of the form
at the two, finite, end points
$a$ and
$b$. The functions
$p$ and
$q$, which are realvalued, are defined by
COEFFN.
For the theoretical basis of the numerical method to be valid, the following conditions should hold on the coefficient functions:
(a) 
$p\left(x\right)$ must be nonzero and must not change sign throughout the closed interval $\left[a,b\right]$; 
(b) 
$\frac{\partial q}{\partial \lambda}$ must not change sign and must be nonzero throughout the open interval $\left(a,b\right)$ and for all relevant values of $\lambda $, and must not be identically zero as $x$ varies, for any relevant value $\lambda $; and, 
(c) 
$p$ and $q$ should (as functions of $x$) have continuous derivatives, preferably up to the fourthorder, on $\left[a,b\right]$. The differential equation code used will integrate through mild discontinuities, but probably with severely reduced efficiency. Therefore, if $p$ and $q$ violate this condition, D02KDF should be used. 
The eigenvalue
$\stackrel{~}{\lambda}$ is determined by a shooting method based on a Pruefer transformation of the differential equations. Providing certain assumptions are met, the computed value of
$\stackrel{~}{\lambda}$ will be correct to within a mixed absolute/relative error specified by
TOL. D02KAF is a driver routine for the more complicated routine
D02KDF whose specification provides more details of the techniques used.
A good account of the theory of Sturm–Liouville systems, with some description of Pruefer transformations, is given in Chapter X of
Birkhoff and Rota (1962). An introduction to the use of Pruefer transformations for the numerical solution of eigenvalue problems arising from physics and chemistry is given in
Bailey (1966).
4 References
Bailey P B (1966) Sturm–Liouville eigenvalues via a phase function SIAM J. Appl. Math. 14 242–249
Birkhoff G and Rota G C (1962) Ordinary Differential Equations Ginn & Co., Boston and New York
5 Parameters
 1: XL – REAL (KIND=nag_wp)Input
 2: XR – REAL (KIND=nag_wp)Input
On entry: the left and righthand end points $a$ and $b$ respectively, of the interval of definition of the problem.
Constraint:
${\mathbf{XL}}<{\mathbf{XR}}$.
 3: COEFFN – SUBROUTINE, supplied by the user.External Procedure
COEFFN must compute the values of the coefficient functions
$p\left(x\right)$ and
$q\left(x;\lambda \right)$ for given values of
$x$ and
$\lambda $.
Section 3 states the conditions which
$p$ and
$q$ must satisfy.
The specification of
COEFFN is:
INTEGER 
JINT 
REAL (KIND=nag_wp) 
P, Q, DQDL, X, ELAM 

 1: P – REAL (KIND=nag_wp)Output
On exit: the value of $p\left(x\right)$ for the current value of $x$.
 2: Q – REAL (KIND=nag_wp)Output
On exit: the value of $q\left(x;\lambda \right)$ for the current value of $x$ and the current trial value of $\lambda $.
 3: DQDL – REAL (KIND=nag_wp)Output
On exit: the value of
$\frac{\partial q}{\partial \lambda}\left(x;\lambda \right)$ for the current value of
$x$ and the current trial value of
$\lambda $. However
DQDL is only used in error estimation and, in the rare cases where it may be difficult to evaluate, an approximation (say to within
$20\%$) will suffice.
 4: X – REAL (KIND=nag_wp)Input
On entry: the current value of $x$.
 5: ELAM – REAL (KIND=nag_wp)Input
On entry: the current trial value of the eigenvalue parameter $\lambda $.
 6: JINT – INTEGERInput
This parameter is included for compatibility with the more complex routine
D02KDF (which is called by D02KAF).
On entry: need not be set.
COEFFN must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which D02KAF is called. Parameters denoted as
Input must
not be changed by this procedure.
 4: BCOND($3$,$2$) – REAL (KIND=nag_wp) arrayInput/Output
On entry:
${\mathbf{BCOND}}\left(1,1\right)$ and
${\mathbf{BCOND}}\left(2,1\right)$ must contain the numbers
${a}_{1}$,
${a}_{2}$ specifying the lefthand boundary condition in the form
where
$\left{a}_{2}\right+\left{a}_{1}p\left(a\right)\right\ne 0$.
${\mathbf{BCOND}}\left(1,2\right)$ and
${\mathbf{BCOND}}\left(2,2\right)$ must contain
${b}_{1}$,
${b}_{2}$ such that
where
$\left{b}_{2}\right+\left{b}_{1}p\left(b\right)\right\ne 0$.
Note the occurrence of $p\left(a\right)$, $p\left(b\right)$ in these formulae.
On exit:
${\mathbf{BCOND}}\left(3,1\right)$ and
${\mathbf{BCOND}}\left(3,2\right)$ hold values
${\sigma}_{l},{\sigma}_{r}$ estimating the sensitivity of the computed eigenvalue to changes in the boundary conditions. These values should only be of interest if the boundary conditions are, in some sense, an approximation to some ‘true’ boundary conditions. For example, if the range [
XL,
XR] should really be
$\left[0,\infty \right]$ but instead
XR has been given a large value and the boundary conditions at infinity applied at
XR, then the sensitivity parameter
${\sigma}_{r}$ may be of interest. Refer to
Section 8.5 in D02KDF, for the actual meaning of
${\sigma}_{r}$ and
${\sigma}_{l}$.
 5: K – INTEGERInput
On entry:
$k$, the index of the required eigenvalue when the eigenvalues are ordered
Constraint:
${\mathbf{K}}\ge 0$.
 6: TOL – REAL (KIND=nag_wp)Input
On entry: the tolerance parameter which determines the accuracy of the computed eigenvalue. The error estimate held in
DELAM on exit satisfies the mixed absolute/relative error test
where
ELAM is the final estimate of the eigenvalue.
DELAM is usually somewhat smaller than the righthand side of
(1) but not several orders of magnitude smaller.
Constraint:
${\mathbf{TOL}}>0.0$.
 7: ELAM – REAL (KIND=nag_wp)Input/Output
On entry: an initial estimate of the eigenvalue $\stackrel{~}{\lambda}$.
On exit: the final computed estimate, whether or not an error occurred.
 8: DELAM – REAL (KIND=nag_wp)Input/Output
On entry: an indication of the scale of the problem in the
$\lambda $direction.
DELAM holds the initial ‘search step’ (positive or negative). Its value is not critical, but the first two trial evaluations are made at
ELAM and
${\mathbf{ELAM}}+{\mathbf{DELAM}}$, so the routine will work most efficiently if the eigenvalue lies between these values. A reasonable choice (if a closer bound is not known) is about half the distance between adjacent eigenvalues in the neighbourhood of the one sought. In practice, there will often be a problem, similar to the one in hand but with known eigenvalues, which will help one to choose initial values for
ELAM and
DELAM.
If ${\mathbf{DELAM}}=0.0$ on entry, it is given the default value of $0.25\times \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1.0,\left{\mathbf{ELAM}}\right\right)$.
On exit: if
${\mathbf{IFAIL}}={\mathbf{0}}$,
DELAM holds an estimate of the absolute error in the computed eigenvalue, that is
$\left\stackrel{~}{\lambda}{\mathbf{ELAM}}\right\simeq {\mathbf{DELAM}}$, where
$\stackrel{~}{\lambda}$ is the true eigenvalue.
If
${\mathbf{IFAIL}}\ne {\mathbf{0}}$,
DELAM may hold an estimate of the error, or its initial value, depending on the value of
IFAIL. See
Section 6 for further details.
 9: MONIT – SUBROUTINE, supplied by the NAG Library or the user.External Procedure
MONIT is called by D02KAF at the end of each iteration for
$\lambda $ and allows you to monitor the course of the computation by printing out the parameters (see
Section 9 for an example).
If no monitoring is required, the dummy (sub)program D02KAY may be used. (D02KAY is included in the NAG Library.)
The specification of
MONIT is:
INTEGER 
NIT, IFLAG 
REAL (KIND=nag_wp) 
ELAM, FINFO(15) 

 1: NIT – INTEGERInput
On entry: 15 minus the number of iterations used so far in the search for $\stackrel{~}{\lambda}$. (Up to $15$ iterations are permitted.)
 2: IFLAG – INTEGERInput
On entry: describes what phase the computation is in.
 ${\mathbf{IFLAG}}<0$
 An error occurred in the computation at this iteration; an error exit from D02KAF will follow.
 ${\mathbf{IFLAG}}=1$
 The routine is trying to bracket the eigenvalue $\stackrel{~}{\lambda}$.
 ${\mathbf{IFLAG}}=2$
 The routine is converging to the eigenvalue $\stackrel{~}{\lambda}$ (having already bracketed it).
Normally, the iteration will terminate after a sequence of iterates with ${\mathbf{IFLAG}}=2$, but occasionally the bracket on $\stackrel{~}{\lambda}$ thus determined will not be sufficiently small and the iteration will be repeated with tighter accuracy control.
 3: ELAM – REAL (KIND=nag_wp)Input
On entry: the current trial value of $\stackrel{~}{\lambda}$.
 4: FINFO($15$) – REAL (KIND=nag_wp) arrayInput
On entry: information about the behaviour of the shooting method, and diagnostic information in the case of errors. It should not normally be printed in full if no error has occurred (that is, if
${\mathbf{IFLAG}}\ge 0$), though the first few components may be of interest to you. In case of an error (
${\mathbf{IFLAG}}<0$) all the components of
FINFO should be printed.
The contents of
FINFO are as follows:
 ${\mathbf{FINFO}}\left(1\right)$
 The current value of the ‘missdistance’ or ‘residual’ function $f\left(\lambda \right)$ on which the shooting method is based. $f\left(\stackrel{~}{\lambda}\right)=0$ in theory. This is set to zero if ${\mathbf{IFLAG}}<0$.
 ${\mathbf{FINFO}}\left(2\right)$
 An estimate of the quantity $\partial \lambda $ defined as follows. Consider the perturbation in the missdistance $f\left(\lambda \right)$ that would result if the local error in the solution of the differential equation were always positive and equal to its maximum permitted value. Then $\partial \lambda $ is the perturbation in $\lambda $ that would have the same effect on $f\left(\lambda \right)$. Thus, at the zero of $f\left(\lambda \right),\left\partial \lambda \right$ is an approximate bound on the perturbation of the zero (that is the eigenvalue) caused by errors in numerical solution. If $\partial \lambda $ is very large then it is possible that there has been a programming error in COEFFN such that $q$ is independent of $\lambda $. If this is the case, an error exit with ${\mathbf{IFAIL}}={\mathbf{5}}$ should follow. ${\mathbf{FINFO}}\left(2\right)$ is set to zero if ${\mathbf{IFLAG}}<0$.
 ${\mathbf{FINFO}}\left(3\right)$
 The number of internal iterations, using the same value of $\lambda $ and tighter accuracy tolerances, needed to bring the accuracy (that is, the value of $\partial \lambda $) to an acceptable value. Its value should normally be $1.0$, and should almost never exceed $2.0$.
 ${\mathbf{FINFO}}\left(4\right)$
 The number of calls to COEFFN at this iteration.
 ${\mathbf{FINFO}}\left(5\right)$
 The number of successful steps taken by the internal differential equation solver at this iteration. A step is successful if it is used to advance the integration.
 ${\mathbf{FINFO}}\left(6\right)$
 The number of unsuccessful steps used by the internal integrator at this iteration.
 ${\mathbf{FINFO}}\left(7\right)$
 The number of successful steps at the maximum step size taken by the internal integrator at this iteration.
 ${\mathbf{FINFO}}\left(8\right)$
 Not used.
 ${\mathbf{FINFO}}\left(9\right)$ to ${\mathbf{FINFO}}\left(15\right)$
 Set to zero, unless ${\mathbf{IFLAG}}<0$ in which case they hold the following values describing the point of failure:
 ${\mathbf{FINFO}}\left(9\right)$
 1 or $2$ depending on whether integration was in a forward or backward direction at the time of failure.
 ${\mathbf{FINFO}}\left(10\right)$
 The value of the independent variable, $x$, the point at which the error occurred.
 ${\mathbf{FINFO}}\left(11\right)$, ${\mathbf{FINFO}}\left(12\right)$, ${\mathbf{FINFO}}\left(13\right)$
 The current values of the Pruefer dependent variables $\beta $, $\varphi $ and $\rho $ respectively. See Section 3 in D02KEF for a description of these variables.
 ${\mathbf{FINFO}}\left(14\right)$
 The localerror tolerance being used by the internal integrator at the point of failure.
 ${\mathbf{FINFO}}\left(15\right)$
 The last integration mesh point.
MONIT must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which D02KAF is called. Parameters denoted as
Input must
not be changed by this procedure.
 10: IFAIL – INTEGERInput/Output

On entry:
IFAIL must be set to
$0$,
$1\text{ or}1$. If you are unfamiliar with this parameter you should refer to
Section 3.3 in the Essential Introduction for details.
For environments where it might be inappropriate to halt program execution when an error is detected, the value
$1\text{ or}1$ is recommended. If the output of error messages is undesirable, then the value
$1$ is recommended. Otherwise, if you are not familiar with this parameter, the recommended value is
$0$.
When the value $\mathbf{1}\text{ or}\mathbf{1}$ is used it is essential to test the value of IFAIL on exit.
On exit:
${\mathbf{IFAIL}}={\mathbf{0}}$ unless the routine detects an error or a warning has been flagged (see
Section 6).
6 Error Indicators and Warnings
If on entry
${\mathbf{IFAIL}}={\mathbf{0}}$ or
${{\mathbf{1}}}$, explanatory error messages are output on the current error message unit (as defined by
X04AAF).
Errors or warnings detected by the routine:
 ${\mathbf{IFAIL}}=1$
On entry,  ${\mathbf{K}}<0$, 
or  ${\mathbf{TOL}}\le 0.0$. 
 ${\mathbf{IFAIL}}=2$
On entry,  ${a}_{1}=p\left(a\right){a}_{2}=0$, 
or  ${b}_{1}=p\left(b\right){b}_{2}=0$, 
(the array
BCOND has been set up incorrectly).
 ${\mathbf{IFAIL}}=3$

At some point between
XL and
XR the value of
$p\left(x\right)$ computed by
COEFFN became zero or changed sign. See the last call of
MONIT for details.
 ${\mathbf{IFAIL}}=4$
After $15$ iterations the eigenvalue had not been found to the required accuracy.
 ${\mathbf{IFAIL}}=5$
The ‘bracketing’ phase (with
IFLAG of the
MONIT equal to
$1$) failed to bracket the eigenvalue within ten iterations. This is caused by an error in formulating the problem (for example,
$q$ is independent of
$\lambda $), or by very poor initial estimates of
ELAM and
DELAM.
On exit,
ELAM and
${\mathbf{ELAM}}+{\mathbf{DELAM}}$ give the end points of the interval within which no eigenvalue was located by the routine.
 ${\mathbf{IFAIL}}=6$

To obtain the desired accuracy the local error tolerance was set so small at the start of some subinterval that the differential equation solver could not choose an initial step size large enough to make significant progress. See the last call of
MONIT for diagnostics.
 ${\mathbf{IFAIL}}=7$
At some point the step size in the differential equation solver was reduced to a value too small to make significant progress (for the same reasons as with
${\mathbf{IFAIL}}={\mathbf{6}}$). This could be due to pathological behaviour of
$p\left(x\right)$ and
$q\left(x;\lambda \right)$ or to an unreasonable accuracy requirement or to the current value of
$\lambda $ making the equations ‘stiff’. See the last call of
MONIT for details.
 ${\mathbf{IFAIL}}=8$
TOL is too small for the problem being solved and the
machine precision being used. The local value of
ELAM should be a very good approximation to the eigenvalue
$\stackrel{~}{\lambda}$.
 ${\mathbf{IFAIL}}=9$
C05AZF, called by D02KAF, has terminated with the error exit corresponding to a pole of the matching function. This error exit should not occur, but if it does, try solving the problem again with a smaller value for
TOL.
Note: error exits with
${\mathbf{IFAIL}}={\mathbf{2}}$,
${\mathbf{3}}$,
${\mathbf{6}}$,
${\mathbf{7}}$ or
${\mathbf{9}}$ are caused by the inability to set up or solve the differential equation at some iteration and will be immediately preceded by a call of
MONIT giving diagnostic information.
 ${\mathbf{IFAIL}}=10$
 ${\mathbf{IFAIL}}=11$
 ${\mathbf{IFAIL}}=12$

A serious error has occurred in an internal call to an interpolation routine. Check all (sub)program calls and array dimensions. Seek expert help.
7 Accuracy
The absolute accuracy of the computed eigenvalue is usually within a mixed absolute/relative bound defined by
TOL (as defined above).
The time taken by D02KAF depends on the complexity of the coefficient functions, whether they or their derivatives are rapidly changing, the tolerance demanded, and how many iterations are needed to obtain convergence. The amount of work per iteration is roughly doubled when
TOL is divided by
$16$. To make the most economical use of the routine, one should try to obtain good initial values for
ELAM and
DELAM.
See
Section 8 in D02KDF for a discussion of the technique used.
9 Example
This example finds the fourth eigenvalue of Mathieu's equations
with boundary conditions
and
$q=5$. We use a starting value
${\mathbf{ELAM}}=15.0$ and a step
${\mathbf{DELAM}}=4.0$. We illustrate the effect of varying
TOL by choosing
${\mathbf{TOL}}=\text{1.0E\u22125}$ and
$\text{1.0E\u22126}$ (note the change in the output value of the error estimate
DELAM). The range of integration and initial estimates are read from a data file.
9.1 Program Text
Program Text (d02kafe.f90)
9.2 Program Data
Program Data (d02kafe.d)
9.3 Program Results
Program Results (d02kafe.r)