NAG Toolbox: nag_ode_dae_dassl_gen (d02ne)

Purpose

Syntax

Description

.
.
.
% 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

Parameters

Compulsory Input Parameters

1:     $\mathrm{t}$ – double scalar

On initial entry: the initial value of the independent variable, $t$.

2:     $\mathrm{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: $\mathbf{tout}\ne \mathbf{t}$.

3:     $\mathrm{y}\left(\mathbf{neq}\right)$ – double array

On initial entry: the vector of initial values of the dependent variables $y$.

4:     $\mathrm{ydot}\left(\mathbf{neq}\right)$ – double array

On initial entry: ydot must contain approximations to the time derivatives $y\text{'}$ of the vector $y$ evaluated at the initial value of the independent variable.

5:     $\mathrm{rtol}\left(:\right)$ – 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 $\mathbf{neq}$ if $\mathbf{itol}=\mathit{true}$ and at least $1$ if $\mathbf{itol}=\mathit{false}$

The relative local error tolerance.

Constraint: $\mathbf{rtol}\left(\mathit{i}\right)\ge 0.0$, for $\mathit{i}=1,2,\dots ,n$

where $n=\mathbf{neq}$ when $\mathbf{itol}=\mathit{true}$ and $n=1$ otherwise.

6:     $\mathrm{atol}\left(:\right)$ – 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 $\mathbf{neq}$ if $\mathbf{itol}=\mathit{true}$ and at least $1$ if $\mathbf{itol}=\mathit{false}$

The absolute local error tolerance.

Constraint: $\mathbf{atol}\left(\mathit{i}\right)\ge 0.0$, for $\mathit{i}=1,2,\dots ,n$

where $n=\mathbf{neq}$ when $\mathbf{itol}=\mathit{true}$ and $n=1$ otherwise.

7:     $\mathrm{itask}$ – int64int32nag_int scalar

On initial entry: need not be set.

8:     $\mathrm{res}$ – function handle or string containing name of m-file

res must evaluate the residual

$$R=F\left(t,y,y^{\prime }\right)\text{.}$$

[r, ires, user] = res(neq, t, y, ydot, ires, user)

Input Parameters

1:     $\mathrm{neq}$ – int64int32nag_int scalar

The number of differential-algebraic equations being solved.

2:     $\mathrm{t}$ – double scalar

$t$, the current value of the independent variable.

3:     $\mathrm{y}\left(\mathbf{neq}\right)$ – double array

$y_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,\mathbf{neq}$, the current solution component.

4:     $\mathrm{ydot}\left(\mathbf{neq}\right)$ – double array

The derivative of the solution at the current point $t$.

5:     $\mathrm{ires}$ – int64int32nag_int scalar

Is always equal to zero.

6:     $\mathrm{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:     $\mathrm{r}\left(\mathbf{neq}\right)$ – double array

$\mathbf{r}\left(\mathit{i}\right)$ must contain the $\mathit{i}$th component of $R$, for $\mathit{i}=1,2,\dots ,\mathbf{neq}$ where

$$R=F\left(\mathbf{t},\mathbf{y},\mathbf{ydot}\right)\text{.}$$

2:     $\mathrm{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 $\mathbf{ifail}=\mathbf{23}$.

3:     $\mathrm{user}$ – Any MATLAB object

9:     $\mathrm{jac}$ – function handle or string containing name of m-file

Evaluates the matrix of partial derivatives, $J$, where

$$J_{ij}=\frac{\partial F_i}{\partial y_j}+\mathbf{cj}\times \frac{\partial F_i}{\partial {y^{\prime }}_j}\text{,\hspace{1em}}i,j=1,2,\dots ,\mathbf{neq}\text{.}$$

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:     $\mathrm{neq}$ – int64int32nag_int scalar

The number of differential-algebraic equations being solved.

2:     $\mathrm{t}$ – double scalar

$t$, the current value of the independent variable.

3:     $\mathrm{y}\left(\mathbf{neq}\right)$ – double array

$y_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,\mathbf{neq}$, the current solution component.

4:     $\mathrm{ydot}\left(\mathbf{neq}\right)$ – double array

The derivative of the solution at the current point $t$.

5:     $\mathrm{pd}\left(:\right)$ – double array

pd is preset to zero before the call to jac.

6:     $\mathrm{cj}$ – double scalar

cj is a scalar constant which will be defined in nag_ode_dae_dassl_gen (d02ne).

7:     $\mathrm{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:     $\mathrm{pd}\left(:\right)$ – double array

If the Jacobian is full then $\mathbf{pd}\left(\left(\mathit{j}-1\right)\times \mathbf{neq}+\mathit{i}\right)=J_{\mathit{i}\mathit{j}}$, for $\mathit{i}=1,2,\dots ,\mathbf{neq}$ and $\mathit{j}=1,2,\dots ,\mathbf{neq}$; if the Jacobian is banded then $\mathbf{pd}\left(\left(j-1\right)\times \left(2\mathbf{ml}+\mathbf{mu}+1\right)+\mathbf{ml}+\mathbf{mu}+i-j+1\right)=J_{ij}$, for $\max \left(1,j-\mathbf{mu}\right)\le i\le \min \left(n,j+\mathbf{ml}\right)$; (see also in nag_lapack_dgbtrf (f07bd)).

2:     $\mathrm{user}$ – Any MATLAB object

10:   $\mathrm{icom}\left(50+\mathbf{neq}\right)$ – 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:

$\mathbf{icom}\left(22\right)$

The order of the method to be attempted on the next step.

$\mathbf{icom}\left(23\right)$

The order of the method used on the last step.

$\mathbf{icom}\left(26\right)$

The number of steps taken so far.

$\mathbf{icom}\left(27\right)$

The number of calls to res so far.

$\mathbf{icom}\left(28\right)$

The number of evaluations of the matrix of partial derivatives needed so far.

$\mathbf{icom}\left(29\right)$

The total number of error test failures so far.

$\mathbf{icom}\left(30\right)$

The total number of convergence test failures so far.

11:   $\mathrm{com}\left(\mathit{lcom}\right)$ – double array

lcom, the dimension of the array, must satisfy the constraint $\mathit{lcom}\ge 40+\left(\mathit{maxorder}+4\right)\times \mathbf{neq}+\mathbf{neq}\times p+q$ where $\mathit{maxorder}$ is the maximum order that can be used by the integration method (see maxord in nag_ode_dae_dassl_setup (d02mw)); $p=\mathbf{neq}$ when the Jacobian is full and $p=\left(2\times \mathbf{ml}+\mathbf{mu}+1\right)$ when the Jacobian is banded; and, $q=\left(\mathbf{neq}/\left(\mathbf{ml}+\mathbf{mu}+1\right)\right)+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:

$\mathbf{com}\left(3\right)$

The step size to be attempted on the next step.

$\mathbf{com}\left(4\right)$

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 ($\mathbf{itask}=3$).

Optional Input Parameters

1:     $\mathrm{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: $\mathbf{neq}\ge 1$.

2:     $\mathrm{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:     $\mathrm{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:     $\mathrm{y}\left(\mathbf{neq}\right)$ – 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=\mathbf{tout}$).

3:     $\mathrm{ydot}\left(\mathbf{neq}\right)$ – double array

The time derivatives $y\text{'}$ of the vector $y$ at the last integration point.

4:     $\mathrm{rtol}\left(:\right)$ – 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 $\mathbf{neq}$ if $\mathbf{itol}=\mathit{true}$ and at least $1$ if $\mathbf{itol}=\mathit{false}$

rtol remains unchanged unless nag_ode_dae_dassl_gen (d02ne) exits with $\mathbf{ifail}=\mathbf{16}$ in which case the values may have been increased to values estimated to be appropriate for continuing the integration.

5:     $\mathrm{atol}\left(:\right)$ – 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 $\mathbf{neq}$ if $\mathbf{itol}=\mathit{true}$ and at least $1$ if $\mathbf{itol}=\mathit{false}$

atol remains unchanged unless nag_ode_dae_dassl_gen (d02ne) exits with $\mathbf{ifail}=\mathbf{16}$ in which case the values may have been increased to values estimated to be appropriate for continuing the integration.

6:     $\mathrm{itask}$ – int64int32nag_int scalar

The task performed by the integrator on successful completion or an indicator that a problem occurred during integration.

$\mathbf{itask}=2$

The integration to tout was successfully completed ($\mathbf{t}=\mathbf{tout}$) by stepping exactly to tout.

$\mathbf{itask}=3$

The integration to tout was successfully completed ($\mathbf{t}=\mathbf{tout}$) by stepping past tout. y and ydot are obtained by interpolation.

$\mathbf{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:     $\mathrm{icom}\left(50+\mathbf{neq}\right)$ – int64int32nag_int array

8:     $\mathrm{com}\left(\mathit{lcom}\right)$ – double array

9:     $\mathrm{user}$ – Any MATLAB object

10:   $\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}=3$

Constraint: $\mathbf{tout}\ne \mathbf{t}$.

tout is behind t in the direction of $h$.

tout is too close to t to start integration.

   $\mathbf{ifail}=6$

Some element of rtol is less than zero.

   $\mathbf{ifail}=7$

Some element of atol is less than zero.

   $\mathbf{ifail}=8$

A previous call to this function returned with $\mathbf{itask}=\_$ and no appropriate action was taken.

   $\mathbf{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.

   $\mathbf{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.

   $\mathbf{ifail}=13$

com array is of insufficient length; length required $=\_$; actual length $\mathit{lcom}=\_$.

   $\mathbf{ifail}=14$

All elements of rtol and atol are zero.

   $\mathbf{ifail}=15$

Maximum number of steps taken on this call before reaching tout.

   $\mathbf{ifail}=16$

Too much accuracy requested for precision of machine.

   $\mathbf{ifail}=17$

A solution component has become zero when a purely relative tolerance (zero absolute tolerance) was selected for that component.

   $\mathbf{ifail}=18$

The error test failed repeatedly with $\left|h\right|=\mathit{hmin}$.

   $\mathbf{ifail}=19$

The corrector repeatedly failed to converge with $\left|h\right|=\mathit{hmin}$.

   $\mathbf{ifail}=20$

The iteration matrix is singular.

   $\mathbf{ifail}=21$

The corrector could not converge and the error test failed repeatedly.

   $\mathbf{ifail}=22$

ires was set to $-1$ during a call to res and could not be resolved.

   $\mathbf{ifail}=23$

ires was set to $-2$ during a call to res.

   $\mathbf{ifail}=24$

The initial ydot could not be computed.

   $\mathbf{ifail}=25$

Repeated occurrences of input constraint violations have been detected.

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

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