NAG Library Sub-chapter Introduction

D02M–N – Integrators for Stiff Ordinary Differential Systems

This sub-chapter contains the specifications of the integrators, the setup routines and diagnostic routines which have been developed from the SPRINT package, Berzins and Furzeland (1985), and from the DASSL package, Brenan et al. (1996).

The SPRINT explicit integrators D02NBF, D02NCF and D02NDF are designed for solving stiff systems of explicitly defined ordinary differential equations,

The SPRINT implicit integrators D02NGF, D02NHF and D02NJF are designed for solving stiff systems of implicitly defined ordinary differential equations,

The DASSL integrator D02NEF is designed for solving systems of the form, $F\left(t,y,y\prime \right)=0$. These formulations permits solution of differential/algebraic systems (DAEs). The facilities provided are essentially those of the explicit solvers.

The SPRINT integrator routines have almost identical calling sequences but each is designed to solve a problem where the Jacobian is of a particular structure: full matrix (D02NBF and D02NGF), banded matrix (D02NCF and D02NHF) or sparse matrix (D02NDF and D02NJF). Each of these structures has associated with it a linear algebra setup routine: D02NSF, D02NTF and D02NUF respectively. A linear algebra setup routine must be called before the first call to the appropriate integrator. These linear algebra setup routines check various parameters of the corresponding integrator routine and set certain parameters for the linear algebra computations. A routine, D02NXF, is supplied which permits extraction of diagnostic information after a call to either of the sparse linear algebra solvers D02NDF and D02NJF.

With the SPRINT integrators are also associated three integrator setup routines D02MVF, D02NVF and D02NWF, one of which must be called before the first call to any integrator routine. They provide input to the Backward Differentiation Formulae (BDF), the Blend Formulae and the special fixed leading coefficient BDF codes respectively. On return from an integrator, if it is feasible to continue the integration, D02NZF may be called to reset various integration parameters. It is often of considerable interest to determine statistics concerning the integration process. D02NYF is provided with this aim in mind. It should prove especially useful to those who wish to integrate many similar problems as it provides suitable values for many of the input parameters and indications of the difficulties encountered when solving the problem.

Hence, the general form of a program calling one of the integrator routines D02NBF, D02NCF, D02NDF, D02NGF, D02NHF or D02NJF will be

declarations 
   .  
   .  
call linear algebra setup routine 
call integrator setup routine 
call integrator 
call integrator diagnostic routine (if required)  
call linear algebra diagnostic routine (if appropriate and if required)  
   .  
   .  
STOP  
END

The DASSL integrator, D02NEF, has an associated setup routine D02MWF which must be called first. On return from the integrator, if it is feasible to continue the integration, the associated continuation call routine is D02MCF may be called to rest various integration parameters. The structure of the Jacobian is assumed to be full unless D02NPF is called following a call to the setup routine to specify that the Jacobian is banded and to supply its bandwidths.

The required calling sequence for different Jacobian structures and system types is represented diagrammatically in Figure 1.

The integrators D02NMF and D02NNF are reverse communication routines designed for solving explicit and implicit stiff ordinary differential systems respectively. You are warned that you should use these routines only when the integrators mentioned above are inadequate for their application. For example, if it is difficult to write one or more of the user-supplied routine FCN (RESID) or JAC (or MONITR) or if the integrators are to be embedded in a package, it may be advisable to consider these routines.

Since these routines use reverse communication you do not need to define any routines with a prescribed argument list. This makes them especially suitable for large scale computations where encapsulation of the definition of the differential system or its Jacobian matrix in a prescribed routine form may be particularly difficult to achieve.

D02NMF is the reverse communication counterpart of the forward communication routines D02NBF, D02NCF and D02NDF whereas D02NNF is the reverse communication counterpart of the forward communication routines D02NGF, D02NHF and D02NJF. When using these reverse communication routines it is necessary to call the same linear algebra and integrator setup routines as for the forward communication counterpart. All the other continuation and interrogation routines available for use with the forward communication routines are also available to you when calling the reverse communication routines.

There is also a routine, D02NRF, to tell you how to supply the Jacobian when the sparse linear algebra option is being employed with either of D02NMF and D02NNF. Hence, the general form of a program calling one of the integrator routines D02NMF or D02NNF will be

     declarations 
     call linear algebra setup routine 
     call integrator setup routine 
     IREVCM = 0  
1000 call integrator( ..., IREVCM, ...)  
     IF (IREVCM.GT.0) THEN  
     
        evaluate residual and Jacobian (including a call to D02NRF if 
        sparse linear algebra is being used), call the MONITR routine etc.  
    
        GO TO 1000  
     ENDIF  
 
     call integrator diagnostic routine (if required)  
     call linear algebra diagnostic routine (if appropriate and if required)  

     STOP 
     END

The required calling sequence in the case of reverse communication, is represented diagramatically in Figure 2.

In the example programs for the eight SPRINT integrators D02NBF, D02NCF, D02NDF, D02NGF, D02NHF, D02NJF, D02NMF and D02NNF we attempt to illustrate the various options available. Many of these options are available in all the routines and you are invited to scan all the example programs for illustrations of their use. In each case we use as an example the stiff Robertson problem despite the fact that it is not a sensible choice to use either the banded or the sparse linear algebra for this problem. Their use here serves for illustration of the techniques involved. For the implicit integrators D02NGF, D02NHF and D02NJF we write the Robertson problem in residual form, as an implicit differential system and as a differential/algebraic system respectively. Here we are exploiting the fact that $a+b+c$ is constant and hence one of the equations may be replaced by ${\left(a+b+c\right)}^{\prime }=0.0$ or $a+b+c=1.0$ (for our particular choice of initial conditions). For the reverse communication routines D02NMF and D02NNF our examples are intended only to illustrate the reverse communication technique.

The DASSL integrator D02NEF can solve DAEs of the fully implicit form $F\left(t,y,y\prime \right)=0$ and therefore has increased functionality over the SPRINT integrators. Additionally D02NEF can be used to solve difficult algebraic problems by continuation; for example, the nonlinear algebraic problem can be solved by integrating solutions of where the solution to $f\left(x\right)+g\left(x\right)=0$ is known. The solution of this type of problem is illustrated in Section 9 in D02NEF.

Berzins M and Furzeland R M (1985) A user's manual for SPRINT – A versatile software package for solving systems of algebraic, ordinary and partial differential equations: Part 1 – Algebraic and ordinary differential equations Report TNER.85.085 Shell Research Limited

Brenan K, Campbell S and Petzold L (1996) Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations SIAM, Philadelphia