NAG Toolbox: nag_ode_dae_dassl_setup (d02mw)

Purpose

Syntax

Description

References

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.

See nag_ode_dae_dassl_gen (d02ne).

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)\times \mathbf{neq}+{\mathbf{neq}}^2$;
• if the system Jacobian is banded, $\mathbf{lcom}\ge 40+\left(\mathbf{maxord}+4\right)\times \mathbf{neq}+\left(2\times \mathbf{ml}+\mathbf{mu}+1\right)\times \mathbf{neq}+2\times \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).

Optional Input Parameters

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

   $\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$

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

   $\mathbf{ifail}=-399$

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

   $\mathbf{ifail}=-999$

Dynamic memory allocation failed.

Accuracy

Further Comments

Example

Open in the MATLAB editor: d02mw_example

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;
itask = int64(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);

fprintf('\nThe integrator returned with itask = %d\n', itask);
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