d02njf 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
It is designed specifically for the case where the resulting Jacobian is a sparse 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
d02njf is given below. It calls
the sparse matrix linear algebra setup routine
d02nuf,
the Backward Differentiation Formula (BDF) integrator setup
routine
d02nvf, its diagnostic counterpart
d02nyf,
and
the sparse matrix linear algebra diagnostic routine
d02nxf.
! Declarations
External resid, jac, monitr
.
.
.
ifail = 0
Call d02nvf(...,ifail)
Call d02nuf(neq, neqmax, jceval, nwkjac, ia, nia, ja, nja, &
jacpvt, njcpvt, sens, u, eta, lblock, isplit, &
rwork, ifail)
ifail = 1
Call d02njf(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 d02nxf(...)
Call d02nyf(...)
.
.
.
Stop
End
The linear algebra setup routine
d02nuf and
one of
the integrator setup
routines,
d02mvf,
d02nvf or
d02nwf,
must be called prior to the call of
d02njf.
Either or both of the integrator diagnostic routine
d02nyf, or the
sparse matrix linear algebra diagnostic routine
d02nxf, may be called after the call to
d02njf.
There is also a routine,
d02nzf, designed to permit you to change step size on a continuation call to
d02njf without restarting the integration process.
If on entry
${\mathbf{ifail}}=0$ or
$1$, explanatory error messages are output on the current error message unit (as defined by
x04aaf).
 ${\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$ and
${\mathbf{njcpvt}}=\u2329\mathit{\text{value}}\u232a$ in
d02nuf.
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$ and
${\mathbf{nwkjac}}=\u2329\mathit{\text{value}}\u232a$ in
d02nuf.
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$.
ires was set to an illegal value during initialization.
${\mathbf{ires}}=\u2329\mathit{\text{value}}\u232a$ at time
$\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$.
monitr appears to have overwritten the solution vector.
Further integration will not be attempted.
monitr set
${\mathbf{imon}}=\u2329\mathit{\text{value}}\u232a$.
Constraint:
$2\le {\mathbf{imon}}\le 4$.
monitr set
${\mathbf{inln}}=\u2329\mathit{\text{value}}\u232a$.
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}}=\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,
${\mathbf{ires}}=\u2329\mathit{\text{value}}\u232a$ and
$\frac{dy}{dt}=0.0$ for all elements.
Check the evaluation of the residual for this value of
ires.
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$.
The initial stepsize, $h=\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$.
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).
d02njf 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.
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.
Since numerical stability and memory are often conflicting requirements when solving ordinary differential systems where the Jacobian matrix is sparse we provide a diagnostic routine,
d02nxf, whose aim is to inform you how much memory is required to solve the problem and to give you some indicators of numerical stability.
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).
This example solves the wellknown stiff Robertson problem written as a mixed differential/algebraic system in implicit form
exploiting the fact that, from the initial conditions
$a=1.0$ and
$b=c=0.0$, we know that
$a+b+c=1$ for all time. We integrate over the range
$\left[0,10.0\right]$ with vector relative error control and scalar absolute error control (
${\mathbf{itol}}=3$) and using the BDF integrator (setup routine
d02nvf)
and a modified Newton method.
The Jacobian is evaluated, in turn, using the 'A' (Analytical) and 'F' (Full information) options.
We provide a monitor routine to terminate the integration when the value of the component a falls below
$0.9$.