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_reg_finite (d02ka)

Purpose

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

Syntax

[bcond, elam, delam, ifail] = d02ka(xl, xr, coeffn, bcond, k, tol, elam, delam, monit)
[bcond, elam, delam, ifail] = nag_ode_sl2_reg_finite(xl, xr, coeffn, bcond, k, tol, elam, delam, monit)

Description

nag_ode_sl2_reg_finite (d02ka) 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 boundary conditions of the form
a2y(a) = a1p(a)y(a)
a2y(a)=a1p(a)y(a)
b2y(b) = b1p(b)y(b)
b2y(b)=b1p(b)y(b)
at the two, finite, end points aa and bb. The functions pp and qq, which are real-valued, are defined by coeffn.
For the theoretical basis of the numerical method to be valid, the following conditions should hold on the coefficient functions:
(a) p(x)p(x) must be nonzero and must not change sign throughout the closed interval [a,b][a,b];
(b) (q)/(λ) q λ  must not change sign and must be nonzero throughout the open interval (a,b)(a,b) and for all relevant values of λλ, and must not be identically zero as xx varies, for any relevant value λλ; and,
(c) pp and qq should (as functions of xx) have continuous derivatives, preferably up to the fourth-order, on [a,b][a,b]. The differential equation code used will integrate through mild discontinuities, but probably with severely reduced efficiency. Therefore, if pp and qq violate this condition, nag_ode_sl2_breaks_vals (d02kd) should be used.
The eigenvalue λ̃λ~ is determined by a shooting method based on a Pruefer transformation of the differential equations. Providing certain assumptions are met, the computed value of λ̃λ~ will be correct to within a mixed absolute/relative error specified by tol. nag_ode_sl2_reg_finite (d02ka) is a driver function for the more complicated function nag_ode_sl2_breaks_vals (d02kd) whose specification provides more details of the techniques used.
A good account of the theory of Sturm–Liouville systems, with some description of Pruefer transformations, is given in Chapter X of Birkhoff and Rota (1962). An introduction to the use of Pruefer transformations for the numerical solution of eigenvalue problems arising from physics and chemistry is given in Bailey (1966).

References

Bailey P B (1966) Sturm–Liouville eigenvalues via a phase function SIAM J. Appl. Math. 14 242–249
Birkhoff G and Rota G C (1962) Ordinary Differential Equations Ginn & Co., Boston and New York

Parameters

Compulsory Input Parameters

1:     xl – double scalar
2:     xr – double scalar
The left- and right-hand end points aa and bb respectively, of the interval of definition of the problem.
Constraint: xl < xrxl<xr.
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.
[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
This parameter is included for compatibility with the more complex function nag_ode_sl2_breaks_vals (d02kd) (which is called by nag_ode_sl2_reg_finite (d02ka)).
Need not be set.

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:     bcond(33,22) – double array
bcond(1,1)bcond11 and bcond(2,1)bcond21 must contain the numbers a1a1, a2a2 specifying the left-hand boundary condition in the form
a2 y(a) = a1 p(a) y(a)
a2 y(a) = a1 p(a) y(a)
where |a2| + |a1p(a)|0|a2|+|a1p(a)|0.
bcond(1,2)bcond12 and bcond(2,2)bcond22 must contain b1b1, b2b2 such that
b2 y(b) = b1 p(b) y(b)
b2 y(b) = b1 p(b) y(b)
where |b2| + |b1p(b)|0|b2|+|b1p(b)|0.
Note the occurrence of p(a)p(a), p(b)p(b) in these formulae.
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 about half the distance between adjacent eigenvalues in the neighbourhood of the one sought. In practice, there will often be a problem, similar to the one in hand but with known eigenvalues, which will help one to choose initial values for elam and delam.
If delam = 0.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:     monit – function handle or string containing name of m-file
monit is called by nag_ode_sl2_reg_finite (d02ka) at the end of each iteration for λλ 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
15 minus the number of iterations used so far in the search for λ̃λ~. (Up to 1515 iterations are permitted.)
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_reg_finite (d02ka) 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).
Normally, the iteration will terminate after a sequence of iterates with iflag = 2iflag=2, but occasionally the bracket on λ̃λ~ thus determined will not be sufficiently small and the iteration will be repeated with tighter accuracy control.
3:     elam – 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 iflag0iflag0), 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. f(λ̃) = 0f(λ~)=0 in theory. This 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
1 or 22 depending on whether integration was in a forward or backward direction at the time of failure.
finfo(10)finfo10
The value of the independent variable, xx, the point at which the error occurred.
finfo(11)finfo11, finfo(12)finfo12, finfo(13)finfo13
The current values of the Pruefer dependent variables ββ, φϕ and ρρ respectively. See Section [Description] in (d02ke) for a description of these variables.
finfo(14)finfo14
The local-error tolerance being used by the internal integrator at the point of failure.
finfo(15)finfo15
The last integration mesh point.

Optional Input Parameters

None.

Input Parameters Omitted from the MATLAB Interface

None.

Output Parameters

1:     bcond(33,22) – double array
bcond(3,1)bcond31 and bcond(3,2)bcond32 hold values σl,σrσl,σr estimating the sensitivity of the computed eigenvalue to changes in the boundary conditions. These values should only be of interest if the boundary conditions are, in some sense, an approximation to some ‘true’ boundary conditions. For example, if the range [xl, xr] should really be [0,][0,] but instead xr has been given a large value and the boundary conditions at infinity applied at xr, then the sensitivity parameter σrσr may be of interest. Refer to Section [Examples of Boundary Conditions at Singular Points] in (d02kd), for the actual meaning of σrσr and σlσl.
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, where λ̃λ~ is the true eigenvalue.
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:     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
On entry,k < 0k<0,
ortol0.0tol0.0.
  ifail = 2ifail=2
On entry,a1 = p(a)a2 = 0a1=p(a)a2=0,
orb1 = p(b)b2 = 0b1=p(b)b2=0,
(the array bcond has been set up incorrectly).
  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
After 1515 iterations the eigenvalue had not been found to the required accuracy.
  ifail = 5ifail=5
The ‘bracketing’ phase (with 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
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 = 7ifail=7
At some point 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 = 6ifail=6). 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 = 8ifail=8
tol is too small for the problem being solved and the machine precision being used. The local value of elam should be a very good approximation to the eigenvalue λ̃λ~.
  ifail = 9ifail=9
nag_roots_contfn_brent_rcomm (c05az), called by nag_ode_sl2_reg_finite (d02ka), has terminated with the error exit corresponding to a pole of the matching function. This error exit should not occur, but if it does, try solving the problem again with a smaller value for tol.
Note: error exits with ifail = 2ifail=2, 33, 66, 77 or 99 are caused by the inability to set up or solve the differential equation at some iteration and will be immediately preceded by a call of monit giving diagnostic information.
  ifail = 10ifail=10
  ifail = 11ifail=11
  ifail = 12ifail=12
A serious error has occurred in an internal call to an interpolation function. Check all (sub)program calls and array dimensions. Seek expert help.

Accuracy

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

Further Comments

The time taken by nag_ode_sl2_reg_finite (d02ka) 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.
See Section [Further Comments] in (d02kd) for a discussion of the technique used.

Example

function nag_ode_sl2_reg_finite_example
xl = 0;
xr = 3.141592653589793;
bcond = [1, 1;
     0, 0;
     0, -0.04564094560555453];
k = int64(4);
tol = 1e-05;
elam = 15;
delam = 4;
[bcondOut, elamOut, delamOut, ifail] = ...
    nag_ode_sl2_reg_finite(xl, xr, @coeffn, bcond, k, tol, elam, delam, @monit)

function [p, q, dqdl] = coeffn(x, elam, jint)
  p = 1;
  dqdl = 1;
  q = elam - 10*cos(2.0*x);
function monit(nit, iflag, elam, finfo)
  if (nit == 14)
    fprintf('\nOutput from Monit\n');
  end

  fprintf('%2d   %d   %6.3f %+6.3f %+6.3g %+6.3f %+6.3f\n', ...
     nit, iflag, elam, finfo(1), finfo(2), finfo(3), finfo(4));
 

Output from Monit
14   1   15.000 -0.322 -0.000108 +1.000 +206.000
13   1   15.000 -0.322 -5.67e-05 +2.000 +234.000
12   1   19.000 +0.257 -6.69e-05 +1.000 +226.000
11   2   17.225 +0.018 -6.75e-05 +1.000 +226.000
10   2   17.097 +0.000 -6.43e-05 +1.000 +226.000
 9   2   17.097 -0.000 -6.41e-05 +1.000 +226.000
 8   2   17.097 -0.000 -6.41e-05 +1.000 +226.000

bcondOut =

    1.0000    1.0000
         0         0
   -0.9064    0.9064


elamOut =

   17.0966


delamOut =

   1.4960e-04


ifail =

                    0


function d02ka_example
xl = 0;
xr = 3.141592653589793;
bcond = [1, 1;
     0, 0;
     0, -0.04564094560555453];
k = int64(4);
tol = 1e-05;
elam = 15;
delam = 4;
[bcondOut, elamOut, delamOut, ifail] = ...
    d02ka(xl, xr, @coeffn, bcond, k, tol, elam, delam, @monit)

function [p, q, dqdl] = coeffn(x, elam, jint)
  p = 1;
  dqdl = 1;
  q = elam - 10*cos(2.0*x);
function monit(nit, iflag, elam, finfo)
  if (nit == 14)
    fprintf('\nOutput from Monit\n');
  end

  fprintf('%2d   %d   %6.3f %+6.3f %+6.3g %+6.3f %+6.3f\n', ...
     nit, iflag, elam, finfo(1), finfo(2), finfo(3), finfo(4));
 

Output from Monit
14   1   15.000 -0.322 -0.000108 +1.000 +206.000
13   1   15.000 -0.322 -5.67e-05 +2.000 +234.000
12   1   19.000 +0.257 -6.69e-05 +1.000 +226.000
11   2   17.225 +0.018 -6.75e-05 +1.000 +226.000
10   2   17.097 +0.000 -6.43e-05 +1.000 +226.000
 9   2   17.097 -0.000 -6.41e-05 +1.000 +226.000
 8   2   17.097 -0.000 -6.41e-05 +1.000 +226.000

bcondOut =

    1.0000    1.0000
         0         0
   -0.9064    0.9064


elamOut =

   17.0966


delamOut =

   1.4960e-04


ifail =

                    0



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