NAG Toolbox: nag_ode_sl2_breaks_funs (d02ke)

Purpose

Syntax

Description

References

Parameters

Compulsory Input Parameters

1:     $\mathrm{xpoint}\left(\mathbf{m}\right)$ – double array

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

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

with the following meanings:

(a)	$x_1$ and $x_m$ are the left- and right-hand end points, $a$ and $b$, of the domain of definition of the Sturm–Liouville system if these are finite. If either $a$ or $b$ is infinite, the corresponding value $x_1$ or $x_m$ may be a more-or-less arbitrarily ‘large’ number of appropriate sign.
(b)	$x_2$ and $x_{m-1}$ are the Boundary Matching Points (BMPs), that is the points at which the left and right boundary conditions computed in bdyval are imposed. If the left-hand end point is a regular point then you should set $x_2=x_1$ $\left(=a\right)$, 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.
(c)	The remaining $m-4$ points $x_3,\dots ,x_{m-2}$, if any, define ‘break-points’ which divide the interval $\left[x_2,x_{m-1}\right]$ into $m-3$ sub-intervals $$i_1=\left[x_2,x_3\right],\dots ,i_{m-3}=\left[x_{m-2},x_{m-1}\right]\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\left(x;\lambda \right)$ 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\left(x\right)$ or $q\left(x;\lambda \right)$, or in their derivatives up to the fifth-order, should appear as break-points. Examples are given in Further Comments and Example. xpoint determines the position of the Shooting Matching Point (SMP), as explained in The Position of the Shooting Matching Point .

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

2:     $\mathrm{match}$ – int64int32nag_int scalar

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 ). If match is set to a value outside the range $\left[2,m-1\right]$ then a default value is taken, corresponding to the break-point nearest the centre of the interval $\left[\mathbf{xpoint}\left(2\right),\mathbf{xpoint}\left(m-1\right)\right]$.

3:     $\mathrm{coeffn}$ – function handle or string containing name of m-file

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

[p, q, dqdl] = coeffn(x, elam, jint)

Input Parameters

1:     $\mathrm{x}$ – double scalar

The current value of $x$.

2:     $\mathrm{elam}$ – double scalar

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

3:     $\mathrm{jint}$ – int64int32nag_int scalar

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

Output Parameters

1:     $\mathrm{p}$ – double scalar

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

2:     $\mathrm{q}$ – double scalar

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

3:     $\mathrm{dqdl}$ – double scalar

The value of $\frac{\partial q}{\partial \lambda }\left(x;\lambda \right)$ for the current value of $x$ and the current trial value of $\lambda$. However dqdl is only used in error estimation and, in the rare cases where it may be difficult to evaluate, an approximation (say to within $20\%$) will suffice.

4:     $\mathrm{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\left(x\right)$ and $p\left(x\right)y^{\prime }\left(x\right)$ which are consistent with the boundary conditions at the end points; only the ratio of the values matters. Here $x$ is a given point (xl or xr) equal to, or close to, the end point.

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

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

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

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

[yl, yr] = bdyval(xl, xr, elam)

Input Parameters

1:     $\mathrm{xl}$ – double scalar

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

2:     $\mathrm{xr}$ – double scalar

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

3:     $\mathrm{elam}$ – double scalar

The current trial value of $\lambda$.

Output Parameters

1:     $\mathrm{yl}\left(3\right)$ – double array

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

2:     $\mathrm{yr}\left(3\right)$ – double array

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

5:     $\mathrm{k}$ – int64int32nag_int scalar

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

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

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

6:     $\mathrm{tol}$ – double scalar

The tolerance argument which determines the accuracy of the computed eigenvalue. The error estimate held in delam on exit satisfies the mixed absolute/relative error test

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

(1)

where elam is the final estimate of the eigenvalue. delam is usually somewhat smaller than the right-hand side of (1) but not several orders of magnitude smaller.

Constraint: $\mathbf{tol}>0.0$.

7:     $\mathrm{elam}$ – double scalar

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

8:     $\mathrm{delam}$ – double scalar

An indication of the scale of the problem in the $\lambda$-direction. delam holds the initial ‘search step’ (positive or negative). Its value is not critical, but the first two trial evaluations are made at elam and $\mathbf{elam}+\mathbf{delam}$, so the 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 $\mathbf{delam}=0.0$ on entry, it is given the default value of $0.25\times \max \left(1.0,\left|\mathbf{elam}\right|\right)$.

9:     $\mathrm{hmax}\left(2,\mathbf{m}\right)$ – double array

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

It is recommended that $\mathbf{hmax}\left(1,j\right)$ be set to zero unless the coefficient functions $p$ and $q$ have features (such as a narrow peak) within the $j$th sub-interval that could be ‘missed’ if a long step were taken. In such a case $\mathbf{hmax}\left(1,j\right)$ 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.

10:   $\mathrm{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 arguments (see 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:     $\mathrm{nit}$ – int64int32nag_int scalar

The current value of the argument maxit of nag_ode_sl2_breaks_funs (d02ke), this is decreased by one at each iteration.

2:     $\mathrm{iflag}$ – int64int32nag_int scalar

Describes what phase the computation is in.

$\mathbf{iflag}<0$

An error occurred in the computation at this iteration; an error exit from nag_ode_sl2_breaks_funs (d02ke) with $\mathbf{ifail}=-\mathbf{iflag}$ will follow.

$\mathbf{iflag}=1$

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

$\mathbf{iflag}=2$

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

3:     $\mathrm{elam}$ – double scalar

The current trial value of $\lambda$.

4:     $\mathrm{finfo}\left(15\right)$ – 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 $\mathbf{iflag}>0$), though the first few components may be of interest to you. In case of an error ($\mathbf{iflag}<0$) all the components of finfo should be printed.

The contents of finfo are as follows:

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

The current value of the ‘miss-distance’ or ‘residual’ function $f\left(\lambda \right)$ on which the shooting method is based. (See General Description of the Algorithm for further information.) $\mathbf{finfo}\left(1\right)$ is set to zero if $\mathbf{iflag}<0$.

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

An estimate of the quantity $\partial \lambda$ defined as follows. Consider the perturbation in the miss-distance $f\left(\lambda \right)$ that would result if the local error in the solution of the differential equation were always positive and equal to its maximum permitted value. Then $\partial \lambda$ is the perturbation in $\lambda$ that would have the same effect on $f\left(\lambda \right)$. Thus, at the zero of $f\left(\lambda \right),\left|\partial \lambda \right|$ is an approximate bound on the perturbation of the zero (that is the eigenvalue) caused by errors in numerical solution. If $\partial \lambda$ is very large then it is possible that there has been a programming error in coeffn such that $q$ is independent of $\lambda$. If this is the case, an error exit with $\mathbf{ifail}=\mathbf{5}$ should follow. $\mathbf{finfo}\left(2\right)$ is set to zero if $\mathbf{iflag}<0$.

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

The number of internal iterations, using the same value of $\lambda$ and tighter accuracy tolerances, needed to bring the accuracy (that is, the value of $\partial \lambda$) to an acceptable value. Its value should normally be $1.0$, and should almost never exceed $2.0$.

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

The number of calls to coeffn at this iteration.

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

The number of successful steps taken by the internal differential equation solver at this iteration. A step is successful if it is used to advance the integration.

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

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

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

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

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

Not used.

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

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

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

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

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

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

$\mathbf{finfo}\left(11\right)$, $\mathbf{finfo}\left(12\right)$, $\mathbf{finfo}\left(13\right)$

The current values of the Pruefer dependent variables $\beta$, $\varphi$ and $\rho$ respectively. These are set to zero in case of an error in bdyval.

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

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

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

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

11:   $\mathrm{report}$ – function handle or string containing name of m-file

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

report(x, v, jint)

Input Parameters

1:     $\mathrm{x}$ – double scalar

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

2:     $\mathrm{v}\left(3\right)$ – double array

$\mathbf{v}\left(1\right)$, $\mathbf{v}\left(2\right)$, $\mathbf{v}\left(3\right)$ hold the current values of the Pruefer variables $\beta$, $\varphi$, $\rho$ respectively.

3:     $\mathrm{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=\mathbf{xpoint}\left(2\right)$) jint is set to $0$ and at the extreme right-hand end point (when $x=x_r=\mathbf{xpoint}\left(m-1\right)$) jint is set to $m-2$.

Optional Input Parameters

1:     $\mathrm{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: $\mathbf{m}\ge 4$.

2:     $\mathrm{maxit}$ – int64int32nag_int scalar

Default: $0$

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

3:     $\mathrm{maxfun}$ – int64int32nag_int scalar

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

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

Default: $0$

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

Output Parameters

1:     $\mathrm{match}$ – int64int32nag_int scalar

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

2:     $\mathrm{elam}$ – double scalar

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

3:     $\mathrm{delam}$ – double scalar

If $\mathbf{ifail}=\mathbf{0}$, delam holds an estimate of the absolute error in the computed eigenvalue, that is $\left|\tilde{\lambda }-\mathbf{elam}\right|\simeq \mathbf{delam}$. (In 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 $\mathbf{ifail}\ne \mathbf{0}$, delam may hold an estimate of the error, or its initial value, depending on the value of ifail. See Error Indicators and Warnings for further details.

4:     $\mathrm{hmax}\left(2,\mathbf{m}\right)$ – double array

$\mathbf{hmax}\left(1,m-1\right)$ and $\mathbf{hmax}\left(1,m\right)$ contain the sensitivity coefficients ${\sigma }_l,{\sigma }_r$, described in The Sensitivity Parameters and . Other entries contain diagnostic output in the case of an error exit (see Error Indicators and Warnings).

5:     $\mathrm{maxit}$ – int64int32nag_int scalar

Default: $0$

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

6:     $\mathrm{ifail}$ – int64int32nag_int scalar

$\mathbf{ifail}=\mathbf{0}$ unless the function detects an error (see Error Indicators and Warnings).

Error Indicators and Warnings

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

   $\mathbf{ifail}=1$

A argument error. All arguments (except ifail) are left unchanged. The reason for the error is shown by the value of $\mathbf{hmax}\left(2,1\right)$ as follows:

$\mathbf{hmax}\left(2,1\right)=1$:	$\mathbf{m}<4$;
$\mathbf{hmax}\left(2,1\right)=2$:	$\mathbf{k}<0$;
$\mathbf{hmax}\left(2,1\right)=3$:	$\mathbf{tol}\le 0.0$;
$\mathbf{hmax}\left(2,1\right)=4$:	$\mathbf{xpoint}\left(1\right)$ to $\mathbf{xpoint}\left(m\right)$ are not in ascending order. $\mathbf{hmax}\left(2,2\right)$ gives the position $i$ in xpoint where this was detected.

   $\mathbf{ifail}=2$

At some call to bdyval, invalid values were returned, that is, either $\mathbf{yl}\left(1\right)=\mathbf{yl}\left(2\right)=0.0$, or $\mathbf{yr}\left(1\right)=\mathbf{yr}\left(2\right)=0.0$ (a programming error in bdyval). See the last call of monit for details.

This error exit will also occur if $p\left(x\right)$ is zero at the point where the boundary condition is imposed. Probably bdyval was called with xl equal to a singular end point $a$ or with xr equal to a singular end point $b$.

   $\mathbf{ifail}=3$

At some point between xl and xr the value of $p\left(x\right)$ computed by coeffn became zero or changed sign. See the last call of monit for details.

   $\mathbf{ifail}=4$

$\mathbf{maxit}>0$ on entry, and after maxit iterations the eigenvalue had not been found to the required accuracy.

   $\mathbf{ifail}=5$

The ‘bracketing’ phase (with argument iflag of the monit equal to $1$) failed to bracket the eigenvalue within ten iterations. This is caused by an error in formulating the problem (for example, $q$ is independent of $\lambda$), or by very poor initial estimates of elam and delam.

On exit, elam and $\mathbf{elam}+\mathbf{delam}$ give the end points of the interval within which no eigenvalue was located by the function.

   $\mathbf{ifail}=6$

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

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

   $\mathbf{ifail}=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 $\mathbf{ifail}=\mathbf{7}$). This could be due to pathological behaviour of $p\left(x\right)$ and $q\left(x;\lambda \right)$ or to an unreasonable accuracy requirement or to the current value of $\lambda$ making the equations ‘stiff’. See the last call of monit for details.

W  $\mathbf{ifail}=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.

   $\mathbf{ifail}=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\left(\lambda \right)$. This error exit should not occur, but if it does, try solving the problem again with a smaller value for tol.

   $\mathbf{ifail}=11$

   $\mathbf{ifail}=12$

A serious error has occurred in an internal call. Check all (sub)program calls and array dimensions. Seek expert help.

   $\mathbf{ifail}=-99$

An unexpected error has been triggered by this routine. Please contact NAG.

   $\mathbf{ifail}=-399$

Your licence key may have expired or may not have been installed correctly.

   $\mathbf{ifail}=-999$

Dynamic memory allocation failed.

Accuracy

Further Comments

Timing

General Description of the Algorithm

nag_ode_sl2_breaks_funs (d02ke) defines a miss-distance $f\left(\lambda \right)$ based on the rotation about the origin of the point $\mathbf{p}\left(x\right)=\left(p\left(x\right)y^{\prime }\left(x\right),y\left(x\right)\right)$ in the Phase Plane as the solution proceeds from $a$ to $b$. The boundary conditions define the ray (i.e., two-sided line through the origin) on which $p\left(x\right)$ should start, and the ray on which it should finish. The eigenvalue $k$ defines the total number of half-turns it should make. Numerical solution is actually done by ‘shooting forward’ from $x=a$ and ‘shooting backward’ from $x=b$ to a matching point $x=c$. Then $f\left(\lambda \right)$ is taken as the angle between the rays to the two resulting points $P_a\left(c\right)$ and $P_b\left(c\right)$. A relative scaling of the $p{y}^{\prime }$ and $y$ axes, based on the behaviour of the coefficient functions $p$ and $q$, is used to improve the numerical behaviour.

Figure 1

The Position of the Shooting Matching Point c

Note that the shooting method used by the code integrates first from the left-hand end $x_l$, then from the right-hand end $x_r$, to meet at the matching point $c$ in the middle. This will of course be reflected in printed or graphical output. The diagram shows a possible sequence of nine mesh points ${\tau }_1$ through ${\tau }_9$ in the order in which they appear, assuming there are just two sub-intervals (so $m=5$).

Figure 2

The choice of matching point $c$ is open. If we choose $c=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.0$ are given schematically in Figure 3.

Figure 3

If we choose $c$ as the maximum of the hump in $q\left(x;\lambda \right)$ (see item (a) above) we instead obtain the accurate results given in Figure 4

Figure 4

Examples of Coding the coeffn

function [p, q, dqdl] = coeffn(x, elam, jint)
  global nu;
  ...
  p = x;
  q = elam*x - nu*nu/x
  dqdl = x;

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

  function [p, q, dqdl] = coeffn(x, elam, jint)
  ...
  p = 1;
  dqdl = 1 - 2*exp(-100*x*x);
  q = elam * dqdl;

Examples of Boundary Conditions at Singular Points

function [yl, yr] = bdyval(xl, xr, elam)
  yl(1) = xl;
  yl(2) = 2;
  yr(1) = 1;
  yr(2) = -sqrt(xr-elam);

The Sensitivity Parameters σl and σr

‘Missed Zeros’

Example

Open in the MATLAB editor: d02ke_example

function d02ke_example


fprintf('d02ke example results\n\n');

% 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);

% report outputs intermediate results.
[match, elam, delam, hmax, maxit, ifail] = ...
d02ke( ...
       xpoint, int64(match), @coeffn, @bdyval, ...
       int64(k), tol, elam, delam, hmax, 'd02kay', @report);

% Output final results.
fprintf('\n\n');
fprintf('Eigenvalue number, k = %d\n',k);
fprintf(' estimated value     = %6.3f\n',elam);
fprintf(' estimated error     = %9.2e\n',delam);
fprintf(' sensitivity to left  boundary = %6.3f\n',hmax(1,m-1));
fprintf(' sensitivity to right boundary = %6.3f\n\n', hmax(1,m));
fig1 = figure;
display_plot(xkeep,pkeep,ykeep)


function [yl, yr] = bdyval(xl, xr, elam)
  % Define the boundary conditions.
  yl(1) = xl;
  yl(2) = 2;
  yl(3) = 0;
  yr(1) = 1;
  yr(2) = -sqrt(xr-elam);
  yr(3) = 0;

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');
  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.
  ncall = ncall+1;
  ykeep(ncall,1) = y;
  xkeep(ncall,1) = x;
  pkeep(ncall,1) = pyp;

function display_plot(xkeep, pkeep, ykeep)

  % 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'});
  % Label the axes.
  xlabel('x'); ylabel('Eigenfunction and Derivative');
  % Add a legend.
  legend('y','p(x)y''','Location','Best');

d02ke example results


 A singular problem


Eigenvalue number, k = 11
 estimated value     = 14.946
 estimated error     =  9.60e-04
 sensitivity to left  boundary = -0.015
 sensitivity to right boundary =  0.000