NAG Library Routine Document
d02nmf (ivp_stiff_exp_revcom)
1
Purpose
d02nmf is a reverse communication routine for integrating stiff systems of explicit ordinary differential equations.
2
Specification
Fortran Interface
Subroutine d02nmf ( 
neq, ldysav, t, tout, y, ydot, rwork, rtol, atol, itol, inform, ysav, sdysav, wkjac, nwkjac, jacpvt, njcpvt, imon, inln, ires, irevcm, itask, itrace, ifail) 
Integer, Intent (In)  ::  neq, ldysav, itol, sdysav, nwkjac, njcpvt, itask, itrace  Integer, Intent (Inout)  ::  inform(23), jacpvt(njcpvt), imon, inln, ires, irevcm, ifail  Real (Kind=nag_wp), Intent (In)  ::  tout, rtol(*), atol(*)  Real (Kind=nag_wp), Intent (Inout)  ::  t, y(neq), ydot(neq), rwork(50+4*neq), ysav(ldysav,sdysav), wkjac(nwkjac) 

C Header Interface
#include <nagmk26.h>
void 
d02nmf_ (const Integer *neq, const Integer *ldysav, double *t, const double *tout, double y[], double ydot[], double rwork[], const double rtol[], const double atol[], const Integer *itol, Integer inform[], double ysav[], const Integer *sdysav, double wkjac[], const Integer *nwkjac, Integer jacpvt[], const Integer *njcpvt, Integer *imon, Integer *inln, Integer *ires, Integer *irevcm, const Integer *itask, const Integer *itrace, Integer *ifail) 

3
Description
d02nmf is a general purpose routine for integrating the initial value problem for a stiff system of explicit ordinary differential equations,
An outline of a typical calling program is given below:
! Declarations
Call linear algebra setup routine
Call integrator setup routine
irevcm=0
1000 Call d02nmf(neq, ldysav, t, tout, y, ydot, rwork, rtol, &
atol, itol, inform, ysave, sdysav, wkjac, nwkjac, &
jacpvt, njcpvt, imon, inln, ires, irevcm, itask, &
itrace, ifail)
If (irevcm.gt.0) Then
If (irevcm. eq. 8) Then
supply the Jacobian matrix (i)
Else If(irevcm.eq.9) Then
perform monitoring tasks requested by the user (ii)
Else If(irecvm.eq.1.or.irevcm.ge.3.and.irevcm.le.5) Then
evaluate the derivative (iii)
Else If(irevcm.eq.10) Then
indicates an unsuccessful step
Endif
Go To 1000
Endif
! post processing (optional linear algebra diagnostic call
! (sparse case only), optional integrator diagnostic call)
Stop
End
There are three major operations that may be required of the calling (sub)program on an intermeditate return (${\mathbf{irevcm}}\ne 0$) from d02nmf; these are denoted (i), (ii) and (iii) above.
The following sections describe in greater detail exactly what is required of each of these operations.
(i) Supply the Jacobian Matrix
You need only provide this facility if the argument
${\mathbf{jceval}}=\text{'A'}$ (or
${\mathbf{jceval}}=\text{'F'}$ if using sparse matrix linear algebra) in a call to the linear algebra setup routine (see
jceval in
d02nsf). If the Jacobian matrix is to be evaluated numerically by the integrator, then the remainder of section (i) can be ignored.
We must define the system of nonlinear equations which is solved internally by the integrator. The time derivative,
${y}^{\prime}$, has the form
where
$h$ is the current step size and
$d$ is an argument that depends on the integration method in use. The vector
$y$ is the current solution and the vector
$z$ depends on information from previous time steps. This means that
$\frac{d}{d{y}^{\prime}}\left(\text{}\right)=\left(hd\right)\frac{d}{dy}\left(\text{}\right)$.
The system of nonlinear equations that is solved has the form
but is solved in the form
where the function
$r$ is defined by
It is the Jacobian matrix
$\frac{\partial r}{\partial y}$ that you must supply as follows:
where
$t$,
$h$ and
$d$ are located in
${\mathbf{rwork}}\left(19\right)$,
${\mathbf{rwork}}\left(16\right)$ and
${\mathbf{rwork}}\left(20\right)$ respectively and the array
y contains the current values of the dependent variables. Only the nonzero elements of the Jacobian need be set, since the locations where it is to be stored are preset to zero.
Hereafter in this document this operation will be referred to as JAC.
(ii) Perform Tasks Requested by You
This operation is essentially a monitoring function and additionally provides the opportunity of changing the current values of
y, HNEXT (the step size that the integrator proposes to take on the next step), HMIN (the minimum step size to be taken on the next step), and HMAX (the maximum step size to be taken on the next step). The scaled local error at the end of a timestep may be obtained by calling real function
d02zaf as follows:
ifail = 1
errloc = d02zaf(neq,rwork(51+neq),rwork(51),ifail)
! CHECK IFAIL BEFORE PROCEEDING
The following gives details of the location within the array
rwork of variables that may be of interest to you:
Variable 
Specification 
Location 
tcurr 
the current value of the independent variable 
${\mathbf{rwork}}\left(19\right)$ 
hlast 
last step size successfully used by the integrator 
${\mathbf{rwork}}\left(15\right)$ 
hnext 
step size that the integrator proposes to take on the next step 
${\mathbf{rwork}}\left(16\right)$ 
hmin 
minimum step size to be taken on the next step 
${\mathbf{rwork}}\left(17\right)$ 
hmax 
maximum step size to be taken on the next step 
${\mathbf{rwork}}\left(18\right)$ 
nqu 
the order of the integrator used on the last step 
${\mathbf{rwork}}\left(10\right)$ 
You are advised to consult the description of
monitr in
d02nbf for details on what optional input can be made.
If
y is changed, then
imon must be set to
$2$ before return to
d02nmf. If either of the values of HMIN or HMAX are changed, then
imon must be set
$\text{}\ge 3$ before return to
d02nmf. If HNEXT is changed, then
imon must be set to
$4$ before return to
d02nmf.
In addition you can force
d02nmf to evaluate the residual vector
by setting
${\mathbf{imon}}=0$ and
${\mathbf{inln}}=3$ and then returning to
d02nmf; on return to this monitoring operation the residual vector will be stored in
${\mathbf{rwork}}\left(50+2\times {\mathbf{neq}}+\mathit{i}\right)$, for
$\mathit{i}=1,2,\dots ,{\mathbf{neq}}$.
Hereafter in this document this operation will be referred to as MONITR.
(iii) Evaluate the Derivative
This operation must evaluate the derivative vector for the explicit ordinary differential equation system defined by
where
$t$ is located in
${\mathbf{rwork}}\left(19\right)$.
Hereafter in this document this operation will be referred to as FCN.
4
References
5
Arguments
Note: this routine uses
reverse communication. Its use involves an initial entry, intermediate exits and reentries, and a final exit, as indicated by the argument
irevcm. Between intermediate exits and reentries,
all arguments other than ydot, rwork, wkjac, imon, inln and ires must remain unchanged.
 1: $\mathbf{neq}$ – IntegerInput

On initial entry: the number of differential equations to be solved.
Constraint:
${\mathbf{neq}}\ge 1$.
 2: $\mathbf{ldysav}$ – IntegerInput

On initial entry: an upper bound on the maximum number of differential equations to be solved during the integration.
Constraint:
${\mathbf{ldysav}}\ge {\mathbf{neq}}$.
 3: $\mathbf{t}$ – Real (Kind=nag_wp)Input/Output

On initial entry:
$t$, the value of the independent variable. The input value of
t is used only on the first call as the initial point of the integration.
On final exit: the value at which the computed solution
$y$ is returned (usually at
tout).
 4: $\mathbf{tout}$ – Real (Kind=nag_wp)Input

On initial entry: the next value of
$t$ at which a computed solution is desired. For the initial
$t$, the input value of
tout is used to determine the direction of integration. Integration is permitted in either direction (see also
itask).
Constraint:
${\mathbf{tout}}\ne {\mathbf{t}}$.
 5: $\mathbf{y}\left({\mathbf{neq}}\right)$ – Real (Kind=nag_wp) arrayInput/Output

On initial entry: the values of the dependent variables (solution). On the first call the first
neq elements of
$y$ must contain the vector of initial values.
On final exit: the computed solution vector evaluated at
t (usually
$t={\mathbf{tout}}$).
 6: $\mathbf{ydot}\left({\mathbf{neq}}\right)$ – Real (Kind=nag_wp) arrayInput/Output

On intermediate reentry: must be set to the derivatives as defined under the description of
irevcm.
On final exit: the time derivatives ${y}^{\prime}$ of the vector $y$ at the last integration point.
 7: $\mathbf{rwork}\left(50+4\times {\mathbf{neq}}\right)$ – Real (Kind=nag_wp) arrayCommunication Array

On initial entry: must be the same array as used by one of the method setup routines
d02mvf,
d02nvf or
d02nwf, and by one of the storage setup routines
d02ntf,
d02nuf or
d02nvf. The contents of
rwork must not be changed between any call to a setup routine and the first call to
d02nmf.
On intermediate reentry: elements of
rwork must be set to quantities as defined under the description of
irevcm.
On intermediate exit:
contains information for JAC, FCN and MONITR operations as described in
Section 3 and the argument
irevcm.
 8: $\mathbf{rtol}\left(*\right)$ – Real (Kind=nag_wp) arrayInput

Note: the dimension of the array
rtol
must be at least
$1$ if
${\mathbf{itol}}=1$ or
$2$, and at least
${\mathbf{neq}}$ otherwise.
On initial entry: the relative local error tolerance.
Constraint:
${\mathbf{rtol}}\left(i\right)\ge 0.0$ for all relevant
$i$ (see
itol).
 9: $\mathbf{atol}\left(*\right)$ – Real (Kind=nag_wp) arrayInput

Note: the dimension of the array
atol
must be at least
$1$ if
${\mathbf{itol}}=1$ or
$3$, and at least
${\mathbf{neq}}$ otherwise.
On initial entry: the absolute local error tolerance.
Constraint:
${\mathbf{atol}}\left(i\right)\ge 0.0$ for all relevant
$i$ (see
itol).
 10: $\mathbf{itol}$ – IntegerInput

On initial entry: a value to indicate the form of the local error test.
itol indicates to
d02nmf whether to interpret either or both of
rtol or
atol as a vector or a scalar. The error test to be satisfied is
$\Vert {e}_{i}/{w}_{i}\Vert <1.0$, where
${w}_{i}$ is defined as follows:
itol  rtol  atol  ${w}_{i}$ 
1  scalar  scalar  ${\mathbf{rtol}}\left(1\right)\times \left{y}_{i}\right+{\mathbf{atol}}\left(1\right)$ 
2  scalar  vector  ${\mathbf{rtol}}\left(1\right)\times \left{y}_{i}\right+{\mathbf{atol}}\left(i\right)$ 
3  vector  scalar  ${\mathbf{rtol}}\left(i\right)\times \left{y}_{i}\right+{\mathbf{atol}}\left(1\right)$ 
4  vector  vector  ${\mathbf{rtol}}\left(i\right)\times \left{y}_{i}\right+{\mathbf{atol}}\left(i\right)$ 
${e}_{i}$ is an estimate of the local error in ${y}_{i}$, computed internally, and the choice of norm to be used is defined by a previous call to an integrator setup routine.
Constraint:
${\mathbf{itol}}=1$, $2$, $3$ or $4$.
 11: $\mathbf{inform}\left(23\right)$ – Integer arrayCommunication Array

 12: $\mathbf{ysav}\left({\mathbf{ldysav}},{\mathbf{sdysav}}\right)$ – Real (Kind=nag_wp) arrayCommunication Array
 13: $\mathbf{sdysav}$ – IntegerInput

On initial entry: the second dimension of the array
ysav as declared in the (sub)program from which
d02nmf is called. An appropriate value for
sdysav is described in the specifications of the integrator setup routines
d02nvf and
d02nwf. This value must be the same as that supplied to the integrator setup routine.
 14: $\mathbf{wkjac}\left({\mathbf{nwkjac}}\right)$ – Real (Kind=nag_wp) arrayInput/Output

On intermediate reentry: elements of the Jacobian as defined under the description of
irevcm. If a numerical Jacobian was requested then
wkjac is used for workspace.
On intermediate exit:
the Jacobian is overwritten.
 15: $\mathbf{nwkjac}$ – IntegerInput

On initial entry: the dimension of the array
wkjac as declared in the (sub)program from which
d02nmf is called. The actual size depends on the linear algebra method used. An appropriate value for
nwkjac is described in the specifications of the linear algebra setup routines
d02nsf,
d02ntf and
d02nuf for full, banded and sparse matrix linear algebra respectively. This value must be the same as that supplied to the linear algebra setup routine.
 16: $\mathbf{jacpvt}\left({\mathbf{njcpvt}}\right)$ – Integer arrayCommunication Array
 17: $\mathbf{njcpvt}$ – IntegerInput

On initial entry: the dimension of the array
jacpvt as declared in the (sub)program from which
d02nmf is called. The actual size depends on the linear algebra method used. An appropriate value for
njcpvt is described in the specifications of the linear algebra setup routines
d02ntf and
d02nuf for banded and sparse matrix linear algebra respectively. This value must be the same as that supplied to the linear algebra setup routine. When full matrix linear algebra is chosen, the array
jacpvt is not used and hence
njcpvt should be set to
$1$.
 18: $\mathbf{imon}$ – IntegerInput/Output

On intermediate exit:
used to pass information between
d02nmf and the MONITR operation (see
Section 3). With
${\mathbf{irevcm}}=9$,
imon contains a flag indicating under what circumstances the return from
d02nmf occurred:
 ${\mathbf{imon}}=2$
 Exit from d02nmf after ${\mathbf{ires}}=4$ caused an early termination (this facility could be used to locate discontinuities).
 ${\mathbf{imon}}=1$
 The current step failed repeatedly.
 ${\mathbf{imon}}=0$
 Exit from d02nmf after a call to the internal nonlinear equation solver.
 ${\mathbf{imon}}=1$
 The current step was successful.
On intermediate reentry: may be reset to determine subsequent action in
d02nmf.
 ${\mathbf{imon}}=2$
 Integration is to be halted. A return will be made from d02nmf to the calling (sub)program with ${\mathbf{ifail}}={\mathbf{12}}$.
 ${\mathbf{imon}}=1$
 Allow d02nmf to continue with its own internal strategy. The integrator will try up to three restarts unless ${\mathbf{imon}}\ne 1$.
 ${\mathbf{imon}}=0$
 Return to the internal nonlinear equation solver, where the action taken is determined by the value of inln.
 ${\mathbf{imon}}=1$
 Normal exit to d02nmf to continue integration.
 ${\mathbf{imon}}=2$
 Restart the integration at the current time point. The integrator will restart from order $1$ when this option is used. The solution y, provided by the monitr operation (see Section 3), will be used for the initial conditions.
 ${\mathbf{imon}}=3$
 Try to continue with the same step size and order as was to be used before entering the monitr operation (see Section 3). hmin and hmax may be altered if desired.
 ${\mathbf{imon}}=4$
 Continue the integration but using a new value of hnext and possibly new values of hmin and hmax.
 19: $\mathbf{inln}$ – IntegerInput/Output

On intermediate reentry: with
${\mathbf{imon}}=0$ and
${\mathbf{irevcm}}=9$,
inln specifies the action to be taken by the internal nonlinear equation solver. By setting
${\mathbf{inln}}=3$ and returning to
d02nmf, the residual vector is evaluated and placed in
${\mathbf{rwork}}\left(50+2\times {\mathbf{neq}}+\mathit{i}\right)$, for
$\mathit{i}=1,2,\dots ,{\mathbf{neq}}$ and then the
monitr operation (see
Section 3) is invoked again. At present this is the only option available:
inln must not be set to any other value.
On intermediate exit:
contains a flag indicating the action to be taken, if any, by the internal nonlinear equation solver.
 20: $\mathbf{ires}$ – IntegerInput/Output

On intermediate exit:
with
${\mathbf{irevcm}}=1$,
$2$,
$3$,
$4$ or
$5$,
ires contains the value
$1$.
On intermediate reentry: should be unchanged unless one of the following actions is required of
d02nmf in which case
ires should be set accordingly.
 ${\mathbf{ires}}=2$
 Indicates to d02nmf that control should be passed back immediately to the calling (sub)program with the error indicator set to ${\mathbf{ifail}}={\mathbf{11}}$.
 ${\mathbf{ires}}=3$
 Indicates to d02nmf that an error condition has occurred in the solution vector, its time derivative or in the value of $t$. The integrator will use a smaller time step to try to avoid this condition. If this is not possible d02nmf returns to the calling (sub)program with the error indicator set to ${\mathbf{ifail}}={\mathbf{7}}$.
 ${\mathbf{ires}}=4$
 Indicates to d02nmf to stop its current operation and to enter the monitr operation (see Section 3) immediately.
 21: $\mathbf{irevcm}$ – IntegerInput/Output

On initial entry: must contain $0$.
On intermediate reentry: should remain unchanged.
On intermediate exit:
indicates what action you must take before reentering. The possible exit values of
irevcm are
$1$,
$3$,
$4$,
$5$,
$8$,
$9$ or
$10$, which should be interpreted as follows:
 ${\mathbf{irevcm}}=1$, $3$, $4$ and $5$
 Indicates that an FCN operation (see Section 3) is required: ${y}^{\prime}=g\left(t,y\right)$ must be supplied, where
${\mathbf{y}}\left(\mathit{i}\right)$ is located in ${y}_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,{\mathbf{neq}}$.
For ${\mathbf{irevcm}}=1$ or $3$,
${y}_{\mathit{i}}^{\prime}$ should be placed in location ${\mathbf{rwork}}\left(50+2\times {\mathbf{neq}}+\mathit{i}\right)$, for $\mathit{i}=1,2,\dots ,{\mathbf{neq}}$.
For ${\mathbf{irevcm}}=4$,
${y}_{\mathit{i}}^{\prime}$ should be placed in location ${\mathbf{rwork}}\left(50+{\mathbf{neq}}+\mathit{i}\right)$, for $\mathit{i}=1,2,\dots ,{\mathbf{neq}}$.
For ${\mathbf{irevcm}}=5$,
${y}_{\mathit{i}}^{\prime}$ should be placed in location ${\mathbf{ydot}}\left(\mathit{i}\right)$, for $\mathit{i}=1,2,\dots ,{\mathbf{neq}}$.
 ${\mathbf{irevcm}}=8$
 Indicates that a JAC operation (see Section 3) is required: the Jacobian matrix must be supplied.
If full matrix linear algebra is being used, the $\left(i,j\right)$th element of the Jacobian must be stored in ${\mathbf{wkjac}}\left(\left(j1\right)\times {\mathbf{neq}}+i\right)$.
If banded matrix linear algebra is being used then the $\left(i,j\right)$th element of the Jacobian
must be stored in ${\mathbf{wkjac}}\left(\left(i1\right)\times {m}_{B}+k\right)$, where ${m}_{B}={m}_{L}+{m}_{U}+1$ and
$k=\mathrm{min}\phantom{\rule{0.125em}{0ex}}\left({m}_{L}i+1,0\right)+j$; here ${m}_{L}$ and ${m}_{U}$ are the number of subdiagonals and superdiagonals, respectively, in the band.
If sparse matrix linear algebra is being used then
d02nrf must be called to determine which column of the Jacobian is required and where it should be stored.
Call d02nrf(j, iplace, inform)
will return in
j the number of the column of the Jacobian that is required and will set
${\mathbf{iplace}}=1$ or
$2$. If
${\mathbf{iplace}}=1$, the
$\left(i,j\right)$th element of the Jacobian must be stored in
${\mathbf{rwork}}\left(50+2\times {\mathbf{neq}}+i\right)$; otherwise it must be stored in
${\mathbf{rwork}}\left(50+{\mathbf{neq}}+i\right)$.
 ${\mathbf{irevcm}}=9$
 Indicates that a MONITR operation (see Section 3) can be performed.
 ${\mathbf{irevcm}}=10$
 Indicates that the current step was not successful, due to error test failure or convergence test failure. The only information supplied to you on this return is the current value of the independent variable $t$, located in ${\mathbf{rwork}}\left(19\right)$. No values must be changed before reentering d02nmf; this facility enables you to determine the number of unsuccessful steps.
On final exit:
${\mathbf{irevcm}}=0$ indicated the userspecified task has been completed or an error has been encountered (see the descriptions for
itask and
ifail).
Constraint:
${\mathbf{irevcm}}=0$, $1$, $3$, $4$, $5$, $8$, $9$ or $10$.
Note: any values you return to d02nmf as part of the reverse communication procedure should not include floatingpoint NaN (Not a Number) or infinity values, since these are not handled by d02nmf. If your code does inadvertently return any NaNs or infinities, d02nmf is likely to produce unexpected results.
 22: $\mathbf{itask}$ – IntegerInput

On initial entry: the task to be performed by the integrator.
 ${\mathbf{itask}}=1$
 Normal computation of output values of $y\left(t\right)$ at $t={\mathbf{tout}}$ (by overshooting and interpolating).
 ${\mathbf{itask}}=2$
 Take one step only and return.
 ${\mathbf{itask}}=3$
 Stop at the first internal integration point at or beyond $t={\mathbf{tout}}$ and return.
 ${\mathbf{itask}}=4$
 Normal computation of output values of $y\left(t\right)$ at $t={\mathbf{tout}}$ but without overshooting $t={\mathbf{tcrit}}$. tcrit must be specified as an option in one of the integrator setup routines before the first call to the integrator, or specified in the optional input routine before a continuation call. tcrit (e.g., see d02nvf) may be equal to or beyond tout, but not before it in the direction of integration.
 ${\mathbf{itask}}=5$
 Take one step only and return, without passing tcrit (e.g., see d02nvf). tcrit must be specified under ${\mathbf{itask}}=4$.
Constraint:
$1\le {\mathbf{itask}}\le 5$.
 23: $\mathbf{itrace}$ – IntegerInput

On initial entry: the level of output that is printed by the integrator.
itrace may take the value
$1$,
$0$,
$1$,
$2$ or
$3$.
 ${\mathbf{itrace}}<1$
 $1$ is assumed and similarly if ${\mathbf{itrace}}>3$, $3$ is assumed.
 ${\mathbf{itrace}}=1$
 No output is generated.
 ${\mathbf{itrace}}=0$
 Only warning messages are printed on the current error message unit (see x04aaf).
 ${\mathbf{itrace}}>0$
 Warning messages are printed as above, and on the current advisory message unit (see x04abf) output is generated which details Jacobian entries, the nonlinear iteration and the time integration. The advisory messages are given in greater detail the larger the value of itrace.
 24: $\mathbf{ifail}$ – IntegerInput/Output

On initial 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, 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 final 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$

Either the integrator setup routine has not been called prior to the first call of this routine, or a communication array has become corrupted.
Either the linear algebra setup routine has not been called prior to the first call of this routine, or a communication array has become corrupted.
Either the routine was entered on a continuation call without a prior call of this routine, or a communication array has become corrupted.
Either the value of
njcpvt is not the same as the value supplied to the setup routine or a communication array has become corrupted.
${\mathbf{njcpvt}}=\u2329\mathit{\text{value}}\u232a$,
njcpvt (setup)
$=\u2329\mathit{\text{value}}\u232a$.
Either the value of
nwkjac is not the same as the value supplied to the setup routine or a communication array has become corrupted.
${\mathbf{nwkjac}}=\u2329\mathit{\text{value}}\u232a$,
nwkjac (setup)
$=\u2329\mathit{\text{value}}\u232a$.
Either the value of
sdysav is not the same as the value supplied to the setup routine or a communication array has become corrupted.
${\mathbf{sdysav}}=\u2329\mathit{\text{value}}\u232a$,
sdysav (setup)
$=\u2329\mathit{\text{value}}\u232a$.
Failure during internal time interpolation.
tout and the current time are too close.
${\mathbf{itask}}=\u2329\mathit{\text{value}}\u232a$ and
${\mathbf{tout}}=\u2329\mathit{\text{value}}\u232a$ and the current time is
$\u2329\mathit{\text{value}}\u232a$.
${\mathbf{itask}}=3$ and
tout is more than an integration step behind the current time.
${\mathbf{tout}}=\u2329\mathit{\text{value}}\u232a$, current time minus step size:
$\u2329\mathit{\text{value}}\u232a$.
On entry, an illegal (negative) maximum number of steps was provided in a prior call to a setup routine. ${\mathbf{maxstp}}=\u2329\mathit{\text{value}}\u232a$.
On entry, an illegal (negative) maximum stepsize was provided in a prior call to a setup routine. ${\mathbf{hmax}}=\u2329\mathit{\text{value}}\u232a$.
On entry, an illegal (negative) minimum stepsize was provided in a prior call to a setup routine. ${\mathbf{hmin}}=\u2329\mathit{\text{value}}\u232a$.
On entry, ${\mathbf{atol}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{atol}}\ge 0.0$.
On entry,
irevcm is invalid:
${\mathbf{irevcm}}=\u2329\mathit{\text{value}}\u232a$.
On entry, ${\mathbf{itask}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: $1\le {\mathbf{itask}}\le 5$.
On entry,
${\mathbf{itask}}=4$ or
$5$ and
tcrit is before the current time in the direction of integration.
${\mathbf{itask}}=\u2329\mathit{\text{value}}\u232a$,
${\mathbf{tcrit}}=\u2329\mathit{\text{value}}\u232a$ and the current time is
$\u2329\mathit{\text{value}}\u232a$.
On entry,
${\mathbf{itask}}=4$ or
$5$ and
tcrit is before
tout in the direction of integration.
${\mathbf{itask}}=\u2329\mathit{\text{value}}\u232a$,
${\mathbf{tcrit}}=\u2329\mathit{\text{value}}\u232a$ and
${\mathbf{tout}}=\u2329\mathit{\text{value}}\u232a$.
On entry, ${\mathbf{itol}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: $1\le {\mathbf{itol}}\le 4$.
On entry, ${\mathbf{neq}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{neq}}\ge 1$.
On entry, ${\mathbf{neq}}=\u2329\mathit{\text{value}}\u232a$ and ${\mathbf{ldysav}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{neq}}\le {\mathbf{ldysav}}$.
On entry, ${\mathbf{rtol}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{rtol}}\ge 0.0$.
On entry,
tout is less than
t with respect to the direction of integration given by the sign of
h0 in a prior call to a setup routine.
${\mathbf{tout}}=\u2329\mathit{\text{value}}\u232a$,
${\mathbf{t}}=\u2329\mathit{\text{value}}\u232a$ and
${\mathbf{h0}}=\u2329\mathit{\text{value}}\u232a$.
On entry,
tout is too close to
t to start integration.
${\mathbf{tout}}=\u2329\mathit{\text{value}}\u232a$ and
${\mathbf{t}}=\u2329\mathit{\text{value}}\u232a$.
On reentry, ${\mathbf{imon}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: $2\le {\mathbf{imon}}\le 4$.
On reentry, ${\mathbf{inln}}\ne 3$ for the case ${\mathbf{irevcm}}=9$ and ${\mathbf{imon}}=0$.
${\mathbf{inln}}=\u2329\mathit{\text{value}}\u232a$.
On reentry, the solution vector appears to have been overwritten.
Further integration will not be attempted.
The initial stepsize, $\u2329\mathit{\text{value}}\u232a$, is too small.
Weight number
$i=\u2329\mathit{\text{value}}\u232a$ used in the local error test is too small. Check the values of
rtol and
atol.
${\mathbf{atol}}\left(i\right)$ and
${\mathbf{y}}\left(i\right)$ may both be zero.
Weight
$i=\u2329\mathit{\text{value}}\u232a$.
 ${\mathbf{ifail}}=2$

At time $\u2329\mathit{\text{value}}\u232a$ the maximum number of allowed steps on this call was taken before reaching the next output point ${\mathbf{tout}}=\u2329\mathit{\text{value}}\u232a$.
Maximum number of steps $=\u2329\mathit{\text{value}}\u232a$.
 ${\mathbf{ifail}}=3$

Too much accuracy requested for precision of the machine at time $\u2329\mathit{\text{value}}\u232a$.
The tolerances should be checked; the requested accuracy should be reduced by a factor of at least $\u2329\mathit{\text{value}}\u232a$.
With the given values of
rtol and
atol no further progress can be made across the integration range from the current point
t. The components
${\mathbf{y}}\left(1\right),{\mathbf{y}}\left(2\right),\dots ,{\mathbf{y}}\left({\mathbf{neq}}\right)$ contain the computed values of the solution at the current point
t.
 ${\mathbf{ifail}}=4$

There were repeated errortest failures on an attempted step, before completing the requested task, but the integration was successful as far as
t. The problem may have a singularity, or the local error requirements may be inappropriate.
 ${\mathbf{ifail}}=5$

There were repeated convergencetest failures on an attempted step, before completing the requested task, but the integration was successful as far as
t. This may be caused by an inaccurate Jacobian matrix or one which is incorrectly computed.
 ${\mathbf{ifail}}=6$

At time
$\u2329\mathit{\text{value}}\u232a$, error weight
$\u2329\mathit{\text{value}}\u232a$ became zero. Check the values of
atol,
rtol and
itol supplied.
 ${\mathbf{ifail}}=7$

The derivative evaluation procedure set the error flag ${\mathbf{ires}}=3$ repeatedly despite attempts by the integrator to avoid this.
 ${\mathbf{ifail}}=9$

A singular Jacobian has been encountered. You should check the problem formulation and Jacobian calculation.
 ${\mathbf{ifail}}=10$

Larger integer workspace required.
Provided: $\u2329\mathit{\text{value}}\u232a$; required: $\u2329\mathit{\text{value}}\u232a$.
Not enough integer store provided for sparse matrix solver.
Units of store needed: $\u2329\mathit{\text{value}}\u232a$. Amount provided: $\u2329\mathit{\text{value}}\u232a$.
Not enough real store provided for sparse matrix solver.
Units of store needed: $\u2329\mathit{\text{value}}\u232a$. Amount provided: $\u2329\mathit{\text{value}}\u232a$.
Workspace error occurred when trying to form the Jacobian matrix in calculating the initial values of the solution and its time derivative.
 ${\mathbf{ifail}}=11$

On reentry, ${\mathbf{ires}}=2$, which signals that the integration should terminate. ${\mathbf{ires}}=\u2329\mathit{\text{value}}\u232a$ at time $\u2329\mathit{\text{value}}\u232a$.
 ${\mathbf{ifail}}=12$

A return was forced by setting
${\mathbf{imon}}=2$, but the integration was successful as far as
t.
 ${\mathbf{ifail}}=13$

The requested task has been completed, but it is estimated that a small change in
rtol and
atol is unlikely to produce any change in the computed solution. (This ONLY applies when you are NOT operating in one step mode; that is, when
${\mathbf{itask}}\ne 2$ or
$5$.)
 ${\mathbf{ifail}}=14$

On entry, too much accuracy requested for precision of the machine at the start of problem. The tolerances should be checked; the requested accuracy should be reduced by a factor of at least $\u2329\mathit{\text{value}}\u232a$.
 ${\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
The accuracy of the numerical solution may be controlled by a careful choice of the arguments
rtol and
atol, and to a much lesser extent by the choice of norm. You are advised to use scalar error control unless the components of the solution are expected to be poorly scaled. For the type of decaying solution typical of many stiff problems, relative error control with a small absolute error threshold will be most appropriate (that is, you are advised to choose
${\mathbf{itol}}=1$ with
${\mathbf{atol}}\left(1\right)$ small but positive).
8
Parallelism and Performance
d02nmf 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.
d02nmf is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.
d02nmf 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 cost of computing a solution depends critically on the size of the differential system and to a lesser extent on the degree of stiffness of the problem; also on the type of linear algebra being used. For further details see
Section 9 of the documents for
d02nbf (full matrix),
d02ncf (banded matrix) or
d02ndf (sparse matrix).
In general, you are advised to choose the backward differentiation formula option (setup routine
d02nvf) but if efficiency is of great importance and especially if it is suspected that
$\frac{\partial g}{\partial y}$ has complex eigenvalues near the imaginary axis for some part of the integration, you should try the BLEND option (setup routine
d02nwf).
10
Example
This example solves the wellknown stiff Robertson problem
over the range
$\left[0,10\right]$ with initial conditions
$a=1.0$ and
$b=c=0.0$ and with scalar error control (
${\mathbf{itol}}=1$). The integration proceeds until
${\mathbf{tout}}=10.0$ is passed, providing
${C}^{1}$ interpolation at intervals of
$2.0$ through a MONITR operation. The integration method used is the BDF method (setup routine
d02nvf) with a modified Newton method. The Jacobian is a full matrix, which is specified using the setup routine
d02nsf; this Jacobian is to be calculated numerically.
10.1
Program Text
Program Text (d02nmfe.f90)
10.2
Program Data
Program Data (d02nmfe.d)
10.3
Program Results
Program Results (d02nmfe.r)