hide long namesshow long names
hide short namesshow short names
Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

NAG Toolbox: nag_ode_dae_dassl_setup (d02mw)


    1  Purpose
    2  Syntax
    7  Accuracy
    9  Example


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.


[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)


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.




Compulsory Input Parameters

1:     neq int64int32nag_int scalar
The number of differential-algebraic equations to be solved.
Constraint: neq1.
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 means that the maximum order used will be 5.
Constraint: 1maxord.
3:     jceval – string (length ≥ 1)
Specifies the technique to be used to compute the Jacobian.
The Jacobian is to be evaluated numerically by the integrator.
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' or 'A'.
4:     hmax – double scalar
The maximum absolute step size to be allowed. Set hmax=0.0 if this option is not required.
Constraint: hmax0.0.
5:     h0 – double scalar
The step size to be attempted on the first step. Set 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.
rtol and atol are single element vectors.
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 or 1.
7:     lcom int64int32nag_int scalar
The dimension of the array com.
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 ;
  • if the system Jacobian is banded, lcom 40 + maxord+4 × neq + 2×ml+mu+1 × neq + 2 × neq / ml + mu + 1 + 1 .
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


Output Parameters

1:     icomlicom int64int32nag_int array
Used to communicate details of the task to be carried out to the integration function nag_ode_dae_dassl_gen (d02ne).
2:     comlcom – 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
ifail=0 unless the function detects an error (see Error Indicators and Warnings).

Error Indicators and Warnings

Errors or warnings detected by the function:
Constraint: neq1.
Constraint: maxord1.
Constraint: jceval='N' or 'A'.
Constraint: hmax0.0.
Constraint: itol=0 or 1.
Constraint: licom50+neq.
An unexpected error has been triggered by this routine. Please contact NAG.
Your licence key may have expired or may not have been installed correctly.
Dynamic memory allocation failed.


Not applicable.

Further Comments



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, y, u, v and using x2+y2-1=0, the corresponding DAE system includes the differential equations and the algebraic equation in λ:
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;
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] = ...
         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)];
  pd(4:neq:neq^2) = [0,-y(5),0,-cj,-y(2)];
  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

PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

© The Numerical Algorithms Group Ltd, Oxford, UK. 2009–2015