NAG FL Interface
d02kdf (sl2_​breaks_​vals)

1 Purpose

2 Specification

3 Description

4 References

5 Arguments

1: $\mathbf{xpoint}\left(\mathbf{m}\right)$ – Real (Kind=nag_wp) array Input

On entry: the points where the boundary conditions computed by bdyval are to be imposed, and also any break-points, i.e., $\mathbf{xpoint}\left(1\right)$ to $\mathbf{xpoint}\left(m\right)$ must contain values $x_1,\dots ,x_m$ such that

$$x_1\le x_2<x_3<\cdots <x_{m-1}\le x_m$$

with the following meanings:

1. (a)$x_1$ and $x_m$ are the left- and right-hand end points, $a$ and $b$, of the domain of definition of the Sturm–Liouville system if these are finite. If either $a$ or $b$ is infinite, the corresponding value $x_1$ or $x_m$ may be a more-or-less arbitrarily ‘large’ number of appropriate sign.
2. (b)$x_2$ and $x_{m-1}$ are the Boundary Matching Points (BMPs), that is the points at which the left and right boundary conditions computed in bdyval are imposed.
   If the left-hand end point is a regular point then you should set $x_2=x_1$ $(=a)$, while if it is a singular point you must set $x_2>x_1$. Similarly $x_{m-1}=x_m$ ($=b$) if the right-hand end point is regular, and $x_{m-1}<x_m$ if it is singular.
3. (c)The remaining $m-4$ points $x_3,\dots ,x_{m-2}$, if any, define ‘break-points’ which divide the interval $[x_2,x_{m-1}]$ into $m-3$ sub-intervals

   $$i_1=[x_2,x_3],\dots ,i_{m-3}=[x_{m-2},x_{m-1}]\text{.}$$

   Numerical integration of the differential equation is stopped and restarted at each break-point. In simple cases no break-points are needed. However, if $p\left(x\right)$ or $q(x;\lambda )$ are given by different formulae in different parts of the interval, integration is more efficient if the range is broken up by break-points in the appropriate way. Similarly points where any jumps occur in $p\left(x\right)$ or $q(x;\lambda )$, or in their derivatives up to the fifth-order, should appear as break-points.
   Examples are given in Section 9. xpoint determines the position of the Shooting Matching Point (SMP), as explained in Section 9.3.

Constraint: $\mathbf{xpoint}\left(1\right)\le \mathbf{xpoint}\left(2\right)<\cdots <\mathbf{xpoint}\left(\mathbf{m}-1\right)\le \mathbf{xpoint}\left(\mathbf{m}\right)$.

2: $\mathbf{m}$ – Integer Input

On entry: the number of points in the array xpoint.

Constraint: $\mathbf{m}\ge 4$.

3: $\mathbf{coeffn}$ – Subroutine, supplied by the user. External Procedure

coeffn must compute the values of the coefficient functions $p\left(x\right)$ and $q(x;\lambda )$ for given values of $x$ and $\lambda$. Section 3 states the conditions which $p$ and $q$ must satisfy. See Section 9.4 for examples.

The specification of coeffn is:

Fortran Interface copy

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 copy

void  coeffn (double *p, double *q, double *dqdl, const double *x, const double *elam, const Integer *jint)

C++ Header Interface copy

#include <nag.h> extern "C" { void  coeffn (double &p, double &q, double &dqdl, const double &x, const double &elam, const Integer &jint) }

1: $\mathbf{p}$ – Real (Kind=nag_wp) Output

On exit: the value of $p\left(x\right)$ for the current value of $x$.

2: $\mathbf{q}$ – Real (Kind=nag_wp) Output

On exit: the value of $q(x;\lambda )$ for the current value of $x$ and the current trial value of $\lambda$.

3: $\mathbf{dqdl}$ – Real (Kind=nag_wp) Output

On exit: the value of $\frac{\partial q}{\partial \lambda }(x;\lambda )$ 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: $\mathbf{x}$ – Real (Kind=nag_wp) Input

On entry: the current value of $x$.

5: $\mathbf{elam}$ – Real (Kind=nag_wp) Input

On entry: the current trial value of the eigenvalue argument $\lambda$.

6: $\mathbf{jint}$ – Integer Input

On entry: the index $j$ of the sub-interval $i_j$ (see specification of xpoint) in which $x$ lies.

coeffn must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which d02kdf 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 d02kdf. If your code inadvertently does return any NaNs or infinities, d02kdf is likely to produce unexpected results.

4: $\mathbf{bdyval}$ – Subroutine, supplied by the user. External Procedure

bdyval must define the boundary conditions. For each end point, bdyval must return (in yl or yr) values of $y\left(x\right)$ and $p\left(x\right)y^{\prime }\left(x\right)$ which are consistent with the boundary conditions at the end points; only the ratio of the values matters. Here $x$ is a given point (xl or xr) equal to, or close to, the end point.

For a regular end point ($a$, say), $x=a$, a boundary condition of the form

$$c_{1}y\left(a\right)+c_2{y}^{\prime }\left(a\right)=0$$

can be handled by returning constant values in yl, e.g., $\mathbf{yl}\left(1\right)=c_2$ and $\mathbf{yl}\left(2\right)=-c_{1}p\left(a\right)$.

For a singular end point however, $\mathbf{yl}\left(1\right)$ and $\mathbf{yl}\left(2\right)$ will in general be functions of xl and elam, and $\mathbf{yr}\left(1\right)$ and $\mathbf{yr}\left(2\right)$ functions of xr and elam, usually derived analytically from a power-series or asymptotic expansion. Examples are given in Section 9.5.

The specification of bdyval is:

Fortran Interface copy

Subroutine bdyval ( xl, xr, elam, yl, yr) Real (Kind=nag_wp), Intent (In) :: xl, xr, elam Real (Kind=nag_wp), Intent (Out) :: yl(3), yr(3)

C Header Interface copy

void  bdyval (const double *xl, const double *xr, const double *elam, double yl[], double yr[])

C++ Header Interface copy

#include <nag.h> extern "C" { void  bdyval (const double &xl, const double &xr, const double &elam, double yl[], double yr[]) }

1: $\mathbf{xl}$ – Real (Kind=nag_wp) Input

On entry: if $a$ is a regular end point of the system (so that $a=x_1=x_2$), xl contains $a$. If $a$ is a singular point (so that $a\le x_1<x_2$), xl contains a point $x$ such that $x_1<x\le x_2$.

2: $\mathbf{xr}$ – Real (Kind=nag_wp) Input

On entry: if $b$ is a regular end point of the system (so that $x_{m-1}=x_m=b$), xr contains $b$. If $b$ is a singular point (so that $x_{m-1}<x_m\le b$), xr contains a point $x$ such that $x_{m-1}\le x<x_m$.

3: $\mathbf{elam}$ – Real (Kind=nag_wp) Input

On entry: the current trial value of $\lambda$.

4: $\mathbf{yl}\left(3\right)$ – Real (Kind=nag_wp) array Output

On exit: $\mathbf{yl}\left(1\right)$ and $\mathbf{yl}\left(2\right)$ should contain values of $y\left(x\right)$ and $p\left(x\right)y^{\prime }\left(x\right)$ respectively (not both zero) which are consistent with the boundary condition at the left-hand end point, given by $x=\mathbf{xl}$. $\mathbf{yl}\left(3\right)$ should not be set.

5: $\mathbf{yr}\left(3\right)$ – Real (Kind=nag_wp) array Output

On exit: $\mathbf{yr}\left(1\right)$ and $\mathbf{yr}\left(2\right)$ should contain values of $y\left(x\right)$ and $p\left(x\right)y^{\prime }\left(x\right)$ respectively (not both zero) which are consistent with the boundary condition at the right-hand end point, given by $x=\mathbf{xr}$. $\mathbf{yr}\left(3\right)$ should not be set.

bdyval must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which d02kdf is called. Arguments denoted as Input must not be changed by this procedure.

Note: bdyval should not return floating-point NaN (Not a Number) or infinity values, since these are not handled by d02kdf. If your code inadvertently does return any NaNs or infinities, d02kdf is likely to produce unexpected results.

5: $\mathbf{k}$ – Integer Input

On entry: $k$, the index of the required eigenvalue when the eigenvalues are ordered

$${\lambda }_0<{\lambda }_1<{\lambda }_2<\cdots <{\lambda }_k<\cdots \text{.}$$

Constraint: $\mathbf{k}\ge 0$.

6: $\mathbf{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

$$\mathbf{delam}\le \mathbf{tol}\times \max (1.0,\left|\mathbf{elam}\right|)\text{,}$$

(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: $\mathbf{tol}>0.0$.

7: $\mathbf{elam}$ – Real (Kind=nag_wp) Input/Output

On entry: an initial estimate of the eigenvalue $\tilde{\lambda }$.

On exit: the final computed estimate, whether or not an error occurred.

8: $\mathbf{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 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 \max (1.0,\left|\mathbf{elam}\right|)$.

On exit: if $\mathbf{ifail}=\mathbf{0}$, delam holds an estimate of the absolute error in the computed eigenvalue, that is $|\tilde{\lambda }-\mathbf{elam}|\simeq \mathbf{delam}$. (In Section 9.2 we discuss the assumptions under which this is true.) The true error is rarely more than twice, or less than a tenth, of the estimated error.

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: $\mathbf{hmax}(2,\mathbf{m})$ – Real (Kind=nag_wp) array Input/Output

On entry: $\mathbf{hmax}(1,\mathit{j})$ should contain a maximum step size to be used by the differential equation code in the $\mathit{j}$th sub-interval ${\mathit{i}}_{\mathit{j}}$ (as described in the specification of argument xpoint), for $\mathit{j}=1,2,\dots ,m-3$. If it is zero the routine generates a maximum step size internally.

It is recommended that $\mathbf{hmax}(1,j)$ be set to zero unless the coefficient functions $p$ and $q$ have features (such as a narrow peak) within the $j$th sub-interval that could be ‘missed’ if a long step were taken. In such a case $\mathbf{hmax}(1,j)$ should be set to about half the distance over which the feature should be observed. Too small a value will increase the computing time for the routine. See Section 9 for further suggestions.

The rest of the array is used as workspace.

On exit: $\mathbf{hmax}(1,m-1)$ and $\mathbf{hmax}(1,m)$ contain the sensitivity coefficients ${\sigma }_l,{\sigma }_r$, described in Section 9.6. Other entries contain diagnostic output in the case of an error exit (see Section 6).

10: $\mathbf{maxit}$ – Integer Input/Output

On entry: a bound on $n_r$, the number of root-finding iterations allowed, that is the number of trial values of $\lambda$ that are used. If $\mathbf{maxit}\le 0$, no such bound is assumed. (See also maxfun.)

Suggested value: $\mathbf{maxit}=0$.

On exit: will have been decreased by the number of iterations actually performed, whether or not it was positive on entry.

11: $\mathbf{maxfun}$ – Integer Input

On entry: a bound on $n_f$, the number of calls to coeffn made in any one root-finding iteration. If $\mathbf{maxfun}\le 0$, no such bound is assumed.

Suggested value: $\mathbf{maxfun}=0$.

maxfun and maxit may be used to limit the computational cost of a call to d02kdf, which is roughly proportional to $n_r\times n_f$.

12: $\mathbf{monit}$ – Subroutine, supplied by the NAG Library or the user. External Procedure

monit is called by d02kdf at the end of each root-finding iteration and allows you to monitor the course of the computation by printing out the arguments.

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 copy

Subroutine monit ( nit, iflag, elam, finfo) Integer, Intent (In) :: nit, iflag Real (Kind=nag_wp), Intent (In) :: elam, finfo(15)

C Header Interface copy

void  monit (const Integer *nit, const Integer *iflag, const double *elam, const double finfo[])

C++ Header Interface copy

#include <nag.h> extern "C" { void  monit (const Integer &nit, const Integer &iflag, const double &elam, const double finfo[]) }

1: $\mathbf{nit}$ – Integer Input

On entry: the current value of the argument maxit of d02kdf, this is decreased by one at each iteration.

2: $\mathbf{iflag}$ – Integer Input

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 d02kdf with $\mathbf{ifail}=-\mathbf{iflag}$ will follow.

$\mathbf{iflag}=1$

The routine is trying to bracket the eigenvalue $\tilde{\lambda }$.

$\mathbf{iflag}=2$

The routine is converging to the eigenvalue $\tilde{\lambda }$ (having already bracketed it).

3: $\mathbf{elam}$ – Real (Kind=nag_wp) Input

On entry: the current trial value of $\lambda$.

4: $\mathbf{finfo}\left(15\right)$ – 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 $\mathbf{iflag}>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 ‘miss-distance’ or ‘residual’ function $f\left(\lambda \right)$ on which the shooting method is based. (See Section 9.2 for further information.) $\mathbf{finfo}\left(1\right)$ 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 miss-distance $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)$

The index of the sub-interval where failure occurred, in the range $1$ to $m-3$. In case of an error in bdyval, it is set to $0$ or $m-2$ depending on whether the left or right boundary condition caused the error.

$\mathbf{finfo}\left(10\right)$

The value of the independent variable, $x$, the point at which the error occurred. In case of an error in bdyval, it is set to the value of xl or xr as appropriate (see the specification of bdyval).

$\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. These are set to zero in case of an error in bdyval. (See d02kef for a description of these variables.)

$\mathbf{finfo}\left(14\right)$

The local-error tolerance being used by the internal integrator at the point of failure. This is set to zero in the case of an error in bdyval.

$\mathbf{finfo}\left(15\right)$

The last integration mesh point. This is set to zero in the case of an error in bdyval.

monit must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which d02kdf is called. Arguments denoted as Input must not be changed by this procedure.

13: $\mathbf{ifail}$ – Integer Input/Output

On entry: ifail must be set to $0$, $\mathrm{-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 $\mathrm{-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 $\mathrm{-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 $-\mathbf{1}$ 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

$\mathbf{ifail}=1$

On entry, $\mathbf{k}=⟨\mathit{\text{value}}⟩$.
Constraint: $\mathbf{k}\ge 0$.

On entry, $\mathbf{m}=⟨\mathit{\text{value}}⟩$.
Constraint: $\mathbf{m}\ge 4$.

On entry, the sequence of boundary and break points is not monotonic.
For $i=⟨\mathit{\text{value}}⟩$, point $i-1=⟨\mathit{\text{value}}⟩$ and point $i=⟨\mathit{\text{value}}⟩$.

On entry, $\mathbf{tol}=⟨\mathit{\text{value}}⟩$.
Constraint: $\mathbf{tol}>0.0$.

$\mathbf{ifail}=2$

The constraint on the specification of left- or right-hand boundary condition has been violated.

At some call to bdyval, invalid values were returned, that is, either $\mathbf{yl}\left(1\right)=\mathbf{yl}\left(2\right)=0.0$, or $\mathbf{yr}\left(1\right)=\mathbf{yr}\left(2\right)=0.0$ (a programming error in bdyval). See the last call of monit for details.This error exit will also occur if $p\left(x\right)$ is zero at the point where the boundary condition is imposed. Probably bdyval was called with xl equal to a singular end point $a$ or with xr equal to a singular end point $b$.

$\mathbf{ifail}=3$

The function $p\left(x\right)$ became zero or changed sign in the interval $[a,b]$.

$\mathbf{ifail}=4$

After $⟨\mathit{\text{value}}⟩$ iterations the eigenvalue had not been found to the required accuracy.

$\mathbf{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 $\lambda$), or by very poor initial estimates for the eigenvalue and search step.

On exit, $\mathbf{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$

The last iteration was terminated because more than $⟨\mathit{\text{value}}⟩$ calls to evaluate $p$, $q$ and its derivative were performed.

$\mathbf{ifail}=7$

Too small a step size was required at the start of a sub-interval to obtain the desired accuracy.

$\mathbf{ifail}=8$

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 $\mathbf{ifail}=\mathbf{7}$). This could be due to pathological behaviour of $p\left(x\right)$ and $q(x;\lambda )$ 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}=9$

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.

$\mathbf{ifail}=10$

An internal error has occurred corresponding to a pole of the matching function. Try solving the problem again with a smaller tolerance.

$\mathbf{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.

$\mathbf{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.

$\mathbf{ifail}=-999$

Dynamic memory allocation failed.

See Section 9 in the Introduction to the NAG Library FL Interface for further information.

7 Accuracy

8 Parallelism and Performance

9 Further Comments

9.1 Timing

9.2 General Description of the Algorithm

9.3 The Position of the Shooting Matching Point $\mathit{c}$

9.4 Examples of Coding the coeffn

Subroutine coeffn(p,q,dqdl,x,elam,jint)
...
p = x
q = elam*x - nu*nu/x
dqdl = x
Return
End

Subroutine coeffn(p,q,dqdl,x,elam,jint)
...
If (jint.eq.2) Then
  q = elam + x*x - 10.0E0
Else
  q = elam + 6.0E0/abs(x)
Endif
p = 1.0E0
dqdl = 1.0E0
Return
End

Subroutine coeffn(p,q,dqdl,x,elam,jint)
...
p = 1.0E0
dqdl = 1.0E0 - 2.0E0*exp(-100.0E0*x*x)
q = elam*dqdl
Return
End

9.5 Examples of Boundary Conditions at Singular Points

Subroutine bdyval(xl,xr,elam,yl,yr)
Real xl, xr, elam, yl(3), yr(3)
yl(1) = xl
yl(2) = 2.0E0
yr(1) = 1.0E0
yr(2) = -sqrt(xr-elam)
Return
End

9.6 The Sensitivity Parameters ${\mathit{\sigma }}_{\mathit{l}}$ and ${\mathit{\sigma }}_{\mathit{r}}$

9.7 ‘Missed Zeros’

10 Example

10.1 Program Text

10.2 Program Data

10.3 Program Results