D02NDF is a direct communication routine for integrating stiff systems of explicit ordinary differential equations when the Jacobian is a sparse matrix.
D02NDF is a general purpose routine for integrating the initial value problem for a stiff system of explicit ordinary differential equations,
It is designed specifically for the case where the Jacobian
$\frac{\partial g}{\partial y}$ is a sparse matrix.
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 D02NDF 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 linear algebra diagnostic routine
D02NXF.
! Declarations
EXTERNAL FCN, 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 D02NDF(NEQ, NEQMAX, T, TOUT, Y, YDOT, RWORK, RTOL, &
ATOL, ITOL, INFORM, FCN, YSAVE, NY2DIM, JAC, &
WKJAC,NWKJAC, JACPVT, NJCPVT, MONITR, 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,
D02NVF or
D02NWF, must be called prior to the call of D02NDF. 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 D02NDF. There is also a routine,
D02NZF, designed to permit you to change step size on a continuation call to D02NDF 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).
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).
D02NDF 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.
D02NDF is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.
D02NDF 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.
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 indication of numerical stability.
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).
This example solves the well-known stiff Robertson problem
over the range
$\left[0,10.0\right]$ with initial conditions
$a=1.0$ and
$b=c=0.0$ using scalar error control (
${\mathbf{ITOL}}=1$). The solution is computed up to
$10.0$ by overshooting and interpolating (
${\mathbf{ITASK}}=1$) and the intermediate solution computed on an equispaced mesh through
MONITR. The integration algorithm used is the BDF method (setup routine
D02NVF) and a modified Newton method is also used. The use of the 'N' (Numerical) and 'S' (Structural) options are illustrated in turn for calculating the Jacobian.