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_setup (d02mw)

## Purpose

nag_ode_dae_dassl_setup (d02mw) is a setup function which must be called prior to the integrator nag_ode_dae_dassl_gen (d02ne), if the DASSL implementation of Backward Differentiation Formulae (BDF) is to be used.

## Syntax

[icom, com, ifail] = d02mw(neq, maxord, jceval, hmax, h0, itol, lcom)
[icom, com, ifail] = nag_ode_dae_dassl_setup(neq, maxord, jceval, hmax, h0, itol, lcom)

## Description

This integrator setup function must be called before the first call to the integrator nag_ode_dae_dassl_gen (d02ne). This setup function nag_ode_dae_dassl_setup (d02mw) permits you to define options for the DASSL integrator, such as: whether the Jacobian is to be provided or is to be approximated numerically by the integrator; the initial and maximum step-sizes for the integration; whether relative and absolute tolerances are system wide or per system equation; and the maximum order of BDF method permitted.

None.

## Parameters

### Compulsory Input Parameters

1:     neq – int64int32nag_int scalar
The number of differential-algebraic equations to be solved.
Constraint: neq1${\mathbf{neq}}\ge 1$.
2:     maxord – int64int32nag_int scalar
The maximum order to be used for the BDF method. Orders up to 5th order are available; setting maxord > 5${\mathbf{maxord}}>5$ means that the maximum order used will be 5$5$.
Constraint: 1maxord$1\le {\mathbf{maxord}}$.
3:     jceval – string (length ≥ 1)
Specifies the technique to be used to compute the Jacobian.
jceval = 'N'${\mathbf{jceval}}=\text{'N'}$
The Jacobian is to be evaluated numerically by the integrator.
jceval = 'A'${\mathbf{jceval}}=\text{'A'}$
You must supply a function to evaluate the Jacobian on a call to the integrator.
Only the first character of the actual paramater jceval is passed to nag_ode_dae_dassl_setup (d02mw); hence it is permissible for the actual argument to be more descriptive, e.g., ‘Numerical’ or ‘Analytical’, on a call to nag_ode_dae_dassl_setup (d02mw).
Constraint: jceval = 'N'${\mathbf{jceval}}=\text{'N'}$ or 'A'$\text{'A'}$.
4:     hmax – double scalar
The maximum absolute step size to be allowed. Set hmax = 0.0${\mathbf{hmax}}=0.0$ if this option is not required.
Constraint: hmax0.0${\mathbf{hmax}}\ge 0.0$.
5:     h0 – double scalar
The step size to be attempted on the first step. Set h0 = 0.0${\mathbf{h0}}=0.0$ if the initial step size is calculated internally.
6:     itol – int64int32nag_int scalar
A value to indicate the form of the local error test.
itol = 0${\mathbf{itol}}=0$
rtol and atol are single element vectors.
itol = 1${\mathbf{itol}}=1$
rtol and atol are vectors. This should be chosen if you want to apply different tolerances to each equation in the system.
Note: the tolerances must either both be single element vectors or both be vectors of length neq.
Constraint: itol = 0${\mathbf{itol}}=0$ or 1$1$.
7:     lcom – int64int32nag_int scalar
The dimension of the array com as declared in the (sub)program from which nag_ode_dae_dassl_setup (d02mw) is called.
Constraints:
the array com must be large enough for the requirements of nag_ode_dae_dassl_gen (d02ne). That is:
• if the system Jacobian is dense, lcom 40 + (maxord + 4) × neq + neq2 ${\mathbf{lcom}}\ge 40+\left({\mathbf{maxord}}+4\right)×{\mathbf{neq}}+{{\mathbf{neq}}}^{2}$;
• if the system Jacobian is banded,   lcom 40 + (maxord + 4) × neq + (2 × ml + mu + 1) × neq + 2 × (neq / (ml + mu + 1) + 1) $\phantom{\rule{0ex}{0ex}}{\mathbf{lcom}}\ge 40+\left({\mathbf{maxord}}+4\right)×{\mathbf{neq}}+\left(2×{\mathbf{ml}}+{\mathbf{mu}}+1\right)×{\mathbf{neq}}+2×\left({\mathbf{neq}}/\left({\mathbf{ml}}+{\mathbf{mu}}+1\right)+1\right)$.
Here ml and mu are the lower and upper bandwidths respectively that are to be specified in a subsequent call to nag_ode_dae_dassl_linalg (d02np).

None.

licom

### Output Parameters

1:     icom(licom) – int64int32nag_int array
licomneq + 50$\mathit{licom}\ge {\mathbf{neq}}+50$.
Used to communicate details of the task to be carried out to the integration function nag_ode_dae_dassl_gen (d02ne).
2:     com(lcom) – double array
Used to communicate problem parameters to the integration function nag_ode_dae_dassl_gen (d02ne). This must be the same communication array as the array com supplied to nag_ode_dae_dassl_gen (d02ne). In particular, the values of hmax and h0 are contained in com.
3:     ifail – int64int32nag_int scalar
${\mathrm{ifail}}={\mathbf{0}}$ unless the function detects an error (see [Error Indicators and Warnings]).

## Error Indicators and Warnings

Errors or warnings detected by the function:
ifail = 1${\mathbf{ifail}}=1$
 On entry, neq < 1${\mathbf{neq}}<1$.
ifail = 2${\mathbf{ifail}}=2$
 On entry, maxord < 1${\mathbf{maxord}}<1$, or maxord > 5${\mathbf{maxord}}>5$.
ifail = 3${\mathbf{ifail}}=3$
 On entry, jceval ≠ 'N'${\mathbf{jceval}}\ne \text{'N'}$ or 'A'$\text{'A'}$.
ifail = 4${\mathbf{ifail}}=4$
 On entry, hmax < 0.0${\mathbf{hmax}}<0.0$.
ifail = 6${\mathbf{ifail}}=6$
 On entry, itol ≠ 0${\mathbf{itol}}\ne 0$ or 1$1$.
ifail = 8${\mathbf{ifail}}=8$
 On entry, licom < neq + 50$\mathit{licom}<{\mathbf{neq}}+50$.

Not applicable.

None.

## Example

```function nag_ode_dae_dassl_setup_example
neq = 3;
maxord = 5;
mu = 2;
ml = 1;
lcom = 40+(maxord+4)*neq+(2*ml+mu+1)*neq+2*(neq/(ml+mu+1)+1)
itol = int64(1);
rtol = [1e-3; 1e-3; 1e-3];
atol = [1e-6; 1e-6; 1e-6];
ydot = zeros(neq,1);
% Set initial values
y    = [1; 0; 0];
% Initialize the problem, specifying that the Jacobian is to be
% evaluated analytically using the provided routine jac.
jceval = 'Analytic';
hmax = 0;
ho = 0;
t = 0;
tout = 0.02;
[icom, com, ifail] = ...
nag_ode_dae_dassl_setup(int64(neq), int64(maxord), jceval, hmax, ho, itol, int64(lcom));

% Specify that the Jacobian is banded
if ifail == 0
[icom, ifail] = nag_ode_dae_dassl_linalg(int64(neq), int64(ml), int64(mu), icom);
end

% 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(1)        y(2)        y(3)   \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] = ...
nag_ode_dae_dassl_gen(t, tout, y, ydot, rtol, atol, itask, @res, @jac, ...
icom, com, 'user', user);
fprintf(' %8.4f   %12.6f %12.6f %12.6f\n', t, y);
tout = tout + 0.02;
icom = nag_ode_dae_dassl_cont(icom);
end
end

function [pd, user] = jac(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] = res(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);
```
```

lcom =

85.5000

t            y(1)        y(2)        y(3)
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

```
```function d02mw_example
neq = 3;
maxord = 5;
mu = 2;
ml = 1;
lcom = 40+(maxord+4)*neq+(2*ml+mu+1)*neq+2*(neq/(ml+mu+1)+1)
itol = int64(1);
rtol = [1e-3; 1e-3; 1e-3];
atol = [1e-6; 1e-6; 1e-6];
ydot = zeros(neq,1);
% Set initial values
y    = [1; 0; 0];
% Initialize the problem, specifying that the Jacobian is to be
% evaluated analytically using the provided routine jac.
jceval = 'Analytic';
hmax = 0;
ho = 0;
t = 0;
tout = 0.02;
[icom, com, ifail] = d02mw(int64(neq), int64(maxord), jceval, hmax, ho, itol, int64(lcom));

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

% 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(1)        y(2)        y(3)   \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, @res, @jac, ...
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 [pd, user] = jac(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] = res(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);
```
```

lcom =

85.5000

t            y(1)        y(2)        y(3)
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