NAG FL Interface
d02qff (ivp_adams_roots)
1
Purpose
d02qff is a routine for integrating a nonstiff system of firstorder ordinary differential equations using a variableorder variablestep Adams' method. A rootfinding facility is provided.
2
Specification
Fortran Interface
Subroutine d02qff ( 
fcn, neqf, t, y, tout, g, neqg, root, rwork, lrwork, iwork, liwork, ifail) 
Integer, Intent (In) 
:: 
neqf, neqg, lrwork, liwork 
Integer, Intent (Inout) 
:: 
iwork(liwork), ifail 
Real (Kind=nag_wp), External 
:: 
g 
Real (Kind=nag_wp), Intent (In) 
:: 
tout 
Real (Kind=nag_wp), Intent (Inout) 
:: 
t, y(neqf), rwork(lrwork) 
Logical, Intent (Out) 
:: 
root 
External 
:: 
fcn 

C Header Interface
#include <nag.h>
void 
d02qff_ ( void (NAG_CALL *fcn)(const Integer *neqf, const double *x, const double y[], double f[]), const Integer *neqf, double *t, double y[], const double *tout, double (NAG_CALL *g)(const Integer *neqf, const double *x, const double y[], const double yp[], const Integer *k), const Integer *neqg, logical *root, double rwork[], const Integer *lrwork, Integer iwork[], const Integer *liwork, Integer *ifail) 

C++ Header Interface
#include <nag.h> extern "C" {
void 
d02qff_ ( void (NAG_CALL *fcn)(const Integer &neqf, const double &x, const double y[], double f[]), const Integer &neqf, double &t, double y[], const double &tout, double (NAG_CALL *g)(const Integer &neqf, const double &x, const double y[], const double yp[], const Integer &k), const Integer &neqg, logical &root, double rwork[], const Integer &lrwork, Integer iwork[], const Integer &liwork, Integer &ifail) 
}

The routine may be called by the names d02qff or nagf_ode_ivp_adams_roots.
3
Description
Given the initial values
$x,{y}_{1},{y}_{2},\dots ,{y}_{{\mathbf{neqf}}}$ d02qff integrates a nonstiff system of firstorder differential equations of the type
from
$x={\mathbf{t}}$ to
$x={\mathbf{tout}}$ using a variableorder variablestep Adams' method. The system is defined by
fcn, which evaluates
${f}_{i}$ in terms of
$x$ and
${y}_{1},{y}_{2},\dots ,{y}_{{\mathbf{neqf}}}$, and
${y}_{1},{y}_{2},\dots ,{y}_{{\mathbf{neqf}}}$ are supplied at
$x={\mathbf{t}}$. The routine is capable of finding roots (values of
$x$) of prescribed event functions of the form
(See
d02qwf for the specification of
neqg.)
Each
${g}_{j}$ is considered to be independent of the others so that roots are sought of each
${g}_{j}$ individually. The root reported by the routine will be the first root encountered by any
${g}_{j}$. Two techniques for determining the presence of a root in an integration step are available: the sophisticated method described in
Watts (1985) and a simplified method whereby sign changes in each
${g}_{j}$ are looked for at the ends of each integration step. The event functions are defined by
g supplied by you which evaluates
${g}_{j}$ in terms of
$x,{y}_{1},\dots ,{y}_{{\mathbf{neqf}}}$ and
${y}_{1}^{\prime},\dots ,{y}_{{\mathbf{neqf}}}^{\prime}$. In onestep mode the routine returns an approximation to the solution at each integration point. In interval mode this value is returned at the end of the integration range. If a root is detected this approximation is given at the root. You select the mode of operation, the error control, the rootfinding technique and various optional inputs by a prior call to the setup routine
d02qwf.
For a description of the practical implementation of an Adams' formula see
Shampine and Gordon (1975) and
Shampine and Watts (1979).
4
References
Shampine L F and Gordon M K (1975) Computer Solution of Ordinary Differential Equations – The Initial Value Problem W H Freeman & Co., San Francisco
Shampine L F and Watts H A (1979) DEPAC – design of a user oriented package of ODE solvers Report SAND792374 Sandia National Laboratory
Watts H A (1985) RDEAM – An Adams ODE code with root solving capability Report SAND851595 Sandia National Laboratory
5
Arguments

1:
$\mathbf{fcn}$ – Subroutine, supplied by the user.
External Procedure

fcn must evaluate the functions
${f}_{i}$ (that is the first derivatives
${y}_{i}^{\prime}$) for given values of its arguments
$x$,
${y}_{1},{y}_{2},\dots ,{y}_{{\mathbf{neqf}}}$.
The specification of
fcn is:
Fortran Interface
Integer, Intent (In) 
:: 
neqf 
Real (Kind=nag_wp), Intent (In) 
:: 
x, y(neqf) 
Real (Kind=nag_wp), Intent (Out) 
:: 
f(neqf) 

C Header Interface
void 
fcn_ (const Integer *neqf, const double *x, const double y[], double f[]) 

C++ Header Interface
#include <nag.h> extern "C" {
void 
fcn_ (const Integer &neqf, const double &x, const double y[], double f[]) 
}


1:
$\mathbf{neqf}$ – Integer
Input

On entry: the number of differential equations.

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

On entry: the current value of the argument $x$.

3:
$\mathbf{y}\left({\mathbf{neqf}}\right)$ – Real (Kind=nag_wp) array
Input

On entry: ${y}_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,{\mathbf{neqf}}$, the current value of the argument.

4:
$\mathbf{f}\left({\mathbf{neqf}}\right)$ – Real (Kind=nag_wp) array
Output

On exit: the value of
${f}_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,{\mathbf{neqf}}$.
fcn must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which
d02qff 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
d02qff. If your code inadvertently
does return any NaNs or infinities,
d02qff is likely to produce unexpected results.

2:
$\mathbf{neqf}$ – Integer
Input

On entry: the number of firstorder ordinary differential equations to be solved by
d02qff. It must contain the same value as the argument
neqf used in a prior call to
d02qwf.
Constraint:
${\mathbf{neqf}}\ge 1$.

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

On entry: after a call to
d02qwf with
${\mathbf{statef}}=\text{'S'}$ (i.e., an initial entry),
t must be set to the initial value of the independent variable
$x$.
On exit: the value of
$x$ at which
$y$ has been computed. This may be an intermediate output point, a root,
tout or a point at which an error has occurred. If the integration is to be continued, possibly with a new value for
tout,
t must not be changed.

4:
$\mathbf{y}\left({\mathbf{neqf}}\right)$ – Real (Kind=nag_wp) array
Input/Output

On entry: the initial values of the solution ${y}_{1},{y}_{2},\dots ,{y}_{{\mathbf{neqf}}}$.
On exit: the computed values of the solution at the exit value of
t. If the integration is to be continued, possibly with a new value for
tout, these values must not be changed.

5:
$\mathbf{tout}$ – Real (Kind=nag_wp)
Input

On entry: the next value of
$x$ at which a computed solution is required. For the initial
t, the input value of
tout is used to determine the direction of integration. Integration is permitted in either direction. If
${\mathbf{tout}}={\mathbf{t}}$ on exit,
tout must be reset beyond
t in the direction of integration, before any continuation call.

6:
$\mathbf{g}$ – real (Kind=nag_wp) Function, supplied by the user.
External Procedure

g must evaluate a given component of
$g\left(x,y,{y}^{\prime}\right)$ at a specified point.
If rootfinding is not required the actual argument for
g must be the dummy routine
d02qfz. (
d02qfz is included in the NAG Library.)
The specification of
g is:
Fortran Interface
Real (Kind=nag_wp) 
:: 
g 
Integer, Intent (In) 
:: 
neqf, k 
Real (Kind=nag_wp), Intent (In) 
:: 
x, y(neqf), yp(neqf) 

C Header Interface
double 
g_ (const Integer *neqf, const double *x, const double y[], const double yp[], const Integer *k) 

C++ Header Interface
#include <nag.h> extern "C" {
double 
g_ (const Integer &neqf, const double &x, const double y[], const double yp[], const Integer &k) 
}


1:
$\mathbf{neqf}$ – Integer
Input

On entry: the number of differential equations being solved.

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

On entry: the current value of the independent variable.

3:
$\mathbf{y}\left({\mathbf{neqf}}\right)$ – Real (Kind=nag_wp) array
Input

On entry: the current values of the dependent variables.

4:
$\mathbf{yp}\left({\mathbf{neqf}}\right)$ – Real (Kind=nag_wp) array
Input

On entry: the current values of the derivatives of the dependent variables.

5:
$\mathbf{k}$ – Integer
Input

On entry: the component of $g$ which must be evaluated.
g must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which
d02qff is called. Arguments denoted as
Input must
not be changed by this procedure.
Note: g should not return floatingpoint NaN (Not a Number) or infinity values, since these are not handled by
d02qff. If your code inadvertently
does return any NaNs or infinities,
d02qff is likely to produce unexpected results.

7:
$\mathbf{neqg}$ – Integer
Input

On entry: the number of event functions which you are defining for rootfinding. If rootfinding is not required the value for
neqg must be
$\text{}\le 0$. Otherwise it must be the same argument
neqg used in the prior call to
d02qwf.

8:
$\mathbf{root}$ – Logical
Output

On exit: if rootfinding was required (
${\mathbf{neqg}}>0$ on entry),
root specifies whether or not the output value of the argument
t is a root of one of the event functions. If
${\mathbf{root}}=\mathrm{.FALSE.}$, no root was detected, whereas
${\mathbf{root}}=\mathrm{.TRUE.}$ indicates a root and you should make a call to
d02qyf for further information.
If rootfinding was not required (${\mathbf{neqg}}=0$ on entry), then on exit ${\mathbf{root}}=\mathrm{.FALSE.}$.

9:
$\mathbf{rwork}\left({\mathbf{lrwork}}\right)$ – Real (Kind=nag_wp) array
Communication Array

This
must be the same argument
rwork as supplied to
d02qwf. It is used to pass information from
d02qwf to
d02qff, and from
d02qff to
d02qxf,
d02qyf and
d02qzf. Therefore the contents of this array
must not be changed before the call to
d02qff or calling any of the routines
d02qxf,
d02qyf and
d02qzf.

10:
$\mathbf{lrwork}$ – Integer
Input

On entry: the dimension of the array
rwork as declared in the (sub)program from which
d02qff is called.
This must be the same argument
lrwork as supplied to
d02qwf.

11:
$\mathbf{iwork}\left({\mathbf{liwork}}\right)$ – Integer array
Communication Array

This
must be the same argument
iwork as supplied to
d02qwf. It is used to pass information from
d02qwf to
d02qff, and from
d02qff to
d02qxf,
d02qyf and
d02qzf. Therefore the contents of this array
must not be changed before the call to
d02qff or calling any of the routines
d02qxf,
d02qyf and
d02qzf.

12:
$\mathbf{liwork}$ – Integer
Input

On entry: the dimension of the array
iwork as declared in the (sub)program from which
d02qff is called.
This must be the same argument
liwork as supplied to
d02qwf.

13:
$\mathbf{ifail}$ – Integer
Input/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 4 in the Introduction to the NAG Library FL Interface 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, because for this routine the values of the output arguments may be useful even if
${\mathbf{ifail}}\ne {\mathbf{0}}$ on exit, the recommended value is
$1$.
When the value $\mathbf{1}\text{or}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$

On entry,
${\mathbf{liwork}}=\u2329\mathit{\text{value}}\u232a$ and
${\mathbf{liwork}}=\u2329\mathit{\text{value}}\u232a$ in
d02qwf.
Constraint:
${\mathbf{liwork}}={\mathbf{liwork}}$ in
d02qwf.
On entry,
${\mathbf{lrwork}}=\u2329\mathit{\text{value}}\u232a$ and
${\mathbf{lrwork}}=\u2329\mathit{\text{value}}\u232a$ in
d02qwf.
Constraint:
${\mathbf{lrwork}}={\mathbf{lrwork}}$ in
d02qwf.
On entry,
${\mathbf{neqf}}=\u2329\mathit{\text{value}}\u232a$ and
${\mathbf{neqf}}=\u2329\mathit{\text{value}}\u232a$ in
d02qwf.
Constraint:
${\mathbf{neqf}}={\mathbf{neqf}}$ in
d02qwf.
On entry,
${\mathbf{neqg}}=\u2329\mathit{\text{value}}\u232a$ and
${\mathbf{neqg}}=0$ in
d02qwf.
Constraint: if
${\mathbf{neqg}}>0$ then
${\mathbf{neqg}}>0$ in
d02qwf.
On entry,
${\mathbf{neqg}}=\u2329\mathit{\text{value}}\u232a$ and
${\mathbf{neqg}}=0$ in
d02qwf.
Constraint: if
${\mathbf{neqg}}\le 0$ then
${\mathbf{neqg}}={\mathbf{neqg}}$ in
d02qwf.
On entry,
${\mathbf{neqg}}=\u2329\mathit{\text{value}}\u232a$ and
${\mathbf{neqg}}=\u2329\mathit{\text{value}}\u232a$ in
d02qwf.
Constraint: if
${\mathbf{neqg}}\le 0$ then
${\mathbf{neqg}}\le 0$ in
d02qwf.
On entry,
${\mathbf{tout}}=\u2329\mathit{\text{value}}\u232a$ but
${\mathbf{crit}}=\mathrm{.TRUE.}$ in
d02qwf.
Integration cannot be attempted beyond
${\mathbf{tcrit}}=\u2329\mathit{\text{value}}\u232a$.
On entry, ${\mathbf{tout}}={\mathbf{t}}=\u2329\mathit{\text{value}}\u232a$.
The call to setup routine
d02qwf produced an error.
The setup routine
d02qwf has not been called.
The value of
t has been changed from
$\u2329\mathit{\text{value}}\u232a$ to
$\u2329\mathit{\text{value}}\u232a$.
This is not permitted on a continuation call.
The value of
tout,
$\u2329\mathit{\text{value}}\u232a$, indicates a change in the integration direction. This is not permitted on a continuation call.
 ${\mathbf{ifail}}=2$

The maximum number of steps has been attempted.
If integration is to be continued then the routine may be called again and a further
maxstp in
d02qwf steps will be attempted.
 ${\mathbf{ifail}}=3$

The error tolerances are too stringent.
 ${\mathbf{ifail}}=4$

Error weight
$i$ has become zero during the integration.
${\mathbf{atol}}\left(i\right)=0.0$ in
d02qwf but
${\mathbf{y}}\left(i\right)$ is now
$0.0$. Integration successful as far as
t.
 ${\mathbf{ifail}}=5$

The problem appears to be stiff.
(See the
D02 Chapter Introduction for a discussion of the term ‘stiff’). Although it is inefficient to use this integrator to solve stiff problems, integration may be continued by resetting
ifail and calling the routine again.
 ${\mathbf{ifail}}=6$

A change in sign of an event function has been detected but the rootfinding process appears to have converged to a singular point of
t rather than a root. Integration may be continued by calling the routine again.
 ${\mathbf{ifail}}=7$

Two successive errors detected at the current value of
t,
$\u2329\mathit{\text{value}}\u232a$.
 ${\mathbf{ifail}}=99$
An unexpected error has been triggered by this routine. Please
contact
NAG.
See
Section 7 in the Introduction to the NAG Library FL Interface for further information.
 ${\mathbf{ifail}}=399$
Your licence key may have expired or may not have been installed correctly.
See
Section 8 in the Introduction to the NAG Library FL Interface for further information.
 ${\mathbf{ifail}}=999$
Dynamic memory allocation failed.
See
Section 9 in the Introduction to the NAG Library FL Interface for further information.
7
Accuracy
The accuracy of integration is determined by the arguments
vectol,
rtol and
atol in a prior call to
d02qwf. Note that only the local error at each step is controlled by these arguments. The error estimates obtained are not strict bounds but are usually reliable over one step. Over a number of steps the overall error may accumulate in various ways, depending on the properties of the differential equation system. The code is designed so that a reduction in the tolerances should lead to an approximately proportional reduction in the error. You are strongly recommended to call
d02qff with more than one set of tolerances and to compare the results obtained to estimate their accuracy.
The accuracy obtained depends on the type of error test used. If the solution oscillates around zero a relative error test should be avoided, whereas if the solution is exponentially increasing an absolute error test should not be used. If different accuracies are required for different components of the solution then a componentwise error test should be used. For a description of the error test see the specifications of the arguments
vectol,
atol and
rtol in the routine document for
d02qwf.
The accuracy of any roots located will depend on the accuracy of integration and may also be restricted by the numerical properties of $g\left(x,y,{y}^{\prime}\right)$. When evaluating $g$ you should try to write the code so that unnecessary cancellation errors will be avoided.
8
Parallelism and Performance
d02qff is not thread safe and should not be called from a multithreaded user program. Please see
Section 1 in FL Interface Multithreading for more information on thread safety.
d02qff is not threaded in any implementation.
If
d02qff fails with
${\mathbf{ifail}}={\mathbf{3}}$ then the combination of
atol and
rtol may be so small that a solution cannot be obtained, in which case the routine should be called again with larger values for
rtol and/or
atol (see
d02qwf). If the accuracy requested is really needed then you should consider whether there is a more fundamental difficulty. For example:

(a)in the region of a singularity the solution components will usually be of a large magnitude. d02qff could be used in onestep mode to monitor the size of the solution with the aim of trapping the solution before the singularity. In any case numerical integration cannot be continued through a singularity, and analytical treatment may be necessary;

(b)for ‘stiff’ equations, where the solution contains rapidly decaying components, the routine will require a very small step size to preserve stability. This will usually be exhibited by excessive computing time and sometimes an error exit with ${\mathbf{ifail}}={\mathbf{3}}$, but usually an error exit with ${\mathbf{ifail}}={\mathbf{2}}$ or ${\mathbf{5}}$. The Adams' methods are not efficient in such cases and you should consider using a routine from the Subchapter D02MN. A high proportion of failed steps
(see argument nfail) may indicate stiffness but there may be other reasons for this phenomenon.
d02qff can be used for producing results at short intervals (for example, for graph plotting); you should set
${\mathbf{crit}}=\mathrm{.TRUE.}$ and
tcrit to the last output point required in a prior call to
d02qwf and then set
tout appropriately for each output point in turn in the call to
d02qff.
10
Example
This example solves the equation
reposed as
over the range
$\left[0,10.0\right]$ with initial conditions
${y}_{1}=0.0$ and
${y}_{2}=1.0$ using vector error control (
${\mathbf{vectol}}=\mathrm{.TRUE.}$) and computation of the solution at
${\mathbf{tout}}=10.0$ with
${\mathbf{tcrit}}=10.0$ (
${\mathbf{crit}}=\mathrm{.TRUE.}$). Also, we use
d02qff to locate the positions where
${y}_{1}=0.0$ or where the first component has a turningpoint, that is
${y}_{1}^{\prime}=0.0$.
10.1
Program Text
10.2
Program Data
10.3
Program Results