# NAG FL Interfaced02nhf (ivp_​stiff_​imp_​bandjac)

## ▸▿ Contents

Settings help

FL Name Style:

FL Specification Language:

## 1Purpose

d02nhf is a direct communication routine for integrating stiff systems of implicit ordinary differential equations coupled with algebraic equations when the Jacobian is a banded matrix.

## 2Specification

Fortran Interface
 Subroutine d02nhf ( neq, t, tout, y, ydot, rtol, atol, itol, ysav, jac,
 Integer, Intent (In) :: neq, ldysav, itol, sdysav, nwkjac, njcpvt, itask, itrace Integer, Intent (Inout) :: inform(23), jacpvt(njcpvt), ifail Real (Kind=nag_wp), Intent (In) :: rtol(*), atol(*) Real (Kind=nag_wp), Intent (Inout) :: t, tout, y(neq), ydot(neq), rwork(50+4*neq), ysav(ldysav,sdysav), wkjac(nwkjac) Logical, Intent (Inout) :: lderiv(2) External :: resid, jac, monitr
C Header Interface
#include <nag.h>
 void d02nhf_ (const Integer *neq, const Integer *ldysav, double *t, double *tout, double y[], double ydot[], double rwork[], const double rtol[], const double atol[], const Integer *itol, Integer inform[], void (NAG_CALL *resid)(const Integer *neq, const double *t, const double y[], const double ydot[], double r[], Integer *ires),double ysav[], const Integer *sdysav, void (NAG_CALL *jac)(const Integer *neq, const double *t, const double y[], const double ydot[], const double *h, const double *d, const Integer *ml, const Integer *mu, double p[]),double wkjac[], const Integer *nwkjac, Integer jacpvt[], const Integer *njcpvt, void (NAG_CALL *monitr)(const Integer *neq, const Integer *ldysav, const double *t, const double *hlast, double *hnext, double y[], const double ydot[], const double ysav[], const double r[], const double acor[], Integer *imon, Integer *inln, double *hmin, double *hmax, const Integer *nqu),logical lderiv[], const Integer *itask, const Integer *itrace, Integer *ifail)
The routine may be called by the names d02nhf or nagf_ode_ivp_stiff_imp_bandjac.

## 3Description

d02nhf is a general purpose routine for integrating the initial value problem for a stiff system of implicit ordinary differential equations coupled with algebraic equations, written in the form
 $A(t,y)y′=g(t,y).$
It is designed specifically for the case where the resulting Jacobian is a banded matrix (see the description of jac).
Both interval and step oriented modes of operation are available and also modes designed to permit intermediate output within an interval oriented mode.
An outline of a typical calling program for d02nhf is given below. It calls the banded matrix linear algebra setup routine d02ntf, the Backward Differentiation Formula (BDF) integrator setup routine d02nvf, and its diagnostic counterpart d02nyf.
```!     Declarations

External resid, jac, monitr
.
.
.
ifail = 0
Call d02nvf(...,ifail)
Call d02ntf(neq, neqmax, jceval, ml, mu, nwkjac, njcpvt, &
rwork, ifail)
ifail = -1
Call d02nhf(neq, neqmax, t, tout, y, ydot, rwork, rtol,  &
atol, itol, inform, resid, ysave, ny2dim,    &
jac, wkjac, nwkjac, jacpvt, njcpvt, monitr,  &
lderiv, itask, itrace, ifail)
If (ifail.eq.1 .or. ifail.ge.14) Stop
ifail = 0
Call d02nyf(...)
.
.
.
Stop
End```
The linear algebra setup routine d02ntf and one of the integrator setup routines, d02mvf, d02nvf or d02nwf, must be called prior to the call of d02nhf. The integrator diagnostic routine d02nyf may be called after the call to d02nhf. There is also a routine, d02nzf, designed to permit you to change step size on a continuation call to d02nhf without restarting the integration process.

## 4References

See the D02M–N Sub-chapter Introduction.

## 5Arguments

1: $\mathbf{neq}$Integer Input
On entry: the number of differential equations to be solved.
Constraint: ${\mathbf{neq}}\ge 1$.
2: $\mathbf{ldysav}$Integer Input
On entry: a bound on the maximum number of equations to be solved during the integration.
Constraint: ${\mathbf{ldysav}}\ge {\mathbf{neq}}$.
3: $\mathbf{t}$Real (Kind=nag_wp) Input/Output
On 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 exit: the value at which the computed solution $y$ is returned (usually at tout).
4: $\mathbf{tout}$Real (Kind=nag_wp) Input/Output
On 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}}$.
On exit: normally unchanged. However, when ${\mathbf{itask}}=6$, tout contains the value of t at which initial values have been computed without performing any integration. See descriptions of itask and lderiv.
5: $\mathbf{y}\left({\mathbf{neq}}\right)$Real (Kind=nag_wp) array Input/Output
On 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 exit: the computed solution vector, evaluated at t (usually ${\mathbf{t}}={\mathbf{tout}}$).
6: $\mathbf{ydot}\left({\mathbf{neq}}\right)$Real (Kind=nag_wp) array Input/Output
On entry: if ${\mathbf{lderiv}}\left(1\right)=\mathrm{.TRUE.}$, ydot must contain approximations to the time derivatives ${y}^{\prime }$ of the vector $y$.
If ${\mathbf{lderiv}}\left(1\right)=\mathrm{.FALSE.}$, ydot need not be set on entry.
On exit: the time derivatives ${y}^{\prime }$ of the vector $y$ at the last integration point.
7: $\mathbf{rwork}\left(50+4×{\mathbf{neq}}\right)$Real (Kind=nag_wp) array Communication Array
8: $\mathbf{rtol}\left(*\right)$Real (Kind=nag_wp) array Input
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 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) array Input
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 entry: the absolute local error tolerance.
Constraint: ${\mathbf{atol}}\left(i\right)\ge 0.0$ for all relevant $i$ (see itol).
10: $\mathbf{itol}$Integer Input
On entry: a value to indicate the form of the local error test. itol indicates to d02nhf whether to interpret either or both of rtol or atol as a vector or a scalar. The error test to be satisfied is $‖{e}_{i}/{w}_{i}‖<1.0$, where ${w}_{i}$ is defined as follows:
itol rtol atol ${\mathbit{w}}_{\mathbit{i}}$
1 scalar scalar ${\mathbf{rtol}}\left(1\right)×|{y}_{i}|+{\mathbf{atol}}\left(1\right)$
2 scalar vector ${\mathbf{rtol}}\left(1\right)×|{y}_{i}|+{\mathbf{atol}}\left(i\right)$
3 vector scalar ${\mathbf{rtol}}\left(i\right)×|{y}_{i}|+{\mathbf{atol}}\left(1\right)$
4 vector vector ${\mathbf{rtol}}\left(i\right)×|{y}_{i}|+{\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 array Communication Array
12: $\mathbf{resid}$Subroutine, supplied by the user. External Procedure
resid must evaluate the residual
 $r=g(t,y)-A(t,y)y′$
in one case and
 $r=-A(t,y)y′$
in another.
The specification of resid is:
Fortran Interface
 Subroutine resid ( neq, t, y, ydot, r, ires)
 Integer, Intent (In) :: neq Integer, Intent (Inout) :: ires Real (Kind=nag_wp), Intent (In) :: t, y(neq), ydot(neq) Real (Kind=nag_wp), Intent (Out) :: r(neq)
C Header Interface
 void resid (const Integer *neq, const double *t, const double y[], const double ydot[], double r[], Integer *ires)
1: $\mathbf{neq}$Integer Input
On entry: the number of equations being solved.
2: $\mathbf{t}$Real (Kind=nag_wp) Input
On entry: $t$, the current value of the independent variable.
3: $\mathbf{y}\left({\mathbf{neq}}\right)$Real (Kind=nag_wp) array Input
On entry: the value of ${y}_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,{\mathbf{neq}}$.
4: $\mathbf{ydot}\left({\mathbf{neq}}\right)$Real (Kind=nag_wp) array Input
On entry: the value of ${y}_{\mathit{i}}^{\prime }$, for $\mathit{i}=1,2,\dots ,{\mathbf{neq}}$, at $t$.
5: $\mathbf{r}\left({\mathbf{neq}}\right)$Real (Kind=nag_wp) array Output
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
 $r=g(t,y)-A(t,y)y′$ (1)
or
 $r=-A(t,y)y′$ (2)
and where the definition of $r$ is determined by the input value of ires.
6: $\mathbf{ires}$Integer Input/Output
On entry: the form of the residual that must be returned in array r.
${\mathbf{ires}}=-1$
The residual defined in equation (2) must be returned.
${\mathbf{ires}}=1$
The residual defined in equation (1) must be returned.
On exit: should be unchanged unless one of the following actions is required of the integrator, in which case ires should be set accordingly.
${\mathbf{ires}}=2$
Indicates to the integrator 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 the integrator 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, the integrator returns to the calling (sub)program with the error indicator set to ${\mathbf{ifail}}={\mathbf{7}}$.
${\mathbf{ires}}=4$
Indicates to the integrator to stop its current operation and to enter monitr immediately with argument ${\mathbf{imon}}=-2$.
resid must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which d02nhf is called. Arguments denoted as Input must not be changed by this procedure.
Note: resid should not return floating-point NaN (Not a Number) or infinity values, since these are not handled by d02nhf. If your code inadvertently does return any NaNs or infinities, d02nhf is likely to produce unexpected results.
13: $\mathbf{ysav}\left({\mathbf{ldysav}},{\mathbf{sdysav}}\right)$Real (Kind=nag_wp) array Communication Array
14: $\mathbf{sdysav}$Integer Input
On entry: the second dimension of the array ysav as declared in the (sub)program from which d02nhf is called. An appropriate value for sdysav is described in the specifications of the integrator setup routines d02mvf, d02nvf and d02nwf. This value must be the same as that supplied to the integrator setup routine.
15: $\mathbf{jac}$Subroutine, supplied by the NAG Library or the user. External Procedure
jac must evaluate the Jacobian of the system. If this option is not required, the actual argument for jac must be the dummy routine d02nhz. (d02nhz is included in the NAG Library.) You must indicate to the integrator whether this option is to be used by setting the argument jceval appropriately in a call to the banded linear algebra setup routine d02ntf.
First we must define the system of nonlinear equations which is solved internally by the integrator. The time derivative, ${y}^{\prime }$, generated internally, has the form
 $y′ = (y-z) / (hd) ,$
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
 $A (t,y) y′ - g (t,y) = 0$
but it is solved in the form
 $r (t,y) = 0 ,$
where $r$ is the function defined by
 $r (t,y) = (hd) ( A (t,y) (y-z) / (hd) - g (t,y) ) .$
It is the Jacobian matrix $\frac{\partial r}{\partial y}$ that you must supply in jac as follows:
 $∂ri ∂yj = aij (t,y) + (hd) ∂ ∂yj ( ∑ k=1 neq aik(t,y)yk′-gi(t,y)) .$
The specification of jac is:
Fortran Interface
 Subroutine jac ( neq, t, y, ydot, h, d, ml, mu, p)
 Integer, Intent (In) :: neq, ml, mu Real (Kind=nag_wp), Intent (In) :: t, y(neq), ydot(neq), h, d Real (Kind=nag_wp), Intent (Inout) :: p(ml+mu+1,neq)
C Header Interface
 void jac (const Integer *neq, const double *t, const double y[], const double ydot[], const double *h, const double *d, const Integer *ml, const Integer *mu, double p[])
1: $\mathbf{neq}$Integer Input
On entry: the number of equations being solved.
2: $\mathbf{t}$Real (Kind=nag_wp) Input
On entry: $t$, the current value of the independent variable.
3: $\mathbf{y}\left({\mathbf{neq}}\right)$Real (Kind=nag_wp) array Input
On entry: ${y}_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,{\mathbf{neq}}$, the current solution component.
4: $\mathbf{ydot}\left({\mathbf{neq}}\right)$Real (Kind=nag_wp) array Input
On entry: the derivative of the solution at the current point $t$.
5: $\mathbf{h}$Real (Kind=nag_wp) Input
On entry: the current step size.
6: $\mathbf{d}$Real (Kind=nag_wp) Input
On entry: the argument $d$ which depends on the integration method.
7: $\mathbf{ml}$Integer Input
8: $\mathbf{mu}$Integer Input
On entry: the number of subdiagonals and superdiagonals respectively in the band.
9: $\mathbf{p}\left({\mathbf{ml}}+{\mathbf{mu}}+1,{\mathbf{neq}}\right)$Real (Kind=nag_wp) array Input/Output
On entry: is set to zero.
On exit: elements of the Jacobian matrix $\frac{\partial r}{\partial y}$ stored as specified by the following pseudocode:
```   Do i = 1, neq
j1 = max(i-ml,1)
j2 = min(i+mu,neq)
Do j = j1, j2
k = min(ml+1-i,0)+j
p(k,i) = δr/δy(i,j)
End Do
End Do```
See also f07bdf.
Only nonzero elements of this array need be set, since it is preset to zero before the call to jac.
jac must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which d02nhf is called. Arguments denoted as Input must not be changed by this procedure.
Note: jac should not return floating-point NaN (Not a Number) or infinity values, since these are not handled by d02nhf. If your code inadvertently does return any NaNs or infinities, d02nhf is likely to produce unexpected results.
16: $\mathbf{wkjac}\left({\mathbf{nwkjac}}\right)$Real (Kind=nag_wp) array Communication Array
17: $\mathbf{nwkjac}$Integer Input
On entry: the dimension of the array wkjac as declared in the (sub)program from which d02nhf is called. This value must be the same as that supplied to the linear algebra setup routine d02ntf.
Constraint: ${\mathbf{nwkjac}}\ge \left(2{m}_{L}+{m}_{U}+1\right)×{\mathbf{neq}}$, where ${m}_{L}$ and ${m}_{U}$ are the number of subdiagonals and superdiagonals respectively in the band, defined by a call to d02ntf.
18: $\mathbf{jacpvt}\left({\mathbf{njcpvt}}\right)$Integer array Communication Array
19: $\mathbf{njcpvt}$Integer Input
On entry: the size of the array jacpvt. This value must be the same as that supplied to the linear algebra setup routine d02ntf.
Constraint: ${\mathbf{njcpvt}}\ge {\mathbf{neq}}$.
20: $\mathbf{monitr}$Subroutine, supplied by the NAG Library or the user. External Procedure
monitr performs tasks requested by you. If this option is not required, the actual argument for monitr must be the dummy routine d02nby. (d02nby is included in the NAG Library.)
The specification of monitr is:
Fortran Interface
 Subroutine monitr ( neq, t, y, ydot, ysav, r, acor, imon, inln, hmin, hmax, nqu)
 Integer, Intent (In) :: neq, ldysav, nqu Integer, Intent (Inout) :: imon Integer, Intent (Out) :: inln Real (Kind=nag_wp), Intent (In) :: t, hlast, ydot(neq), ysav(ldysav,*), r(neq), acor(neq,2) Real (Kind=nag_wp), Intent (Inout) :: hnext, y(neq), hmin, hmax
C Header Interface
 void monitr (const Integer *neq, const Integer *ldysav, const double *t, const double *hlast, double *hnext, double y[], const double ydot[], const double ysav[], const double r[], const double acor[], Integer *imon, Integer *inln, double *hmin, double *hmax, const Integer *nqu)
1: $\mathbf{neq}$Integer Input
On entry: the number of equations being solved.
2: $\mathbf{ldysav}$Integer Input
On entry: an upper bound on the number of equations to be solved.
3: $\mathbf{t}$Real (Kind=nag_wp) Input
On entry: the current value of the independent variable.
4: $\mathbf{hlast}$Real (Kind=nag_wp) Input
On entry: the last step size successfully used by the integrator.
5: $\mathbf{hnext}$Real (Kind=nag_wp) Input/Output
On entry: the step size that the integrator proposes to take on the next step.
On exit: the next step size to be used. If this is different from the input value, imon must be set to $4$.
6: $\mathbf{y}\left({\mathbf{neq}}\right)$Real (Kind=nag_wp) array Input/Output
On entry: $y$, the values of the dependent variables evaluated at $t$.
On exit: these values must not be changed unless imon is set to $2$.
7: $\mathbf{ydot}\left({\mathbf{neq}}\right)$Real (Kind=nag_wp) array Input
On entry: the time derivatives ${y}^{\prime }$ of the vector $y$.
8: $\mathbf{ysav}\left({\mathbf{ldysav}},*\right)$Real (Kind=nag_wp) array Communication Array
Note: the second dimension of ysav is sdysav as in the call of d02nhf.
On entry: workspace to enable you to carry out interpolation using either of the routines d02xjf or d02xkf.
9: $\mathbf{r}\left({\mathbf{neq}}\right)$Real (Kind=nag_wp) array Input
On entry: if ${\mathbf{imon}}=0$ and ${\mathbf{inln}}=3$, the first neq elements contain the residual vector $A\left(t,y\right){y}^{\prime }-g\left(t,y\right)$.
10: $\mathbf{acor}\left({\mathbf{neq}},2\right)$Real (Kind=nag_wp) array Input
On entry: with ${\mathbf{imon}}=1$, ${\mathbf{acor}}\left(i,1\right)$ contains the weight used for the $i$th equation when the norm is evaluated, and ${\mathbf{acor}}\left(i,2\right)$ contains the estimated local error for the $i$th equation. The scaled local error at the end of a timestep may be obtained by calling d02zaf as follows:
```ifail = 1
errloc = d02zaf(neq, acor(1,2), acor(1,1), ifail)
! CHECK IFAIL BEFORE PROCEEDING```
11: $\mathbf{imon}$Integer Input/Output
On entry: a flag indicating under what circumstances monitr was called:
${\mathbf{imon}}=-2$
Entry from the integrator after ${\mathbf{ires}}=4$ (set in resid) caused an early termination (this facility could be used to locate discontinuities).
${\mathbf{imon}}=-1$
The current step failed repeatedly.
${\mathbf{imon}}=0$
Entry after a call to the internal nonlinear equation solver (see inln).
${\mathbf{imon}}=1$
The current step was successful.
On exit: may be reset to determine subsequent action in d02nhf.
${\mathbf{imon}}=-2$
Integration is to be halted. A return will be made from the integrator to the calling (sub)program with ${\mathbf{ifail}}={\mathbf{12}}$.
${\mathbf{imon}}=-1$
Allow the integrator to continue with its own internal strategy. The integrator will try up to three restarts unless ${\mathbf{imon}}\ne -1$ on exit.
${\mathbf{imon}}=0$
Return to the internal nonlinear equation solver, where the action taken is determined by the value of inln (see inln).
${\mathbf{imon}}=1$
Normal exit to the integrator 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 monitr, 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 the call to monitr. 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.
12: $\mathbf{inln}$Integer Output
On exit: the action to be taken by the internal nonlinear equation solver when monitr is exited with ${\mathbf{imon}}=0$. By setting ${\mathbf{inln}}=3$ and returning to the integrator, the residual vector is evaluated and placed in the array r, and then monitr is called again. At present this is the only option available: inln must not be set to any other value.
13: $\mathbf{hmin}$Real (Kind=nag_wp) Input/Output
On entry: the minimum step size to be taken on the next step.
On exit: the minimum step size to be used. If this is different from the input value, imon must be set to $3$ or $4$.
14: $\mathbf{hmax}$Real (Kind=nag_wp) Input/Output
On entry: the maximum step size to be taken on the next step.
On exit: the maximum step size to be used. If this is different from the input value, imon must be set to $3$ or $4$. If hmax is set to zero, no limit is assumed.
15: $\mathbf{nqu}$Integer Input
On entry: the order of the integrator used on the last step. This is supplied to enable you to carry out interpolation using either of the routines d02xjf or d02xkf.
monitr must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which d02nhf is called. Arguments denoted as Input must not be changed by this procedure.
Note: monitr should not return floating-point NaN (Not a Number) or infinity values, since these are not handled by d02nhf. If your code inadvertently does return any NaNs or infinities, d02nhf is likely to produce unexpected results.
21: $\mathbf{lderiv}\left(2\right)$Logical array Input/Output
On entry: ${\mathbf{lderiv}}\left(1\right)$ must be set to .TRUE. if you have supplied both an initial $y$ and an initial ${y}^{\prime }$. ${\mathbf{lderiv}}\left(1\right)$ must be set to .FALSE. if only the initial $y$ has been supplied.
${\mathbf{lderiv}}\left(2\right)$ must be set to .TRUE. if the integrator is to use a modified Newton method to evaluate the initial $y$ and ${y}^{\prime }$. Note that $y$ and ${y}^{\prime }$, if supplied, are used as initial estimates. This method involves taking a small step at the start of the integration, and if ${\mathbf{itask}}=6$ on entry, t and tout will be set to the result of taking this small step. ${\mathbf{lderiv}}\left(2\right)$ must be set to .FALSE. if the integrator is to use functional iteration to evaluate the initial $y$ and ${y}^{\prime }$, and if this fails a modified Newton method will then be attempted. ${\mathbf{lderiv}}\left(2\right)=\mathrm{.TRUE.}$ is recommended if there are implicit equations or the initial $y$ and ${y}^{\prime }$ are zero.
On exit: ${\mathbf{lderiv}}\left(1\right)$ is normally unchanged. However if ${\mathbf{itask}}=6$ and internal initialization was successful then ${\mathbf{lderiv}}\left(1\right)=\mathrm{.TRUE.}$.
${\mathbf{lderiv}}\left(2\right)=\mathrm{.TRUE.}$, if implicit equations were detected. Otherwise ${\mathbf{lderiv}}\left(2\right)=\mathrm{.FALSE.}$.
22: $\mathbf{itask}$Integer Input
On 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 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. tcrit must be specified as under ${\mathbf{itask}}=4$.
${\mathbf{itask}}=6$
The integrator will solve for the initial values of $y$ and ${y}^{\prime }$ only and then return to the calling (sub)program without doing the integration. This option can be used to check the initial values of $y$ and ${y}^{\prime }$. Functional iteration or a ‘small’ backward Euler method used in conjunction with a damped Newton iteration is used to calculate these values (see lderiv). Note that if a backward Euler step is used then the value of $t$ will have been advanced a short distance from the initial point.
Note:  if d02nhf is recalled with a different value of itask (and tout altered), the initialization procedure is repeated, possibly leading to different initial conditions.
Constraint: $1\le {\mathbf{itask}}\le 6$.
23: $\mathbf{itrace}$Integer Input
On 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}$Integer Input/Output
On entry: ifail must be set to $0$, $-1$ or $1$ to set behaviour on detection of an error; these values have no effect when no error is detected.
A value of $0$ causes the printing of an error message and program execution will be halted; otherwise program execution continues. A value of $-1$ means that an error message is printed while a value of $1$ means that it is not.
If halting is not appropriate, the value $-1$ or $1$ is recommended. If message printing is undesirable, then the value $1$ is recommended. Otherwise, the value $-1$ is recommended since useful values can be provided in some output arguments even when ${\mathbf{ifail}}\ne {\mathbf{0}}$ on exit. When the value $-\mathbf{1}$ 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).

## 6Error 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:
Note: in some cases d02nhf may return useful information.
${\mathbf{ifail}}=1$
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}}=⟨\mathit{\text{value}}⟩$ and ${\mathbf{njcpvt}}=⟨\mathit{\text{value}}⟩$ in d02ntf.
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}}=⟨\mathit{\text{value}}⟩$ and ${\mathbf{nwkjac}}=⟨\mathit{\text{value}}⟩$ in d02ntf.
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}}=⟨\mathit{\text{value}}⟩$, sdysav (setup) $=⟨\mathit{\text{value}}⟩$.
Failure during internal time interpolation. tout and the current time are too close.
${\mathbf{itask}}=⟨\mathit{\text{value}}⟩$ and ${\mathbf{tout}}=⟨\mathit{\text{value}}⟩$ and the current time is $⟨\mathit{\text{value}}⟩$.
ires was set to an illegal value during initialization.
${\mathbf{ires}}=⟨\mathit{\text{value}}⟩$ at time $⟨\mathit{\text{value}}⟩$.
${\mathbf{itask}}=3$ and tout is more than an integration step behind the current time.
${\mathbf{tout}}=⟨\mathit{\text{value}}⟩$, current time minus step size: $⟨\mathit{\text{value}}⟩$.
monitr appears to have overwritten the solution vector.
Further integration will not be attempted.
monitr set ${\mathbf{imon}}=⟨\mathit{\text{value}}⟩$.
Constraint: $-2\le {\mathbf{imon}}\le 4$.
monitr set ${\mathbf{inln}}=⟨\mathit{\text{value}}⟩$.
Constraint: ${\mathbf{inln}}=3$.
On entry, an illegal (negative) maximum number of steps was provided in a prior call to a setup routine. ${\mathbf{maxstp}}=⟨\mathit{\text{value}}⟩$.
On entry, an illegal (negative) maximum stepsize was provided in a prior call to a setup routine. ${\mathbf{hmax}}=⟨\mathit{\text{value}}⟩$.
On entry, an illegal (negative) minimum stepsize was provided in a prior call to a setup routine. ${\mathbf{hmin}}=⟨\mathit{\text{value}}⟩$.
On entry, ${\mathbf{atol}}=⟨\mathit{\text{value}}⟩$.
Constraint: ${\mathbf{atol}}\ge 0.0$.
On entry, ${\mathbf{ires}}=⟨\mathit{\text{value}}⟩$ and $\frac{dy}{dt}=0.0$ for all elements.
Check the evaluation of the residual for this value of ires.
On entry, ${\mathbf{itask}}=⟨\mathit{\text{value}}⟩$.
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}}=⟨\mathit{\text{value}}⟩$, ${\mathbf{tcrit}}=⟨\mathit{\text{value}}⟩$ and the current time is $⟨\mathit{\text{value}}⟩$.
On entry, ${\mathbf{itask}}=4$ or $5$ and tcrit is before tout in the direction of integration.
${\mathbf{itask}}=⟨\mathit{\text{value}}⟩$, ${\mathbf{tcrit}}=⟨\mathit{\text{value}}⟩$ and ${\mathbf{tout}}=⟨\mathit{\text{value}}⟩$.
On entry, ${\mathbf{itol}}=⟨\mathit{\text{value}}⟩$.
Constraint: $1\le {\mathbf{itol}}\le 4$.
On entry, ${\mathbf{neq}}=⟨\mathit{\text{value}}⟩$.
Constraint: ${\mathbf{neq}}\ge 1$.
On entry, ${\mathbf{neq}}=⟨\mathit{\text{value}}⟩$ and ${\mathbf{ldysav}}=⟨\mathit{\text{value}}⟩$.
Constraint: ${\mathbf{neq}}\le {\mathbf{ldysav}}$.
On entry, ${\mathbf{rtol}}=⟨\mathit{\text{value}}⟩$.
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}}=⟨\mathit{\text{value}}⟩$, ${\mathbf{t}}=⟨\mathit{\text{value}}⟩$ and ${\mathbf{h0}}=⟨\mathit{\text{value}}⟩$.
On entry, tout is too close to t to start integration.
${\mathbf{tout}}=⟨\mathit{\text{value}}⟩$ and ${\mathbf{t}}=⟨\mathit{\text{value}}⟩$.
The initial stepsize, $⟨\mathit{\text{value}}⟩$, is too small.
Weight number $i=⟨\mathit{\text{value}}⟩$ 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=⟨\mathit{\text{value}}⟩$.
${\mathbf{ifail}}=2$
At time $⟨\mathit{\text{value}}⟩$ the maximum number of allowed steps on this call was taken before reaching the next output point ${\mathbf{tout}}=⟨\mathit{\text{value}}⟩$.
Maximum number of steps $=⟨\mathit{\text{value}}⟩$.
${\mathbf{ifail}}=3$
Too much accuracy requested for precision of the machine at time $⟨\mathit{\text{value}}⟩$.
The tolerances should be checked; the requested accuracy should be reduced by a factor of at least $⟨\mathit{\text{value}}⟩$.
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 error-test 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 convergence-test 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 $⟨\mathit{\text{value}}⟩$, error weight $⟨\mathit{\text{value}}⟩$ became zero. Check the values of atol, rtol and itol supplied.
${\mathbf{ifail}}=7$
resid set ${\mathbf{ires}}=3$, which signals that an error condition has occurred in the solution vector, its time derivative or in the value of $t$. It was not possible to remove this condition.
${\mathbf{ires}}=⟨\mathit{\text{value}}⟩$ at $t=⟨\mathit{\text{value}}⟩$.
${\mathbf{ifail}}=8$
Attempt was made to reduce the step size to a value less than the minimum step size during the calculation of initial values.
Minimum stepsize: $⟨\mathit{\text{value}}⟩$.
Nonlinear solver failed to converge using a damped Newton method to solve for initial values.
Damping factor: $⟨\mathit{\text{value}}⟩$; convergence rate: $⟨\mathit{\text{value}}⟩$.
The residual routine returned an error when calculating the initial values of the solution and its time derivative.
The user problem has one or more inconsistencies between the ${\mathbf{ires}}=1$ and ${\mathbf{ires}}=-1$ parts. Integration will not be attempted.
Workspace error occurred when trying to form the Jacobian matrix in calculating the initial values of the solution and its time derivative.
${\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: $⟨\mathit{\text{value}}⟩$; required: $⟨\mathit{\text{value}}⟩$.
Not enough integer store provided for internal storage pointers.
Units of store needed: $⟨\mathit{\text{value}}⟩$. Amount provided: $⟨\mathit{\text{value}}⟩$.
Not enough integer store provided for sparse matrix solver.
Units of store needed: $⟨\mathit{\text{value}}⟩$. Amount provided: $⟨\mathit{\text{value}}⟩$.
Not enough real store provided for sparse matrix solver.
Units of store needed: $⟨\mathit{\text{value}}⟩$. Amount provided: $⟨\mathit{\text{value}}⟩$.
${\mathbf{ifail}}=11$
resid set ${\mathbf{ires}}=2$, which signals that the integration should terminate. ${\mathbf{ires}}=⟨\mathit{\text{value}}⟩$ at time $⟨\mathit{\text{value}}⟩$.
${\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 $⟨\mathit{\text{value}}⟩$.
${\mathbf{ifail}}=15$
Either the banded matrix linear algebra setup routine was not called first or a communication array has become corrupted.
${\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.

## 7Accuracy

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

## 8Parallelism and Performance

Background information to multithreading can be found in the Multithreading documentation.
d02nhf 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.
d02nhf is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.
d02nhf 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 implementation-specific information.

## 9Further Comments

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 d02nhf the cost is proportional to ${\mathbf{neq}}×{\left({\mathbf{ml}}+{\mathbf{mu}}+1\right)}^{2}$, though for problems which are only mildly nonlinear the cost may be dominated by factors proportional to ${\mathbf{neq}}×\left({\mathbf{ml}}+{\mathbf{mu}}+1\right)$ except for very large problems.
In general, you are advised to choose the BDF option (setup routine d02nvf) but if efficiency is of great importance and especially if it is suspected that $\frac{\partial }{\partial y}\left({A}^{-1}g\right)$ has complex eigenvalues near the imaginary axis for some part of the integration, you should try the BLEND option (setup routine d02nwf).

## 10Example

This example solves the well-known stiff Robertson problem written as an implicit differential system and in implicit form
 $r1 = a′ + b′ + c′ r2 = 0.04a - 1.0E4bc - 3.0E7⁢b2 - b′ r3 = 3.0E7⁢b2 - c′$
exploiting the fact that we can show that ${\left(a+b+c\right)}^{\prime }=0$ for all time. Integration is over the range $\left[0,10.0\right]$ with initial conditions $a=1.0$ and $b=c=0.0$ using scalar relative error control and vector absolute error control (${\mathbf{itol}}=2$). We integrate using a BDF method (setup routine d02nvf) and a modified Newton method. The Jacobian is calculated numerically and we employ a default monitor, dummy routine d02nby. We perform a normal integration (${\mathbf{itask}}=1$) to obtain the value at ${\mathbf{tout}}=10.0$ by integrating past tout and interpolating. We also illustrate the use of ${\mathbf{itask}}=6$ to calculate initial values of $y$ and ${y}^{\prime }$ and then return without integrating further.

### 10.1Program Text

Program Text (d02nhfe.f90)

### 10.2Program Data

Program Data (d02nhfe.d)

### 10.3Program Results

Program Results (d02nhfe.r)