D02NHF is a forward communication routine for integrating stiff systems of implicit ordinary differential equations coupled with algebraic equations when the Jacobian is a banded matrix.
SUBROUTINE D02NHF ( |
NEQ, LDYSAV, T, TOUT, Y, YDOT, RWORK, RTOL, ATOL, ITOL, INFORM, RESID, YSAV, SDYSAV, JAC, WKJAC, NWKJAC, JACPVT, NJCPVT, MONITR, LDERIV, ITASK, ITRACE, IFAIL) |
INTEGER |
NEQ, LDYSAV, ITOL, INFORM(23), SDYSAV, NWKJAC, JACPVT(NJCPVT), NJCPVT, ITASK, ITRACE, IFAIL |
REAL (KIND=nag_wp) |
T, TOUT, Y(NEQ), YDOT(NEQ), RWORK(50+4*NEQ), RTOL(*), ATOL(*), YSAV(LDYSAV,SDYSAV), WKJAC(NWKJAC) |
LOGICAL |
LDERIV(2) |
EXTERNAL |
RESID, JAC, MONITR |
|
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
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.
If on entry
or
, 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 parameters
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
with
small but positive).
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.
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 , though for problems which are only mildly nonlinear the cost may be dominated by factors proportional to 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
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 written as an implicit differential system and in implicit form
exploiting the fact that we can show that
for all time. Integration is over the range
with initial conditions
and
using scalar relative error control and vector absolute error control (
). 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 (
) to obtain the value at
by integrating past
TOUT and interpolating. We also illustrate the use of
to calculate initial values of
and
and then return without integrating further.