NAG Library Function Document

nag_dae_ivp_dassl_gen (d02nec)

1  Purpose

2  Specification

3  Description

/* declarations */
   EXTERN res, jac
        .
        .
        .
/* Initialize the integrator */
   nag_dae_ivp_dassl_setup(...);
/* Is the Jacobian matrix banded? */
   if (banded) {nag_dae_ivp_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_dae_ivp_dassl_gen (...);
     if (itask != 1)
        tout = MIN(tout+dt,tend);
        /* Print the solution */
   }
       .
       .
       .

4  References

5  Arguments

1:     neq – IntegerInput

On entry: the number of differential-algebraic equations to be solved.

Constraint: $\mathbf{neq}\ge 1$.

2:     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:     tout – doubleInput

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:     y[neq] – doubleInput/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:     ydot[neq] – doubleInput/Output

On initial entry: ydot must contain approximations to the time derivatives $y\text{'}$ of the vector $y$ evaluated at the initial value of the independent variable.

On exit: the time derivatives $y\text{'}$ of the vector $y$ at the last integration point.

6:     rtol[$\mathit{dim}$] – doubleInput/Output

Note: the dimension, dim, of the array rtol depends on the value of vector_tol as set in nag_dae_ivp_dassl_setup (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 nag_dae_ivp_dassl_gen (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:     atol[$\mathit{dim}$] – doubleInput/Output

Note: the dimension, dim, of the array atol depends on the value of vector_tol as set in nag_dae_ivp_dassl_setup (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 nag_dae_ivp_dassl_gen (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:     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 nag_dae_ivp_dassl_gen (d02nec).

9:     res – function, supplied by the userExternal Function

res must evaluate the residual

$$R=F\left(t,y,y^{\prime }\right)\text{.}$$

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:     neq – IntegerInput

On entry: the number of differential-algebraic equations being solved.

2:     t – doubleInput

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

3:     y[neq] – const doubleInput

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

4:     ydot[neq] – const doubleInput

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

5:     r[neq] – doubleOutput

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

$$R=F\left(\mathbf{t},\mathbf{y},\mathbf{ydot}\right)\text{.}$$

6:     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$; nag_dae_ivp_dassl_gen (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 nag_dae_ivp_dassl_gen (d02nec) to exit with $\mathbf{fail}\mathbf{.}\mathbf{code}=$ NE_RES_FLAG.

7:     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 nag_dae_ivp_dassl_gen (d02nec) you may allocate memory and initialize these pointers with various quantities for use by res when called from nag_dae_ivp_dassl_gen (d02nec) (see Section 3.2.1.1 in the Essential Introduction).

10:   jac – function, supplied by the userExternal Function

Evaluates the matrix of partial derivatives, $J$, where

$$J_{ij}=\frac{\partial F_i}{\partial y_j}+\mathbf{cj}\times \frac{\partial F_i}{\partial {y^{\prime }}_j}\text{,\hspace{1em}}i,j=1,2,\dots ,\mathbf{neq}\text{.}$$

If this option is not required, the actual argument for jac may 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 nag_dae_ivp_dassl_setup (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:     neq – IntegerInput

On entry: the number of differential-algebraic equations being solved.

2:     t – doubleInput

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

3:     y[neq] – const doubleInput

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

4:     ydot[neq] – const doubleInput

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

5:     pd[$\mathit{dim}$] – doubleInput/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 nag_dae_ivp_dassl_linalg (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(j-1\right)\times \left(2\mathbf{ml}+\mathbf{mu}+1\right)+\mathbf{ml}+\mathbf{mu}+i-j\right]=J_{ij}$, for $\max \left(1,j-\mathbf{mu}\right)\le i\le \min \left(n,j+\mathbf{ml}\right)$; (see also in nag_dgbtrf (f07bdc)).

6:     cj – doubleInput

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

7:     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 nag_dae_ivp_dassl_gen (d02nec) you may allocate memory and initialize these pointers with various quantities for use by jac when called from nag_dae_ivp_dassl_gen (d02nec) (see Section 3.2.1.1 in the Essential Introduction).

11:   icom[$50+\mathbf{neq}$] – IntegerCommunication 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:   com[lcom] – doubleCommunication 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:   lcom – IntegerInput

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 nag_dae_ivp_dassl_setup (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:   comm – Nag_Comm *Communication Structure

The NAG communication argument (see Section 3.2.1.1 in the Essential Introduction).

15:   fail – NagError *Input/Output

The NAG error argument (see Section 3.6 in the Essential Introduction).

6  Error Indicators and Warnings

NE_ALLOC_FAIL

Dynamic memory allocation failed.

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 $⟨\mathit{\text{value}}⟩$ had an illegal value.

NE_CONV_CONT

The corrector could not converge and the error test failed repeatedly. $\mathbf{t}=⟨\mathit{\text{value}}⟩$. Stepsize $h=⟨\mathit{\text{value}}⟩$.

The corrector repeatedly failed to converge with $\left|h\right|=\mathit{hmin}$. $\mathbf{t}=⟨\mathit{\text{value}}⟩$. Stepsize $h=⟨\mathit{\text{value}}⟩$.

NE_CONV_JACOBG

The iteration matrix is singular. $\mathbf{t}=⟨\mathit{\text{value}}⟩$. Stepsize $h=⟨\mathit{\text{value}}⟩$.

NE_CONV_ROUNDOFF

The error test failed repeatedly with $\left|h\right|=\mathit{hmin}$. $\mathbf{t}=⟨\mathit{\text{value}}⟩$. Stepsize $h=⟨\mathit{\text{value}}⟩$.

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}=⟨\mathit{\text{value}}⟩$ and no appropriate action was taken.

NE_INT_2

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

NE_INT_ARG_LT

On entry, $\mathbf{neq}=⟨\mathit{\text{value}}⟩$.
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.

NE_MAX_STEP

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

NE_ODE_TOL

A solution component has become zero when a purely relative tolerance (zero absolute tolerance) was selected for that component. $\mathbf{t}=⟨\mathit{\text{value}}⟩$, $\mathbf{y}\left[\mathit{I}-1\right]=⟨\mathit{\text{value}}⟩$ for component $\mathit{I}=⟨\mathit{\text{value}}⟩$.

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}=⟨\mathit{\text{value}}⟩$, $R=⟨\mathit{\text{value}}⟩$.

NE_REAL_2

tout is behind t in the direction of $h$: $\mathbf{tout}-\mathbf{t}=⟨\mathit{\text{value}}⟩$, $h=⟨\mathit{\text{value}}⟩$.

tout is too close to t to start integration: $\mathbf{tout}-\mathbf{t}=⟨\mathit{\text{value}}⟩$: $\mathit{hmin}=⟨\mathit{\text{value}}⟩$.

NE_REAL_ARG_EQ

On entry, $\mathbf{t}=⟨\mathit{\text{value}}⟩$.
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}=⟨\mathit{\text{value}}⟩$. Stepsize $h=⟨\mathit{\text{value}}⟩$.

ires was set to $-2$ during a call to res. $\mathbf{t}=⟨\mathit{\text{value}}⟩$. Stepsize $=⟨\mathit{\text{value}}⟩$.

Repeated occurrences of input constraint violations have been detected. This could result in a potential infinite loop: $\mathbf{itask}=⟨\mathit{\text{value}}⟩$. Current violation corresponds to exit with $\mathbf{fail}\mathbf{.}\mathbf{code}=⟨\mathit{\text{value}}⟩$.

NE_SINGULAR_POINT

The initial ydot could not be computed. $\mathbf{t}=⟨\mathit{\text{value}}⟩$. Stepsize $h=⟨\mathit{\text{value}}⟩$.

7  Accuracy

8  Parallelism and Performance

9  Further Comments

10  Example

10.1  Program Text

Program Text (d02nece.c)

10.2  Program Data

10.3  Program Results

Program Results (d02nece.r)