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_linalg (d02np)


    1  Purpose
    2  Syntax
    7  Accuracy
    9  Example


nag_ode_dae_dassl_linalg (d02np) is a setup function which you must call prior to nag_ode_dae_dassl_gen (d02ne) and after a call to nag_ode_dae_dassl_setup (d02mw), if the Jacobian is to be considered as having a banded structure.


[icom, ifail] = d02np(neq, ml, mu, icom, 'licom', licom)
[icom, ifail] = nag_ode_dae_dassl_linalg(neq, ml, mu, icom, 'licom', licom)


A call to nag_ode_dae_dassl_linalg (d02np) specifies that the Jacobian to be used is banded in structure. If nag_ode_dae_dassl_linalg (d02np) is not called before a call to nag_ode_dae_dassl_gen (d02ne) then the Jacobian is assumed to be full.




Compulsory Input Parameters

1:     neq int64int32nag_int scalar
The number of differential-algebraic equations to be solved.
Constraint: 1neq.
2:     ml int64int32nag_int scalar
mL, the number of subdiagonals in the band.
Constraint: 0mlneq-1.
3:     mu int64int32nag_int scalar
mU, the number of superdiagonals in the band.
Constraint: 0muneq-1.
4:     icomlicom int64int32nag_int array
icom is used to communicate details of the integration from nag_ode_dae_dassl_setup (d02mw) and details of the banded structure of the Jacobian to nag_ode_dae_dassl_gen (d02ne).

Optional Input Parameters

1:     licom int64int32nag_int scalar
Default: the dimension of the array icom.
The dimension of the array icom.
Constraint: licom50+neq.

Output Parameters

1:     icomlicom int64int32nag_int array
2:     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: mlneq-1.
Constraint: ml0.
Constraint: muneq-1.
Constraint: mu0.
Either the initialization function has not been called prior to the first call of this function or the communication array has become corrupted.
On entry, licom is too small.
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



See Example in nag_ode_dae_dassl_gen (d02ne) and nag_ode_dae_dassl_setup (d02mw).
function d02np_example

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

% 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(1)        y(2)        y(3)   \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, @res, @jac, icom, com, 'user', user);
    fprintf(' %8.4f   %12.6f %12.6f %12.6f\n', t, y);
    tout = tout + 0.02;
    icom = d02mc(icom);

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

function [pd, user] = jac(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] = res(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);
d02np example results

    t            y(1)        y(2)        y(3)   
   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

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