NAG FL Interface
d02kaf (sl2_​reg_​finite)

1 Purpose

d02kaf finds a specified eigenvalue of a regular second-order Sturm–Liouville system defined on a finite range, using a Pruefer transformation and a shooting method.

2 Specification

Fortran Interface
Subroutine d02kaf ( xl, xr, coeffn, bcond, k, tol, elam, delam, monit, ifail)
Integer, Intent (In) :: k
Integer, Intent (Inout) :: ifail
Real (Kind=nag_wp), Intent (In) :: xl, xr, tol
Real (Kind=nag_wp), Intent (Inout) :: bcond(3,2), elam, delam
External :: coeffn, monit
C Header Interface
#include <nag.h>
void  d02kaf_ (const double *xl, const double *xr,
void (NAG_CALL *coeffn)(double *p, double *q, double *dqdl, const double *x, const double *elam, const Integer *jint),
double bcond[], const Integer *k, const double *tol, double *elam, double *delam,
void (NAG_CALL *monit)(const Integer *nit, const Integer *iflag, const double *elam, const double finfo[]),
Integer *ifail)
The routine may be called by the names d02kaf or nagf_ode_sl2_reg_finite.

3 Description

d02kaf finds a specified eigenvalue λ~ of a Sturm–Liouville system defined by a self-adjoint differential equation of the second-order
px y + q x;λ y=0 ,   a<x<b ,  
together with boundary conditions of the form
a2ya=a1paya  
b2yb=b1pbyb  
at the two, finite, end points a and b. The functions p and q, which are real-valued, are defined by coeffn.
For the theoretical basis of the numerical method to be valid, the following conditions should hold on the coefficient functions:
  1. (a)px must be nonzero and must not change sign throughout the closed interval a,b;
  2. (b) q λ must not change sign and must be nonzero throughout the open interval a,b and for all relevant values of λ, and must not be identically zero as x varies, for any relevant value λ; and,
  3. (c)p and q should (as functions of x) have continuous derivatives, preferably up to the fourth-order, on a,b. 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 λ~ is determined by a shooting method based on a Pruefer transformation of the differential equations. Providing certain assumptions are met, the computed value of λ~ 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 Arguments

1: xl Real (Kind=nag_wp) Input
2: xr Real (Kind=nag_wp) Input
On entry: the left- and right-hand end points a and b respectively, of the interval of definition of the problem.
Constraint: xl<xr.
3: coeffn Subroutine, supplied by the user. External Procedure
coeffn must compute the values of the coefficient functions px and qx;λ for given values of x and λ. Section 3 states the conditions which p and q must satisfy.
The specification of coeffn is:
Fortran Interface
Subroutine coeffn ( p, q, dqdl, x, elam, jint)
Integer, Intent (In) :: jint
Real (Kind=nag_wp), Intent (In) :: x, elam
Real (Kind=nag_wp), Intent (Out) :: p, q, dqdl
C Header Interface
void  coeffn_ (double *p, double *q, double *dqdl, const double *x, const double *elam, const Integer *jint)
1: p Real (Kind=nag_wp) Output
On exit: the value of px for the current value of x.
2: q Real (Kind=nag_wp) Output
On exit: the value of qx;λ for the current value of x and the current trial value of λ.
3: dqdl Real (Kind=nag_wp) Output
On exit: the value of q λ x;λ for the current value of x and the current trial value of λ. 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 argument λ.
6: jint Integer Input
This argument 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. Arguments denoted as Input must not be changed by this procedure.
Note: coeffn should not return floating-point NaN (Not a Number) or infinity values, since these are not handled by d02kaf. If your code inadvertently does return any NaNs or infinities, d02kaf is likely to produce unexpected results.
4: bcond32 Real (Kind=nag_wp) array Input/Output
On entry: bcond11 and bcond21 must contain the numbers a1, a2 specifying the left-hand boundary condition in the form
a2 ya = a1 pa ya  
where a2+a1pa0.
bcond12 and bcond22 must contain b1, b2 such that
b2 yb = b1 pb yb  
where b2+b1pb0.
Note the occurrence of pa, pb in these formulae.
On exit: bcond31 and bcond32 hold values σl,σ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 0, but instead xr has been given a large value and the boundary conditions at infinity applied at xr, the sensitivity argument σr may be of interest. Refer to Section 9.5 in d02kdf, for the actual meaning of σr and σl.
5: k Integer Input
On entry: k, the index of the required eigenvalue when the eigenvalues are ordered
λ0 < λ1 < λ2 < < λk < .  
Constraint: k0.
6: tol Real (Kind=nag_wp) Input
On entry: the tolerance argument which determines the accuracy of the computed eigenvalue. The error estimate held in delam on exit satisfies the mixed absolute/relative error test
delamtol×max1.0,elam, (1)
where elam is the final estimate of the eigenvalue. delam is usually somewhat smaller than the right-hand side of (1) but not several orders of magnitude smaller.
Constraint: tol>0.0.
7: elam Real (Kind=nag_wp) Input/Output
On entry: an initial estimate of the eigenvalue λ~.
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 λ-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 elam+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 delam=0.0 on entry, it is given the default value of 0.25×max1.0,elam.
On exit: if ifail=0, delam holds an estimate of the absolute error in the computed eigenvalue, that is λ~-elamdelam, where λ~ is the true eigenvalue.
If ifail0, 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 λ and allows you to monitor the course of the computation by printing out the arguments (see Section 10 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:
Fortran Interface
Subroutine monit ( nit, iflag, elam, finfo)
Integer, Intent (In) :: nit, iflag
Real (Kind=nag_wp), Intent (In) :: elam, finfo(15)
C Header Interface
void  monit_ (const Integer *nit, const Integer *iflag, const double *elam, const double finfo[])
1: nit Integer Input
On entry: 15 minus the number of iterations used so far in the search for λ~. (Up to 15 iterations are permitted.)
2: iflag Integer Input
On entry: describes what phase the computation is in.
iflag<0
An error occurred in the computation at this iteration; an error exit from d02kaf will follow.
iflag=1
The routine is trying to bracket the eigenvalue λ~.
iflag=2
The routine is converging to the eigenvalue λ~ (having already bracketed it).
Normally, the iteration will terminate after a sequence of iterates with iflag=2, but occasionally the bracket on λ~ 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 λ~.
4: finfo15 Real (Kind=nag_wp) array Input
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 iflag0), though the first few components may be of interest to you. In case of an error (iflag<0) all the components of finfo should be printed.
The contents of finfo are as follows:
finfo1
The current value of the ‘miss-distance’ or ‘residual’ function fλ on which the shooting method is based. fλ~=0 in theory. This is set to zero if iflag<0.
finfo2
An estimate of the quantity λ defined as follows. Consider the perturbation in the miss-distance fλ 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 λ is the perturbation in λ that would have the same effect on fλ. Thus, at the zero of fλ,λ is an approximate bound on the perturbation of the zero (that is the eigenvalue) caused by errors in numerical solution. If λ is very large then it is possible that there has been a programming error in coeffn such that q is independent of λ. If this is the case, an error exit with ifail=5 should follow. finfo2 is set to zero if iflag<0.
finfo3
The number of internal iterations, using the same value of λ and tighter accuracy tolerances, needed to bring the accuracy (that is, the value of λ) to an acceptable value. Its value should normally be 1.0, and should almost never exceed 2.0.
finfo4
The number of calls to coeffn at this iteration.
finfo5
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.
finfo6
The number of unsuccessful steps used by the internal integrator at this iteration.
finfo7
The number of successful steps at the maximum step size taken by the internal integrator at this iteration.
finfo8
Not used.
finfo9 to finfo15
Set to zero, unless iflag<0 in which case they hold the following values describing the point of failure:
finfo9
1 or 2 depending on whether integration was in a forward or backward direction at the time of failure.
finfo10
The value of the independent variable, x, the point at which the error occurred.
finfo11, finfo12, finfo13
The current values of the Pruefer dependent variables β, ϕ and ρ respectively. See Section 3 in d02kef for a description of these variables.
finfo14
The local-error tolerance being used by the internal integrator at the point of failure.
finfo15
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. Arguments denoted as Input must not be changed by this procedure.
10: ifail Integer Input/Output
On entry: ifail must be set to 0, -1 or 1 to set behaviour on detection of an error; these values have no effect when no error is detected.
A value of 0 causes the printing of an error message and program execution will be halted; otherwise program execution continues. A value of -1 means that an error message is printed while a value of 1 means that it is not.
If halting is not appropriate, the value -1 or 1 is recommended. If message printing is undesirable, then the value 1 is recommended. Otherwise, the value 0 is recommended. When the value -1 or 1 is used it is essential to test the value of ifail on exit.
On exit: ifail=0 unless the routine detects an error or a warning has been flagged (see Section 6).

6 Error Indicators and Warnings

If on entry ifail=0 or -1, explanatory error messages are output on the current error message unit (as defined by x04aaf).
Errors or warnings detected by the routine:
ifail=1
On entry, k=value.
Constraint: k0.
On entry, tol=value.
Constraint: tol>0.0.
ifail=2
The constraint on the left- or right-hand boundary condition has been violated.
ifail=3
The function px became zero or changed sign in the interval a,b.
ifail=4
After 15 iterations the eigenvalue had not been found to the required accuracy.
ifail=5
The bracketing phase failed to bracket the eigenvalue within ten iterations. This may be due to an error in formulating the problem (for example, q is independent of λ), or by very poor initial estimates for the eigenvalue and search step.
On exit, elam and elam+delam give the end points of the interval within which no eigenvalue was located by the routine.
ifail=6
Too small a step size was required at the start of a sub-interval to obtain the desired accuracy.
ifail=7
Too small a step size was required during solution in a sub-interval to obtain the desired accuracy.
At some point inside a sub-interval 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 ifail=7). This could be due to pathological behaviour of px and qx;λ or to an unreasonable accuracy requirement or to the current value of λ making the equations ‘stiff’. See the last call of monit for details.
ifail=8
The tolerance set is too small for the problem being solved and the machine precision being used. The local eigenvalue estimate returned is likely to be a very good approximation.
ifail=9
An internal error has occurred corresponding to a pole of the matching function. Try solving the problem again with a smaller tolerance.
ifail=-99
An unexpected error has been triggered by this routine. Please contact NAG.
See Section 7 in the Introduction to the NAG Library FL Interface for further information.
ifail=-399
Your licence key may have expired or may not have been installed correctly.
See Section 8 in the Introduction to the NAG Library FL Interface for further information.
ifail=-999
Dynamic memory allocation failed.
See Section 9 in the Introduction to the NAG Library FL Interface for further information.

7 Accuracy

The absolute accuracy of the computed eigenvalue is usually within a mixed absolute/relative bound defined by tol (as defined above).

8 Parallelism and Performance

d02kaf is not threaded in any implementation.

9 Further Comments

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 9 in d02kdf for a discussion of the technique used.

10 Example

This example finds the fourth eigenvalue of Mathieu's equations
y+λ-2qcos2xy=0  
with boundary conditions
y0=yπ =0  
and q=5. We use a starting value elam=15.0 and a step delam=4.0. We illustrate the effect of varying tol by choosing tol=1.0E−5 and 1.0E−6 (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.

10.1 Program Text

Program Text (d02kafe.f90)

10.2 Program Data

Program Data (d02kafe.d)

10.3 Program Results

Program Results (d02kafe.r)