NAG Library Routine Document
d02saf (bvp_shoot_genpar_algeq)
1
Purpose
d02saf solves a twopoint boundary value problem for a system of firstorder ordinary differential equations with boundary conditions, combined with additional algebraic equations. It uses initial value techniques and a modified Newton iteration in a shooting and matching method.
2
Specification
Fortran Interface
Subroutine d02saf ( 
p, m, n, n1, pe, pf, e, dp, npoint, swp, ldswp, icount, range, bc, fcn, eqn, constr, ymax, monit, prsol, w, ldw, sdw, ifail) 
Integer, Intent (In)  ::  m, n, n1, npoint, ldswp, icount, ldw, sdw  Integer, Intent (Inout)  ::  ifail  Real (Kind=nag_wp), Intent (In)  ::  pe(m), e(n), ymax  Real (Kind=nag_wp), Intent (Inout)  ::  p(m), pf(m), dp(m), swp(ldswp,6), w(ldw,sdw)  Logical, External  ::  constr  External  ::  range, bc, fcn, eqn, monit, prsol 

C Header Interface
#include nagmk26.h
void 
d02saf_ (double p[], const Integer *m, const Integer *n, const Integer *n1, const double pe[], double pf[], const double e[], double dp[], const Integer *npoint, double swp[], const Integer *ldswp, const Integer *icount, void (NAG_CALL *range)(double x[], const Integer *npoint, const double p[], const Integer *m), void (NAG_CALL *bc)(double g1[], double g2[], const double p[], const Integer *m, const Integer *n), void (NAG_CALL *fcn)(const double *x, const double y[], double f[], const Integer *n, const double p[], const Integer *m, const Integer *i), void (NAG_CALL *eqn)(double e[], const Integer *q, const double p[], const Integer *m), logical (NAG_CALL *constr)(const double p[], const Integer *m), const double *ymax, void (NAG_CALL *monit)(const Integer *istate, const Integer *iflag, const Integer *ifail1, const double p[], const Integer *m, const double f[], const double *pnorm, const double *pnorm1, const double *eps, const double d[]), void (NAG_CALL *prsol)(double *z, const double y[], const Integer *n), double w[], const Integer *ldw, const Integer *sdw, Integer *ifail) 

3
Description
d02saf solves a twopoint boundary value problem for a system of
$n$ firstorder ordinary differential equations with separated boundary conditions by determining certain unknown arguments
${p}_{1},{p}_{2},\dots ,{p}_{m}$. (There may also be additional algebraic equations to be solved in the determination of the arguments and, if so, these equations are defined by
eqn.) The arguments may be, but need not be, boundary values; they may include eigenvalues, arguments in the coefficients of the differential equations, coefficients in series expansions or asymptotic expansions for boundary values, the length of the range of definition of the system of differential equations, etc.
It is assumed that we have a system of
$n$ differential equations of the form
where
$p={\left({p}_{1},{p}_{2},\dots ,{p}_{m}\right)}^{\mathrm{T}}$ is the vector of arguments, and that the derivative
$f$ is evaluated by
fcn. Also,
${n}_{1}$ of the equations are assumed to depend on
$p$. For
${n}_{1}<n$ the
$n{n}_{1}$ equations of the system are not involved in the matching process. These are the driving equations; they should be independent of
$p$ and of the solution of the other
${n}_{1}$ equations. In numbering the equations in
fcn and
bc the driving equations must be put
first (as they naturally occur in most applications). The range of definition [
$a,b$] of the differential equations is defined by
range and may depend on the arguments
${p}_{1},{p}_{2},\dots ,{p}_{m}$ (that is, on
$p$).
range must define the points
${x}_{1},{x}_{2},\dots ,{x}_{{\mathbf{npoint}}}$,
${\mathbf{npoint}}\ge 2$, which must satisfy
(or a similar relationship with all the inequalities reversed).
If
${\mathbf{npoint}}>2$ the points
${x}_{1},{x}_{2},\dots ,{x}_{{\mathbf{npoint}}}$ can be used to break up the range of definition. Integration is restarted at each of these points. This means that the differential equations
(1) can be defined differently in each subinterval
$\left[{x}_{\mathit{i}},{x}_{\mathit{i}+1}\right]$, for
$\mathit{i}=1,2,\dots ,{\mathbf{npoint}}1$. Also, since initial and maximum integration step sizes can be supplied on each subinterval (via the array
swp), you can indicate parts of the range
$\left[a,b\right]$ where the solution
$y\left(x\right)$ may be difficult to obtain accurately and can take appropriate action.
The boundary conditions may also depend on the arguments and are applied at
$a={x}_{1}$ and
$b={x}_{{\mathbf{npoint}}}$. They are defined (in
bc) in the form
The boundary value problem is solved by determining the unknown arguments
$p$ by a shooting and matching technique. The differential equations are always integrated from
$a$ to
$b$ with initial values
$y\left(a\right)={g}_{1}\left(p\right)$. The solution vector thus obtained at
$x=b$ is subtracted from the vector
${g}_{2}\left(p\right)$ to give the
${n}_{1}$ residuals
${r}_{1}\left(p\right)$, ignoring the first
$n{n}_{1}$, driving equations. Because the direction of integration is always from
$a$ to
$b$, it is unnecessary, in
bc, to supply values for the first
$n{n}_{1}$ boundary values at
$b$, that is the first
$n{n}_{1}$ components of
${g}_{2}$ in
(3). For
${n}_{1}<m$ then
${r}_{1}\left(p\right)$. Together with the
$m{n}_{1}$ equations defined by
eqn,
these give a vector of residuals
$r$, which at the solution,
$p$, must satisfy
These equations are solved by a pseudoNewton iteration which uses a modified singular value decomposition of
$J=\frac{\partial r}{\partial p}$ when solving the linear equations which arise. The Jacobian
$J$ used in Newton's method is obtained by numerical differentiation. The arguments at each Newton iteration are accepted only if the norm
${\Vert {D}^{1}{\stackrel{~}{J}}^{+}r\Vert}_{2}$ is much reduced from its previous value. Here
${\stackrel{~}{J}}^{+}$ is the pseudoinverse, calculated from the singular value decomposition, of a modified version of the Jacobian
$J$ (
${J}^{+}$ is actually the inverse of the Jacobian in wellconditioned cases).
$D$ is a diagonal matrix with
where
pf is an array of floor values.
See
Deuflhard (1974) for further details of the variants of Newton's method used,
Gay (1976) for the modification of the singular value decomposition and
Gladwell (1979) for an overview of the method used.
Two facilities are provided to prevent the pseudoNewton iteration running into difficulty. First, you are permitted to specify constraints on the values of the arguments
$p$ via a
constr. These constraints are only used to prevent the Newton iteration using values for
$p$ which would violate them; that is, they are not used to determine the values of
$p$. Secondly, you are permitted to specify a maximum value
${y}_{\mathrm{max}}$ for
${\Vert y\left(x\right)\Vert}_{\infty}$ at all points in the range
$\left[a,b\right]$. It is intended that this facility be used to prevent machine ‘overflow’ in the integrations of equation
(1) due to poor choices of the arguments
$p$ which might arise during the Newton iteration. When using this facility, it is presumed that you have an estimate of the likely size of
${\Vert y\left(x\right)\Vert}_{\infty}$ at all points
$x\in \left[a,b\right]$.
${y}_{\mathrm{max}}$ should then be chosen rather larger (say by a factor of
$10$) than this estimate.
You are strongly advised to supply a
monit (or to call the ‘default’ routine d02hbx, see
monit) to monitor the progress of the pseudoNewton iteration. You can output the solution of the problem
$y\left(x\right)$ by supplying a suitable
prsol (an example is given in
Section 10 of a routine designed to output the solution at equally spaced points).
d02saf is designed to try all possible options before admitting failure and returning to you. Provided the routine can start the Newton iteration from the initial point $p$ it will exhaust all the options available to it (though you can override this by specifying a maximum number of iterations to be taken). The fact that all its options have been exhausted is the only error exit from the iteration. Other error exits are possible, however, whilst setting up the Newton iteration and when computing the final solution.
If you require more background information about the solution of boundary value problems by shooting methods you are recommended to read the appropriate chapters of
Hall and Watt (1976), and for a detailed description of
d02saf
Gladwell (1979) is recommended.
4
References
Deuflhard P (1974) A modified Newton method for the solution of illconditioned systems of nonlinear equations with application to multiple shooting Numer. Math. 22 289–315
Gay D (1976) On modifying singular values to solve possibly singular systems of nonlinear equations Working Paper 125 Computer Research Centre, National Bureau for Economics and Management Science, Cambridge, MA
Gladwell I (1979) The development of the boundary value codes in the ordinary differential equations chapter of the NAG Library Codes for Boundary Value Problems in Ordinary Differential Equations. Lecture Notes in Computer Science (eds B Childs, M Scott, J W Daniel, E Denman and P Nelson) 76 Springer–Verlag
Hall G and Watt J M (ed.) (1976) Modern Numerical Methods for Ordinary Differential Equations Clarendon Press, Oxford
5
Arguments
 1: $\mathbf{p}\left({\mathbf{m}}\right)$ – Real (Kind=nag_wp) arrayInput/Output

On entry: ${\mathbf{p}}\left(\mathit{i}\right)$ must be set to an estimate of the $\mathit{i}$th argument, ${p}_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,m$.
On exit: the corrected value for the $i$th argument, unless an error has occurred, when it contains the last calculated value of the argument.
 2: $\mathbf{m}$ – IntegerInput

On entry: $m$, the number of arguments.
Constraint:
${\mathbf{m}}>0$.
 3: $\mathbf{n}$ – IntegerInput

On entry: $n$, the total number of differential equations.
Constraint:
${\mathbf{n}}>0$.
 4: $\mathbf{n1}$ – IntegerInput

On entry:
${n}_{1}$, the number of differential equations active in the matching process. The active equations must be placed last in the numbering in
fcn and
bc. The
first ${\mathbf{n}}{\mathbf{n1}}$ equations are used as the driving equations.
Constraint:
${\mathbf{n1}}\le {\mathbf{n}}$, ${\mathbf{n1}}\le {\mathbf{m}}$ and ${\mathbf{n1}}>0$.
 5: $\mathbf{pe}\left({\mathbf{m}}\right)$ – Real (Kind=nag_wp) arrayInput

On entry:
${\mathbf{pe}}\left(\mathit{i}\right)$, for
$\mathit{i}=1,2,\dots ,m$, must be set to a positive value for use in the convergence test in the
$i$th argument
${p}_{i}$. See the description of
pf for further details.
Constraint:
${\mathbf{pe}}\left(\mathit{i}\right)>0.0$, for $\mathit{i}=1,2,\dots ,m$.
 6: $\mathbf{pf}\left({\mathbf{m}}\right)$ – Real (Kind=nag_wp) arrayInput/Output

On entry:
${\mathbf{pf}}\left(\mathit{i}\right)$, for
$\mathit{i}=1,2,\dots ,m$, should be set to a ‘floor’ value in the convergence test on the
$i$th argument
${p}_{i}$. If
${\mathbf{pf}}\left(i\right)\le 0.0$ on entry then it is set to the small positive value
$\sqrt{\epsilon}$ (where
$\epsilon $ may in most cases be considered to be
machine precision); otherwise it is used unchanged.
The Newton iteration is presumed to have converged if a full Newton step is taken (
${\mathbf{istate}}=1$ in the specification of
monit), the singular values of the Jacobian are not being significantly perturbed (also see
monit) and if the Newton correction
${C}_{i}$ satisfies
where
${p}_{i}$ is the current value of the
$i$th argument. The values
${\mathbf{pf}}\left(i\right)$ are also used in determining the Newton iterates as discussed in
Section 3, see equation
(6).
On exit: the values actually used.
 7: $\mathbf{e}\left({\mathbf{n}}\right)$ – Real (Kind=nag_wp) arrayInput

On entry: values for use in controlling the local error in the integration of the differential equations. If
${\mathit{err}}_{\mathit{i}}$ is an estimate of the local error in
${y}_{\mathit{i}}$, for
$\mathit{i}=1,2,\dots ,n$,
where
$\epsilon $ may in most cases be considered to be
machine precision.
Suggested value:
${\mathbf{e}}\left(i\right)={10}^{5}$.
Constraint:
${\mathbf{e}}\left(\mathit{i}\right)>0.0$, for $\mathit{i}=1,2,\dots ,{\mathbf{n}}$.
 8: $\mathbf{dp}\left({\mathbf{m}}\right)$ – Real (Kind=nag_wp) arrayInput/Output

On entry: a value to be used in perturbing the argument
${p}_{i}$ in the numerical differentiation to estimate the Jacobian used in Newton's method. If
${\mathbf{dp}}\left(i\right)=0.0$ on entry, an estimate is made internally by setting
where
${p}_{i}$ is the initial value of the argument supplied by you and
$\epsilon $ may in most cases be considered to be
machine precision. The estimate of the Jacobian,
$J$, is made using forward differences, that is for each
$\mathit{i}$, for
$\mathit{i}=1,2,\dots ,m$,
${p}_{\mathit{i}}$ is perturbed to
${p}_{\mathit{i}}+{\mathbf{dp}}\left(\mathit{i}\right)$ and the
$\mathit{i}$th column of
$J$ is estimated as
where the other components of
$p$ are unchanged (see
(3) for the notation used). If this fails to produce a Jacobian with significant columns, backward differences are tried by perturbing
${p}_{\mathit{i}}$ to
${p}_{\mathit{i}}{\mathbf{dp}}\left(\mathit{i}\right)$ and if this also fails then central differences are used with
${p}_{\mathit{i}}$ perturbed to
${p}_{\mathit{i}}+10.0\times {\mathbf{dp}}\left(\mathit{i}\right)$. If this also fails then the calculation of the Jacobian is abandoned. If the Jacobian has not previously been calculated then an error exit is taken. If an earlier estimate of the Jacobian is available then the current argument set,
${p}_{\mathit{i}}$, for
$\mathit{i}=1,2,\dots ,m$, is abandoned in favour of the last argument set from which useful progress was made and the singular values of the Jacobian used at the point are modified before proceeding with the Newton iteration. You are recommended to use the default value
${\mathbf{dp}}\left(i\right)=0.0$ unless you have prior knowledge of a better choice. If any of the perturbations described are likely to lead to an unfortunate set of argument values then you should use the LOGICAL FUNCTION
constr to prevent such perturbations (all changes of arguments are checked by a call to
constr).
On exit: the values actually used.
 9: $\mathbf{npoint}$ – IntegerInput

On entry: 2 plus the number of breakpoints in the range of definition of the system of differential equations
(1).
Constraint:
${\mathbf{npoint}}\ge 2$.
 10: $\mathbf{swp}\left({\mathbf{ldswp}},6\right)$ – Real (Kind=nag_wp) arrayInput/Output

On entry:
${\mathbf{swp}}\left(i,1\right)$ must contain an estimate for an initial step size for integration across the
$i$th subinterval
$\left[{\mathbf{x}}\left(\mathit{i}\right),{\mathbf{x}}\left(\mathit{i}+1\right)\right]$, for
$\mathit{i}=1,2,\dots ,{\mathbf{npoint}}1$, (see
range).
${\mathbf{swp}}\left(i,1\right)$ should have the same sign as
${\mathbf{x}}\left(i+1\right){\mathbf{x}}\left(i\right)$ if it is nonzero. If
${\mathbf{swp}}\left(i,1\right)=0.0$, on entry, a default value for the initial step size is calculated internally. This is the recommended mode of entry.
${\mathbf{swp}}\left(i,3\right)$ must contain a lower bound for the modulus of the step size on the
$\mathit{i}$th subinterval $\left[{\mathbf{x}}\left(\mathit{i}\right),{\mathbf{x}}\left(\mathit{i}+1\right)\right]$, for $\mathit{i}=1,2,\dots ,{\mathbf{npoint}}1$. If ${\mathbf{swp}}\left(i,3\right)=0.0$ on entry, a very small default value is used. By setting ${\mathbf{swp}}\left(i,3\right)>0.0$ but smaller than the expected step sizes (assuming you have some insight into the likely step sizes) expensive integrations with arguments $p$ far from the solution can be avoided.
${\mathbf{swp}}\left(\mathit{i},2\right)$ must contain an upper bound on the modulus of the step size to be used in the integration on $\left[{\mathbf{x}}\left(\mathit{i}\right),{\mathbf{x}}\left(\mathit{i}+1\right)\right]$, for $\mathit{i}=1,2,\dots ,{\mathbf{npoint}}1$. If ${\mathbf{swp}}\left(i,2\right)=0.0$ on entry no bound is assumed. This is the recommended mode of entry unless the solution is expected to have important features which might be ‘missed’ in the integration if the step size were permitted to be chosen freely.
On exit:
${\mathbf{swp}}\left(\mathit{i},1\right)$ contains the initial step size used on the last integration on
$\left[{\mathbf{x}}\left(\mathit{i}\right),{\mathbf{x}}\left(\mathit{i}+1\right)\right]$, for
$\mathit{i}=1,2,\dots ,{\mathbf{npoint}}1$, (excluding integrations during the calculation of the Jacobian).
${\mathbf{swp}}\left(\mathit{i},2\right)$, for $\mathit{i}=1,2,\dots ,{\mathbf{npoint}}1$, is usually unchanged. If the maximum step size ${\mathbf{swp}}\left(i,2\right)$ is so small or the length of the range $\left[{\mathbf{x}}\left(i\right),{\mathbf{x}}\left(i+1\right)\right]$ is so short that on the last integration the step size was not controlled in the main by the size of the error tolerances ${\mathbf{e}}\left(i\right)$ but by these other factors, ${\mathbf{swp}}\left({\mathbf{npoint}},2\right)$ is set to the floatingpoint value of $i$ if the problem last occurred in $\left[{\mathbf{x}}\left(i\right),{\mathbf{x}}\left(i+1\right)\right]$. Any results obtained when this value is returned as nonzero should be viewed with caution.
${\mathbf{swp}}\left(\mathit{i},3\right)$, for $\mathit{i}=1,2,\dots ,{\mathbf{npoint}}1$, are unchanged.
If an error exit with
${\mathbf{ifail}}={\mathbf{4}}$,
${\mathbf{5}}$ or
${\mathbf{6}}$ (see
Section 6) occurs on the integration made from
${\mathbf{x}}\left(i\right)$ to
${\mathbf{x}}\left(i+1\right)$ the floatingpoint value of
$i$ is returned in
${\mathbf{swp}}\left({\mathbf{npoint}},1\right)$. The actual point
$x\in \left[{\mathbf{x}}\left(i\right),{\mathbf{x}}\left(i+1\right)\right]$ where the error occurred is returned in
${\mathbf{swp}}\left(1,5\right)$ (see also the specification of
w). The floatingpoint value of
npoint is returned in
${\mathbf{swp}}\left({\mathbf{npoint}},1\right)$ if the error exit is caused by a call to
bc.
If an error exit occurs when estimating the Jacobian matrix (
${\mathbf{ifail}}={\mathbf{7}}$,
${\mathbf{8}}$,
${\mathbf{9}}$,
${\mathbf{10}}$,
${\mathbf{11}}$ or
${\mathbf{12}}$, see
Section 6) and if argument
${p}_{i}$ was the cause of the failure then on exit
${\mathbf{swp}}\left({\mathbf{npoint}},1\right)$ contains the floatingpoint value of
$i$.
${\mathbf{swp}}\left(\mathit{i},4\right)$ contains the point ${\mathbf{x}}\left(\mathit{i}\right)$, for $\mathit{i}=1,2,\dots ,{\mathbf{npoint}}$, used at the solution $p$ or at the final values of $p$ if an error occurred.
swp is also partly used as workspace.
 11: $\mathbf{ldswp}$ – IntegerInput

On entry: the first dimension of the array
swp as declared in the (sub)program from which
d02saf is called.
Constraint:
${\mathbf{ldswp}}\ge {\mathbf{npoint}}$.
 12: $\mathbf{icount}$ – IntegerInput

On entry: an upper bound on the number of Newton iterations. If ${\mathbf{icount}}=0$ on entry, no check on the number of iterations is made (this is the recommended mode of entry).
Constraint:
${\mathbf{icount}}\ge 0$.
 13: $\mathbf{range}$ – Subroutine, supplied by the user.External Procedure

range must specify the breakpoints
${x}_{\mathit{i}}$, for
$\mathit{i}=1,2,\dots ,{\mathbf{npoint}}$, which may depend on the arguments
${p}_{\mathit{j}}$, for
$\mathit{j}=1,2,\dots ,m$.
The specification of
range is:
Fortran Interface
Integer, Intent (In)  ::  npoint, m  Real (Kind=nag_wp), Intent (In)  ::  p(m)  Real (Kind=nag_wp), Intent (Out)  ::  x(npoint) 

C Header Interface
#include nagmk26.h
void 
range (double x[], const Integer *npoint, const double p[], const Integer *m) 

 1: $\mathbf{x}\left({\mathbf{npoint}}\right)$ – Real (Kind=nag_wp) arrayOutput

On exit: the
$\mathit{i}$th breakpoint, for
$\mathit{i}=1,2,\dots ,{\mathbf{npoint}}$. The sequence
$\left({\mathbf{x}}\left(i\right)\right)$ must be strictly monotonic, that is either
or
 2: $\mathbf{npoint}$ – IntegerInput

On entry: $2$ plus the number of breakpoints in $\left(a,b\right)$.
 3: $\mathbf{p}\left({\mathbf{m}}\right)$ – Real (Kind=nag_wp) arrayInput

On entry: the current estimate of the
$\mathit{i}$th argument, for $\mathit{i}=1,2,\dots ,m$.
 4: $\mathbf{m}$ – IntegerInput

On entry: $m$, the number of arguments.
range must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which
d02saf is called. Arguments denoted as
Input must
not be changed by this procedure.
Note: range should not return floatingpoint NaN (Not a Number) or infinity values, since these are not handled by
d02saf. If your code inadvertently
does return any NaNs or infinities,
d02saf is likely to produce unexpected results.
 14: $\mathbf{bc}$ – Subroutine, supplied by the user.External Procedure

bc must place in
g1 and
g2 the boundary conditions at
$a$ and
$b$ respectively.
The specification of
bc is:
Fortran Interface
Integer, Intent (In)  ::  m, n  Real (Kind=nag_wp), Intent (In)  ::  p(m)  Real (Kind=nag_wp), Intent (Out)  ::  g1(n), g2(n) 

C Header Interface
#include nagmk26.h
void 
bc (double g1[], double g2[], const double p[], const Integer *m, const Integer *n) 

 1: $\mathbf{g1}\left({\mathbf{n}}\right)$ – Real (Kind=nag_wp) arrayOutput

On exit: the value of
${y}_{\mathit{i}}\left(a\right)$, for $\mathit{i}=1,2,\dots ,n$, (where this may be a known value or a function of the parameters
${p}_{\mathit{j}}$, for $\mathit{j}=1,2,\dots ,m$).
 2: $\mathbf{g2}\left({\mathbf{n}}\right)$ – Real (Kind=nag_wp) arrayOutput

On exit: the value of
${y}_{\mathit{i}}\left(b\right)$, for
$\mathit{i}=1,2,\dots ,n$, (where these may be known values or functions of the parameters
${p}_{\mathit{j}}$, for
$\mathit{j}=1,2,\dots ,m$). If
$n>{n}_{1}$, so that there are some driving equations, the first
$n{n}_{1}$ values of
g2 need not be set since they are never used.
 3: $\mathbf{p}\left({\mathbf{m}}\right)$ – Real (Kind=nag_wp) arrayInput

On entry: an estimate of the
$\mathit{i}$th argument, ${p}_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,m$.
 4: $\mathbf{m}$ – IntegerInput

On entry: $m$, the number of arguments.
 5: $\mathbf{n}$ – IntegerInput

On entry: $n$, the number of differential equations.
bc must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which
d02saf is called. Arguments denoted as
Input must
not be changed by this procedure.
Note: bc should not return floatingpoint NaN (Not a Number) or infinity values, since these are not handled by
d02saf. If your code inadvertently
does return any NaNs or infinities,
d02saf is likely to produce unexpected results.
 15: $\mathbf{fcn}$ – Subroutine, supplied by the user.External Procedure

fcn must evaluate the functions
${f}_{\mathit{i}}$ (i.e., the derivatives
${y}_{\mathit{i}}^{\prime}$), for
$\mathit{i}=1,2,\dots ,n$.
The specification of
fcn is:
Fortran Interface
Subroutine fcn ( 
x, y, f, n, p, m, i) 
Integer, Intent (In)  ::  n, m, i  Real (Kind=nag_wp), Intent (In)  ::  x, y(n), p(m)  Real (Kind=nag_wp), Intent (Out)  ::  f(n) 

C Header Interface
#include nagmk26.h
void 
fcn (const double *x, const double y[], double f[], const Integer *n, const double p[], const Integer *m, const Integer *i) 

 1: $\mathbf{x}$ – Real (Kind=nag_wp)Input

On entry: $x$, the value of the argument.
 2: $\mathbf{y}\left({\mathbf{n}}\right)$ – Real (Kind=nag_wp) arrayInput

On entry: ${y}_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,\mathit{n}$, the value of the argument.
 3: $\mathbf{f}\left({\mathbf{n}}\right)$ – Real (Kind=nag_wp) arrayOutput

On exit: the derivative of
${y}_{\mathit{i}}$, for
$\mathit{i}=1,2,\dots ,n$, evaluated at
$x$.
${\mathbf{f}}\left(i\right)$ may depend upon the parameters
${p}_{\mathit{j}}$, for
$\mathit{j}=1,2,\dots ,m$. If there are any driving equations (see
Section 3) then these must be numbered first in the ordering of the components of
f.
 4: $\mathbf{n}$ – IntegerInput

On entry: $\mathit{n}$, the number of equations.
 5: $\mathbf{p}\left({\mathbf{m}}\right)$ – Real (Kind=nag_wp) arrayInput

On entry: the current estimate of the
$\mathit{i}$th argument ${p}_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,m$.
 6: $\mathbf{m}$ – IntegerInput

On entry: $m$, the number of arguments.
 7: $\mathbf{i}$ – IntegerInput

On entry: specifies the subinterval $\left[{x}_{i},{x}_{i+1}\right]$ on which the derivatives are to be evaluated.
fcn must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which
d02saf is called. Arguments denoted as
Input must
not be changed by this procedure.
Note: fcn should not return floatingpoint NaN (Not a Number) or infinity values, since these are not handled by
d02saf. If your code inadvertently
does return any NaNs or infinities,
d02saf is likely to produce unexpected results.
 16: $\mathbf{eqn}$ – Subroutine, supplied by the NAG Library or the user.External Procedure

eqn is used to describe the additional algebraic equations to be solved in the determination of the parameters,
${p}_{\mathit{i}}$, for
$\mathit{i}=1,2,\dots ,m$. If there are no additional algebraic equations (i.e.,
$m={n}_{1}$) then
eqn is never called and the dummy routine d02hbz should be used as the actual argument.
The specification of
eqn is:
Fortran Interface
Subroutine eqn ( 
e, q, p, m) 
Integer, Intent (In)  ::  q, m  Real (Kind=nag_wp), Intent (In)  ::  p(m)  Real (Kind=nag_wp), Intent (Out)  ::  e(q) 

C Header Interface
#include nagmk26.h
void 
eqn (double e[], const Integer *q, const double p[], const Integer *m) 

 1: $\mathbf{e}\left({\mathbf{q}}\right)$ – Real (Kind=nag_wp) arrayOutput

On exit: the vector of residuals, ${r}_{2}\left(p\right)$, that is the amount by which the current estimates of the arguments fail to satisfy the algebraic equations.
 2: $\mathbf{q}$ – IntegerInput

On entry: the number of algebraic equations, $m{n}_{1}$.
 3: $\mathbf{p}\left({\mathbf{m}}\right)$ – Real (Kind=nag_wp) arrayInput

On entry: the current estimate of the
$\mathit{i}$th argument ${p}_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,m$.
 4: $\mathbf{m}$ – IntegerInput

On entry: $m$, the number of arguments.
eqn must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which
d02saf is called. Arguments denoted as
Input must
not be changed by this procedure.
Note: eqn should not return floatingpoint NaN (Not a Number) or infinity values, since these are not handled by
d02saf. If your code inadvertently
does return any NaNs or infinities,
d02saf is likely to produce unexpected results.
 17: $\mathbf{constr}$ – Logical Function, supplied by the user.External Procedure

constr is used to prevent the pseudoNewton iteration running into difficulty.
constr should return the value .TRUE. if the constraints are satisfied by the parameters
${p}_{1},{p}_{2},\dots ,{p}_{m}$. Otherwise
constr should return the value .FALSE.. Usually the dummy function d02hby, which returns the value .TRUE. at all times, will suffice and in the first instance this is recommended as the actual argument.
The specification of
constr is:
Fortran Interface
Logical  ::  constr  Integer, Intent (In)  ::  m  Real (Kind=nag_wp), Intent (In)  ::  p(m) 

C Header Interface
#include nagmk26.h
Nag_Boolean 
constr (const double p[], const Integer *m) 

 1: $\mathbf{p}\left({\mathbf{m}}\right)$ – Real (Kind=nag_wp) arrayInput

On entry: an estimate of the
$\mathit{i}$th argument, ${p}_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,m$.
 2: $\mathbf{m}$ – IntegerInput

On entry: $m$, the number of arguments.
constr must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which
d02saf is called. Arguments denoted as
Input must
not be changed by this procedure.
 18: $\mathbf{ymax}$ – Real (Kind=nag_wp)Input

On entry: a nonnegative value which is used as a bound on all values ${\Vert y\left(x\right)\Vert}_{\infty}$ where $y\left(x\right)$ is the solution at any point $x$ between ${\mathbf{x}}\left(1\right)$ and ${\mathbf{x}}\left({\mathbf{npoint}}\right)$ for the current arguments ${p}_{1},{p}_{2},\dots ,{p}_{m}$. If this bound is exceeded the integration is terminated and the current arguments are rejected. Such a rejection will result in an error exit if it prevents the initial residual or Jacobian, or the final solution, being calculated. If ${\mathbf{ymax}}=0$ on entry, no bound on the solution $y$ is used; that is the integrations proceed without any checking on the size of ${\Vert y\Vert}_{\infty}$.
 19: $\mathbf{monit}$ – Subroutine, supplied by the NAG Library or the user.External Procedure

monit enables you to monitor the values of various quantities during the calculation.
It
is called by
d02saf after every calculation of the norm
${\Vert {{\mathbf{d}}}^{1}{\stackrel{~}{J}}^{+}r\Vert}_{2}$ which determines the strategy of the Newton method, every time there is an internal error exit leading to a change of strategy, and before an error exit when calculating the initial Jacobian.
Usually the routine d02hbx will be adequate and you are advised to use this as the actual argument for
monit in the first instance. (In this case a call to
x04abf must be made before the call of
d02saf.) If no monitoring is required, the dummy routine d02sas may be used.
The specification of
monit is:
Fortran Interface
Integer, Intent (In)  ::  istate, iflag, ifail1, m  Real (Kind=nag_wp), Intent (In)  ::  p(m), f(m), pnorm, pnorm1, eps, d(m) 

C Header Interface
#include nagmk26.h
void 
monit (const Integer *istate, const Integer *iflag, const Integer *ifail1, const double p[], const Integer *m, const double f[], const double *pnorm, const double *pnorm1, const double *eps, const double d[]) 

 1: $\mathbf{istate}$ – IntegerInput

On entry: the state of the Newton iteration.
 ${\mathbf{istate}}=0$
 The calculation of the residual, Jacobian and ${\Vert {{\mathbf{d}}}^{1}{\stackrel{~}{J}}^{+}r\Vert}_{2}$ are taking place.
 ${\mathbf{istate}}=1$ to $5$
 During the Newton iteration a factor of ${2}^{\left({\mathbf{istate}}+1\right)}$ of the Newton step is being used to try to reduce the norm.
 ${\mathbf{istate}}=6$
 The current Newton step has been rejected and the Jacobian is being recalculated.
 ${\mathbf{istate}}=6$ to $1$
 An internal error exit has caused the rejection of the current set of argument values, $p$. ${\mathbf{istate}}$ is the value which istate would have taken if the error had not occurred.
 ${\mathbf{istate}}=7$
 An internal error exit has occurred when calculating the initial Jacobian.
 2: $\mathbf{iflag}$ – IntegerInput

On entry: whether or not the Jacobian being used has been calculated at the beginning of the current iteration. If the Jacobian has been updated then ${\mathbf{iflag}}=1$; otherwise ${\mathbf{iflag}}=2$. The Jacobian is only calculated when convergence to the current argument values has been slow.
 3: $\mathbf{ifail1}$ – IntegerInput

On entry: if
$6\le {\mathbf{istate}}\le 1$,
ifail1 specifies the
ifail error number that would be produced were control returned to you.
ifail1 is unspecified for values of
istate outside this range.
 4: $\mathbf{p}\left({\mathbf{m}}\right)$ – Real (Kind=nag_wp) arrayInput

On entry: the current estimate of the
$\mathit{i}$th argument ${p}_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,m$.
 5: $\mathbf{m}$ – IntegerInput

On entry: $m$, the number of arguments.
 6: $\mathbf{f}\left({\mathbf{m}}\right)$ – Real (Kind=nag_wp) arrayInput

On entry:
$r$, the residual corresponding to the current argument values, provided
$1\le {\mathbf{istate}}\le 5$ or
${\mathbf{istate}}=7$.
f is unspecified for other values of
istate.
 7: $\mathbf{pnorm}$ – Real (Kind=nag_wp)Input

On entry: a quantity against which all reductions in norm are currently measured.
 8: $\mathbf{pnorm1}$ – Real (Kind=nag_wp)Input

On entry:
$p$, the norm of the current arguments. It is set for
$1\le {\mathbf{istate}}\le 5$ and is undefined for other values of
istate.
 9: $\mathbf{eps}$ – Real (Kind=nag_wp)Input

On entry: gives some indication of the convergence rate. It is the current singular value modification factor (see
Gay (1976)). It is zero initially and whenever convergence is proceeding steadily.
eps is
${\epsilon}^{3/8}$ or greater (where
$\epsilon $ may in most cases be considered
machine precision) when the singular values of
$J$ are approximately zero or when convergence is not being achieved. The larger the value of
eps the worse the convergence rate. When
eps becomes too large the Newton iteration is terminated.
 10: $\mathbf{d}\left({\mathbf{m}}\right)$ – Real (Kind=nag_wp) arrayInput

On entry:
$J$, the singular values of the current modified Jacobian matrix. If
${\mathbf{d}}\left(m\right)$ is small relative to
${\mathbf{d}}\left(1\right)$ for a number of Jacobians corresponding to different argument values then the computed results should be viewed with suspicion. It could be that the matching equations do not depend significantly on some argument (which could be due to a programming error in
fcn,
bc,
range or
eqn). Alternatively, the system of differential equations may be very illconditioned when viewed as an initial value problem, in which case
d02saf is unsuitable. This may also be indicated by some singular values being very large. These values of
${\mathbf{d}}\left(\mathit{i}\right)$, for
$\mathit{i}=1,2,\dots ,m$, should not be changed.
monit must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which
d02saf is called. Arguments denoted as
Input must
not be changed by this procedure.
Note: monit should not return floatingpoint NaN (Not a Number) or infinity values, since these are not handled by
d02saf. If your code inadvertently
does return any NaNs or infinities,
d02saf is likely to produce unexpected results.
 20: $\mathbf{prsol}$ – Subroutine, supplied by the NAG Library or the user.External Procedure

prsol can be used to obtain values of the solution
$y$ at a selected point
$z$ by integration across the final range
$\left[{\mathbf{x}}\left(1\right),{\mathbf{x}}\left({\mathbf{npoint}}\right)\right]$. If no output is required d02hbw can be used as the actual argument.
The specification of
prsol is:
Fortran Interface
Subroutine prsol ( 
z, y, n) 
Integer, Intent (In)  ::  n  Real (Kind=nag_wp), Intent (In)  ::  y(n)  Real (Kind=nag_wp), Intent (Inout)  ::  z 

C Header Interface
#include nagmk26.h
void 
prsol (double *z, const double y[], const Integer *n) 

 1: $\mathbf{z}$ – Real (Kind=nag_wp)Input/Output

On entry: contains
${x}_{1}$ on the first call. On subsequent calls
z contains its previous output value.
On exit: the next point at which output is required. The new point must be nearer
${\mathbf{x}}\left({\mathbf{npoint}}\right)$ than the old.
If
z is set to a point outside
$\left[{\mathbf{x}}\left(1\right),{\mathbf{x}}\left({\mathbf{npoint}}\right)\right]$ the process stops and control returns from
d02saf to the (sub)program from which
d02saf is called. Otherwise the next call to
prsol is made by
d02saf at the point
z, with solution values
${y}_{1},{y}_{2},\dots ,{y}_{n}$ at
z contained in
y. If
z is set to
${\mathbf{x}}\left({\mathbf{npoint}}\right)$ exactly, the final call to
prsol is made with
${y}_{1},{y}_{2},\dots ,{y}_{n}$ as values of the solution at
${\mathbf{x}}\left({\mathbf{npoint}}\right)$ produced by the integration. In general the solution values obtained at
${\mathbf{x}}\left({\mathbf{npoint}}\right)$ from
prsol will differ from the values obtained at this point by a call to
bc. The difference between the two solutions is the residual
$r$. You are reminded that the points
${\mathbf{x}}\left(1\right),{\mathbf{x}}\left(2\right),\dots ,{\mathbf{x}}\left({\mathbf{npoint}}\right)$ are available in the locations
${\mathbf{swp}}\left(1,4\right),{\mathbf{swp}}\left(2,4\right),\dots ,{\mathbf{swp}}\left({\mathbf{npoint}},4\right)$ at all times.
 2: $\mathbf{y}\left({\mathbf{n}}\right)$ – Real (Kind=nag_wp) arrayInput

On entry: the solution value
${y}_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,n$, at $z$.
 3: $\mathbf{n}$ – IntegerInput

On entry: $\mathit{n}$, the total number of differential equations.
prsol must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which
d02saf is called. Arguments denoted as
Input must
not be changed by this procedure.
Note: prsol should not return floatingpoint NaN (Not a Number) or infinity values, since these are not handled by
d02saf. If your code inadvertently
does return any NaNs or infinities,
d02saf is likely to produce unexpected results.
 21: $\mathbf{w}\left({\mathbf{ldw}},{\mathbf{sdw}}\right)$ – Real (Kind=nag_wp) arrayOutput

On exit: in the case of an error exit of the type where the point of failure is returned in
${\mathbf{swp}}\left(1,5\right)$, the solution at this point of failure is returned in
${\mathbf{w}}\left(\mathit{i},1\right)$, for
$\mathit{i}=1,2,\dots ,n$.
Otherwise
w is used for workspace.
 22: $\mathbf{ldw}$ – IntegerInput

On entry: the first dimension of the array
w as declared in the (sub)program from which
d02saf is called.
Constraint:
${\mathbf{ldw}}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left({\mathbf{n}},{\mathbf{m}}\right)$.
 23: $\mathbf{sdw}$ – IntegerInput

On entry: the second dimension of the array
w as declared in the (sub)program from which
d02saf is called.
Constraint:
${\mathbf{sdw}}\ge 3\times {\mathbf{m}}+12+\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(11,{\mathbf{m}}\right)$.
 24: $\mathbf{ifail}$ – IntegerInput/Output

On entry:
ifail must be set to
$0$,
$1\text{ or}1$. If you are unfamiliar with this argument you should refer to
Section 3.4 in How to Use the NAG Library and its Documentation for details.
For environments where it might be inappropriate to halt program execution when an error is detected, the value
$1\text{ or}1$ is recommended. If the output of error messages is undesirable, then the value
$1$ is recommended. Otherwise, if you are not familiar with this argument, the recommended value is
$0$.
When the value $\mathbf{1}\text{ or}\mathbf{1}$ is used it is essential to test the value of ifail on exit.
On exit:
${\mathbf{ifail}}={\mathbf{0}}$ unless the routine detects an error or a warning has been flagged (see
Section 6).
6
Error Indicators and Warnings
If on entry
${\mathbf{ifail}}=0$ or
$1$, explanatory error messages are output on the current error message unit (as defined by
x04aaf).
Errors or warnings detected by the routine:
 ${\mathbf{ifail}}=1$

One or more of the arguments
n,
n1,
m,
ldswp,
npoint,
icount,
ldw,
sdw,
e,
pe or
ymax is incorrectly set.
 ${\mathbf{ifail}}=2$

The constraints have been violated by the initial arguments.
 ${\mathbf{ifail}}=3$

The condition
${\mathbf{x}}\left(1\right)<{\mathbf{x}}\left(2\right)<\cdots <{\mathbf{x}}\left({\mathbf{npoint}}\right)$ (or
${\mathbf{x}}\left(1\right)>{\mathbf{x}}\left(2\right)>\cdots >{\mathbf{x}}\left({\mathbf{npoint}}\right)$) has been violated on a call to
range with the initial arguments.
 ${\mathbf{ifail}}=4$

In the integration from ${\mathbf{x}}\left(1\right)$ to ${\mathbf{x}}\left({\mathbf{npoint}}\right)$ with the initial or the final arguments, the step size was reduced too far for the integration to proceed. Consider reversing the order of the points ${\mathbf{x}}\left(1\right),{\mathbf{x}}\left(2\right),\dots ,{\mathbf{x}}\left({\mathbf{npoint}}\right)$. If this error exit still results, it is likely that d02saf is not a suitable method for solving the problem, or the initial choice of arguments is very poor, or the accuracy requirement specified by ${\mathbf{e}}\left(\mathit{i}\right)$, for $\mathit{i}=1,2,\dots ,n$, is too stringent.
 ${\mathbf{ifail}}=5$

In the integration from ${\mathbf{x}}\left(1\right)$ to ${\mathbf{x}}\left({\mathbf{npoint}}\right)$ with the initial or final arguments, an initial step could not be found to start the integration on one of the intervals ${\mathbf{x}}\left(i\right)$ to ${\mathbf{x}}\left(i+1\right)$. Consider reversing the order of the points. If this error exit still results it is likely that d02saf is not a suitable routine for solving the problem, or the initial choice of arguments is very poor, or the accuracy requirement specified by ${\mathbf{e}}\left(\mathit{i}\right)$, for $\mathit{i}=1,2,\dots ,n$, is much too stringent.
 ${\mathbf{ifail}}=6$

In the integration from
${\mathbf{x}}\left(1\right)$ to
${\mathbf{x}}\left({\mathbf{npoint}}\right)$ with the initial or final arguments, the solution exceeded
ymax in magnitude (when
${\mathbf{ymax}}>0.0$). It is likely that the initial choice of arguments was very poor or
ymax was incorrectly set.
Note: on an error with
${\mathbf{ifail}}={\mathbf{4}}$,
${\mathbf{5}}$ or
${\mathbf{6}}$ with the initial arguments, the interval in which failure occurs is contained in
${\mathbf{swp}}\left({\mathbf{npoint}},1\right)$. If a
monit similar to the one in
Section 10 is being used then it is a simple matter to distinguish between errors using the initial and final arguments. None of the error exits
${\mathbf{ifail}}={\mathbf{4}}$,
${\mathbf{5}}$ or
${\mathbf{6}}$ should occur on the
final integration (when computing the solution) as this integration has already been performed previously with exactly the same arguments
${p}_{\mathit{i}}$, for
$\mathit{i}=1,2,\dots ,m$. Seek expert help if this error occurs.
 ${\mathbf{ifail}}=7$

On calculating the initial approximation to the Jacobian, the constraints were violated.
 ${\mathbf{ifail}}=8$

On perturbing the arguments when calculating the initial approximation to the Jacobian, the condition ${\mathbf{x}}\left(1\right)<{\mathbf{x}}\left(2\right)<\cdots <{\mathbf{x}}\left({\mathbf{npoint}}\right)$ (or ${\mathbf{x}}\left(1\right)>{\mathbf{x}}\left(2\right)>\cdots >{\mathbf{x}}\left({\mathbf{npoint}}\right)$) is violated.
 ${\mathbf{ifail}}=9$

On calculating the initial approximation to the Jacobian, the integration step size was reduced too far to make further progress (see ${\mathbf{ifail}}={\mathbf{4}}$).
 ${\mathbf{ifail}}=10$

On calculating the initial approximation to the Jacobian, the initial integration step size on some interval was too small (see ${\mathbf{ifail}}={\mathbf{5}}$).
 ${\mathbf{ifail}}=11$

On calculating the initial approximation to the Jacobian, the solution of the system of differential equations exceeded
ymax in magnitude (when
${\mathbf{ymax}}>0.0$).
Note: all the error exits
${\mathbf{ifail}}={\mathbf{7}}$,
${\mathbf{8}}$,
${\mathbf{9}}$,
${\mathbf{10}}$ and
${\mathbf{11}}$ can be treated by reducing the size of some or all the elements of
dp.
 ${\mathbf{ifail}}=12$

On calculating the initial approximation to the Jacobian, a column of the Jacobian is found to be insignificant. This could be due to an element ${\mathbf{dp}}\left(i\right)$ being too small (but nonzero) or the solution having no dependence on one of the arguments (a programming error).
Note: on an error exit with ${\mathbf{ifail}}={\mathbf{7}}$, ${\mathbf{8}}$, ${\mathbf{9}}$, ${\mathbf{10}}$, ${\mathbf{11}}$ or ${\mathbf{12}}$, if a perturbation of the argument ${p}_{i}$ is the cause of the error then ${\mathbf{swp}}\left({\mathbf{npoint}},1\right)$ will contain the floatingpoint value of $i$.
 ${\mathbf{ifail}}=13$

After calculating the initial approximation to the Jacobian, the calculation of its singular value decomposition failed. It is likely that the error will never occur as it is usually associated with the Jacobian having multiple singular values. To remedy the error it should only be necessary to change the initial arguments. If the error persists it is likely that the problem has not been correctly formulated.
 ${\mathbf{ifail}}=14$

The Newton iteration has failed to converge after exercising all its options. You are strongly recommended to monitor the progress of the iteration via
monit. There are many possible reasons for the iteration not converging. Amongst the most likely are:
(a) 
there is no solution; 
(b) 
the initial arguments are too far away from the correct arguments; 
(c) 
the problem is too illconditioned as an initial value problem for Newton's method to choose suitable corrections; 
(d) 
the accuracy requirements for convergence are too restrictive, that is some of the components of pe (and maybe pf) are too small – in this case the final value of this norm output via monit will usually be very small; or 
(e) 
the initial arguments are so close to the solution arguments $p$ that the Newton iteration cannot find improved arguments. The norm output by monit should be very small. 
 ${\mathbf{ifail}}=15$

The number of iterations permitted by
icount has been exceeded (in the case when
${\mathbf{icount}}>0$ on entry).
 ${\mathbf{ifail}}=16$
 ${\mathbf{ifail}}=17$
 ${\mathbf{ifail}}=18$
 ${\mathbf{ifail}}=19$

These indicate that there has been a serious error in an internal call. Check all subroutine calls and array dimensions. Seek expert help.
 ${\mathbf{ifail}}=99$
An unexpected error has been triggered by this routine. Please
contact
NAG.
See
Section 3.9 in How to Use the NAG Library and its Documentation for further information.
 ${\mathbf{ifail}}=399$
Your licence key may have expired or may not have been installed correctly.
See
Section 3.8 in How to Use the NAG Library and its Documentation for further information.
 ${\mathbf{ifail}}=999$
Dynamic memory allocation failed.
See
Section 3.7 in How to Use the NAG Library and its Documentation for further information.
7
Accuracy
If the iteration converges, the accuracy to which the unknown arguments are determined is usually close to that specified by you. The accuracy of the solution (output via
prsol) depends on the error tolerances
${\mathbf{e}}\left(\mathit{i}\right)$, for
$\mathit{i}=1,2,\dots ,n$. You are strongly recommended to vary all tolerances to check the accuracy of the arguments
$p$ and the solution
$y$.
8
Parallelism and Performance
d02saf is not thread safe and should not be called from a multithreaded user program. Please see
Section 3.12.1 in How to Use the NAG Library and its Documentation for more information on thread safety.
d02saf makes calls to BLAS and/or LAPACK routines, which may be threaded within the vendor library used by this implementation. Consult the documentation for the vendor library for further information.
Please consult the
X06 Chapter Introduction for information on how to control and interrogate the OpenMP environment used within this routine. Please also consult the
Users' Note for your implementation for any additional implementationspecific information.
The time taken by
d02saf depends on the complexity of the system of differential equations and on the number of iterations required. In practice, the integration of the differential system
(1) is usually by far the most costly process involved. The computing time for integrating the differential equations can sometimes depend critically on the quality of the initial estimates for the arguments
$p$. If it seems that too much computing time is required and, in particular, if the values of the residuals (output in
monit) are much larger than expected given your knowledge of the expected solution, then the coding of
fcn,
eqn,
range and
bc should be checked for errors. If no errors can be found then an independent attempt should be made to improve the initial estimates
$p$.
In the case of an error exit in the integration of the differential system indicated by
${\mathbf{ifail}}={\mathbf{4}}$,
${\mathbf{5}}$,
${\mathbf{9}}$ or
${\mathbf{10}}$ you are strongly recommended to perform trial integrations with
d02pff to determine the effects of changes of the local error tolerances and of changes to the initial choice of the arguments
${p}_{\mathit{i}}$, for
$\mathit{i}=1,2,\dots ,m$, (that is the initial choice of
$p$).
It is possible that by following the advice given in
Section 6 an error exit with
${\mathbf{ifail}}={\mathbf{7}}$,
${\mathbf{8}}$,
${\mathbf{9}}$,
${\mathbf{10}}$ or
${\mathbf{11}}$ might be followed by one with
${\mathbf{ifail}}={\mathbf{12}}$ (or viceversa) where the advice given is the opposite. If you are unable to refine the choice of
${\mathbf{dp}}\left(\mathit{i}\right)$, for
$\mathit{i}=1,2,\dots ,n$, such that both these types of exits are avoided then the problem should be rescaled if possible or the method must be abandoned.
The choice of the ‘floor’ values
${\mathbf{pf}}\left(\mathit{i}\right)$, for $\mathit{i}=1,2,\dots ,m$, may be critical in the convergence of the Newton iteration. For each value $i$, the initial choice of ${p}_{i}$ and the choice of ${\mathbf{pf}}\left(i\right)$ should not both be very small unless it is expected that the final argument ${p}_{i}$ will be very small and that it should be determined accurately in a relative sense.
For many problems it is critical that a good initial estimate be found for the arguments
$p$ or the iteration will not converge or may even break down with an error exit. There are many mathematical techniques which obtain good initial estimates for
$p$ in simple cases but which may fail to produce useful estimates in harder cases. If no such technique is available it is recommended that you try a continuation (homotopy) technique preferably based on a physical argument (e.g., the Reynolds or Prandtl number is often a suitable continuation argument). In a continuation method a sequence of problems is solved, one for each choice of the continuation argument, starting with the problem of interest. At each stage the arguments
$p$ calculated at earlier stages are used to compute a good initial estimate for the arguments at the current stage (see
Hall and Watt (1976) for more details).
10
Example
This example intends to illustrate the use of the breakpoint and equation solving facilities of
d02saf. Most of the facilities which are common to
d02saf and
d02hbf are illustrated in the example in the specification of
d02hbf (which should also be consulted).
The program solves a projectile problem in two media determining the position of change of media, ${p}_{3}$, and the gravity and viscosity in the second medium (${p}_{2}$ represents gravity and ${p}_{4}$ represents viscosity).
10.1
Program Text
Program Text (d02safe.f90)
10.2
Program Data
Program Data (d02safe.d)
10.3
Program Results
Program Results (d02safe.r)