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:     $\mathrm{neq}$int64int32nag_int scalar
The number of differential-algebraic equations to be solved.
Constraint: ${\mathbf{neq}}\ge 1$.
2:     $\mathrm{maxord}$int64int32nag_int scalar
The maximum order to be used for the BDF method. Orders up to 5th order are available; setting ${\mathbf{maxord}}>5$ means that the maximum order used will be $5$.
Constraint: $1\le {\mathbf{maxord}}$.
3:     $\mathrm{jceval}$ – string (length ≥ 1)
Specifies the technique to be used to compute the Jacobian.
${\mathbf{jceval}}=\text{'N'}$
The Jacobian is to be evaluated numerically by the integrator.
${\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: ${\mathbf{jceval}}=\text{'N'}$ or $\text{'A'}$.
4:     $\mathrm{hmax}$ – double scalar
The maximum absolute step size to be allowed. Set ${\mathbf{hmax}}=0.0$ if this option is not required.
Constraint: ${\mathbf{hmax}}\ge 0.0$.
5:     $\mathrm{h0}$ – double scalar
The step size to be attempted on the first step. Set ${\mathbf{h0}}=0.0$ if the initial step size is calculated internally.
6:     $\mathrm{itol}$int64int32nag_int scalar
A value to indicate the form of the local error test.
${\mathbf{itol}}=0$
rtol and atol are single element vectors.
${\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: ${\mathbf{itol}}=0$ or $1$.
7:     $\mathrm{lcom}$int64int32nag_int scalar
The dimension of the array com.
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, ${\mathbf{lcom}}\ge 40+\left({\mathbf{maxord}}+4\right)×{\mathbf{neq}}+{{\mathbf{neq}}}^{2}$;
• if the system Jacobian is banded, $\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.

### Output Parameters

1:     $\mathrm{icom}\left(\mathit{licom}\right)$int64int32nag_int array
$\mathit{licom}={\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:     $\mathrm{com}\left({\mathbf{lcom}}\right)$ – 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:     $\mathrm{ifail}$int64int32nag_int scalar
${\mathbf{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:
${\mathbf{ifail}}=1$
Constraint: ${\mathbf{neq}}\ge 1$.
${\mathbf{ifail}}=2$
Constraint: ${\mathbf{maxord}}\ge 1$.
${\mathbf{ifail}}=3$
Constraint: ${\mathbf{jceval}}=\text{'N'}$ or $\text{'A'}$.
${\mathbf{ifail}}=4$
Constraint: ${\mathbf{hmax}}\ge 0.0$.
${\mathbf{ifail}}=6$
Constraint: ${\mathbf{itol}}=0$ or $1$.
${\mathbf{ifail}}=8$
Constraint: $\mathit{licom}\ge 50+{\mathbf{neq}}$.
${\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.

Not applicable.

None.

## Example

This example solves the plane pendulum problem, defined by the following equations:
 $x′ = u y′ = v u′ = -λx v′ = -λy-1 x2+y2 = 1.$
Differentiating the algebraic constraint once, a new algebraic constraint is obtained
 $xu+yv=0 .$
Differentiating the algebraic constraint one more time, substituting for ${x}^{\prime }$, ${y}^{\prime }$, ${u}^{\prime }$, ${v}^{\prime }$ and using ${x}^{2}+{y}^{2}-1=0$, the corresponding DAE system includes the differential equations and the algebraic equation in $\lambda$:
 $u2 + v2 - λ - y = 0 .$
We solve the reformulated DAE system
 $y1′ = y3 y2′ = y4 y3′ = -y5×y1 y4′ = -y5×y2-1 y32 + y42 - y5 - y2 = 0.$
For our experiments, we take consistent initial values
 $y10 = 1 , ​ y20 = 0 , ​ y30 = 0 , ​ y40 = 1 ​ and ​ y50 = 1$
at $t=0$.
```function d02mw_example

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

% Initialize the problem, specifying that the Jacobian is to be
% evaluated analytically using the provided routine jac.
neq    = int64(5);
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(1:neq) = 1e-8;
atol = rtol;
y    = [1  0  0  1  1];
ydot = zeros(neq,1);
t    = 0;
tout = 1;

fprintf('\n    t            y\n');
fprintf('%6.4f   %10.6f %10.6f %10.6f %10.6f %10.6f\n', t, y);

[t, y, ydot, rtol, atol, itask, icom, com, user, ifail] = ...
d02ne(...
t, tout, y, ydot, rtol, atol, itask, @res, @jac, icom, com);
fprintf('%6.4f   %10.6f %10.6f %10.6f %10.6f %10.6f\n', t, y);

g1 = y(1)^2+y(2)^2 - 1;
g2 = y(1)*y(3) + y(2)*y(4);
fprintf('The position-level constraint G1 = %12.4e\n',g1);
fprintf('The velocity-level constraint G2 = %12.4e\n',g2);

function [pd, user] = jac(neq, t, y, ydot, pd, cj, user)
% r1
pd(1:neq:neq^2) = [-cj,0,1,0,0];
% r2
pd(2:neq:neq^2) = [0,-cj,0,1,0];
% r3
pd(3:neq:neq^2) = [-y(5),0,-cj,0,-y(1)];
%r4
pd(4:neq:neq^2) = [0,-y(5),0,-cj,-y(2)];
%r5
pd(5:neq:neq^2) = [0,-1,2*y(3),2*y(4),-1];

function [r, ires, user] = res(neq, t, y, ydot, ires, user)
r = zeros(neq, 1);
r(1) = y(3) - ydot(1);
r(2) = y(4) - ydot(2);
r(3) = -y(5)*y(1)     - ydot(3);
r(4) = -y(5)*y(2) - 1 - ydot(4);
r(5) = y(3)^2 + y(4)^2 - y(5) - y(2);
```
```d02mw example results

t            y
0.0000     1.000000   0.000000   0.000000   1.000000   1.000000
1.0000     0.867349   0.497701  -0.033748   0.058813  -0.493103

The integrator returned with itask = 3
The position-level constraint G1 =  -8.5802e-09
The velocity-level constraint G2 =  -3.0051e-08
```