naginterfaces.library.ode.sl2_breaks_vals¶
- naginterfaces.library.ode.sl2_breaks_vals(xpoint, coeffn, bdyval, k, tol, elam, delam, hmax, maxit=0, maxfun=0, monit=None, data=None)[source]¶
sl2_breaks_vals
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. Provision is made for discontinuities in the coefficient functions or their derivatives.For full information please refer to the NAG Library document for d02kd
https://www.nag.com/numeric/nl/nagdoc_29.2/flhtml/d02/d02kdf.html
- Parameters
- xpointfloat, array-like, shape
The points where the boundary conditions computed by are to be imposed, and also any break-points, i.e., to must contain values such that
with the following meanings:
and are the left- and right-hand end points, and , of the domain of definition of the Sturm–Liouville system if these are finite. If either or is infinite, the corresponding value or may be a more-or-less arbitrarily ‘large’ number of appropriate sign.
and are the Boundary Matching Points (BMPs), that is the points at which the left and right boundary conditions computed in are imposed.
If the left-hand end point is a regular point then you should set , while if it is a singular point you must set .
Similarly () if the right-hand end point is regular, and if it is singular.
The remaining points , if any, define ‘break-points’ which divide the interval into sub-intervals
Numerical integration of the differential equation is stopped and restarted at each break-point. In simple cases no break-points are needed. However, if or 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 or , or in their derivatives up to the fifth-order, should appear as break-points.
Examples are given in Further Comments. determines the position of the Shooting Matching Point (SMP), as explained in The Position of the Shooting Matching Point c.
- coeffncallable (p, q, dqdl) = coeffn(x, elam, jint, data=None)
must compute the values of the coefficient functions and for given values of and . Notes states the conditions which and must satisfy.
See Further Comments for examples.
- Parameters
- xfloat
The current value of .
- elamfloat
The current trial value of the eigenvalue argument .
- jintint
The index of the sub-interval (see specification of ) in which lies.
- dataarbitrary, optional, modifiable in place
User-communication data for callback functions.
- Returns
- pfloat
The value of for the current value of .
- qfloat
The value of for the current value of and the current trial value of .
- dqdlfloat
The value of for the current value of and the current trial value of . However is only used in error estimation and, in the rare cases where it may be difficult to evaluate, an approximation (say to within ) will suffice.
- bdyvalcallable (yl, yr) = bdyval(xl, xr, elam, data=None)
must define the boundary conditions.
For each end point, must return (in or ) values of and which are consistent with the boundary conditions at the end points; only the ratio of the values matters.
Here is a given point ( or ) equal to, or close to, the end point.
For a regular end point (, say), , a boundary condition of the form
can be handled by returning constant values in , e.g., and .
For a singular end point however, and will in general be functions of and , and and functions of and , usually derived analytically from a power-series or asymptotic expansion.
Examples are given in Further Comments.
- Parameters
- xlfloat
If is a regular end point of the system (so that ), contains . If is a singular point (so that ), contains a point such that .
- xrfloat
If is a regular end point of the system (so that ), contains . If is a singular point (so that ), contains a point such that .
- elamfloat
The current trial value of .
- dataarbitrary, optional, modifiable in place
User-communication data for callback functions.
- Returns
- ylfloat, array-like, shape
and should contain values of and respectively (not both zero) which are consistent with the boundary condition at the left-hand end point, given by . should not be set.
- yrfloat, array-like, shape
and should contain values of and respectively (not both zero) which are consistent with the boundary condition at the right-hand end point, given by . should not be set.
- kint
, the index of the required eigenvalue when the eigenvalues are ordered
- tolfloat
The tolerance argument which determines the accuracy of the computed eigenvalue. The error estimate held in on exit satisfies the mixed absolute/relative error test
where is the final estimate of the eigenvalue. 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 .
- delamfloat
An indication of the scale of the problem in the -direction. holds the initial ‘search step’ (positive or negative). Its value is not critical, but the first two trial evaluations are made at and , 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 and .
If on entry, it is given the default value of .
- hmaxfloat, array-like, shape
should contain a maximum step size to be used by the differential equation code in the th sub-interval (as described in the specification of argument ), for . If it is zero the function generates a maximum step size internally.
It is recommended that be set to zero unless the coefficient functions and have features (such as a narrow peak) within the th sub-interval that could be ‘missed’ if a long step were taken.
In such a case 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.
- maxitint, optional
A bound on , the number of root-finding iterations allowed, that is the number of trial values of that are used. If , no such bound is assumed. (See also .)
- maxfunint, optional
A bound on , the number of calls to made in any one root-finding iteration. If , 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.
is called by
sl2_breaks_vals
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 of
sl2_breaks_vals
, this is decreased by one at each iteration.- iflagint
Describes what phase the computation is in.
An error occurred in the computation at this iteration; an error exit from
sl2_breaks_vals
with will follow.The function is trying to bracket the eigenvalue .
The function is converging to the eigenvalue (having already bracketed it).
- elamfloat
The current trial value of .
- finfofloat, ndarray, shape
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 ), though the first few components may be of interest to you. In case of an error () all the components of should be printed.
The contents of are as follows:
The current value of the ‘miss-distance’ or ‘residual’ function on which the shooting method is based. (See Further Comments for further information.) is set to zero if .
An estimate of the quantity defined as follows. Consider the perturbation in the miss-distance 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 . Thus, at the zero of 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 such that is independent of . If this is the case, an error exit with = 5 should follow. is set to zero if .
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 , and should almost never exceed .
The number of calls to at this iteration.
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.
The number of unsuccessful steps used by the internal integrator at this iteration.
The number of successful steps at the maximum step size taken by the internal integrator at this iteration.
Not used.
to
Set to zero, unless in which case they hold the following values describing the point of failure:
The index of the sub-interval where failure occurred, in the range to . In case of an error in , it is set to or depending on whether the left or right boundary condition caused the error.
The value of the independent variable, , the point at which the error occurred. In case of an error in , it is set to the value of or as appropriate (see the specification of ).
, ,
The current values of the Pruefer dependent variables , and respectively. These are set to zero in case of an error in . (See
sl2_breaks_funs()
for a description of these variables.)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 .
The last integration mesh point. This is set to zero in the case of an error in .
- dataarbitrary, optional, modifiable in place
User-communication data for callback functions.
- dataarbitrary, optional
User-communication data for callback functions.
- Returns
- elamfloat
The final computed estimate, whether or not an error occurred.
- delamfloat
If no exception or warning is raised, holds an estimate of the absolute error in the computed eigenvalue, that is . (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, may hold an estimate of the error, or its initial value, depending on the value of .
See Exceptions for further details.
- hmaxfloat, ndarray, shape
and contain the sensitivity coefficients , 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 )
On entry, the sequence of boundary and break points is not monotonic.
For , point and point .
- (errno )
On entry, .
Constraint: .
- (errno )
On entry, .
Constraint: .
- (errno )
On entry, .
Constraint: .
- (errno )
The constraint on the specification of left- or right-hand boundary condition has been violated.
- (errno )
The function became zero or changed sign in the interval .
- (errno )
The last iteration was terminated because more than calls to evaluate , and its derivative were performed.
- (errno )
Too small a step size was required at the start of a sub-interval to obtain the desired accuracy.
- (errno )
Too small a step size was required during solution in a sub-interval to obtain the desired accuracy.
- Warns
- NagAlgorithmicWarning
- (errno )
After iterations the eigenvalue had not been found to the required accuracy.
- (errno )
The bracketing phase failed to bracket the eigenvalue within ten iterations. This may be due to an error in formulating the problem (for example, is independent of ), or by very poor initial estimates for the eigenvalue and search step.
- (errno )
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.
- (errno )
An internal error has occurred corresponding to a pole of the matching function. Try solving the problem again with a smaller tolerance.
- Notes
No equivalent traditional C interface for this routine exists in the NAG Library.
sl2_breaks_vals
finds a specified eigenvalue of a Sturm–Liouville system defined by a self-adjoint differential equation of the second-ordertogether with appropriate boundary conditions at the two, finite or infinite, end points and . The functions and , which are real-valued, are defined by . The boundary conditions must be defined by , and, in the case of a singularity at or , take the form of an asymptotic formula for the solution near the relevant end point.
For the theoretical basis of the numerical method to be valid, the following conditions should hold on the coefficient functions:
must be nonzero and must not change sign throughout the interval ; and,
must not change sign throughout the interval for all relevant values of , and must not be identically zero as varies, for any .
Points of discontinuity in the functions and or their derivatives are allowed, and should be included as ‘break-points’ in the array .
The eigenvalue is determined by a shooting method based on the Scaled Pruefer form of the differential equation as described in Pryce (1981), with certain modifications. The Pruefer equations are integrated by a special internal function using Merson’s Runge–Kutta formula with automatic control of local error. Providing certain assumptions (see Timing) are met, the computed value of will be correct to within a mixed absolute/relative error specified by .
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))