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_gen (d02ne)

 Contents

    1  Purpose
    2  Syntax
    7  Accuracy
    9  Example

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 (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 Ft,y,ydot=0). The function solves the system from t=t to t=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(...);
  if (tout>=tend || itask<0)
    break
  end
  if (itask ~= 1)
    tout = min(tout + dt, tend);
  end
  % Print the solution
  [...] = d02mc(...);
end
.
.
.

References

None.

Parameters

Compulsory Input Parameters

1:     t – double scalar
On initial entry: the initial value of the independent variable, t.
2:     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: toutt.
3:     yneq – double array
On initial entry: the vector of initial values of the dependent variables y.
4:     ydotneq – double array
On initial entry: ydot must contain approximations to the time derivatives y' of the vector y evaluated at the initial value of the independent variable.
5:     rtol: – 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 neq if itol=true and at least 1 if itol=false
The relative local error tolerance.
Constraint: rtoli0.0, for i=1,2,,n
where n=neq when itol=true and n=1 otherwise.
6:     atol: – 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 neq if itol=true and at least 1 if itol=false
The absolute local error tolerance.
Constraint: atoli0.0, for i=1,2,,n
where n=neq when itol=true and n=1 otherwise.
7:     itask int64int32nag_int scalar
On initial entry: need not be set.
8:     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:     neq int64int32nag_int scalar
The number of differential-algebraic equations being solved.
2:     t – double scalar
t, the current value of the independent variable.
3:     yneq – double array
yi, for i=1,2,,neq, the current solution component.
4:     ydotneq – double array
The derivative of the solution at the current point t.
5:     ires int64int32nag_int scalar
Is always equal to zero.
6:     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:     rneq – double array
ri must contain the ith component of R, for i=1,2,,neq where
R = F t,y,ydot .  
2:     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 ifail=23.
3:     user – Any MATLAB object
9:     jac – function handle or string containing name of m-file
Evaluates the matrix of partial derivatives, J, where
J ij = Fi yj + cj × Fi yj ,  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:     neq int64int32nag_int scalar
The number of differential-algebraic equations being solved.
2:     t – double scalar
t, the current value of the independent variable.
3:     yneq – double array
yi, for i=1,2,,neq, the current solution component.
4:     ydotneq – double array
The derivative of the solution at the current point t.
5:     pd: – double array
pd is preset to zero before the call to jac.
6:     cj – double scalar
cj is a scalar constant which will be defined in nag_ode_dae_dassl_gen (d02ne).
7:     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:     pd: – double array
If the Jacobian is full then pdj-1×neq+i=Jij, for i=1,2,,neq and j=1,2,,neq; if the Jacobian is banded then pdj-1×2ml+mu+1+ml+mu+i-j+1=Jij, for max1,j-muiminn,j+ml; (see also in nag_lapack_dgbtrf (f07bd)).
2:     user – Any MATLAB object
10:   icom50+neq 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:
icom22
The order of the method to be attempted on the next step.
icom23
The order of the method used on the last step.
icom26
The number of steps taken so far.
icom27
The number of calls to res so far.
icom28
The number of evaluations of the matrix of partial derivatives needed so far.
icom29
The total number of error test failures so far.
icom30
The total number of convergence test failures so far.
11:   comlcom – double array
lcom, the dimension of the array, must satisfy the constraint lcom40+maxorder+4×neq+neq×p+q where maxorder is the maximum order that can be used by the integration method (see maxord in nag_ode_dae_dassl_setup (d02mw)); p=neq when the Jacobian is full and p=2×ml+mu+1 when the Jacobian is banded; and, q=neq/ml+mu+1+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:
com3
The step size to be attempted on the next step.
com4
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 (itask=3).

Optional Input Parameters

1:     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: neq1.
2:     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:     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:     yneq – 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=tout).
3:     ydotneq – double array
The time derivatives y' of the vector y at the last integration point.
4:     rtol: – 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 neq if itol=true and at least 1 if itol=false
rtol remains unchanged unless nag_ode_dae_dassl_gen (d02ne) exits with ifail=16 in which case the values may have been increased to values estimated to be appropriate for continuing the integration.
5:     atol: – 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 neq if itol=true and at least 1 if itol=false
atol remains unchanged unless nag_ode_dae_dassl_gen (d02ne) exits with ifail=16 in which case the values may have been increased to values estimated to be appropriate for continuing the integration.
6:     itask int64int32nag_int scalar
The task performed by the integrator on successful completion or an indicator that a problem occurred during integration.
itask=2
The integration to tout was successfully completed (t=tout) by stepping exactly to tout.
itask=3
The integration to tout was successfully completed (t=tout) by stepping past tout. y and ydot are obtained by interpolation.
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:     icom50+neq int64int32nag_int array
8:     comlcom – double array
9:     user – Any MATLAB object
10:   ifail int64int32nag_int scalar
ifail=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:
   ifail=1
Constraint: neq1.
   ifail=3
Constraint: toutt.
tout is behind t in the direction of h.
tout is too close to t to start integration.
   ifail=6
Some element of rtol is less than zero.
   ifail=7
Some element of atol is less than zero.
   ifail=8
A previous call to this function returned with itask=_ and no appropriate action was taken.
   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.
   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.
   ifail=13
com array is of insufficient length; length required =_; actual length lcom=_.
   ifail=14
All elements of rtol and atol are zero.
   ifail=15
Maximum number of steps taken on this call before reaching tout.
   ifail=16
Too much accuracy requested for precision of machine.
   ifail=17
A solution component has become zero when a purely relative tolerance (zero absolute tolerance) was selected for that component.
   ifail=18
The error test failed repeatedly with h=hmin.
   ifail=19
The corrector repeatedly failed to converge with h=hmin.
   ifail=20
The iteration matrix is singular.
   ifail=21
The corrector could not converge and the error test failed repeatedly.
   ifail=22
ires was set to -1 during a call to res and could not be resolved.
   ifail=23
ires was set to -2 during a call to res.
   ifail=24
The initial ydot could not be computed.
   ifail=25
Repeated occurrences of input constraint violations have been detected.
   ifail=-99
An unexpected error has been triggered by this routine. Please contact NAG.
   ifail=-399
Your licence key may have expired or may not have been installed correctly.
   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 itol=0 with atol1 small but positive).

Further Comments

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 neq× ml+ mu+1 2 , while for full systems the cost is proportional to neq3. Note however that for moderately sized problems which are only mildly nonlinear the cost may be dominated by factors proportional to neq×ml+mu+1 and neq2 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.0E7b2 - b r3 = 3.0E7b2 - c  
with initial conditions a=1.0 and b=c=0.0 over the range 0,0.1 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.1eyt=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);

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

  fprintf('\nThe integrator completed task, ITASK = %d\n', itask);

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

  itask = int64(0);
  % 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
  fprintf('\nThe integrator completed task, ITASK = %d\n', itask);

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

The integrator completed task, ITASK = 3

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

The integrator completed task, ITASK = 3

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