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)

Purpose

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.

Syntax

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

Description

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.

References

None.

Parameters

Compulsory Input Parameters

1:     neq – int64int32nag_int scalar
The number of differential-algebraic equations to be solved.
Constraint: 1neq1neq.
2:     ml – int64int32nag_int scalar
mLmL, the number of subdiagonals in the band.
Constraint: 0mlneq10mlneq-1.
3:     mu – int64int32nag_int scalar
mUmU, the number of superdiagonals in the band.
Constraint: 0muneq10muneq-1.
4:     icom(licom) – int64int32nag_int array
licom, the dimension of the array, must satisfy the constraint licom50 + neqlicom50+neq.
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 as declared in the (sub)program from which nag_ode_dae_dassl_linalg (d02np) is called.
Constraint: licom50 + neqlicom50+neq.

Input Parameters Omitted from the MATLAB Interface

None.

Output Parameters

1:     icom(licom) – int64int32nag_int array
2:     ifail – int64int32nag_int scalar
ifail = 0ifail=0 unless the function detects an error (see [Error Indicators and Warnings]).

Error Indicators and Warnings

Errors or warnings detected by the function:
  ifail = 1ifail=1
On entry,neq < 1neq<1.
  ifail = 2ifail=2
On entry,ml < 0ml<0 or ml > neq1ml>neq-1.
  ifail = 3ifail=3
On entry,mu < 0mu<0 or mu > neq1mu>neq-1.
  ifail = 4ifail=4
Either nag_ode_dae_dassl_setup (d02mw) has not been called before this call or the communication array icom has been corrupted.
  ifail = 5ifail=5
On entry,licom < 50 + neqlicom<50+neq.

Accuracy

Not applicable.

Further Comments

None.

Example

function nag_ode_dae_dassl_linalg_example
neq = 3;
maxord = 5;
mu = 2;
ml = 1;
lcom = 40+(maxord+4)*neq+(2*ml+mu+1)*neq+2*(neq/(ml+mu+1)+1)
itol = int64(1);
rtol = [1e-3; 1e-3; 1e-3];
atol = [1e-6; 1e-6; 1e-6];
ydot = zeros(neq,1);
% Set initial values
y    = [1; 0; 0];
% Initialize the problem, specifying that the Jacobian is to be
% evaluated analytically using the provided routine jac.
jceval = 'Analytic';
hmax = 0;
ho = 0;
t = 0;
tout = 0.02;
[icom, com, ifail] = ...
    nag_ode_dae_dassl_setup(int64(neq), int64(maxord), jceval, hmax, ho, itol, int64(lcom));

% Specify that the Jacobian is banded
if ifail == 0
  [icom, ifail] = nag_ode_dae_dassl_linalg(int64(neq), int64(ml), int64(mu), icom);
end

% 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] = ...
      nag_ode_dae_dassl_gen(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 = nag_ode_dae_dassl_cont(icom);
  end
end

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

lcom =

   85.5000


    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

function d02np_example
neq = 3;
maxord = 5;
mu = 2;
ml = 1;
lcom = 40+(maxord+4)*neq+(2*ml+mu+1)*neq+2*(neq/(ml+mu+1)+1)
itol = int64(1);
rtol = [1e-3; 1e-3; 1e-3];
atol = [1e-6; 1e-6; 1e-6];
ydot = zeros(neq,1);
% Set initial values
y    = [1; 0; 0];
% Initialize the problem, specifying that the Jacobian is to be
% evaluated analytically using the provided routine jac.
jceval = 'Analytic';
hmax = 0;
ho = 0;
t = 0;
tout = 0.02;
[icom, com, ifail] = d02mw(int64(neq), int64(maxord), jceval, hmax, ho, itol, int64(lcom));

% Specify that the Jacobian is banded
if ifail == 0
  [icom, ifail] = d02np(int64(neq), int64(ml), int64(mu), icom);
end

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

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

lcom =

   85.5000


    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–2013