D02NDF (PDF version)
D02 Chapter Contents
D02 Chapter Introduction
NAG Library Manual

NAG Library Routine Document

D02NDF

Note:  before using this routine, please read the Users' Note for your implementation to check the interpretation of bold italicised terms and other implementation-dependent details.

+ Contents

    1  Purpose
    7  Accuracy

1  Purpose

D02NDF is a forward communication routine for integrating stiff systems of explicit ordinary differential equations when the Jacobian is a sparse matrix.

2  Specification

SUBROUTINE D02NDF ( NEQ, LDYSAV, T, TOUT, Y, YDOT, RWORK, RTOL, ATOL, ITOL, INFORM, FCN, YSAV, SDYSAV, JAC, WKJAC, NWKJAC, JACPVT, NJCPVT, MONITR, 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)
EXTERNAL  FCN, JAC, MONITR

3  Description

D02NDF is a general purpose routine for integrating the initial value problem for a stiff system of explicit ordinary differential equations,
y = gt,y .
It is designed specifically for the case where the Jacobian g 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.

4  References

See the D02M–N sub-chapter Introduction.

5  Parameters

1:     NEQ – INTEGERInput
On entry: the number of differential equations to be solved.
Constraint: NEQ1.
2:     LDYSAV – INTEGERInput
On entry: a bound on the maximum number of differential equations to be solved during the integration.
Constraint: LDYSAVNEQ.
3:     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:     TOUT – REAL (KIND=nag_wp)Input
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: TOUTT.
5:     Y(NEQ) – REAL (KIND=nag_wp) arrayInput/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 T=TOUT).
6:     YDOT(NEQ) – REAL (KIND=nag_wp) arrayOutput
On exit: the time derivatives y of the vector y at the last integration point.
7:     RWORK(50+4×NEQ) – REAL (KIND=nag_wp) arrayCommunication Array
8:     RTOL(*) – REAL (KIND=nag_wp) arrayInput
Note: the dimension of the array RTOL must be at least 1 if ITOL=1 or 2, and at least NEQ otherwise.
On entry: the relative local error tolerance.
Constraint: RTOLi0.0 for all relevant i (see ITOL).
9:     ATOL(*) – REAL (KIND=nag_wp) arrayInput
Note: the dimension of the array ATOL must be at least 1 if ITOL=1 or 3, and at least NEQ otherwise.
On entry: the absolute local error tolerance.
Constraint: ATOLi0.0 for all relevant i (see ITOL).
10:   ITOL – INTEGERInput
On entry: a value to indicate the form of the local error test. ITOL indicates to D02NDF whether to interpret either or both of RTOL or ATOL as a vector or a scalar. The error test to be satisfied is ei/wi<1.0, where wi is defined as follows:
ITOL RTOL ATOL wi
1 scalar scalar RTOL1×yi+ATOL1
2 scalar vector RTOL1×yi+ATOLi
3 vector scalar RTOLi×yi+ATOL1
4 vector vector RTOLi×yi+ATOLi
ei is an estimate of the local error in yi, computed internally, and the choice of norm to be used is defined by a previous call to an integrator setup routine.
Constraint: ITOL=1, 2, 3 or 4.
11:   INFORM(23) – INTEGER arrayCommunication Array
12:   FCN – SUBROUTINE, supplied by the user.External Procedure
FCN must evaluate the derivative vector for the explicit ordinary differential equation system, defined by y=gt,y.
The specification of FCN is:
SUBROUTINE FCN ( NEQ, T, Y, F, IRES)
INTEGER  NEQ, IRES
REAL (KIND=nag_wp)  T, Y(NEQ), F(NEQ)
1:     NEQ – INTEGERInput
On entry: the number of differential 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: the value of yi, for i=1,2,,NEQ.
4:     F(NEQ) – REAL (KIND=nag_wp) arrayOutput
On exit: the value yi, given by yi=git,y, for i=1,2,,NEQ.
5:     IRES – INTEGERInput/Output
On entry: IRES=1.
On exit: you may set IRES as follows to indicate certain conditions in FCN to the integrator:
IRES=1
Indicates a normal return from FCN, that is IRES has not been altered by you and integration continues.
IRES=2
Indicates to the integrator that control should be passed back immediately to the calling (sub)program with the error indicator set to IFAIL=11.
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 IFAIL=7.
IRES=4
Indicates to the integrator to stop its current operation and to enter MONITR immediately with parameter IMON=-2.
FCN must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which D02NDF is called. Parameters denoted as Input must not be changed by this procedure.
13:   YSAV(LDYSAV,SDYSAV) – REAL (KIND=nag_wp) arrayCommunication Array
14:   SDYSAV – INTEGERInput
On entry: the second dimension of the array YSAV as declared in the (sub)program from which D02NDF is called. An appropriate value for SDYSAV is described in the specification of the integrator setup routines D02NVF and D02NWF. This value must be the same as that supplied to the integrator setup routine.
15:   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 D02NDZ. (D02NDZ 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 sparse linear algebra setup routine D02NUF.
First we must define the system of nonlinear equations which is solved internally by the integrator. The time derivative, y, generated internally, has the form
y = y-z / hd ,
where h is the current step size and d is a parameter 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 d dy ​ ​ = hd d dy ​ ​ . The system of nonlinear equations that is solved has the form
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 y-z / hd - g t,y .
It is the Jacobian matrix r y  that you must supply in JAC as follows:
ri yj =1-hd gi yj , if ​i=j, ri yj = -hd gi yj , otherwise.
The specification of JAC is:
SUBROUTINE JAC ( NEQ, T, Y, H, D, J, PDJ)
INTEGER  NEQ, J
REAL (KIND=nag_wp)  T, Y(NEQ), H, D, PDJ(NEQ)
1:     NEQ – INTEGERInput
On entry: the number of differential 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: yi, for i=1,2,,NEQ, the current solution component.
4:     H – REAL (KIND=nag_wp)Input
On entry: the current step size.
5:     D – REAL (KIND=nag_wp)Input
On entry: the parameter d which depends on the integration method.
6:     J – INTEGERInput
On entry: the column of the Jacobian that JAC must return in the array PDJ.
7:     PDJ(NEQ) – REAL (KIND=nag_wp) arrayInput/Output
On entry: is set to zero.
On exit: PDJi should be set to the i,jth element of the Jacobian, where j is given by J. 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 D02NDF is called. Parameters denoted as Input must not be changed by this procedure.
16:   WKJAC(NWKJAC) – REAL (KIND=nag_wp) arrayCommunication Array
17:   NWKJAC – INTEGERInput
On entry: the dimension of the array WKJAC as declared in the (sub)program from which D02NDF is called. The actual size depends on whether the sparsity structure is supplied or whether it is to be estimated. An appropriate value for NWKJAC is described in the specification of the linear algebra setup routine D02NUF. This value must be the same as that supplied to D02NUF.
18:   JACPVT(NJCPVT) – INTEGER arrayCommunication Array
19:   NJCPVT – INTEGERInput
On entry: the dimension of the array JACPVT as declared in the (sub)program from which D02NDF is called. The actual size depends on whether the sparsity structure is supplied or whether it is to be estimated. An appropriate value for NJCPVT is described in the specification of the linear algebra setup routine D02NUF. This value must be the same as that supplied to D02NUF.
20:   MONITR – SUBROUTINE, supplied by the NAG Library or the user.External Procedure
MONITR performs tasks requested by you. If this option is not required, then the actual argument for MONITR must be the dummy routine D02NBY. (D02NBY is included in the NAG Library.)
The specification of MONITR is:
SUBROUTINE MONITR ( NEQ, LDYSAV, T, HLAST, HNEXT, Y, YDOT, YSAV, R, ACOR, IMON, INLN, HMIN, HMAX, NQU)
INTEGER  NEQ, LDYSAV, IMON, INLN, NQU
REAL (KIND=nag_wp)  T, HLAST, HNEXT, Y(NEQ), YDOT(NEQ), YSAV(LDYSAV,sdysav), R(NEQ), ACOR(NEQ,2), HMIN, HMAX
where sdysav is the numerical value of SDYSAV in the call of D02NDF.
1:     NEQ – INTEGERInput
On entry: the number of differential equations being solved.
2:     LDYSAV – INTEGERInput
On entry: an upper bound on the number of differential equations to be solved.
3:     T – REAL (KIND=nag_wp)Input
On entry: the current value of the independent variable.
4:     HLAST – REAL (KIND=nag_wp)Input
On entry: the last step size successfully used by the integrator.
5:     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, then IMON must be set to 4.
6:     Y(NEQ) – REAL (KIND=nag_wp) arrayInput/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:     YDOT(NEQ) – REAL (KIND=nag_wp) arrayInput
On entry: the time derivatives y of the vector y.
8:     YSAV(LDYSAV,sdysav) – REAL (KIND=nag_wp) arrayInput
On entry: workspace to enable you to carry out interpolation using either of the routines D02XJF or D02XKF.
9:     R(NEQ) – REAL (KIND=nag_wp) arrayInput
On entry: if IMON=0 and INLN=3, the first NEQ elements contain the residual vector, y-gt,y.
10:   ACOR(NEQ,2) – REAL (KIND=nag_wp) arrayInput
On entry: with IMON=1, ACORi1 contains the weight used for the ith equation when the norm is evaluated, and ACORi2 contains the estimated local error for the ith equation. The scaled local error at the end of a timestep may be obtained by calling the real function D02ZAF as follows:
 IFAIL = 1 ERRLOC = D02ZAF(NEQ, ACOR(1,2), ACOR(1,1), IFAIL) ! CHECK IFAIL BEFORE PROCEEDING 
11:   IMON – INTEGERInput/Output
On entry: a flag indicating under what circumstances MONITR was called:
IMON=-2
Entry from the integrator after IRES=4 (set in FCN) caused an early termination (this facility could be used to locate discontinuities).
IMON=-1
The current step failed repeatedly.
IMON=0
Entry after a call to the internal nonlinear equation solver (see INLN).
IMON=1
The current step was successful.
On exit: may be reset to determine subsequent action in D02NDF.
IMON=-2
Integration is to be halted. A return will be made from the integrator to the calling (sub)program with IFAIL=12.
IMON=-1
Allow the integrator to continue with its own internal strategy. The integrator will try up to three restarts unless IMON is set -1 on exit.
IMON=0
Return to the internal nonlinear equation solver, where the action taken is determined by the value of INLN (see INLN).
IMON=1
Normal exit to the integrator to continue integration.
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.
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.
IMON=4
Continue the integration but using a new value of HNEXT and possibly new values of HMIN and HMAX.
12:   INLN – INTEGEROutput
On exit: the action to be taken by the internal nonlinear equation solver when MONITR is exited with IMON=0. By setting 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:   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, then IMON must be set to 3 or 4.
14:   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, then IMON must be set to 3 or 4. If HMAX is set to zero, no limit is assumed.
15:   NQU – INTEGERInput
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 D02NDF is called. Parameters denoted as Input must not be changed by this procedure.
21:   ITASK – INTEGERInput
On entry: the task to be performed by the integrator.
ITASK=1
Normal computation of output values of yt at t=TOUT (by overshooting and interpolating).
ITASK=2
Take one step only and return.
ITASK=3
Stop at the first internal integration point at or beyond t=TOUT and return.
ITASK=4
Normal computation of output values of yt at t=TOUT but without overshooting t=TCRIT (e.g., see D02MVF). 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.
ITASK=5
Take one step only and return, without passing TCRIT (e.g., see D02MVF). TCRIT must be specified as under ITASK=4.
Constraint: ITASK=1, 2, 3, 4 or 5.
22:   ITRACE – INTEGERInput
On entry: the level of output that is printed by the integrator. ITRACE may take the value -1, 0, 1, 2 or 3.
ITRACE<-1
-1 is assumed and similarly if ITRACE>3, then 3 is assumed.
ITRACE=-1
No output is generated.
ITRACE=0
Only warning messages are printed on the current error message unit (see X04AAF).
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.
23:   IFAIL – INTEGERInput/Output
On entry: IFAIL must be set to 0, -1​ 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​ 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 IFAIL0 on exit, the recommended value is -1. When the value -1​ or ​1 is used it is essential to test the value of IFAIL on exit.
On exit: IFAIL=0 unless the routine detects an error or a warning has been flagged (see Section 6).

6  Error Indicators and Warnings

If on entry 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:
IFAIL=1
An illegal input was detected on entry, or after an internal call to MONITR. If ITRACE>-1, then the form of the error will be detailed on the current error message unit (see X04AAF).
IFAIL=2
The maximum number of steps specified has been taken (see the description of optional inputs in the integrator setup routines and the optional input continuation routine, D02NZF).
IFAIL=3
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 Y1,Y2,,YNEQ contain the computed values of the solution at the current point T.
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.
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.
IFAIL=6
Some error weight wi became zero during the integration (see the description of ITOL). Pure relative error control (ATOLi=0.0) was requested on a variable (the ith) which has now vanished. The integration was successful as far as T.
IFAIL=7
FCN set its error flag (IRES=3) continually despite repeated attempts by the integrator to avoid this.
IFAIL=8
Not used for the integrator.
IFAIL=9
A singular Jacobian r y  has been encountered. This error exit is unlikely to be taken when solving explicit ordinary differential equations. You should check the problem formulation and Jacobian calculation.
IFAIL=10
An error occurred during Jacobian formulation or back-substitution (a more detailed error description may be directed to the current error message unit, see X04AAF).
IFAIL=11
FCN signalled the integrator to halt the integration and return (IRES=2). Integration was successful as far as T.
IFAIL=12
MONITR set IMON=-2 and so forced a return but the integration was successful as far as T.
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. (Only applies when you are not operating in one step mode, that is when ITASK2 or 5.)
IFAIL=14
The values of RTOL and ATOL are so small that D02NDF is unable to start the integration.
IFAIL=15
The linear algebra setup routine D02NUF was not called prior to calling D02NDF.

7  Accuracy

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 ITOL=1 with ATOL1 small but positive).

8  Further Comments

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 g y  has complex eigenvalues near the imaginary axis for some part of the integration, you should try the BLEND option (setup routine D02NWF).

9  Example

This example solves the well-known stiff Robertson problem
a = -0.04a+1.0E4bc b = -0.04a-1.0E4bc -3.0E7b2 c = -3.0E7b2
over the range 0,10.0 with initial conditions a=1.0 and b=c=0.0 using scalar error control (ITOL=1). The solution is computed up to 10.0 by overshooting and interpolating (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.

9.1  Program Text

Program Text (d02ndfe.f90)

9.2  Program Data

Program Data (d02ndfe.d)

9.3  Program Results

Program Results (d02ndfe.r)

Produced by GNUPLOT 4.4 patchlevel 0 0 0.2 0.4 0.6 0.8 1 0 2 4 6 8 10 0 5e-06 1e-05 1.5e-05 2e-05 2.5e-05 3e-05 3.5e-05 4e-05 Solution (a,c) Solution (b) x Example Program Stiff Robertson Problem BDF Method, Sparse Numerical Jacobian with Structure Supplied a b c

D02NDF (PDF version)
D02 Chapter Contents
D02 Chapter Introduction
NAG Library Manual

© The Numerical Algorithms Group Ltd, Oxford, UK. 2012