NAG CL Interface
d02nec (dae_dassl_gen)
1
Purpose
d02nec is a function for integrating stiff systems of implicit ordinary differential equations coupled with algebraic equations.
2
Specification
void 
d02nec (Integer neq,
double *t,
double tout,
double y[],
double ydot[],
double rtol[],
double atol[],
Integer *itask,
void 
(*res)(Integer neq,
double t,
const double y[],
const double ydot[],
double r[],
Integer *ires,
Nag_Comm *comm),


void 
(*jac)(Integer neq,
double t,
const double y[],
const double ydot[],
double pd[],
double cj,
Nag_Comm *comm),


Integer icom[],
double com[],
Integer lcom,
Nag_Comm *comm,
NagError *fail) 

The function may be called by the names: d02nec, nag_ode_dae_dassl_gen or nag_dae_ivp_dassl_gen.
3
Description
d02nec is a general purpose function for integrating the initial value problem for a stiff system of implicit ordinary differential equations with coupled algebraic equations written in the form
d02nec uses the DASSL implementation of the Backward Differentiation Formulae (BDF) of orders one to five to solve a system of the above form for
$y$ (
y) and
${y}^{\prime}$ (
ydot). Values for
y and
ydot at the initial time must be given as input. These values must be consistent, (i.e., if
t,
y,
ydot are the given initial values, they must satisfy
$F\left({\mathbf{t}},{\mathbf{y}},{\mathbf{ydot}}\right)=0$). The function solves the system from
$t={\mathbf{t}}$ to
$t={\mathbf{tout}}$.
An outline of a typical calling program for
d02nec is given below. It calls the DASSL implementation of the BDF integrator setup function
d02mwc and the banded matrix setup function
d02npc (if required), and, if the integration needs to proceed, calls
d02mcc before continuing the integration.
/* declarations */
EXTERN res, jac
.
.
.
/* Initialize the integrator */
nag_ode_dae_dassl_setup(...);
/* Is the Jacobian matrix banded? */
if (banded) {nag_ode_dae_dassl_linalg(...);}
/* Set dt to the required temporal resolution */
/* Set tend to the final time */
/* Call the integrator for each temporal value */
itask = 0;
while (tout<tend && itask>1) {
nag_ode_dae_dassl_gen (...);
if (itask != 1)
tout = MIN(tout+dt,tend);
/* Print the solution */
}
.
.
.
4
References
None.
5
Arguments

1:
$\mathbf{neq}$ – Integer
Input

On entry: the number of differentialalgebraic equations to be solved.
Constraint:
${\mathbf{neq}}\ge 1$.

2:
$\mathbf{t}$ – double *
Input/Output

On initial entry: the initial value of the independent variable, $t$.
On intermediate exit:
$t$, the current value of the independent variable.
On final exit: the value of the independent variable at which the computed solution
$y$ is returned (usually at
tout).

3:
$\mathbf{tout}$ – double
Input

On entry: the next value of $t$ at which a computed solution is desired.
On initial entry:
tout is used to determine the direction of integration. Integration is permitted in either direction (see also
itask).
Constraint:
${\mathbf{tout}}\ne {\mathbf{t}}$.

4:
$\mathbf{y}\left[{\mathbf{neq}}\right]$ – double
Input/Output

On initial entry: the vector of initial values of the dependent variables $y$.
On intermediate exit:
the computed solution vector, $y$, evaluated at $t$.
On final exit: the computed solution vector, evaluated at $t$ (usually $t={\mathbf{tout}}$).

5:
$\mathbf{ydot}\left[{\mathbf{neq}}\right]$ – double
Input/Output

On initial entry:
ydot must contain approximations to the time derivatives
${y}^{\prime}$ of the vector
$y$ evaluated at the initial value of the independent variable.
On exit: the time derivatives ${y}^{\prime}$ of the vector $y$ at the last integration point.

6:
$\mathbf{rtol}\left[\mathit{dim}\right]$ – double
Input/Output

Note: the dimension,
dim, of the array
rtol depends on the value of
vector_tol as set in
d02mwc; it
must be at least
 ${\mathbf{neq}}$ when ${\mathbf{vector\_tol}}=\mathrm{Nag\_TRUE}$;
 $1$ when ${\mathbf{vector\_tol}}=\mathrm{Nag\_FALSE}$.
On entry: the relative local error tolerance.
Constraint:
${\mathbf{rtol}}\left[\mathit{i}1\right]\ge 0.0$, for
$\mathit{i}=1,2,\dots ,n$where $n={\mathbf{neq}}$ when ${\mathbf{vector\_tol}}=\mathrm{Nag\_TRUE}$ and $n=1$ otherwise.
On exit:
rtol remains unchanged unless
d02nec exits with
${\mathbf{fail}}\mathbf{.}\mathbf{code}=$ NE_ODE_TOL in which case the values may have been increased to values estimated to be appropriate for continuing the integration.

7:
$\mathbf{atol}\left[\mathit{dim}\right]$ – double
Input/Output

Note: the dimension,
dim, of the array
atol depends on the value of
vector_tol as set in
d02mwc; it
must be at least
 ${\mathbf{neq}}$ when ${\mathbf{vector\_tol}}=\mathrm{Nag\_TRUE}$;
 $1$ when ${\mathbf{vector\_tol}}=\mathrm{Nag\_FALSE}$.
On entry: the absolute local error tolerance.
Constraint:
${\mathbf{atol}}\left[\mathit{i}1\right]\ge 0.0$, for
$\mathit{i}=1,2,\dots ,n$where $n={\mathbf{neq}}$ when ${\mathbf{vector\_tol}}=\mathrm{Nag\_TRUE}$ and $n=1$ otherwise.
On exit:
atol remains unchanged unless
d02nec exits with
${\mathbf{fail}}\mathbf{.}\mathbf{code}=$ NE_ODE_TOL in which case the values may have been increased to values estimated to be appropriate for continuing the integration.

8:
$\mathbf{itask}$ – Integer *
Input/Output

On initial entry: need not be set.
On exit: the task performed by the integrator on successful completion or an indicator that a problem occurred during integration.
 ${\mathbf{itask}}=2$
 The integration to tout was successfully completed (${\mathbf{t}}={\mathbf{tout}}$) by stepping exactly to tout.
 ${\mathbf{itask}}=3$
 The integration to tout was successfully completed (${\mathbf{t}}={\mathbf{tout}}$) by stepping past tout. y and ydot are obtained by interpolation.
 ${\mathbf{itask}}<0$
 Different negative values of itask returned correspond to different failure exits. fail should always be checked in such cases and the corrective action taken where appropriate.
itask must remain
unchanged between calls to
d02nec.

9:
$\mathbf{res}$ – function, supplied by the user
External Function

res must evaluate the residual
The specification of
res is:
void 
res (Integer neq,
double t,
const double y[],
const double ydot[],
double r[],
Integer *ires,
Nag_Comm *comm)



1:
$\mathbf{neq}$ – Integer
Input

On entry: the number of differentialalgebraic equations being solved.

2:
$\mathbf{t}$ – double
Input

On entry: $t$, the current value of the independent variable.

3:
$\mathbf{y}\left[{\mathbf{neq}}\right]$ – const double
Input

On entry: ${y}_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,{\mathbf{neq}}$, the current solution component.

4:
$\mathbf{ydot}\left[{\mathbf{neq}}\right]$ – const double
Input

On entry: the derivative of the solution at the current point $t$.

5:
$\mathbf{r}\left[{\mathbf{neq}}\right]$ – double
Output

On exit:
${\mathbf{r}}\left[\mathit{i}1\right]$ must contain the
$\mathit{i}$th component of
$R$, for
$\mathit{i}=1,2,\dots ,{\mathbf{neq}}$ where

6:
$\mathbf{ires}$ – Integer *
Input/Output

On entry: is always equal to zero.
On exit:
ires should normally be left unchanged. However, if an illegal value of
y is encountered,
ires should be set to
$1$;
d02nec will then attempt to resolve the problem so that illegal values of
y are not encountered.
ires should be set to
$2$ if you wish to return control to the calling function; this will cause
d02nec to exit with
${\mathbf{fail}}\mathbf{.}\mathbf{code}=$ NE_RES_FLAG.

7:
$\mathbf{comm}$ – Nag_Comm *
Pointer to structure of type Nag_Comm; the following members are relevant to
res.
 user – double *
 iuser – Integer *
 p – Pointer
The type Pointer will be
void *. Before calling
d02nec you may allocate memory and initialize these pointers with various quantities for use by
res when called from
d02nec (see
Section 3.1.1 in the Introduction to the NAG Library CL Interface).
Note: res should not return floatingpoint NaN (Not a Number) or infinity values, since these are not handled by
d02nec. If your code inadvertently
does return any NaNs or infinities,
d02nec is likely to produce unexpected results.

10:
$\mathbf{jac}$ – function, supplied by the user
External Function

Evaluates the matrix of partial derivatives,
$J$, where
If this option is not required, the actual argument for
jac must be specified as NULLFN. You must indicate to the integrator whether this option is to be used by setting the argument
jceval appropriately in a call to the setup function
d02mwc.
The specification of
jac is:
void 
jac (Integer neq,
double t,
const double y[],
const double ydot[],
double pd[],
double cj,
Nag_Comm *comm)



1:
$\mathbf{neq}$ – Integer
Input

On entry: the number of differentialalgebraic equations being solved.

2:
$\mathbf{t}$ – double
Input

On entry: $t$, the current value of the independent variable.

3:
$\mathbf{y}\left[{\mathbf{neq}}\right]$ – const double
Input

On entry: ${y}_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,{\mathbf{neq}}$, the current solution component.

4:
$\mathbf{ydot}\left[{\mathbf{neq}}\right]$ – const double
Input

On entry: the derivative of the solution at the current point $t$.

5:
$\mathbf{pd}\left[\mathit{dim}\right]$ – double
Input/Output

Note: the dimension of the array
pd will be
${\mathbf{neq}}\times {\mathbf{neq}}$ when the Jacobian is full and will be
$\left(2\times {\mathbf{ml}}+{\mathbf{mu}}+1\right)\times {\mathbf{neq}}$ when the Jacobian is banded (that is, a prior call to
d02npc has been made).
On entry:
pd is preset to zero before the call to
jac.
On exit: if the Jacobian is full then
${\mathbf{pd}}\left[\left(\mathit{j}1\right)\times {\mathbf{neq}}+\mathit{i}1\right]={J}_{\mathit{i}\mathit{j}}$, for
$\mathit{i}=1,2,\dots ,{\mathbf{neq}}$ and
$\mathit{j}=1,2,\dots ,{\mathbf{neq}}$; if the Jacobian is banded then
${\mathbf{pd}}\left[\left(j1\right)\times \left(2{\mathbf{ml}}+{\mathbf{mu}}+1\right)+{\mathbf{ml}}+{\mathbf{mu}}+ij\right]={J}_{ij}$, for
$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,j{\mathbf{mu}}\right)\le i\le \mathrm{min}\phantom{\rule{0.125em}{0ex}}\left(n,j+{\mathbf{ml}}\right)$; (see also in
f07bdc).

6:
$\mathbf{cj}$ – double
Input

On entry:
cj is a scalar constant which will be defined in
d02nec.

7:
$\mathbf{comm}$ – Nag_Comm *
Pointer to structure of type Nag_Comm; the following members are relevant to
jac.
 user – double *
 iuser – Integer *
 p – Pointer
The type Pointer will be
void *. Before calling
d02nec you may allocate memory and initialize these pointers with various quantities for use by
jac when called from
d02nec (see
Section 3.1.1 in the Introduction to the NAG Library CL Interface).
Note: jac should not return floatingpoint NaN (Not a Number) or infinity values, since these are not handled by
d02nec. If your code inadvertently
does return any NaNs or infinities,
d02nec is likely to produce unexpected results.

11:
$\mathbf{icom}\left[50+{\mathbf{neq}}\right]$ – Integer
Communication Array

icom contains information which is usually of no interest, but is necessary for subsequent calls. However you may find the following useful:
 ${\mathbf{icom}}\left[21\right]$
 The order of the method to be attempted on the next step.
 ${\mathbf{icom}}\left[22\right]$
 The order of the method used on the last step.
 ${\mathbf{icom}}\left[25\right]$
 The number of steps taken so far.
 ${\mathbf{icom}}\left[26\right]$
 The number of calls to res so far.
 ${\mathbf{icom}}\left[27\right]$
 The number of evaluations of the matrix of partial derivatives needed so far.
 ${\mathbf{icom}}\left[28\right]$
 The total number of error test failures so far.
 ${\mathbf{icom}}\left[29\right]$
 The total number of convergence test failures so far.

12:
$\mathbf{com}\left[{\mathbf{lcom}}\right]$ – double
Communication Array

com contains information which is usually of no interest, but is necessary for subsequent calls. However you may find the following useful:
 ${\mathbf{com}}\left[2\right]$
 The step size to be attempted on the next step.
 ${\mathbf{com}}\left[3\right]$
 The current value of the independent variable, i.e., the farthest point integration has reached. This will be different from t only when interpolation has been performed (${\mathbf{itask}}=3$).

13:
$\mathbf{lcom}$ – Integer
Input

On entry: the dimension of the array
com.
Constraint:
${\mathbf{lcom}}\ge 40+\left(\mathit{maxorder}+4\right)\times {\mathbf{neq}}+{\mathbf{neq}}\times p+q$ where
$\mathit{maxorder}$ is the maximum order that can be used by the integration method (see
maxord in
d02mwc);
$p={\mathbf{neq}}$ when the Jacobian is full and
$p=\left(2\times {\mathbf{ml}}+{\mathbf{mu}}+1\right)$ when the Jacobian is banded; and,
$q=\left({\mathbf{neq}}/\left({\mathbf{ml}}+{\mathbf{mu}}+1\right)\right)+1$ when the Jacobian is to be evaluated numerically and
$q=0$ otherwise.

14:
$\mathbf{comm}$ – Nag_Comm *

The NAG communication argument (see
Section 3.1.1 in the Introduction to the NAG Library CL Interface).

15:
$\mathbf{fail}$ – NagError *
Input/Output

The NAG error argument (see
Section 7 in the Introduction to the NAG Library CL Interface).
6
Error Indicators and Warnings
 NE_ALLOC_FAIL

Dynamic memory allocation failed.
See
Section 3.1.2 in the Introduction to the NAG Library CL Interface for further information.
 NE_ARRAY_INPUT

All elements of
rtol and
atol are zero.
Some element of
atol is less than zero.
Some element of
rtol is less than zero.
 NE_BAD_PARAM

On entry, argument $\u2329\mathit{\text{value}}\u232a$ had an illegal value.
 NE_CONV_CONT

The corrector could not converge and the error test failed repeatedly. ${\mathbf{t}}=\u2329\mathit{\text{value}}\u232a$. Stepsize $h=\u2329\mathit{\text{value}}\u232a$.
The corrector repeatedly failed to converge with $\lefth\right=\mathit{hmin}$. ${\mathbf{t}}=\u2329\mathit{\text{value}}\u232a$. Stepsize $h=\u2329\mathit{\text{value}}\u232a$.
 NE_CONV_JACOBG

The iteration matrix is singular. ${\mathbf{t}}=\u2329\mathit{\text{value}}\u232a$. Stepsize $h=\u2329\mathit{\text{value}}\u232a$.
 NE_CONV_ROUNDOFF

The error test failed repeatedly with $\lefth\right=\mathit{hmin}$. ${\mathbf{t}}=\u2329\mathit{\text{value}}\u232a$. Stepsize $h=\u2329\mathit{\text{value}}\u232a$.
 NE_INITIALIZATION

Either the initialization function has not been called prior to the first call of this function or a communication array has become corrupted.
 NE_INT

A previous call to this function returned with ${\mathbf{itask}}=\u2329\mathit{\text{value}}\u232a$ and no appropriate action was taken.
 NE_INT_2

com array is of insufficient length; length required
$\text{}=\u2329\mathit{\text{value}}\u232a$; actual length
${\mathbf{lcom}}=\u2329\mathit{\text{value}}\u232a$.
 NE_INT_ARG_LT

On entry, ${\mathbf{neq}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{neq}}\ge 1$.
 NE_INTERNAL_ERROR

An internal error has occurred in this function. Check the function call and any array sizes. If the call is correct then please contact
NAG for assistance.
See
Section 7.5 in the Introduction to the NAG Library CL Interface for further information.
 NE_MAX_STEP

Maximum number of steps taken on this call before reaching
tout:
${\mathbf{t}}=\u2329\mathit{\text{value}}\u232a$, maximum number of steps
$\text{}=\u2329\mathit{\text{value}}\u232a$.
 NE_NO_LICENCE

Your licence key may have expired or may not have been installed correctly.
See
Section 8 in the Introduction to the NAG Library CL Interface for further information.
 NE_ODE_TOL

A solution component has become zero when a purely relative tolerance (zero absolute tolerance) was selected for that component. ${\mathbf{t}}=\u2329\mathit{\text{value}}\u232a$, ${\mathbf{y}}\left[\mathit{I}1\right]=\u2329\mathit{\text{value}}\u232a$ for component $\mathit{I}=\u2329\mathit{\text{value}}\u232a$.
Too much accuracy requested for precision of machine.
rtol and
atol were increased by scale factor
$R$. Try running again with these scaled tolerances.
${\mathbf{t}}=\u2329\mathit{\text{value}}\u232a$,
$R=\u2329\mathit{\text{value}}\u232a$.
 NE_REAL_2

tout is behind
t in the direction of
$h$:
${\mathbf{tout}}{\mathbf{t}}=\u2329\mathit{\text{value}}\u232a$,
$h=\u2329\mathit{\text{value}}\u232a$.
tout is too close to
t to start integration:
${\mathbf{tout}}{\mathbf{t}}=\u2329\mathit{\text{value}}\u232a$:
$\mathit{hmin}=\u2329\mathit{\text{value}}\u232a$.
 NE_REAL_ARG_EQ

On entry, ${\mathbf{t}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{tout}}\ne {\mathbf{t}}$.
 NE_RES_FLAG

ires was set to
$1$ during a call to
res and could not be resolved.
${\mathbf{t}}=\u2329\mathit{\text{value}}\u232a$. Stepsize
$h=\u2329\mathit{\text{value}}\u232a$.
ires was set to
$2$ during a call to
res.
${\mathbf{t}}=\u2329\mathit{\text{value}}\u232a$. Stepsize
$\text{}=\u2329\mathit{\text{value}}\u232a$.
Repeated occurrences of input constraint violations have been detected. This could result in a potential infinite loop: ${\mathbf{itask}}=\u2329\mathit{\text{value}}\u232a$. Current violation corresponds to exit with ${\mathbf{fail}}\mathbf{.}\mathbf{code}=\u2329\mathit{\text{value}}\u232a$.
 NE_SINGULAR_POINT

The initial
ydot could not be computed.
${\mathbf{t}}=\u2329\mathit{\text{value}}\u232a$. Stepsize
$h=\u2329\mathit{\text{value}}\u232a$.
7
Accuracy
The accuracy of the numerical solution may be controlled by a careful choice of the arguments
rtol and
atol. 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{vector\_tol}}=\mathrm{Nag\_FALSE}$ with
${\mathbf{atol}}\left[0\right]$ small but positive).
8
Parallelism and Performance
d02nec is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.
d02nec 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 function. Please also consult the
Users' Note for your implementation for any additional implementationspecific 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 banded systems the cost is proportional to ${\mathbf{neq}}\times {\left({\mathbf{ml}}+{\mathbf{mu}}+1\right)}^{2}$, while for full systems the cost is proportional to
${{\mathbf{neq}}}^{3}$. Note however that for moderately sized problems which are only mildly nonlinear the cost may be dominated by factors proportional to ${\mathbf{neq}}\times \left({\mathbf{ml}}+{\mathbf{mu}}+1\right)$ and ${{\mathbf{neq}}}^{2}$ respectively.
10
Example
This example solves the wellknown stiff Robertson problem written in implicit form
with initial conditions
$a=1.0$ and
$b=c=0.0$ over the range
$\left[0,0.1\right]$ the BDF method (setup function
d02mwc and
d02npc).
10.1
Program Text
10.2
Program Data
None.
10.3
Program Results