NAG Library Routine Document
D02NEF
1 Purpose
D02NEF is a routine for integrating stiff systems of implicit ordinary differential equations coupled with algebraic equations.
2 Specification
SUBROUTINE D02NEF ( 
NEQ, T, TOUT, Y, YDOT, RTOL, ATOL, ITASK, RES, JAC, ICOM, COM, LCOM, IUSER, RUSER, IFAIL) 
INTEGER 
NEQ, ITASK, ICOM(50+NEQ), LCOM, IUSER(*), IFAIL 
REAL (KIND=nag_wp) 
T, TOUT, Y(NEQ), YDOT(NEQ), RTOL(*), ATOL(*), COM(LCOM), RUSER(*) 
EXTERNAL 
RES, JAC 

3 Description
D02NEF is a general purpose routine for integrating the initial value problem for a stiff system of implicit ordinary differential equations with coupled algebraic equations written in the form
D02NEF uses the DASSL implementation of the Backward Differentiation Formulae (BDF) of orders one to five to solve a system of the above form for
$y$ (
Y) and
${y}^{\prime}$ (
YDOT). Values for
Y and
YDOT at the initial time must be given as input. These values must be consistent, (i.e., if
T,
Y,
YDOT are the given initial values, they must satisfy
$F\left({\mathbf{T}},{\mathbf{Y}},{\mathbf{YDOT}}\right)=0$). The routine solves the system from
$t={\mathbf{T}}$ to
$t={\mathbf{TOUT}}$.
An outline of a typical calling program for D02NEF is given below. It calls the DASSL implementation of the BDF integrator setup routine
D02MWF and the banded matrix setup routine
D02NPF (if required), and, if the integration needs to proceed, calls
D02MCF before continuing the integration.
! Declarations
EXTERNAL RES, JAC
.
.
.
! Initialize the integrator
CALL D02MWF(...)
! Is the Jacobian matrix banded?
IF (BANDED) CALL D02NPF(...)
! Set DT to the required temporal resolution
! Set TEND to the final time
! Call the integrator for each temporal value:
1000 CALL D02NEF(...,RES,JAC,...)
! Continue integration?
IF (TOUT.LT.TEND .AND. ITASK.GE.0) THEN
IF (ITASK.NE.1) TOUT = MIN(TOUT+DT,TEND)
! Print solution
CALL D02MCF(...)
GO TO 1000
ENDIF
.
.
.
4 References
None.
5 Parameters
 1: NEQ – INTEGERInput
On entry: the number of differentialalgebraic equations to be solved.
Constraint:
${\mathbf{NEQ}}\ge 1$.
 2: T – REAL (KIND=nag_wp)Input/Output
On initial entry: the initial value of the independent variable, $t$.
On intermediate exit:
$t$, the current value of the independent variable.
On final exit: the value of the independent variable at which the computed solution
$y$ is returned (usually at
TOUT).
 3: TOUT – REAL (KIND=nag_wp)Input
On entry: the next value of $t$ at which a computed solution is desired.
On initial entry:
TOUT is used to determine the direction of integration. Integration is permitted in either direction (see also
ITASK).
Constraint:
${\mathbf{TOUT}}\ne {\mathbf{T}}$.
 4: Y(NEQ) – REAL (KIND=nag_wp) arrayInput/Output
On initial entry: the vector of initial values of the dependent variables $y$.
On intermediate exit:
the computed solution vector, $y$, evaluated at $t=T$.
On final exit: the computed solution vector, evaluated at $t$ (usually $t={\mathbf{TOUT}}$).
 5: YDOT(NEQ) – REAL (KIND=nag_wp) arrayInput/Output
On initial entry:
YDOT must contain approximations to the time derivatives
$y\text{'}$ of the vector
$y$ evaluated at the initial value of the independent variable.
On exit: the time derivatives $y\text{'}$ of the vector $y$ at the last integration point.
 6: RTOL($*$) – REAL (KIND=nag_wp) arrayInput/Output
Note: the dimension of the array
RTOL depends on the value of
ITOL as set in
D02MWF; it
must be at least
${\mathbf{NEQ}}$ if
${\mathbf{ITOL}}=\mathrm{.TRUE.}$ and at least
$1$ if
${\mathbf{ITOL}}=\mathrm{.FALSE.}$.
On entry: the relative local error tolerance.
Constraint:
${\mathbf{RTOL}}\left(\mathit{i}\right)\ge 0.0$, for
$\mathit{i}=1,2,\dots ,n$where $n={\mathbf{NEQ}}$ when ${\mathbf{ITOL}}=\mathrm{.TRUE.}$ and $n=1$ otherwise.
On exit:
RTOL remains unchanged unless D02NEF exits with
${\mathbf{IFAIL}}={\mathbf{16}}$ in which case the values may have been increased to values estimated to be appropriate for continuing the integration.
 7: ATOL($*$) – REAL (KIND=nag_wp) arrayInput/Output
Note: the dimension of the array
ATOL depends on the value of
ITOL as set in
D02MWF; it
must be at least
${\mathbf{NEQ}}$ if
${\mathbf{ITOL}}=\mathrm{.TRUE.}$ and at least
$1$ if
${\mathbf{ITOL}}=\mathrm{.FALSE.}$.
On entry: the absolute local error tolerance.
Constraint:
${\mathbf{ATOL}}\left(\mathit{i}\right)\ge 0.0$, for
$\mathit{i}=1,2,\dots ,n$where $n={\mathbf{NEQ}}$ when ${\mathbf{ITOL}}=\mathrm{.TRUE.}$ and $n=1$ otherwise.
On exit:
ATOL remains unchanged unless D02NEF exits with
${\mathbf{IFAIL}}={\mathbf{16}}$ in which case the values may have been increased to values estimated to be appropriate for continuing the integration.
 8: ITASK – INTEGERInput/Output
On initial entry: need not be set.
On exit: the task performed by the integrator on successful completion or an indicator that a problem occurred during integration.
 ${\mathbf{ITASK}}=2$
 The integration to TOUT was successfully completed (${\mathbf{T}}={\mathbf{TOUT}}$) by stepping exactly to TOUT.
 ${\mathbf{ITASK}}=3$
 The integration to TOUT was successfully completed (${\mathbf{T}}={\mathbf{TOUT}}$) by stepping past TOUT. Y and YDOT are obtained by interpolation.
 ${\mathbf{ITASK}}<0$
 Different negative values of ITASK returned correspond to different failure exits. IFAIL should always be checked in such cases and the corrective action taken where appropriate.
ITASK must remain
unchanged between calls to D02NEF.
 9: RES – SUBROUTINE, supplied by the user.External Procedure
RES must evaluate the residual
The specification of
RES is:
INTEGER 
NEQ, IRES, IUSER(*) 
REAL (KIND=nag_wp) 
T, Y(NEQ), YDOT(NEQ), R(NEQ), RUSER(*) 

 1: NEQ – INTEGERInput
On entry: the number of differentialalgebraic equations being solved.
 2: T – REAL (KIND=nag_wp)Input
On entry: $t$, the current value of the independent variable.
 3: Y(NEQ) – REAL (KIND=nag_wp) arrayInput
On entry: ${y}_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,{\mathbf{NEQ}}$, the current solution component.
 4: YDOT(NEQ) – REAL (KIND=nag_wp) arrayInput
On entry: the derivative of the solution at the current point $t$.
 5: R(NEQ) – REAL (KIND=nag_wp) arrayOutput
On exit:
${\mathbf{R}}\left(\mathit{i}\right)$ must contain the
$\mathit{i}$th component of
$R$, for
$\mathit{i}=1,2,\dots ,{\mathbf{NEQ}}$ where
 6: IRES – INTEGERInput/Output
On entry: is always equal to zero.
On exit:
IRES should normally be left unchanged. However, if an illegal value of
Y is encountered,
IRES should be set to
$1$; D02NEF will then attempt to resolve the problem so that illegal values of
Y are not encountered.
IRES should be set to
$2$ if you wish to return control to the calling (sub)routine; this will cause D02NEF to exit with
${\mathbf{IFAIL}}={\mathbf{23}}$.
 7: IUSER($*$) – INTEGER arrayUser Workspace
 8: RUSER($*$) – REAL (KIND=nag_wp) arrayUser Workspace

RES is called with the parameters
IUSER and
RUSER as supplied to D02NEF. You are free to use the arrays
IUSER and
RUSER to supply information to
RES as an alternative to using COMMON global variables.
RES must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which D02NEF is called. Parameters denoted as
Input must
not be changed by this procedure.
 10: JAC – SUBROUTINE, supplied by the NAG Library or the user.External Procedure
Evaluates the matrix of partial derivatives,
$J$, where
If this option is not required, the actual argument for
JAC must be the dummy routine D02NEZ. (D02NEZ is included in the NAG Library.) You must indicate to the integrator whether this option is to be used by setting the parameter
JCEVAL appropriately in a call to the setup routine
D02MWF.
The specification of
JAC is:
INTEGER 
NEQ, IUSER(*) 
REAL (KIND=nag_wp) 
T, Y(NEQ), YDOT(NEQ), PD(*), CJ, RUSER(*) 

 1: NEQ – INTEGERInput
On entry: the number of differentialalgebraic equations being solved.
 2: T – REAL (KIND=nag_wp)Input
On entry: $t$, the current value of the independent variable.
 3: Y(NEQ) – REAL (KIND=nag_wp) arrayInput
On entry: ${y}_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,{\mathbf{NEQ}}$, the current solution component.
 4: YDOT(NEQ) – REAL (KIND=nag_wp) arrayInput
On entry: the derivative of the solution at the current point $t$.
 5: PD($*$) – REAL (KIND=nag_wp) arrayInput/Output
On entry:
PD is preset to zero before the call to
JAC.
On exit: if the Jacobian is full then
${\mathbf{PD}}\left(\left(\mathit{i}\mathit{j}\right)\times {\mathbf{NEQ}}+\mathit{i}\right)={J}_{\mathit{i}\mathit{j}}$, for
$\mathit{i}=1,2,\dots ,{\mathbf{NEQ}}$ and
$\mathit{j}=1,2,\dots ,{\mathbf{NEQ}}$; if the Jacobian is banded then
${\mathbf{PD}}\left(\left(j1\right)\times {\mathbf{NEQ}}+{\mathbf{ML}}+{\mathbf{MU}}+ij\right)={J}_{ij}$, for
$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,j{\mathbf{MU}}\right)\le i\le \mathrm{min}\phantom{\rule{0.125em}{0ex}}\left(n,j+{\mathbf{ML}}\right)$; (see also in
F07BDF (DGBTRF)).
 6: CJ – REAL (KIND=nag_wp)Input
On entry:
CJ is a scalar constant which will be defined in D02NEF.
 7: IUSER($*$) – INTEGER arrayUser Workspace
 8: RUSER($*$) – REAL (KIND=nag_wp) arrayUser Workspace

JAC is called with the parameters
IUSER and
RUSER as supplied to D02NEF. You are free to use the arrays
IUSER and
RUSER to supply information to
JAC as an alternative to using COMMON global variables.
JAC must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which D02NEF is called. Parameters denoted as
Input must
not be changed by this procedure.
 11: ICOM($50+{\mathbf{NEQ}}$) – INTEGER arrayCommunication Array
ICOM contains information which is usually of no interest, but is necessary for subsequent calls. However you may find the following useful:
 ${\mathbf{ICOM}}\left(22\right)$
 The order of the method to be attempted on the next step.
 ${\mathbf{ICOM}}\left(23\right)$
 The order of the method used on the last step.
 ${\mathbf{ICOM}}\left(26\right)$
 The number of steps taken so far.
 ${\mathbf{ICOM}}\left(27\right)$
 The number of calls to RES so far.
 ${\mathbf{ICOM}}\left(28\right)$
 The number of evaluations of the matrix of partial derivatives needed so far.
 ${\mathbf{ICOM}}\left(29\right)$
 The total number of error test failures so far.
 ${\mathbf{ICOM}}\left(30\right)$
 The total number of convergence test failures so far.
 12: COM(LCOM) – REAL (KIND=nag_wp) arrayCommunication Array
COM contains information which is usually of no interest, but is necessary for subsequent calls. However you may find the following useful:
 ${\mathbf{COM}}\left(3\right)$
 The step size to be attempted on the next step.
 ${\mathbf{COM}}\left(4\right)$
 The current value of the independent variable, i.e., the farthest point integration has reached. This will be different from T only when interpolation has been performed (${\mathbf{ITASK}}=3$).
 13: LCOM – INTEGERInput
On entry: the dimension of the array
COM as declared in the (sub)program from which D02NEF is called.
Constraint:
${\mathbf{LCOM}}\ge 40+\left(\mathit{maxorder}+4\right)\times {\mathbf{NEQ}}+{\mathbf{NEQ}}\times p+q$ where
$\mathit{maxorder}$ is the maximum order that can be used by the integration method (see
MAXORD in
D02MWF);
$p={\mathbf{NEQ}}$ when the Jacobian is full and
$p=\left(2\times {\mathbf{ML}}+{\mathbf{MU}}+1\right)$ when the Jacobian is banded; and,
$q=\left({\mathbf{NEQ}}/\left({\mathbf{ML}}+{\mathbf{MU}}+1\right)\right)+1$ when the Jacobian is to be evaluated numerically and
$q=0$ otherwise.
 14: IUSER($*$) – INTEGER arrayUser Workspace
 15: RUSER($*$) – REAL (KIND=nag_wp) arrayUser Workspace

IUSER and
RUSER are not used by D02NEF, but are passed directly to
RES and
JAC and may be used to pass information to these routines as an alternative to using COMMON global variables.
 16: IFAIL – INTEGERInput/Output

On entry:
IFAIL must be set to
$0$,
$1\text{ or}1$. If you are unfamiliar with this parameter you should refer to
Section 3.3 in the Essential Introduction 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 parameters 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}}={\mathbf{0}}$ or
${{\mathbf{1}}}$, explanatory error messages are output on the current error message unit (as defined by
X04AAF).
Note: D02NEF may return useful information for one or more of the following detected errors or warnings.
Errors or warnings detected by the routine:
 ${\mathbf{IFAIL}}=1$
On entry,  ${\mathbf{NEQ}}<1$ 
 ${\mathbf{IFAIL}}=3$
On entry,  ${\mathbf{TOUT}}={\mathbf{T}}$, 
or  TOUT is too close to T to start integration, 
or  TOUT is behind T in the direction of H0 (see D02MWF), 
 ${\mathbf{IFAIL}}=6$
On entry,  ${\mathbf{RTOL}}\left(\mathit{i}\right)<0.0$, for $\mathit{i}=1,2,\dots ,n$, where $n={\mathbf{NEQ}}$ when ${\mathbf{ITOL}}=1$ and $n=1$ otherwise, 
or  ${\mathbf{RTOL}}\left(i\right)={\mathbf{ATOL}}\left(i\right)=0.0$ for all relevant $i$. 
 ${\mathbf{IFAIL}}=7$
On entry,  ${\mathbf{ATOL}}\left(\mathit{i}\right)<0.0$, for $\mathit{i}=1,2,\dots ,n$, where $n={\mathbf{NEQ}}$ when ${\mathbf{ITOL}}=1$ and $n=1$ otherwise. 
 ${\mathbf{IFAIL}}=8$
On entry,  a previous call to this routine returned ${\mathbf{ITASK}}<0$ and ${\mathbf{IFAIL}}\ne 0$, but no appropriate action was taken. For example, if a call returns with ${\mathbf{IFAIL}}={\mathbf{15}}$ (${\mathbf{ITASK}}=1$) then a call to D02MCF must be made prior to making a continuation call to D02NEF. 
 ${\mathbf{IFAIL}}=12$
Either the initialization routine
D02MWF has not been called prior to the first call of this routine or one of the communication arrays
ICOM,
COM has become corrupted.
 ${\mathbf{IFAIL}}=13$
On entry,  ${\mathbf{LCOM}}<40+\left(\mathit{maxorder}+4\right)\times {\mathbf{NEQ}}+{\mathbf{NEQ}}\times {\mathbf{NEQ}}$, and the Jacobian is full; 
or  ${\mathbf{LCOM}}<40+\left(\mathit{maxorder}+4\right)\times {\mathbf{NEQ}}+{\mathbf{NEQ}}\times \left(2\times {\mathbf{ML}}+{\mathbf{MU}}+1\right)+q$, and the Jacobian is banded, 
where $\mathit{maxorder}$ is the maximum order of the integration method to be used, and $q=\left({\mathbf{NEQ}}/\left({\mathbf{ML}}+{\mathbf{MU}}+1\right)\right)+1$ when the Jacobian is to be evaluated numerically and $q=0$ otherwise.
 ${\mathbf{IFAIL}}=15$
The maximum number of steps (
$500$) has been taken on this call and the integration has not reached
TOUT. The integration can proceed by calling
D02MCF prior to calling D02NEF again; this will reset the step counter to zero.
 ${\mathbf{IFAIL}}=16$
Too much accuracy was requested for the precision of the machine. On output
RTOL and
ATOL were increased by an appropriate scale factor to prevent this error exit. Try running the problem again with these scaled tolerances.
 ${\mathbf{IFAIL}}=17$
A purely relative tolerance was selected for a given solution component, but that solution component has become zero and a pure relative error test is impossible for this component. Perhaps an absolute tolerance requirement is also necessary for this component.
 ${\mathbf{IFAIL}}=18$
The error test failed repeatedly using the minimum stepsize and so the integration could not proceed. If a (nonzero) minimum stepsize was specified in a call to
D02MWF then either a reduction in this minimum stepsize or in the specified tolerances should be considered.
 ${\mathbf{IFAIL}}=19$
The corrector step to obtain the approximate solution at the next time step repeatedly failed to converge using the minimum stepsize. This may be due to an inconsistency between the system of equations to be solved and initial conditions.
 ${\mathbf{IFAIL}}=20$
The iteration matrix (Jacobian) has become singular. Please check the Jacobian evaluations when this is performed analytically. Also check for invalid solution values in the call to
RES; these should be flagged by returning
${\mathbf{IRES}}=1$. This is probably due to an error in the analytic evaluation of the Jacobian.
 ${\mathbf{IFAIL}}=21$
The corrector step could not converge and the error test failed repeatedly.
 ${\mathbf{IFAIL}}=22$
IRES was set to
$1$ during a call to
RES and the problem could not be resolved.
 ${\mathbf{IFAIL}}=23$
IRES was set to
$2$ during a call to
RES.
 ${\mathbf{IFAIL}}=24$
The initial
YDOT could not be computed. This could happen because the initial approximation to
YDOT was very poor or because no
YDOT exists that is consistent with the initial
$Y$.
 ${\mathbf{IFAIL}}=25$
Repeated occurrences of input constraint violations have been detected. This could result in a potential infinite loop.
7 Accuracy
The accuracy of the numerical solution may be controlled by a careful choice of the parameters
RTOL and
ATOL. 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}}=0$ with
${\mathbf{ATOL}}\left(1\right)$ small but positive).
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. For banded systems the cost is proportional to ${\mathbf{NEQ}}\times {\left({\mathbf{ML}}+{\mathbf{MU}}+1\right)}^{2}$, while for full systems the cost is proportional to
${{\mathbf{NEQ}}}^{3}$. Note however that for moderately sized problems which are only mildly nonlinear the cost may be dominated by factors proportional to ${\mathbf{NEQ}}\times \left({\mathbf{ML}}+{\mathbf{MU}}+1\right)$ and ${{\mathbf{NEQ}}}^{2}$ respectively.
9 Example
For this routine two examples are presented. There is a single example program for D02NEF, with a main program and the code to solve the two example problems given in Example 1 (EX1) and Example 2 (EX2).
Example 1 (EX1)
This example solves the wellknown stiff Robertson problem written in implicit form
with initial conditions
$a=1.0$ and
$b=c=0.0$ over the range
$\left[0,0.1\right]$ the BDF method (setup routine
D02MWF and
D02NPF).
Example 2 (EX2)
This example illustrates the use of D02NEF to solve a simple algebraic problem by continuation. The equation $42y+0.1{e}^{y}t=0$ from $t=0$ (where $y=2$) to $t=1$.
9.1 Program Text
Program Text (d02nefe.f90)
9.2 Program Data
Program Data (d02nefe.d)
9.3 Program Results
Program Results (d02nefe.r)