Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

Chapter Contents
Chapter Introduction
NAG Toolbox

# NAG Toolbox: nag_ode_dae_dassl_gen (d02ne)

## Purpose

nag_ode_dae_dassl_gen (d02ne) is a function for integrating stiff systems of implicit ordinary differential equations coupled with algebraic equations.

## Syntax

[t, y, ydot, rtol, atol, itask, icom, com, user, ifail] = d02ne(t, tout, y, ydot, rtol, atol, itask, res, jac, icom, com, 'neq', neq, 'user', user)
[t, y, ydot, rtol, atol, itask, icom, com, user, ifail] = nag_ode_dae_dassl_gen(t, tout, y, ydot, rtol, atol, itask, res, jac, icom, com, 'neq', neq, 'user', user)

## Description

nag_ode_dae_dassl_gen (d02ne) 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
 $F t,y,y′ = 0 .$
nag_ode_dae_dassl_gen (d02ne) 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 nag_ode_dae_dassl_gen (d02ne) is given below. It calls the DASSL implementation of the BDF integrator setup function nag_ode_dae_dassl_setup (d02mw) and the banded matrix setup function nag_ode_dae_dassl_linalg (d02np) (if required), and, if the integration needs to proceed, calls nag_ode_dae_dassl_cont (d02mc) before continuing the integration.
```.
.
.
% Initialize the integrator
[...] = d02mw(...);
% Is the Jacobian matrix banded?
if (banded)
[...] = d02np(...);
end
% Set DT to the required temporal resolution
% Set TEND to the final time
% Call the integrator for each temporal value
while 1
[tout,...] = d02ne(...);
break
end
tout = min(tout + dt, tend);
end
% Print the solution
[...] = d02mc(...);
end
.
.
.```

None.

## Parameters

### Compulsory Input Parameters

1:     $\mathrm{t}$ – double scalar
On initial entry: the initial value of the independent variable, $t$.
2:     $\mathrm{tout}$ – double scalar
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}}$.
3:     $\mathrm{y}\left({\mathbf{neq}}\right)$ – double array
On initial entry: the vector of initial values of the dependent variables $y$.
4:     $\mathrm{ydot}\left({\mathbf{neq}}\right)$ – double array
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.
5:     $\mathrm{rtol}\left(:\right)$ – double array
The dimension of the array rtol depends on the value of itol as set in nag_ode_dae_dassl_setup (d02mw); it  must be at least ${\mathbf{neq}}$ if ${\mathbf{itol}}=\mathit{true}$ and at least $1$ if ${\mathbf{itol}}=\mathit{false}$
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}}=\mathit{true}$ and $n=1$ otherwise.
6:     $\mathrm{atol}\left(:\right)$ – double array
The dimension of the array atol depends on the value of itol as set in nag_ode_dae_dassl_setup (d02mw); it  must be at least ${\mathbf{neq}}$ if ${\mathbf{itol}}=\mathit{true}$ and at least $1$ if ${\mathbf{itol}}=\mathit{false}$
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}}=\mathit{true}$ and $n=1$ otherwise.
7:     $\mathrm{itask}$int64int32nag_int scalar
On initial entry: need not be set.
8:     $\mathrm{res}$ – function handle or string containing name of m-file
res must evaluate the residual
 $R = F t,y,y′ .$
[r, ires, user] = res(neq, t, y, ydot, ires, user)

Input Parameters

1:     $\mathrm{neq}$int64int32nag_int scalar
The number of differential-algebraic equations being solved.
2:     $\mathrm{t}$ – double scalar
$t$, the current value of the independent variable.
3:     $\mathrm{y}\left({\mathbf{neq}}\right)$ – double array
${y}_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,{\mathbf{neq}}$, the current solution component.
4:     $\mathrm{ydot}\left({\mathbf{neq}}\right)$ – double array
The derivative of the solution at the current point $t$.
5:     $\mathrm{ires}$int64int32nag_int scalar
Is always equal to zero.
6:     $\mathrm{user}$ – Any MATLAB object
res is called from nag_ode_dae_dassl_gen (d02ne) with the object supplied to nag_ode_dae_dassl_gen (d02ne).

Output Parameters

1:     $\mathrm{r}\left({\mathbf{neq}}\right)$ – double array
${\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 t,y,ydot .$
2:     $\mathrm{ires}$int64int32nag_int scalar
ires should normally be left unchanged. However, if an illegal value of y is encountered, ires should be set to $-1$; nag_ode_dae_dassl_gen (d02ne) 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 (sub)routine; this will cause nag_ode_dae_dassl_gen (d02ne) to exit with ${\mathbf{ifail}}={\mathbf{23}}$.
3:     $\mathrm{user}$ – Any MATLAB object
9:     $\mathrm{jac}$ – function handle or string containing name of m-file
Evaluates the matrix of partial derivatives, $J$, where
 $J ij = ∂Fi ∂yj + cj × ∂Fi ∂y′j , i,j=1,2,…,neq .$
If this option is not required, the actual argument for jac must be the string nag_ode_dae_dassl_gen_dummy_jac (d02nez). (nag_ode_dae_dassl_gen_dummy_jac (d02nez) is included in the NAG Toolbox.) 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_ode_dae_dassl_setup (d02mw).
[pd, user] = jac(neq, t, y, ydot, pd, cj, user)

Input Parameters

1:     $\mathrm{neq}$int64int32nag_int scalar
The number of differential-algebraic equations being solved.
2:     $\mathrm{t}$ – double scalar
$t$, the current value of the independent variable.
3:     $\mathrm{y}\left({\mathbf{neq}}\right)$ – double array
${y}_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,{\mathbf{neq}}$, the current solution component.
4:     $\mathrm{ydot}\left({\mathbf{neq}}\right)$ – double array
The derivative of the solution at the current point $t$.
5:     $\mathrm{pd}\left(:\right)$ – double array
pd is preset to zero before the call to jac.
6:     $\mathrm{cj}$ – double scalar
cj is a scalar constant which will be defined in nag_ode_dae_dassl_gen (d02ne).
7:     $\mathrm{user}$ – Any MATLAB object
jac is called from nag_ode_dae_dassl_gen (d02ne) with the object supplied to nag_ode_dae_dassl_gen (d02ne).

Output Parameters

1:     $\mathrm{pd}\left(:\right)$ – double array
If the Jacobian is full then ${\mathbf{pd}}\left(\left(\mathit{j}-1\right)×{\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(\left(j-1\right)×\left(2{\mathbf{ml}}+{\mathbf{mu}}+1\right)+{\mathbf{ml}}+{\mathbf{mu}}+i-j+1\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 nag_lapack_dgbtrf (f07bd)).
2:     $\mathrm{user}$ – Any MATLAB object
10:   $\mathrm{icom}\left(50+{\mathbf{neq}}\right)$int64int32nag_int 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.
11:   $\mathrm{com}\left(\mathit{lcom}\right)$ – double array
lcom, the dimension of the array, must satisfy the constraint $\mathit{lcom}\ge 40+\left(\mathit{maxorder}+4\right)×{\mathbf{neq}}+{\mathbf{neq}}×p+q$ where $\mathit{maxorder}$ is the maximum order that can be used by the integration method (see maxord in nag_ode_dae_dassl_setup (d02mw)); $p={\mathbf{neq}}$ when the Jacobian is full and $p=\left(2×{\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.
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$).

### Optional Input Parameters

1:     $\mathrm{neq}$int64int32nag_int scalar
Default: the dimension of the arrays y, ydot. (An error is raised if these dimensions are not equal.)
The number of differential-algebraic equations to be solved.
Constraint: ${\mathbf{neq}}\ge 1$.
2:     $\mathrm{user}$ – Any MATLAB object
user is not used by nag_ode_dae_dassl_gen (d02ne), but is passed to res and jac. Note that for large objects it may be more efficient to use a global variable which is accessible from the m-files than to use user.

### Output Parameters

1:     $\mathrm{t}$ – double scalar
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).
2:     $\mathrm{y}\left({\mathbf{neq}}\right)$ – double array
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}}$).
3:     $\mathrm{ydot}\left({\mathbf{neq}}\right)$ – double array
The time derivatives $y\text{'}$ of the vector $y$ at the last integration point.
4:     $\mathrm{rtol}\left(:\right)$ – double array
The dimension of the array rtol depends on the value of itol as set in nag_ode_dae_dassl_setup (d02mw); it  will be ${\mathbf{neq}}$ if ${\mathbf{itol}}=\mathit{true}$ and at least $1$ if ${\mathbf{itol}}=\mathit{false}$
rtol remains unchanged unless nag_ode_dae_dassl_gen (d02ne) 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.
5:     $\mathrm{atol}\left(:\right)$ – double array
The dimension of the array atol depends on the value of itol as set in nag_ode_dae_dassl_setup (d02mw); it  will be ${\mathbf{neq}}$ if ${\mathbf{itol}}=\mathit{true}$ and at least $1$ if ${\mathbf{itol}}=\mathit{false}$
atol remains unchanged unless nag_ode_dae_dassl_gen (d02ne) 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.
6:     $\mathrm{itask}$int64int32nag_int scalar
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 nag_ode_dae_dassl_gen (d02ne).
7:     $\mathrm{icom}\left(50+{\mathbf{neq}}\right)$int64int32nag_int array
8:     $\mathrm{com}\left(\mathit{lcom}\right)$ – double array
9:     $\mathrm{user}$ – Any MATLAB object
10:   $\mathrm{ifail}$int64int32nag_int scalar
${\mathbf{ifail}}={\mathbf{0}}$ unless the function detects an error (see Error Indicators and Warnings).

## Error Indicators and Warnings

Note: nag_ode_dae_dassl_gen (d02ne) may return useful information for one or more of the following detected errors or warnings.
Errors or warnings detected by the function:
${\mathbf{ifail}}=1$
Constraint: ${\mathbf{neq}}\ge 1$.
${\mathbf{ifail}}=3$
Constraint: ${\mathbf{tout}}\ne {\mathbf{t}}$.
tout is behind t in the direction of $h$.
tout is too close to t to start integration.
${\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 function returned with ${\mathbf{itask}}=_$ and no appropriate action was taken.
${\mathbf{ifail}}=11$
Either the initialization function has not been called prior to the first call of this function or a communication array has become corrupted.
${\mathbf{ifail}}=12$
Either the initialization function has not been called prior to the first call of this function or a communication array has become corrupted.
${\mathbf{ifail}}=13$
com array is of insufficient length; length required $\text{}=_$; actual length $\mathit{lcom}=_$.
${\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{ifail}}=16$
Too much accuracy requested for precision of machine.
${\mathbf{ifail}}=17$
A solution component has become zero when a purely relative tolerance (zero absolute tolerance) was selected for that component.
${\mathbf{ifail}}=18$
The error test failed repeatedly with $\left|h\right|=\mathit{hmin}$.
${\mathbf{ifail}}=19$
The corrector repeatedly failed to converge with $\left|h\right|=\mathit{hmin}$.
${\mathbf{ifail}}=20$
The iteration matrix is singular.
${\mathbf{ifail}}=21$
The corrector could not converge and the error test failed repeatedly.
${\mathbf{ifail}}=22$
ires was set to $-1$ during a call to res and could not be resolved.
${\mathbf{ifail}}=23$
ires was set to $-2$ during a call to res.
${\mathbf{ifail}}=24$
The initial ydot could not be computed.
${\mathbf{ifail}}=25$
Repeated occurrences of input constraint violations have been detected.
${\mathbf{ifail}}=-99$
${\mathbf{ifail}}=-399$
Your licence key may have expired or may not have been installed correctly.
${\mathbf{ifail}}=-999$
Dynamic memory allocation failed.

## 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{itol}}=0$ with ${\mathbf{atol}}\left(1\right)$ small but positive).

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}}×{\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}}×\left({\mathbf{ml}}+{\mathbf{mu}}+1\right)$ and ${{\mathbf{neq}}}^{2}$ respectively.

## Example

For this function two examples are presented. There is a single example program for nag_ode_dae_dassl_gen (d02ne), with a main program and the code to solve the two example problems given in Example 1 (EX1) and Example 2 (EX2).
Example 1 (EX1)
This example solves the well-known stiff Robertson problem written in implicit form
 $r1 = -0.04a + 1.0E4bc - a′ r2 = 0.04a - 1.0E4bc - 3.0E7⁢b2 - b′ r3 = 3.0E7⁢b2 - c′$
with initial conditions $a=1.0$ and $b=c=0.0$ over the range $\left[0,0.1\right]$ the BDF method (setup function nag_ode_dae_dassl_setup (d02mw) and nag_ode_dae_dassl_linalg (d02np)).
Example 2 (EX2)
This example illustrates the use of nag_ode_dae_dassl_gen (d02ne) to solve a simple algebraic problem by continuation. The equation $4-2y+0.1{e}^{y}t=0$ from $t=0$ (where $y=2$) to $t=1$.
```function d02ne_example

fprintf('d02ne example results\n\n');

fprintf('Example 1\n');
ex1;
fprintf('\nExample 2\n');
ex2;

function ex1
% Stiff Robertson problem

% Initialize the problem, specifying that the Jacobian is to be
% evaluated analytically using the provided routine jac.
neq    = int64(3);
maxord = int64(5);
jceval = 'Analytic';
hmax   = 0;
h0     = 0;
itol   = int64(1);
lcom   = int64(200);
[icom, com, ifail] = d02mw(...
neq, maxord, jceval, hmax, h0, itol, lcom);

% Specify that the Jacobian is banded
mu = int64(2);
ml = int64(1);
[icom, ifail] = d02np(neq, ml, mu, icom);

% Set initial values
rtol = [1e-3; 1e-3; 1e-3];
atol = [1e-6; 1e-6; 1e-6];
y    = [1; 0; 0];
ydot = zeros(neq,1);
t    = 0;
tout = 0.02;

% Use the user parameter to pass the band dimensions through to jac.
% An alternative would be to hard code values for ml and mu in jac.
user = {ml, mu};

fprintf('\n    t              y\n');
fprintf('%8.4f   %12.6f %12.6f %12.6f\n', t, y);

% Obtain the solution at 5 equally spaced values of T.
for j = 1:5
if ifail == 0
[t, y, ydot, rtol, atol, itask, icom, com, user, ifail] = ...
d02ne(...
t, tout, y, ydot, rtol, atol, itask, @res1, @jac1, ...
icom, com, 'user', user);
fprintf('%8.4f   %12.6f %12.6f %12.6f\n', t, y);
tout = tout + 0.02;
icom = d02mc(icom);
end
end

function ex2

% Setup
neq = int64(1);
maxord = int64(5);
jceval = 'Analytic';
hmax   = 0;
h0     = 0;
itol   = int64(1);
lcom   = int64(200);
[icom, com, ifail] = d02mw(...
neq, maxord, jceval, hmax, h0, itol, lcom);

% Set initial values
rtol = [0];
atol = [1e-8];
t      = 0;
tout = 0.2;
y = [2];
ydot = [0];
fprintf('\n    t            y\n');
fprintf('%8.4f   %12.6f\n', t, y);

% Obtain the solution at 5 equally spaced values on continuation path.
for j = 1:5
if ifail == 0
[t, y, ydot, rtol, atol, itask, icom, com, user, ifail] = ...
d02ne(...
t, tout, y, ydot, rtol, atol, itask, @res2, @jac2, ...
icom, com);
fprintf('%8.4f   %12.6f\n', t, y);
tout = tout + 0.2;
icom = d02mc(icom);
end
end

function [pd, user] = jac1(neq, t, y, ydot, pd, cj, user)
ml = user{1};
mu = user{2};

stride = 2*ml+mu+1;
% Main diagonal pdfull(i,i), i=1,neq
md = mu + ml + 1;
pd(md) = -0.04 - cj;
pd(md+stride) = -1.0e4*y(3) - 6.0e7*y(2) - cj;
pd(md+2*stride) = -cj;
% 1 sub-diagonal pdfull(i+1:i), i=1,neq-1
ms = md + 1;
pd(ms) = 0.04;
pd(ms+stride) = 6.0e7*y(2);
% First super-diagonal pdfull(i-1,i), i=2, neq
ms = md - 1;
pd(ms+stride) = 1.0e4*y(3);
pd(ms+2*stride) = -1.0e4*y(2);
% Second super-diagonal pdfull(i-2,i), i=3, neq
ms = md - 2;
pd(ms+2*stride) = 1.0e4*y(2);

function [r, ires, user] = res1(neq, t, y, ydot, ires, user)
r = zeros(neq, 1);
r(1) = -0.04*y(1) + 1.0e4*y(2)*y(3)                   - ydot(1);
r(2) =  0.04*y(1) - 1.0e4*y(2)*y(3) - 3.0e7*y(2)*y(2) - ydot(2);
r(3) =                                3.0e7*y(2)*y(2) - ydot(3);

function [pd, user] = jac2(neq, t, y, ydot, pd, cj, user)
pd(1) = -2*y(1) + 0.1*t*y(1)*exp(y(1));

function [r, ires, user] = res2(neq, t, y, ydot, ires, user)
r(1) = 4 - y(1)^2 + t*0.1*exp(y(1));
```
```d02ne example results

Example 1

t              y
0.0000       1.000000     0.000000     0.000000
0.0200       0.999204     0.000036     0.000760
0.0400       0.998415     0.000036     0.001549
0.0600       0.997631     0.000036     0.002333
0.0800       0.996852     0.000036     0.003112
0.1000       0.996080     0.000036     0.003884

Example 2

t            y
0.0000       2.000000
0.2000       2.038016
0.4000       2.078379
0.6000       2.121462
0.8000       2.167736
1.0000       2.217821