hide long namesshow long names
hide short namesshow short names
Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

NAG Toolbox: nag_ode_sl2_breaks_funs (d02ke)

Purpose

nag_ode_sl2_breaks_funs (d02ke) 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.

Syntax

[match, elam, delam, hmax, maxit, ifail] = d02ke(xpoint, match, coeffn, bdyval, k, tol, elam, delam, hmax, monit, report, 'm', m, 'maxit', maxit, 'maxfun', maxfun)
[match, elam, delam, hmax, maxit, ifail] = nag_ode_sl2_breaks_funs(xpoint, match, coeffn, bdyval, k, tol, elam, delam, hmax, monit, report, 'm', m, 'maxit', maxit, 'maxfun', maxfun)

Description

nag_ode_sl2_breaks_funs (d02ke) has essentially the same purpose as nag_ode_sl2_breaks_vals (d02kd) 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 λ̃λ~ of a Sturm–Liouville system defined by a self-adjoint differential equation of the second-order
(p(x)y) + q (x ; λ) y = 0 ,   a < x < b ,
( p(x) y ) + q (x;λ) y=0 ,   a<x<b ,
together with appropriate boundary conditions at the two, finite or infinite, end points aa and bb. The functions pp and qq, which are real-valued, are defined by coeffn. The boundary conditions must be defined by bdyval, and, in the case of a singularity at aa or bb, take the form of an asymptotic formula for the solution near the relevant end point.
When the final estimate λ = λ̃λ=λ~ of the eigenvalue has been found, the function integrates the differential equation once more with that value of λλ, and with initial conditions chosen so that the integral
b
S = y(x)2(q)/(λ)(x ; λ)dx
a
S=aby (x) 2 q λ (x;λ)dx
is approximately one. When q(x ; λ)q(x;λ) is of the form λw(x) + q(x)λw(x)+q(x), which is the most common case, SS represents the square of the norm of yy induced by the inner product
b
f,g = f(x)g(x)w(x)dx,
a
f,g = ab f(x) g(x) w(x) dx ,
with respect to which the eigenfunctions are mutually orthogonal. This normalization of yy is only approximate, but experience shows that SS generally differs from unity by only one or two per cent.
During this final integration the report is called at each integration mesh point xx. Sufficient information is returned to permit you to compute y(x)y(x) and y(x)y(x) for printing or plotting. For reasons described in Section [General Description of the Algorithm], nag_ode_sl2_breaks_funs (d02ke) passes across to report, not yy and yy, but the Pruefer variables ββ, φϕ and ρρ on which the numerical method is based. Their relationship to yy and yy is given by the equations
p(x) y = sqrt(β) exp(ρ/2) cos(φ/2) ,   y = 1/(sqrt(β)) exp(ρ/2) sin(φ/2) .
p(x) y = β exp(ρ2) cos(ϕ2) ,   y=1β exp(ρ2) sin(ϕ2) .
A specimen report is given in Section [Example] below.
For the theoretical basis of the numerical method to be valid, the following conditions should hold on the coefficient functions:
(a) p(x)p(x) must be nonzero and must not change sign throughout the interval (a,b)(a,b); and,
(b) (q)/(λ) q λ  must not change sign throughout the interval (a,b)(a,b) for all relevant values of λλ, and must not be identically zero as xx varies, for any λλ.
Points of discontinuity in the functions pp and qq 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 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).
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)

Parameters

Compulsory Input Parameters

1:     xpoint(m) – double array
m, the dimension of the array, must satisfy the constraint m4m4.
The points where the boundary conditions computed by bdyval are to be imposed, and also any break points, i.e., xpoint(1)xpoint1 to xpoint(m)xpointm must contain values x1,,xmx1,,xm such that
x1x2 < x3 < < xm1xm
x1x2<x3<<xm-1xm
with the following meanings:
(a) x1x1 and xmxm are the left- and right-hand end points, aa and bb, of the domain of definition of the Sturm–Liouville system if these are finite. If either aa or bb is infinite, the corresponding value x1x1 or xmxm may be a more-or-less arbitrarily ‘large’ number of appropriate sign.
(b) x2x2 and xm1xm-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 x2 = x1x2=x1 ( = a)(=a), while if it is a singular point you must set x2 > x1x2>x1. Similarly xm1 = xmxm-1=xm ( = b=b) if the right-hand end point is regular, and xm1 < xmxm-1<xm if it is singular.
(c) The remaining m4m-4 points x3,,xm2x3,,xm-2, if any, define ‘break points’ which divide the interval [x2,xm1][x2,xm-1] into m3m-3 sub-intervals
i1 = [x2,x3],,im3 = [xm2,xm1].
i1=[x2,x3],,im-3=[xm-2,xm-1].
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(x)p(x) or q(x ; λ)q(x;λ) are given by different formulae in different parts of the interval, then 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(x)p(x) or q(x ; λ)q(x;λ), or in their derivatives up to the fifth-order, should appear as break points.
Examples are given in Sections [Further Comments] and [Example]. xpoint determines the position of the Shooting Matching Point (SMP), as explained in Section [The Position of the Shooting Matching Point ].
Constraint: xpoint(1)xpoint(2) < < xpoint(m1)xpoint(m)xpoint1xpoint2<<xpointm-1xpointm.
2:     match – int64int32nag_int scalar
Must be set to the index of the ‘break point’ to be used as the matching point (see Section [The Position of the Shooting Matching Point ]). If match is set to a value outside the range [2,m1][2,m-1] then a default value is taken, corresponding to the break point nearest the centre of the interval [xpoint(2),xpoint(m1)][xpoint2,xpointm-1].
3:     coeffn – function handle or string containing name of m-file
coeffn must compute the values of the coefficient functions p(x)p(x) and q(x ; λ)q(x;λ) for given values of xx and λλ. Section [Description] states the conditions which pp and qq must satisfy. See Sections [Examples of Coding the coeffn] and [Example] for examples.
[p, q, dqdl] = coeffn(x, elam, jint)

Input Parameters

1:     x – double scalar
The current value of xx.
2:     elam – double scalar
The current trial value of the eigenvalue parameter λλ.
3:     jint – int64int32nag_int scalar
The index jj of the sub-interval ijij (see specification of xpoint) in which xx lies.

Output Parameters

1:     p – double scalar
The value of p(x)p(x) for the current value of xx.
2:     q – double scalar
The value of q(x ; λ)q(x;λ) for the current value of xx and the current trial value of λλ.
3:     dqdl – double scalar
The value of (q)/(λ) (x ; λ) q λ (x;λ) for the current value of xx 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%20%) will suffice.
4:     bdyval – function handle or string containing name of m-file
bdyval must define the boundary conditions. For each end point, bdyval must return (in yl or yr) values of y(x)y(x) and p(x)y(x)p(x)y(x) which are consistent with the boundary conditions at the end points; only the ratio of the values matters. Here xx is a given point (xl or xr) equal to, or close to, the end point.
For a regular end point (aa, say), x = ax=a, a boundary condition of the form
c1y(a) + c2y(a) = 0
c1y(a)+c2y(a)=0
can be handled by returning constant values in yl, e.g., yl(1) = c2yl1=c2 and yl(2) = c1p(a)yl2=-c1p(a).
For a singular end point however, yl(1)yl1 and yl(2)yl2 will in general be functions of xl and elam, and yr(1)yr1 and yr(2)yr2 functions of xr and elam, usually derived analytically from a power-series or asymptotic expansion. Examples are given in Sections [Examples of Boundary Conditions at Singular Points] and [Example].
[yl, yr] = bdyval(xl, xr, elam)

Input Parameters

1:     xl – double scalar
If aa is a regular end point of the system (so that a = x1 = x2a=x1=x2), then xl contains aa. If aa is a singular point (so that ax1 < x2ax1<x2), then xl contains a point xx such that x1 < xx2x1<xx2.
2:     xr – double scalar
If bb is a regular end point of the system (so that xm1 = xm = bxm-1=xm=b), then xr contains bb. If bb is a singular point (so that xm1 < xmbxm-1<xmb), then xr contains a point xx such that xm1x < xmxm-1x<xm.
3:     elam – double scalar
The current trial value of λλ.

Output Parameters

1:     yl(33) – double array
yl(1)yl1 and yl(2)yl2 should contain values of y(x)y(x) and p(x)y(x)p(x)y(x) respectively (not both zero) which are consistent with the boundary condition at the left-hand end point, given by x = xlx=xl. yl(3)yl3 should not be set.
2:     yr(33) – double array
yr(1)yr1 and yr(2)yr2 should contain values of y(x)y(x) and p(x)y(x)p(x)y(x) respectively (not both zero) which are consistent with the boundary condition at the right-hand end point, given by x = xrx=xr. yr(3)yr3 should not be set.
5:     k – int64int32nag_int scalar
kk, the index of the required eigenvalue when the eigenvalues are ordered
λ0 < λ1 < λ2 < < λk < .
λ0 < λ1 < λ2 < < λk < .
Constraint: k0k0.
6:     tol – double scalar
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
delamtol × max (1.0,|elam|),
delamtol×max(1.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.0tol>0.0.
7:     elam – double scalar
An initial estimate of the eigenvalue λ̃λ~.
8:     delam – double scalar
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 + delamelam+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.0delam=0.0 on entry, it is given the default value of 0.25 × max (1.0,|elam|)0.25×max(1.0,|elam|).
9:     hmax(22,m) – double array
hmax(1,j)hmax1j should contain a maximum step size to be used by the differential equation code in the jjth sub-interval ijij (as described in the specification of parameter xpoint), for j = 1,2,,m3j=1,2,,m-3. If it is zero the function generates a maximum step size internally.
It is recommended that hmax(1,j)hmax1j be set to zero unless the coefficient functions pp and qq have features (such as a narrow peak) within the jjth sub-interval that could be ‘missed’ if a long step were taken. In such a case hmax(1,j)hmax1j 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 Section [Further Comments] for further suggestions.
The rest of the array is used as workspace.
10:   monit – function handle or string containing name of m-file
monit is called by nag_ode_sl2_breaks_funs (d02ke) at the end of each root-finding iteration and allows you to monitor the course of the computation by printing out the parameters (see Section [Example] for an example).
If no monitoring is required, the dummy (sub)program nag_ode_sl2_reg_finite_dummy_monit (d02kay) may be used. (nag_ode_sl2_reg_finite_dummy_monit (d02kay) is included in the NAG Toolbox.)
monit(nit, iflag, elam, finfo)

Input Parameters

1:     nit – int64int32nag_int scalar
The current value of the parameter maxit of nag_ode_sl2_breaks_funs (d02ke), this is decreased by one at each iteration.
2:     iflag – int64int32nag_int scalar
Describes what phase the computation is in.
iflag < 0iflag<0
An error occurred in the computation at this iteration; an error exit from nag_ode_sl2_breaks_funs (d02ke) with ifail = iflagifail=-iflag will follow.
iflag = 1iflag=1
The function is trying to bracket the eigenvalue λ̃λ~.
iflag = 2iflag=2
The function is converging to the eigenvalue λ̃λ~ (having already bracketed it).
3:     elam – double scalar
The current trial value of λλ.
4:     finfo(1515) – double array
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 > 0iflag>0), though the first few components may be of interest to you. In case of an error (iflag < 0iflag<0) all the components of finfo should be printed.
The contents of finfo are as follows:
finfo(1)finfo1
The current value of the ‘miss-distance’ or ‘residual’ function f(λ)f(λ) on which the shooting method is based. (See Section [General Description of the Algorithm] for further information.) finfo(1)finfo1 is set to zero if iflag < 0iflag<0.
finfo(2)finfo2
An estimate of the quantity λλ defined as follows. Consider the perturbation in the miss-distance f(λ)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(λ)f(λ). Thus, at the zero of f(λ),|λ|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 qq is independent of λλ. If this is the case, an error exit with ifail = 5ifail=5 should follow. finfo(2)finfo2 is set to zero if iflag < 0iflag<0.
finfo(3)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.01.0, and should almost never exceed 2.02.0.
finfo(4)finfo4
The number of calls to coeffn at this iteration.
finfo(5)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.
finfo(6)finfo6
The number of unsuccessful steps used by the internal integrator at this iteration.
finfo(7)finfo7
The number of successful steps at the maximum step size taken by the internal integrator at this iteration.
finfo(8)finfo8
Not used.
finfo(9)finfo9 to finfo(15)finfo15
Set to zero, unless iflag < 0iflag<0 in which case they hold the following values describing the point of failure:
finfo(9)finfo9
The index of the sub-interval where failure occurred, in the range 11 to m3m-3. In case of an error in bdyval, it is set to 00 or m2m-2 depending on whether the left or right boundary condition caused the error.
finfo(10)finfo10
The value of the independent variable, xx, 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(11)finfo11, finfo(12)finfo12, finfo(13)finfo13
The current values of the Pruefer dependent variables ββ, φϕ and ρρ respectively. These are set to zero in case of an error in bdyval.
finfo(14)finfo14
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(15)finfo15
The last integration mesh point. This is set to zero in the case of an error in bdyval.
11:   report – function handle or string containing name of m-file
report provides the means by which you may compute the eigenfunction y(x)y(x) and its derivative at each integration mesh point xx. (See Section [Further Comments] for an example.)
report(x, v, jint)

Input Parameters

1:     x – double scalar
The current value of the independent variable xx. See Section [The Position of the Shooting Matching Point ] for the order in which values of xx are supplied.
2:     v(33) – double array
v(1)v1, v(2)v2, v(3)v3 hold the current values of the Pruefer variables ββ, φϕ, ρρ respectively.
3:     jint – int64int32nag_int scalar
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(2)x=xpoint2) jint is set to 00 and at the extreme right-hand end point (when x = xr = xpoint(m1)x=xr=xpointm-1) jint is set to m2m-2.

Optional Input Parameters

1:     m – int64int32nag_int scalar
Default: The dimension of the arrays xpoint, hmax. (An error is raised if these dimensions are not equal.)
The number of points in the array xpoint.
Constraint: m4m4.
2:     maxit – int64int32nag_int scalar
A bound on nrnr, the number of root-finding iterations allowed, that is the number of trial values of λλ that are used. If maxit0maxit0, no such bound is assumed. (See also maxfun.)
Default: 00
3:     maxfun – int64int32nag_int scalar
A bound on nfnf, the number of calls to coeffn made in any one root-finding iteration. If maxfun0maxfun0, no such bound is assumed.
Default: 00

Input Parameters Omitted from the MATLAB Interface

None.

Output Parameters

1:     match – int64int32nag_int scalar
The index of the break point actually used as the matching point.
2:     elam – double scalar
The final computed estimate, whether or not an error occurred.
3:     delam – double scalar
If ifail = 0ifail=0, delam holds an estimate of the absolute error in the computed eigenvalue, that is |λ̃elam|delam|λ~-elam|delam. (In Section [General Description of the Algorithm] 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 ifail0ifail0, delam may hold an estimate of the error, or its initial value, depending on the value of ifail. See Section [Error Indicators and Warnings] for further details.
4:     hmax(22,m) – double array
hmax(1,m1)hmax1m-1 and hmax(1,m)hmax1m contain the sensitivity coefficients σl,σrσl,σr, described in Section [The Sensitivity s and ]. Other entries contain diagnostic output in the case of an error exit (see Section [Error Indicators and Warnings]).
5:     maxit – int64int32nag_int scalar
Default: 00
Will have been decreased by the number of iterations actually performed, whether or not it was positive on entry.
6:     ifail – int64int32nag_int scalar
ifail = 0ifail=0 unless the function detects an error (see [Error Indicators and Warnings]).

Error Indicators and Warnings

Errors or warnings detected by the function:

Cases prefixed with W are classified as warnings and do not generate an error of type NAG:error_n. See nag_issue_warnings.

  ifail = 1ifail=1
A parameter error. All parameters (except ifail) are left unchanged. The reason for the error is shown by the value of hmax(2,1)hmax21 as follows:
hmax(2,1) = 1hmax21=1: m < 4m<4;
hmax(2,1) = 2hmax21=2: k < 0k<0;
hmax(2,1) = 3hmax21=3: tol0.0tol0.0;
hmax(2,1) = 4hmax21=4: xpoint(1)xpoint1 to xpoint(m)xpointm are not in ascending order. hmax(2,2)hmax22 gives the position ii in xpoint where this was detected.
  ifail = 2ifail=2
At some call to bdyval, invalid values were returned, that is, either yl(1) = yl(2) = 0.0yl1=yl2=0.0, or yr(1) = yr(2) = 0.0yr1=yr2=0.0 (a programming error in bdyval). See the last call of monit for details.
This error exit will also occur if p(x)p(x) is zero at the point where the boundary condition is imposed. Probably bdyval was called with xl equal to a singular end point aa or with xr equal to a singular end point bb.
  ifail = 3ifail=3
At some point between xl and xr the value of p(x)p(x) computed by coeffn became zero or changed sign. See the last call of monit for details.
  ifail = 4ifail=4
maxit > 0maxit>0 on entry, and after maxit iterations the eigenvalue had not been found to the required accuracy.
  ifail = 5ifail=5
The ‘bracketing’ phase (with parameter iflag of the monit equal to 11) failed to bracket the eigenvalue within ten iterations. This is caused by an error in formulating the problem (for example, qq is independent of λλ), or by very poor initial estimates of elam and delam.
On exit, elam and elam + delamelam+delam give the end points of the interval within which no eigenvalue was located by the function.
  ifail = 6ifail=6
maxfun > 0maxfun>0 on entry, and the last iteration was terminated because more than maxfun calls to coeffn were used. See the last call of monit for details.
  ifail = 7ifail=7
To obtain the desired accuracy the local error tolerance was set so small at the start of some sub-interval 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.
  ifail = 8ifail=8
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 = 7ifail=7). This could be due to pathological behaviour of p(x)p(x) and q(x ; λ)q(x;λ) or to an unreasonable accuracy requirement or to the current value of λλ making the equations ‘stiff’. See the last call of monit for details.
W ifail = 9ifail=9
tol is too small for the problem being solved and the machine precision is being used. The final value of elam should be a very good approximation to the eigenvalue.
  ifail = 10ifail=10
nag_roots_contfn_brent_rcomm (c05az), called by nag_ode_sl2_breaks_funs (d02ke), has terminated with the error exit corresponding to a pole of the residual function f(λ)f(λ). This error exit should not occur, but if it does, try solving the problem again with a smaller value for tol.
  ifail = 11ifail=11
  ifail = 12ifail=12
A serious error has occurred in an internal call. Check all (sub)program calls and array dimensions. Seek expert help.
Note: error exits with ifail = 2ifail=2, 33, 66, 77, 88 or 1111 are caused by being unable to set up or solve the differential equation at some iteration and will be immediately preceded by a call of monit giving diagnostic information. For other errors, diagnostic information is contained in hmax(2,j)hmax2j, for j = 1,2,,mj=1,2,,m, where appropriate.

Accuracy

See the discussion in Section [General Description of the Algorithm].

Further Comments

Timing

The time taken by nag_ode_sl2_breaks_funs (d02ke) 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 1616. To make the most economical use of the function, one should try to obtain good initial values for elam and delam, and, where appropriate, good asymptotic formulae. Also the boundary matching points should not be set unnecessarily close to singular points. The extra time needed to compute the eigenfunction is principally the cost of one additional integration once the eigenvalue has been found.

General Description of the Algorithm

A shooting method, for differential equation problems containing unknown parameters, relies on the construction of a ‘miss-distance function’, which for given trial values of the parameters measures how far the conditions of the problem are from being met. The problem is then reduced to one of finding the values of the parameters for which the miss-distance function is zero, that is to a root-finding process. Shooting methods differ mainly in how the miss-distance is defined.
nag_ode_sl2_breaks_funs (d02ke) defines a miss-distance f(λ)f(λ) based on the rotation about the origin of the point p(x) = (p(x)y(x),y(x))px=(p(x)y(x),y(x)) in the Phase Plane as the solution proceeds from aa to bb. The boundary conditions define the ray (i.e., two-sided line through the origin) on which p(x)p(x) should start, and the ray on which it should finish. The eigenvalue kk defines the total number of half-turns it should make. Numerical solution is actually done by ‘shooting forward’ from x = ax=a and ‘shooting backward’ from x = bx=b to a matching point x = cx=c. Then f(λ)f(λ) is taken as the angle between the rays to the two resulting points Pa(c)Pa(c) and Pb(c)Pb(c). A relative scaling of the pypy and yy axes, based on the behaviour of the coefficient functions pp and qq, is used to improve the numerical behaviour.
Figure 1
Figure 1
The resulting function f(λ)f(λ) is monotonic over < λ < -<λ<, increasing if (q)/(λ) > 0 q λ >0 and decreasing if (q)/(λ) < 0 q λ <0, with a unique zero at the desired eigenvalue λ̃λ~. The function measures f(λ)f(λ) in units of a half-turn. This means that as λλ increases, f(λ)f(λ) varies by about 11 as each eigenvalue is passed. (This feature implies that the values of f(λ)f(λ) at successive iterations – especially in the early stages of the iterative process – can be used with suitable extrapolation or interpolation to help the choice of initial estimates for eigenvalues near to the one currently being found.)
The function actually computes a value for f(λ)f(λ) with errors, arising from the local errors of the differential equation code and from the asymptotic formulae provided by you if singular points are involved. However, the error estimate output in delam is usually fairly realistic, in that the actual error |λ̃elam||λ~-elam| is within an order of magnitude of delam.
We pass the values of ββ, φϕ, ρρ across through report rather than converting them to values of yy, yy inside nag_ode_sl2_breaks_funs (d02ke), for the following reasons. First, there may be cases where auxiliary quantities can be more accurately computed from the Pruefer variables than from yy and yy. Second, in singular problems on an infinite interval yy and yy may underflow towards the end of the range, whereas the Pruefer variables remain well-behaved. Third, with high-order eigenvalues (and therefore highly oscillatory eigenfunctions) the eigenfunction may have a complete oscillation (or more than one oscillation) between two mesh points, so that values of yy and yy at mesh points give a very poor representation of the curve. The probable behaviour of the Pruefer variables in this case is that ββ and ρρ vary slowly whilst φϕ increases quickly: for all three Pruefer variables linear interpolation between the values at adjacent mesh points is probably sufficiently accurate to yield acceptable intermediate values of ββ, φϕ, ρρ (and hence of y,yy,y) for graphical purposes.
Similar considerations apply to the exponentially decaying ‘tails’ of the eigenfunctions that often occur in singular problems. Here φϕ has approximately constant value whilst ρρ increases rapidly in the direction of integration, though the step length is generally fairly small over such a range.
If the solution is output through report at xx values which are too widely spaced, the step length can be controlled by choosing hmax suitably, or, preferably, by reducing tol. Both these choices will lead to more accurate eigenvalues and eigenfunctions but at some computational cost.

The Position of the Shooting Matching Point c

This point is always one of the values xixi in array xpoint. It may be specified using the parameter match. The default value is chosen to be the value of that xixi, 2im12im-1, that lies closest to the middle of the interval [x2,xm1][x2,xm-1]. If there is a tie, the rightmost candidate is chosen. In particular if there are no break points, then c = xm1c=xm-1 ( = x3=x3); that is, the shooting is from left to right in this case. A break point may be inserted purely to move cc to an interior point of the interval, even though the form of the equations does not require it. This often speeds up convergence especially with singular problems.
Note that the shooting method used by the code integrates first from the left-hand end xlxl, then from the right-hand end xrxr, to meet at the matching point cc in the middle. This will of course be reflected in printed or graphical output. The diagram shows a possible sequence of nine mesh points τ1τ1 through τ9τ9 in the order in which they appear, assuming there are just two sub-intervals (so m = 5m=5).
Figure 2
Figure 2
Since the shooting method usually fails to match up the two ‘legs’ of the curve exactly, there is bound to be a jump in yy, or in p(x)yp(x)y or both, at the matching point cc. The code in fact ‘shares’ the discrepancy out so that both yy and p(x)yp(x)y have a jump. A large jump does not imply an inaccurate eigenvalue, but implies either
(a) a badly chosen matching point: if q(x ; λ)q(x;λ) has a ‘humped’ shape, cc should be chosen near the maximum value of qq, especially if qq is negative at the ends of the interval;
(b) an inherently ill-conditioned problem, typically one where another eigenvalue is pathologically close to the one being sought. In this case it is extremely difficult to obtain an accurate eigenfunction.
In Section [Example], we find the 11th eigenvalue and corresponding eigenfunction of the equation
y + (λx2 / x2)y = 0  on  0 < x < ,
y+(λ-x-2/x2)y=0  on  0<x<,
the boundary conditions being that yy should remain bounded as xx tends to 00 and xx tends to . The coding of this problem is discussed in detail in Section [Examples of Boundary Conditions at Singular Points].
The choice of matching point cc is open. If we choose c = 30.0c=30.0 as in nag_ode_sl2_breaks_vals (d02kd) example program we find that the exponentially increasing component of the solution dominates and we get extremely inaccurate values for the eigenfunction (though the eigenvalue is determined accurately). The values of the eigenfunction calculated with c = 30.0c=30.0 are given schematically in Figure 3.
Figure 3
Figure 3
If we choose cc as the maximum of the hump in q(x ; λ)q(x;λ) (see item (a) above) we instead obtain the accurate results given in Figure 4
Figure 4
Figure 4

Examples of Coding the coeffn

Coding coeffn is straightforward except when break points are needed. The examples below show:
(a) a simple case,
(b) a case in which discontinuities in the coefficient functions or their derivatives necessitate break points, and
(c) a case where break points together with the hmax parameter are an efficient way to deal with a coefficient function that is well-behaved except over one short interval.
(Some of these cases are among the examples in Section [Example].)
Example A
The modified Bessel equation
x (xy) + (λx2ν2) y = 0 .
x (xy) + (λx2-ν2) y=0 .
Assuming the interval of solution does not contain the origin and dividing through by xx, we have p(x) = xp(x)=x and q(x ; λ) = λxν2 / xq(x;λ)=λx-ν2/x. The code could be
function [p, q, dqdl] = coeffn(x, elam, jint)
  global nu;
  ...
  p = x;
  q = elam*x - nu*nu/x
  dqdl = x;
where nu (standing for νν) is a double global variable declared in the calling program.
Example B
The Schroedinger equation
y + (λ + q(x))y = 0,
y+(λ+q(x))y=0,
where
q(x) =
{ x2 − 10 (|x| ≤ 4), 6/(|x|) (|x| > 4),
q(x) = { x2- 10 (|x| 4), 6|x| (|x|> 4),
over some interval ‘approximating to (,)(-,)’, say [20,20][-20,20]. Here we need break points at ± 4±4, forming three sub-intervals i1 = [20,4]i1=[-20,-4], i2 = [4,4]i2=[-4,4], i3 = [4,20]i3=[4,20]. The code could be
function [p, q, dqdl] = coeffn(x, elam, jint)
  ...
  p = 1;
  dqdl = 1;
  if (jint == 2)
    q = elam + x*x - 10
  else
    q = elam + 6/abs(x)
  end
The array xpoint would contain the values x1x1, 20.0-20.0, 4.0-4.0, + 4.0+4.0, + 20.0+20.0, x6x6 and mm would be 66. The choice of appropriate values for x1x1 and x6x6 depends on the form of the asymptotic formula computed by bdyval and the technique is discussed in Section [Examples of Boundary Conditions at Singular Points].
Example C
y + λ(12e100x2)y = 0,  10x10.
y+λ(1-2e-100x2)y=0,  -10x10.
Here q(x ; λ)q(x;λ) is nearly constant over the range except for a sharp inverted spike over approximately 0.1x0.1-0.1x0.1. There is a danger that the function will build up to a large step size and ‘step over’ the spike without noticing it. By using break points – say at ± 0.5±0.5 – one can restrict the step size near the spike without impairing the efficiency elsewhere.
The code for coeffn could be
  function [p, q, dqdl] = coeffn(x, elam, jint)
  ...
  p = 1;
  dqdl = 1 - 2*exp(-100*x*x);
  q = elam * dqdl;
xpoint might contain 10.0-10.0, 10.0-10.0, 0.5-0.5, 0.50.5, 10.010.0, 10.010.0 (assuming ± 10±10 are regular points) and mm would be 66. hmax(1,j)hmax1j, for j = 1,2,3j=1,2,3, might contain 0.00.0, 0.10.1 and 0.00.0.

Examples of Boundary Conditions at Singular Points

Quoting from page 243 of Bailey (1966): ‘Usually ... the differential equation has two essentially different types of solution near a singular point, and the boundary condition there merely serves to distinguish one kind from the other. This is the case in all the standard examples of mathematical physics.’
In most cases the behaviour of the ratio p(x)y / yp(x)y/y near the point is quite different for the two types of solution. Essentially what you provide through the bdyval is an approximation to this ratio, valid as xx tends to the singular point (SP).
You must decide (a) how accurate to make this approximation or asymptotic formula, for example how many terms of a series to use, and (b) where to place the boundary matching point (BMP) at which the numerical solution of the differential equation takes over from the asymptotic formula. Taking the BMP closer to the SP will generally improve the accuracy of the asymptotic formula, but will make the computation more expensive as the Pruefer differential equations generally become progressively more ill-behaved as the SP is approached. You are strongly recommended to experiment with placing the BMPs. In many singular problems quite crude asymptotic formulae will do. To help you avoid needlessly accurate formulae, nag_ode_sl2_breaks_funs (d02ke) outputs two ‘sensitivity coefficients’ σl,σrσl,σr which estimate how much the errors at the BMPs affect the computed eigenvalue. They are described in detail in Section [The Sensitivity s and ].
Example of coding bdyval:
The example below illustrates typical situations:
y + (λx2/(x2)) y = 0,  0 < x <
y+(λ-x-2x2) y=0,  0<x<
the boundary conditions being that yy should remain bounded as xx tends to 00 and xx tends to .
At the end x = 0x=0 there is one solution that behaves like x2x2 and another that behaves like x1x-1. For the first of these solutions p(x)y / yp(x)y/y is asymptotically 2 / x2/x while for the second it is asymptotically 1 / x-1/x. Thus the desired ratio is specified by setting
yl(1) = x  and  yl(2) = 2.0.
yl1=x  and  yl2=2.0.
At the end x = x= the equation behaves like Airy's equation shifted through λλ, i.e., like yty = 0y-ty=0 where t = xλt=x-λ, so again there are two types of solution. The solution we require behaves as
exp((2/3)t(3/2)) / root(4,t)
exp(-23t32)/t4
and the other as
exp( + (2/3)t(3/2)) / root(4,t).
exp(+23t32)/t4.
Hence, the desired solution has p(x)y / ysqrt(t)p(x)y/y-t so that we could set yr(1) = 1.0yr1=1.0 and yr(2) = sqrt(xλ)yr2=-x-λ. The complete function might thus be
function [yl, yr] = bdyval(xl, xr, elam)
  yl(1) = xl;
  yl(2) = 2;
  yr(1) = 1;
  yr(2) = -sqrt(xr-elam);
Clearly for this problem it is essential that any value given by nag_ode_sl2_breaks_funs (d02ke) to xr is well to the right of the value of elam, so that you must vary the right-hand BMP with the eigenvalue index kk. One would expect λkλk to be near the kkth zero of the Airy function Ai(x)Ai(x), so there is no problem estimating elam.
More accurate asymptotic formulae are easily found: near x = 0x=0 by the standard Frobenius method, and near x = x= by using standard asymptotics for Ai(x)Ai(x), Ai(x)Ai(x), (see page 448 of Abramowitz and Stegun (1972)).
For example, by the Frobenius method the solution near x = 0x=0 has the expansion
y = x2(c0 + c1x + c2x2 + )
y=x2(c0+c1x+c2x2+)
with
c0 = 1,c1 = 0, c2 = (λ)/10, c3 = (1/18),, cn = (cn 3λ cn 2)/(n(n + 3)).
c0= 1,c1= 0, c2=-λ10, c3=118,, cn=cn- 3-λ cn- 2 n(n+ 3) .
This yields
(p(x)y)/y = (2(2/5)λx2 + )/(x (1λ/10x2 + ) ).
p(x)yy=2-25λx2+ x (1-λ10x2+) .

The Sensitivity Parameters σl and σr

The sensitivity parameters σlσl, σrσr (held in hmax(1,m1)hmax1m-1 and hmax(1,m)hmax1m on output) estimate the effect of errors in the boundary conditions. For sufficiently small errors ΔyΔy, ΔpyΔpy in yy and pypy respectively, the relations
Δλ(y . Δpypy . Δy)lσl
Δλ(y . Δpypy . Δy)rσr
Δλ(y.Δpy-py.Δy)lσl Δλ(y.Δpy-py.Δy)rσr
are satisfied, where the subscripts ll, rr denote errors committed at the left- and right-hand BMPs respectively, and ΔλΔλ denotes the consequent error in the computed eigenvalue.

‘Missed Zeros’

This is a pitfall to beware of at a singular point. If the BMP is chosen so far from the SP that a zero of the desired eigenfunction lies in between them, then the function will fail to ‘notice’ this zero. Since the index of kk of an eigenvalue is the number of zeros of its eigenfunction, the result will be that
(a) the wrong eigenvalue will be computed for the given index kk – in fact some λk + kλk+k will be found where k1k1;
(b) the same index kk can cause convergence to any of several eigenvalues depending on the initial values of elam and delam.
It is up to you to take suitable precautions – for instance by varying the position of the BMPs in the light of knowledge of the asymptotic behaviour of the eigenfunction at different eigenvalues.

Example

This example finds the 11th eigenvalue and eigenfunction of the example of Section [Examples of Boundary Conditions at Singular Points], using the simple asymptotic formulae for the boundary conditions.
Comparison of the results from this example program with the corresponding results from nag_ode_sl2_breaks_vals (d02kd) example program shows that similar output is produced from monit, followed by the eigenfunction values from report, and then a further line of information from monit (corresponding to the integration to find the eigenfunction). Final information is printed within the example program exactly as with nag_ode_sl2_breaks_vals (d02kd).
Note the discrepancy at the matching point cc ( = root(3,4)=43, the maximum of q(x ; λ)q(x;λ), in this case) between the solutions obtained by integrations from left- and right-hand end points.
function nag_ode_sl2_breaks_funs_example
% For communication with display_plot.
global ykeep ncall xkeep pkeep;

% Set up initial values.
xpoint = [0; 0.1; 4d0^(1d0/3d0); 30; 30];
match = 0;
k = 11;
tol = 1d-4;
elam = 14;
delam = 1;
m = 5;
hmax = zeros(2,m);

ncall = 0;
ykeep = zeros(1,1);
xkeep = zeros(1,1);
pkeep = zeros(1,1);

disp('nag_ode_sl2_breaks_funs example program results');

% report outputs intermediate results.
[matchOut, elamOut, delamOut, hmaxOut, maxit, ifail] = ...
    nag_ode_sl2_breaks_funs(xpoint, int64(match), @coeffn, @bdyval, ...
    int64(k), tol, elam, delam, hmax, 'nag_ode_sl2_reg_finite_dummy_monit', @report);
if ifail ~= 0
    % Problems in integration.  Print message and exit.
    error('Warning: nag_ode_sl2_breaks_funs returned with ifail = %1d ',ifail);
end

% Output final results.
fprintf('\nFinal results\n');
fprintf('\nk = %3.0f   elam = %3.3f   delam = %3.2e \n',k,elamOut,delamOut);
fprintf('hmax(1,m-1) = %3.3f    hmax(1,m) = %3.3f\n\n',hmaxOut(1,m-1),hmaxOut(1,m));
fig = figure('Number', 'off');
display_plot(xkeep,pkeep,ykeep)

function [yl, yr] = bdyval(xl, xr, elam)
% Define the boundary conditions.
yl(1) = xl;
yl(2) = 2;
yr(1) = 1;
yr(2) = -sqrt(xr-elam);
function [p, q, dqdl] = coeffn(x, elam, jint)
% Compute the coefficient functions.
p = 1;
dqdl = 1;
q = elam - x - 2/x/x;
function report(x, v, jint)
% For communication with main routine.
global ykeep ncall xkeep pkeep;

% Compute eigenfunction and its dervative at this mesh point.
if (jint == 0)
    fprintf('\n A singular problem\n');
    fprintf('\n Eigenfunction Values\n');
    fprintf('       x            y          pyp\n');
end
sqrtb = sqrt(v(1));
if (0.5*v(3) >= log(nag_machine_real_safe))
    r = exp(0.5*v(3));
else
    r = 0d0;
end
pyp = r*sqrtb*cos(0.5*v(2));
y = r/sqrtb*sin(0.5*v(2));

% Output results and store them for plotting.
fprintf('%10.3f %12.4f %12.4f\n', x, y, pyp);
ncall = ncall+1;
ykeep(ncall,1) = y;
xkeep(ncall,1) = x;
pkeep(ncall,1) = pyp;
function display_plot(xkeep, pkeep, ykeep)
% Formatting for title and axis labels.
titleFmt = {'FontName', 'Helvetica', 'FontWeight', 'Bold', 'FontSize', 14};
labFmt = {'FontName', 'Helvetica', 'FontWeight', 'Bold', 'FontSize', 13};
set(gca, 'FontSize', 13); % for legend, axis tick labels, etc.

% Re-order the arrays before plotting.
n = length(xkeep);
for i = 2:n
    if (xkeep(i-1) > xkeep(i))
        x1 = [xkeep(i-1:n);rot90(xkeep(1:i-2)')];
        p1 = [pkeep(i-1:n);rot90(pkeep(1:i-2)')];
        y1 = [ykeep(i-1:n);rot90(ykeep(1:i-2)')];
        break
    end
end

% Plot the two curves.
plot(x1,y1,'-',x1,p1,'--');
% Add title.
title({'Regular Singular Second-order Sturm-Liouville System',...
    '11th Eigenfunction and Derivative'}, titleFmt{:});
% Label the axes.
xlabel('x', labFmt{:});
ylabel('Eigenfunction and Derivative', labFmt{:});
% Add a legend.
legend('y','p(x)y''','Location','Best');
 
nag_ode_sl2_breaks_funs example program results

 A singular problem

 Eigenfunction Values
       x            y          pyp
     0.100       0.0123       0.2466
     0.168       0.0344       0.3917
     0.216       0.0552       0.4749
     0.312       0.1065       0.5811
     0.408       0.1634       0.5917
     0.578       0.2496       0.3706
     0.725       0.2779      -0.0053
     0.910       0.2268      -0.5395
     1.139       0.0490      -0.9366
     1.454      -0.2160      -0.5684
     1.587      -0.2654      -0.1571
    30.000      -0.0000       0.0000
    29.096      -0.0000       0.0000
    28.629      -0.0000       0.0000
    28.356      -0.0000       0.0000
    28.062      -0.0000       0.0000
    27.713      -0.0000       0.0000
    27.262      -0.0000       0.0000
    26.855      -0.0000       0.0000
    26.432      -0.0000       0.0000
    26.062      -0.0000       0.0000
    25.686      -0.0000       0.0000
    25.301      -0.0000       0.0000
    24.891      -0.0000       0.0000
    24.574      -0.0000       0.0000
    24.249      -0.0000       0.0000
    23.855      -0.0000       0.0000
    23.530      -0.0000       0.0000
    23.157      -0.0000       0.0000
    22.843      -0.0000       0.0000
    22.467      -0.0000       0.0000
    22.140      -0.0000       0.0000
    21.743      -0.0000       0.0000
    21.397      -0.0000       0.0000
    20.979      -0.0000       0.0000
    20.614      -0.0000       0.0000
    20.173      -0.0001       0.0001
    19.786      -0.0001       0.0003
    19.453      -0.0003       0.0006
    19.016      -0.0007       0.0015
    18.601      -0.0017       0.0033
    18.224      -0.0035       0.0066
    17.865      -0.0068       0.0121
    17.503      -0.0127       0.0214
    17.124      -0.0235       0.0371
    16.716      -0.0437       0.0633
    16.248      -0.0829       0.1070
    15.732      -0.1536       0.1681
    15.411      -0.2137       0.2048
    15.079      -0.2863       0.2295
    14.785      -0.3542       0.2285
    14.484      -0.4187       0.1927
    14.237      -0.4592       0.1298
    13.895      -0.4812      -0.0113
    13.519      -0.4378      -0.2267
    13.124      -0.3002      -0.4651
    12.559       0.0249      -0.6286
    12.070       0.2940      -0.4073
    11.605       0.3723       0.1023
    11.133       0.1954       0.6094
    10.652      -0.1371       0.6640
    10.199      -0.3360       0.1374
     9.773      -0.2482      -0.5273
     9.325       0.0680      -0.7601
     8.860       0.3138      -0.1743
     8.441       0.2184       0.5955
     8.005      -0.1129       0.7623
     7.569      -0.3076       0.0199
     7.157      -0.1390      -0.7565
     6.721       0.2052      -0.6183
     6.313       0.2764       0.3190
     5.891      -0.0132       0.8766
     5.480      -0.2783       0.2358
     5.072      -0.1582      -0.7515
     4.649       0.1938      -0.6565
     4.248       0.2499       0.4235
     3.837      -0.0698       0.8915
     3.418      -0.2751      -0.0829
     3.016      -0.0334      -0.9312
     2.607       0.2612      -0.2590
     2.173       0.0886       0.8990
     1.722      -0.2560       0.2925
     1.587      -0.2653      -0.1566

Final results

k =  11   elam = 14.946   delam = 9.60e-04 
hmax(1,m-1) = -0.015    hmax(1,m) = 0.000


function d02ke_example
% For communication with display_plot.
global ykeep ncall xkeep pkeep;

% Set up initial values.
xpoint = [0; 0.1; 4d0^(1d0/3d0); 30; 30];
match = 0;
k = 11;
tol = 1d-4;
elam = 14;
delam = 1;
m = 5;
hmax = zeros(2,m);

ncall = 0;
ykeep = zeros(1,1);
xkeep = zeros(1,1);
pkeep = zeros(1,1);

disp('d02ke example program results');

% report outputs intermediate results.
[matchOut, elamOut, delamOut, hmaxOut, maxit, ifail] = ...
    d02ke(xpoint, int64(match), @coeffn, @bdyval, ...
    int64(k), tol, elam, delam, hmax, 'd02kay', @report);
if ifail ~= 0
    % Problems in integration.  Print message and exit.
    error('Warning: d02ke returned with ifail = %1d ',ifail);
end

% Output final results.
fprintf('\nFinal results\n');
fprintf('\nk = %3.0f   elam = %3.3f   delam = %3.2e \n',k,elamOut,delamOut);
fprintf('hmax(1,m-1) = %3.3f    hmax(1,m) = %3.3f\n\n',hmaxOut(1,m-1),hmaxOut(1,m));
fig = figure('Number', 'off');
display_plot(xkeep,pkeep,ykeep)

function [yl, yr] = bdyval(xl, xr, elam)
% Define the boundary conditions.
yl(1) = xl;
yl(2) = 2;
yr(1) = 1;
yr(2) = -sqrt(xr-elam);
function [p, q, dqdl] = coeffn(x, elam, jint)
% Compute the coefficient functions.
p = 1;
dqdl = 1;
q = elam - x - 2/x/x;
function report(x, v, jint)
% For communication with main routine.
global ykeep ncall xkeep pkeep;

% Compute eigenfunction and its dervative at this mesh point.
if (jint == 0)
    fprintf('\n A singular problem\n');
    fprintf('\n Eigenfunction Values\n');
    fprintf('       x            y          pyp\n');
end
sqrtb = sqrt(v(1));
if (0.5*v(3) >= log(x02am))
    r = exp(0.5*v(3));
else
    r = 0d0;
end
pyp = r*sqrtb*cos(0.5*v(2));
y = r/sqrtb*sin(0.5*v(2));

% Output results and store them for plotting.
fprintf('%10.3f %12.4f %12.4f\n', x, y, pyp);
ncall = ncall+1;
ykeep(ncall,1) = y;
xkeep(ncall,1) = x;
pkeep(ncall,1) = pyp;
function display_plot(xkeep, pkeep, ykeep)
% Formatting for title and axis labels.
titleFmt = {'FontName', 'Helvetica', 'FontWeight', 'Bold', 'FontSize', 14};
labFmt = {'FontName', 'Helvetica', 'FontWeight', 'Bold', 'FontSize', 13};
set(gca, 'FontSize', 13); % for legend, axis tick labels, etc.

% Re-order the arrays before plotting.
n = length(xkeep);
for i = 2:n
    if (xkeep(i-1) > xkeep(i))
        x1 = [xkeep(i-1:n);rot90(xkeep(1:i-2)')];
        p1 = [pkeep(i-1:n);rot90(pkeep(1:i-2)')];
        y1 = [ykeep(i-1:n);rot90(ykeep(1:i-2)')];
        break
    end
end

% Plot the two curves.
plot(x1,y1,'-',x1,p1,'--');
% Add title.
title({'Regular Singular Second-order Sturm-Liouville System',...
    '11th Eigenfunction and Derivative'}, titleFmt{:});
% Label the axes.
xlabel('x', labFmt{:});
ylabel('Eigenfunction and Derivative', labFmt{:});
% Add a legend.
legend('y','p(x)y''','Location','Best');
 
d02ke example program results

 A singular problem

 Eigenfunction Values
       x            y          pyp
     0.100       0.0123       0.2466
     0.168       0.0344       0.3917
     0.216       0.0552       0.4749
     0.312       0.1065       0.5811
     0.408       0.1634       0.5917
     0.578       0.2496       0.3706
     0.725       0.2779      -0.0053
     0.910       0.2268      -0.5395
     1.139       0.0490      -0.9366
     1.454      -0.2160      -0.5684
     1.587      -0.2654      -0.1571
    30.000      -0.0000       0.0000
    29.096      -0.0000       0.0000
    28.629      -0.0000       0.0000
    28.356      -0.0000       0.0000
    28.062      -0.0000       0.0000
    27.713      -0.0000       0.0000
    27.262      -0.0000       0.0000
    26.855      -0.0000       0.0000
    26.432      -0.0000       0.0000
    26.062      -0.0000       0.0000
    25.686      -0.0000       0.0000
    25.301      -0.0000       0.0000
    24.891      -0.0000       0.0000
    24.574      -0.0000       0.0000
    24.249      -0.0000       0.0000
    23.855      -0.0000       0.0000
    23.530      -0.0000       0.0000
    23.157      -0.0000       0.0000
    22.843      -0.0000       0.0000
    22.467      -0.0000       0.0000
    22.140      -0.0000       0.0000
    21.743      -0.0000       0.0000
    21.397      -0.0000       0.0000
    20.979      -0.0000       0.0000
    20.614      -0.0000       0.0000
    20.173      -0.0001       0.0001
    19.786      -0.0001       0.0003
    19.453      -0.0003       0.0006
    19.016      -0.0007       0.0015
    18.601      -0.0017       0.0033
    18.224      -0.0035       0.0066
    17.865      -0.0068       0.0121
    17.503      -0.0127       0.0214
    17.124      -0.0235       0.0371
    16.716      -0.0437       0.0633
    16.248      -0.0829       0.1070
    15.732      -0.1536       0.1681
    15.411      -0.2137       0.2048
    15.079      -0.2863       0.2295
    14.785      -0.3542       0.2285
    14.484      -0.4187       0.1927
    14.237      -0.4592       0.1298
    13.895      -0.4812      -0.0113
    13.519      -0.4378      -0.2267
    13.124      -0.3002      -0.4651
    12.559       0.0249      -0.6286
    12.070       0.2940      -0.4073
    11.605       0.3723       0.1023
    11.133       0.1954       0.6094
    10.652      -0.1371       0.6640
    10.199      -0.3360       0.1374
     9.773      -0.2482      -0.5273
     9.325       0.0680      -0.7601
     8.860       0.3138      -0.1743
     8.441       0.2184       0.5955
     8.005      -0.1129       0.7623
     7.569      -0.3076       0.0199
     7.157      -0.1390      -0.7565
     6.721       0.2052      -0.6183
     6.313       0.2764       0.3190
     5.891      -0.0132       0.8766
     5.480      -0.2783       0.2358
     5.072      -0.1582      -0.7515
     4.649       0.1938      -0.6565
     4.248       0.2499       0.4235
     3.837      -0.0698       0.8915
     3.418      -0.2751      -0.0829
     3.016      -0.0334      -0.9312
     2.607       0.2612      -0.2590
     2.173       0.0886       0.8990
     1.722      -0.2560       0.2925
     1.587      -0.2653      -0.1566

Final results

k =  11   elam = 14.946   delam = 9.60e-04 
hmax(1,m-1) = -0.015    hmax(1,m) = 0.000



PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

© The Numerical Algorithms Group Ltd, Oxford, UK. 2009–2013