NAG FL Interface
d02nef (dae_​dassl_​gen)

1 Purpose

2 Specification

3 Description

!     Declarations

      External res, jac
             .
             .
             .
!     Initialize the integrator
      Call d02mwf(...)
!     Is the Jacobian matrix banded?
      If (banded) Call d02npf(...)

!     Set dt to the required temporal resolution
!     Set tend to the final time
!     Call the integrator for each temporal value:
1000  Call d02nef(...,res,jac,...)
!     Continue integration?
      If (tout.lt.tend .and. itask.ge.0) Then
        If (itask.ne.1) tout = min(tout+dt,tend)
!       Print solution
        Call d02mcf(...)
        Go To 1000
      Endif
             .
             .
             .

4 References

5 Arguments

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

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

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

2: $\mathbf{t}$ – Real (Kind=nag_wp) 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}$ – Real (Kind=nag_wp) 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)$ – Real (Kind=nag_wp) array 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=T$.

On final exit: the computed solution vector, evaluated at $t$ (usually $t=\mathbf{tout}$).

5: $\mathbf{ydot}\left(\mathbf{neq}\right)$ – Real (Kind=nag_wp) array 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(*\right)$ – Real (Kind=nag_wp) array Input/Output

Note: the dimension of the array rtol depends on the value of itol as set in d02mwf; it  must be at least $\mathbf{neq}$ if $\mathbf{itol}=\mathrm{.TRUE.}$ was set in d02mwf and at least $1$ if $\mathbf{itol}=\mathrm{.FALSE.}$ was set in d02mwf.

On entry: the relative local error tolerance.

Constraint: $\mathbf{rtol}\left(\mathit{i}\right)\ge 0.0$, for $\mathit{i}=1,2,\dots ,n$

where $n=\mathbf{neq}$ when $\mathbf{itol}=\mathrm{.TRUE.}$ and $n=1$ otherwise.

On exit: rtol remains unchanged unless d02nef exits with $\mathbf{ifail}=\mathbf{16}$ in which case the values may have been increased to values estimated to be appropriate for continuing the integration.

7: $\mathbf{atol}\left(*\right)$ – Real (Kind=nag_wp) array Input/Output

Note: the dimension of the array atol depends on the value of itol as set in d02mwf; it  must be at least $\mathbf{neq}$ if $\mathbf{itol}=\mathrm{.TRUE.}$ was set in d02mwf and at least $1$ if $\mathbf{itol}=\mathrm{.FALSE.}$ was set in d02mwf.

On entry: the absolute local error tolerance.

Constraint: $\mathbf{atol}\left(\mathit{i}\right)\ge 0.0$, for $\mathit{i}=1,2,\dots ,n$

where $n=\mathbf{neq}$ when $\mathbf{itol}=\mathrm{.TRUE.}$ and $n=1$ otherwise.

On exit: atol remains unchanged unless d02nef exits with $\mathbf{ifail}=\mathbf{16}$ 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. ifail should always be checked in such cases and the corrective action taken where appropriate.

itask must remain unchanged between calls to d02nef.

9: $\mathbf{res}$ – Subroutine, supplied by the user. External Procedure

res must evaluate the residual

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

The specification of res is:

Fortran Interface copy

Subroutine res ( neq, t, y, ydot, r, ires, iuser, ruser) Integer, Intent (In) :: neq Integer, Intent (Inout) :: ires, iuser(*) Real (Kind=nag_wp), Intent (In) :: t, y(neq), ydot(neq) Real (Kind=nag_wp), Intent (Inout) :: ruser(*) Real (Kind=nag_wp), Intent (Out) :: r(neq)

C Header Interface copy

void  res (const Integer *neq, const double *t, const double y[], const double ydot[], double r[], Integer *ires, Integer iuser[], double ruser[])

C++ Header Interface copy

#include <nag.h> extern "C" { void  res (const Integer &neq, const double &t, const double y[], const double ydot[], double r[], Integer &ires, Integer iuser[], double ruser[]) }

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

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

2: $\mathbf{t}$ – Real (Kind=nag_wp) Input

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

3: $\mathbf{y}\left(\mathbf{neq}\right)$ – Real (Kind=nag_wp) array 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)$ – Real (Kind=nag_wp) array Input

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

5: $\mathbf{r}\left(\mathbf{neq}\right)$ – Real (Kind=nag_wp) array Output

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

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

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 $\mathrm{-1}$; d02nef will then attempt to resolve the problem so that illegal values of y are not encountered. ires should be set to $\mathrm{-2}$ if you wish to return control to the calling (sub)routine; this will cause d02nef to exit with $\mathbf{ifail}=\mathbf{23}$.

7: $\mathbf{iuser}\left(*\right)$ – Integer array User Workspace

8: $\mathbf{ruser}\left(*\right)$ – Real (Kind=nag_wp) array User Workspace

res is called with the arguments iuser and ruser as supplied to d02nef. You should use the arrays iuser and ruser to supply information to res.

res must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which d02nef is called. Arguments denoted as Input must not be changed by this procedure.

Note: res should not return floating-point NaN (Not a Number) or infinity values, since these are not handled by d02nef. If your code inadvertently does return any NaNs or infinities, d02nef is likely to produce unexpected results.

10: $\mathbf{jac}$ – Subroutine, supplied by the NAG Library or the user. External Procedure

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 must be the dummy routine d02nez. (d02nez is included in the NAG Library.) 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 routine d02mwf.

The specification of jac is:

Fortran Interface copy

Subroutine jac ( neq, t, y, ydot, pd, cj, iuser, ruser) Integer, Intent (In) :: neq Integer, Intent (Inout) :: iuser(*) Real (Kind=nag_wp), Intent (In) :: t, y(neq), ydot(neq), cj Real (Kind=nag_wp), Intent (Inout) :: pd(*), ruser(*)

C Header Interface copy

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

C++ Header Interface copy

#include <nag.h> extern "C" { void  jac (const Integer &neq, const double &t, const double y[], const double ydot[], double pd[], const double &cj, Integer iuser[], double ruser[]) }

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

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

2: $\mathbf{t}$ – Real (Kind=nag_wp) Input

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

3: $\mathbf{y}\left(\mathbf{neq}\right)$ – Real (Kind=nag_wp) array 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)$ – Real (Kind=nag_wp) array Input

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

5: $\mathbf{pd}\left(*\right)$ – Real (Kind=nag_wp) array Input/Output

On entry: pd is preset to zero before the call to jac.

On exit: if the Jacobian is full then $\mathbf{pd}\left((\mathit{j}-1)\times \mathbf{neq}+\mathit{i}\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((j-1)\times (2\mathbf{ml}+\mathbf{mu}+1)+\mathbf{ml}+\mathbf{mu}+i-j+1\right)=J_{ij}$, for $\max (1,j-\mathbf{mu})\le i\le \min (n,j+\mathbf{ml})$.

6: $\mathbf{cj}$ – Real (Kind=nag_wp) Input

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

7: $\mathbf{iuser}\left(*\right)$ – Integer array User Workspace

8: $\mathbf{ruser}\left(*\right)$ – Real (Kind=nag_wp) array User Workspace

jac is called with the arguments iuser and ruser as supplied to d02nef. You should use the arrays iuser and ruser to supply information to jac.

jac must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which d02nef is called. Arguments denoted as Input must not be changed by this procedure.

Note: jac should not return floating-point NaN (Not a Number) or infinity values, since these are not handled by d02nef. If your code inadvertently does return any NaNs or infinities, d02nef is likely to produce unexpected results.

11: $\mathbf{icom}\left(50+\mathbf{neq}\right)$ – Integer array 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(22\right)$

The order of the method to be attempted on the next step.

$\mathbf{icom}\left(23\right)$

The order of the method used on the last step.

$\mathbf{icom}\left(26\right)$

The number of steps taken so far.

$\mathbf{icom}\left(27\right)$

The number of calls to res so far.

$\mathbf{icom}\left(28\right)$

The number of evaluations of the matrix of partial derivatives needed so far.

$\mathbf{icom}\left(29\right)$

The total number of error test failures so far.

$\mathbf{icom}\left(30\right)$

The total number of convergence test failures so far.

12: $\mathbf{com}\left(\mathbf{lcom}\right)$ – Real (Kind=nag_wp) array 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(3\right)$

The step size to be attempted on the next step.

$\mathbf{com}\left(4\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 as declared in the (sub)program from which d02nef is called.

Constraint: $\mathbf{lcom}\ge 40+(\mathit{maxorder}+4)\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 d02mwf); $p=\mathbf{neq}$ when the Jacobian is full and $p=(2\times \mathbf{ml}+\mathbf{mu}+1)$ when the Jacobian is banded; and, $q=(\mathbf{neq}/(\mathbf{ml}+\mathbf{mu}+1))+1$ when the Jacobian is to be evaluated numerically and $q=0$ otherwise.

14: $\mathbf{iuser}\left(*\right)$ – Integer array User Workspace

15: $\mathbf{ruser}\left(*\right)$ – Real (Kind=nag_wp) array User Workspace

iuser and ruser are not used by d02nef, but are passed directly to res and jac and may be used to pass information to these routines.

16: $\mathbf{ifail}$ – Integer Input/Output

On entry: ifail must be set to $0$, $\mathrm{-1}$ or $1$ to set behaviour on detection of an error; these values have no effect when no error is detected.

A value of $0$ causes the printing of an error message and program execution will be halted; otherwise program execution continues. A value of $\mathrm{-1}$ means that an error message is printed while a value of $1$ means that it is not.

If halting is not appropriate, the value $\mathrm{-1}$ or $1$ is recommended. If message printing is undesirable, then the value $1$ is recommended. Otherwise, the value $\mathrm{-1}$ is recommended since useful values can be provided in some output arguments even when $\mathbf{ifail}\ne \mathbf{0}$ on exit. When the value $-\mathbf{1}$ or $\mathbf{1}$ is used it is essential to test the value of ifail on exit.

On exit: $\mathbf{ifail}=\mathbf{0}$ unless the routine detects an error or a warning has been flagged (see Section 6).

6 Error Indicators and Warnings

$\mathbf{ifail}=1$

On entry, $\mathbf{neq}=⟨\mathit{\text{value}}⟩$.
Constraint: $\mathbf{neq}\ge 1$.

$\mathbf{ifail}=3$

On entry, $\mathbf{t}=⟨\mathit{\text{value}}⟩$.
Constraint: $\mathbf{tout}\ne \mathbf{t}$.

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}}⟩$.

$\mathbf{ifail}=6$

Some element of rtol is less than zero.

$\mathbf{ifail}=7$

Some element of atol is less than zero.

$\mathbf{ifail}=8$

A previous call to this routine returned with $\mathbf{itask}=⟨\mathit{\text{value}}⟩$ and no appropriate action was taken.

$\mathbf{ifail}=11$

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

$\mathbf{ifail}=12$

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

$\mathbf{ifail}=13$

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

$\mathbf{ifail}=14$

All elements of rtol and atol are zero.

$\mathbf{ifail}=15$

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

$\mathbf{ifail}=16$

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}}⟩$.

$\mathbf{ifail}=17$

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}\right)=⟨\mathit{\text{value}}⟩$ for component $\mathit{I}=⟨\mathit{\text{value}}⟩$.

$\mathbf{ifail}=18$

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

$\mathbf{ifail}=19$

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

$\mathbf{ifail}=20$

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

$\mathbf{ifail}=21$

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

$\mathbf{ifail}=22$

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}}⟩$.

$\mathbf{ifail}=23$

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

$\mathbf{ifail}=24$

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

$\mathbf{ifail}=25$

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

$\mathbf{ifail}=-99$

An unexpected error has been triggered by this routine. Please contact NAG.

See Section 7 in the Introduction to the NAG Library FL Interface for further information.

$\mathbf{ifail}=-399$

Your licence key may have expired or may not have been installed correctly.

See Section 8 in the Introduction to the NAG Library FL Interface for further information.

$\mathbf{ifail}=-999$

Dynamic memory allocation failed.

See Section 9 in the Introduction to the NAG Library FL Interface for further information.

7 Accuracy

8 Parallelism and Performance

9 Further Comments

10 Example

10.1 Program Text

10.2 Program Data

10.3 Program Results