naginterfaces.library.ode.sl2_​breaks_​funs

naginterfaces.library.ode.sl2_breaks_funs(xpoint, match, coeffn, bdyval, k, tol, elam, delam, hmax, report, maxit=0, maxfun=0, monit=None, data=None)

sl2_breaks_funs finds a specified eigenvalue of a regular or singular second-order Sturm–Liouville system on a finite or infinite interval, using a Pruefer transformation and a shooting method. It also reports values of the eigenfunction and its derivatives. Provision is made for discontinuities in the coefficient functions or their derivatives.

For full information please refer to the NAG Library document for d02ke

https://support.nag.com/numeric/nl/nagdoc_30.3/flhtml/d02/d02kef.html

Parameters

xpointfloat, array-like, shape $\left(m\right)$

The points where the boundary conditions computed by $bdyval$ are to be imposed, and also any break-points, i.e., $xpoint\left[0\right]$ to $xpoint[m-1]$ 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. $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. $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. 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 Further Comments. $xpoint$ determines the position of the Shooting Matching Point (SMP), as explained in The Position of the Shooting Matching Point c.

matchint

Must be set to the index of the ‘break-point’ to be used as the matching point (see The Position of the Shooting Matching Point c). If $match$ is set to a value outside the range $[2,m-1]$ then a default value is taken, corresponding to the break-point nearest the centre of the interval $[xpoint\left[1\right],xpoint[m-1]]$.

coeffncallable (p, q, dqdl) = coeffn(x, elam, jint, data=None)

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

See Further Comments for examples.

Parameters

xfloat

The current value of $x$.

elamfloat

The current trial value of the eigenvalue argument $\lambda$.

jintint

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

dataarbitrary, optional, modifiable in place

User-communication data for callback functions.

Returns

pfloat

The value of $p\left(x\right)$ for the current value of $x$.

qfloat

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

dqdlfloat

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.

bdyvalcallable (yl, yr) = bdyval(xl, xr, elam, data=None)

$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., $yl\left[0\right]=c_2$ and $yl\left[1\right]=-c_{1}p\left(a\right)$.

For a singular end point however, $yl\left[0\right]$ and $yl\left[1\right]$ will in general be functions of $xl$ and $elam$, and $yr\left[0\right]$ and $yr\left[1\right]$ functions of $xr$ and $elam$, usually derived analytically from a power-series or asymptotic expansion.

Examples are given in Further Comments.

Parameters

xlfloat

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$.

xrfloat

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$.

elamfloat

The current trial value of $\lambda$.

dataarbitrary, optional, modifiable in place

User-communication data for callback functions.

Returns

ylfloat, array-like, shape $\left(3\right)$

$yl\left[0\right]$ and $yl\left[1\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=xl$. $yl\left[2\right]$ should not be set.

yrfloat, array-like, shape $\left(3\right)$

$yr\left[0\right]$ and $yr\left[1\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=xr$. $yr\left[2\right]$ should not be set.

kint

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

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

tolfloat

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

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

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.

elamfloat

An initial estimate of the eigenvalue $\tilde{\lambda }$.

delamfloat

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 $elam+delam$, so the function 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 $delam=0.0$ on entry, it is given the default value of $0.25\times max(1.0,\left|elam\right|)$.

hmaxfloat, array-like, shape $(2,m)$

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

It is recommended that $hmax[0,j-1]$ 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 $hmax[0,j-1]$ 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 function.

See Further Comments for further suggestions.

The rest of the array is used as workspace.

reportcallable report(x, v, jint, data=None)

$report$ provides the means by which you may compute the eigenfunction $y\left(x\right)$ and its derivative at each integration mesh point $x$. (See Further Comments for an example.)

Parameters

xfloat

The current value of the independent variable $x$. See The Position of the Shooting Matching Point c for the order in which values of $x$ are supplied.

vfloat, ndarray, shape $\left(3\right)$

$v\left[0\right]$, $v\left[1\right]$, $v\left[2\right]$ hold the current values of the Pruefer variables $\beta$, $\varphi$, $\rho$ respectively.

jintint

Indicates the sub-interval between break-points in which $x$ lies exactly as for $coeffn$, except that at the extreme left-hand end point (when $x=xpoint\left[1\right]$) $jint$ is set to $0$ and at the extreme right-hand end point (when $x=x_r=xpoint[m-1]$) $jint$ is set to $m-2$.

dataarbitrary, optional, modifiable in place

User-communication data for callback functions.

maxitint, optional

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 $maxit\le 0$, no such bound is assumed. (See also $maxfun$.)

maxfunint, optional

A bound on $n_f$, the number of calls to $coeffn$ made in any one root-finding iteration. If $maxfun\le 0$, no such bound is assumed.

monitNone or callable monit(nit, iflag, elam, finfo, data=None), optional

Note: if this argument is None then a NAG-supplied facility will be used.

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

Parameters

nitint

The current value of the argument $maxit$ of sl2_breaks_funs, this is decreased by one at each iteration.

iflagint

Describes what phase the computation is in.

$iflag<0$

An error occurred in the computation at this iteration; an error exit from sl2_breaks_funs with $\text{errno}=-iflag$ will follow.

$iflag=1$

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

$iflag=2$

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

elamfloat

The current trial value of $\lambda$.

finfofloat, ndarray, shape $\left(15\right)$

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 $iflag>0$), 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:

$finfo\left[0\right]$

The current value of the ‘miss-distance’ or ‘residual’ function $f\left(\lambda \right)$ on which the shooting method is based. (See Further Comments for further information.) $finfo\left[0\right]$ is set to zero if $iflag<0$.

$finfo\left[1\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 $errno$ = 5 should follow. $finfo\left[1\right]$ is set to zero if $iflag<0$.

$finfo\left[2\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$.

$finfo\left[3\right]$

The number of calls to $coeffn$ at this iteration.

$finfo\left[4\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.

$finfo\left[5\right]$

The number of unsuccessful steps used by the internal integrator at this iteration.

$finfo\left[6\right]$

The number of successful steps at the maximum step size taken by the internal integrator at this iteration.

$finfo\left[7\right]$

Not used.

$finfo\left[8\right]$ to $finfo\left[14\right]$

Set to zero, unless $iflag<0$ in which case they hold the following values describing the point of failure:

$finfo\left[8\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.

$finfo\left[9\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$).

$finfo\left[10\right]$, $finfo\left[11\right]$, $finfo\left[12\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$.

$finfo\left[13\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$.

$finfo\left[14\right]$

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

dataarbitrary, optional, modifiable in place

User-communication data for callback functions.

dataarbitrary, optional

User-communication data for callback functions.

Returns

matchint

The index of the break-point actually used as the matching point.

elamfloat

The final computed estimate, whether or not an error occurred.

delamfloat

If no exception or warning is raised, $delam$ holds an estimate of the absolute error in the computed eigenvalue, that is $|\tilde{\lambda }-elam|\simeq delam$. (In Further Comments 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 an exception is raised, $delam$ may hold an estimate of the error, or its initial value, depending on the value of $\text{errno}$.

See Exceptions for further details.

hmaxfloat, ndarray, shape $(2,m)$

$hmax[0,m-2]$ and $hmax[0,m-1]$ contain the sensitivity coefficients ${\sigma }_l,{\sigma }_r$, described in Further Comments. Other entries contain diagnostic output in the case of an error exit (see Exceptions).

maxitint

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

Raises

NagValueError

(errno $1$)

On entry, $m=⟨value⟩$.

Constraint: $m\ge 4$.

(errno $1$)

On entry, $tol=⟨value⟩$.

Constraint: $tol>0.0$.

(errno $1$)

On entry, $k=⟨value⟩$.

Constraint: $k\ge 0$.

(errno $2$)

On entry, the sequence of boundary and break points is not monotonic.

For $i=⟨value⟩$, point $i-1=⟨value⟩$ and point $i=⟨value⟩$.

(errno $2$)

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

(errno $3$)

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

(errno $4$)

After $⟨value⟩$ iterations the eigenvalue had not been found to the required accuracy.

(errno $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.

(errno $6$)

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

(errno $7$)

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

(errno $8$)

Too small a step size was required during solution in a sub-interval to obtain the desired accuracy.

(errno $10$)

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

Warns

NagAlgorithmicWarning

(errno $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.

Notes

No equivalent traditional C interface for this routine exists in the NAG Library.

sl2_breaks_funs has essentially the same purpose as sl2_breaks_vals() with minor modifications to enable values of the eigenfunction to be obtained after convergence to the eigenvalue has been achieved.

It first finds a specified eigenvalue $\tilde{\lambda }$ of a Sturm–Liouville system defined by a self-adjoint differential equation of the second-order

$${\left(p\left(x\right)y^{\prime }\right)}^{\prime }+q(x;\lambda )y=0\text{,}a<x<b\text{,}$$

together with appropriate boundary conditions at the two, finite or infinite, end points $a$ and $b$. The functions $p$ and $q$, which are real-valued, are defined by $coeffn$. The boundary conditions must be defined by $bdyval$, and, in the case of a singularity at $a$ or $b$, take the form of an asymptotic formula for the solution near the relevant end point.

When the final estimate $\lambda =\tilde{\lambda }$ of the eigenvalue has been found, the function integrates the differential equation once more with that value of $\lambda$, and with initial conditions chosen so that the integral

$$S={\int }_a^{b}y{\left(x\right)}^2\frac{\partial q}{\partial \lambda }(x;\lambda )dx$$

is approximately one. When $q(x;\lambda )$ is of the form $\lambda w\left(x\right)+q\left(x\right)$, which is the most common case, $S$ represents the square of the norm of $y$ induced by the inner product

$$⟨f,g⟩={\int }_a^{b}f\left(x\right)g\left(x\right)w\left(x\right)dx\text{,}$$

with respect to which the eigenfunctions are mutually orthogonal. This normalization of $y$ is only approximate, but experience shows that $S$ generally differs from unity by only one or two per cent.

During this final integration the $report$ is called at each integration mesh point $x$. Sufficient information is returned to permit you to compute $y\left(x\right)$ and $y^{\prime }\left(x\right)$ for printing or plotting. For reasons described in Further Comments, sl2_breaks_funs passes across to $report$, not $y$ and $y^{\prime }$, but the Pruefer variables $\beta$, $\varphi$ and $\rho$ on which the numerical method is based. Their relationship to $y$ and $y^{\prime }$ is given by the equations

$$p\left(x\right)y^{\prime }=\sqrt{\beta }exp\left(\frac{\rho }{2}\right)\cos \left(\frac{\varphi }{2}\right)\text{,}y=\frac{1}{\sqrt{\beta }}exp\left(\frac{\rho }{2}\right)\sin \left(\frac{\varphi }{2}\right)\text{.}$$

For the theoretical basis of the numerical method to be valid, the following conditions should hold on the coefficient functions:

1. $p\left(x\right)$ must be nonzero and must not change sign throughout the interval $(a,b)$; and,

2. $\frac{\partial q}{\partial \lambda }$ must not change sign throughout the interval $(a,b)$ for all relevant values of $\lambda$, and must not be identically zero as $x$ varies, for any $\lambda$.

Points of discontinuity in the functions $p$ and $q$ or their derivatives are allowed, and should be included as ‘break-points’ in the array $xpoint$.

A good account of the theory of Sturm–Liouville systems, with some description of Pruefer transformations, is given in Module 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).

The scaled Pruefer method is described in a short note by Pryce and Hargrave (1977) and in some detail in the technical report by Pryce (1981).

References

Abramowitz, M and Stegun, I A, 1972, Handbook of Mathematical Functions, (3rd Edition), Dover Publications

Bailey, P B, 1966, Sturm–Liouville eigenvalues via a phase function, SIAM J. Appl. Math. (14), 242–249

Banks, D O and Kurowski, I, 1968, Computation of eigenvalues of singular Sturm–Liouville Systems, Math. Comput. (22), 304–310

Birkhoff, G and Rota, G C, 1962, Ordinary Differential Equations, Ginn & Co., Boston and New York

Pryce, J D, 1981, Two codes for Sturm–Liouville problems, Technical Report CS-81-01, Department of Computer Science, Bristol University

Pryce, J D and Hargrave, B A, 1977, The scaled Prüfer method for one-parameter and multi-parameter eigenvalue problems in ODEs, IMA Numerical Analysis Newsletter (1(3))