Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

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.

None.

## Parameters

### Compulsory Input Parameters

1:     $\mathrm{neq}$int64int32nag_int scalar
The number of differential-algebraic equations to be solved.
Constraint: $1\le {\mathbf{neq}}$.
2:     $\mathrm{ml}$int64int32nag_int scalar
${m}_{L}$, the number of subdiagonals in the band.
Constraint: $0\le {\mathbf{ml}}\le {\mathbf{neq}}-1$.
3:     $\mathrm{mu}$int64int32nag_int scalar
${m}_{U}$, the number of superdiagonals in the band.
Constraint: $0\le {\mathbf{mu}}\le {\mathbf{neq}}-1$.
4:     $\mathrm{icom}\left({\mathbf{licom}}\right)$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:     $\mathrm{licom}$int64int32nag_int scalar
Default: the dimension of the array icom.
The dimension of the array icom.
Constraint: ${\mathbf{licom}}\ge 50+{\mathbf{neq}}$.

### Output Parameters

1:     $\mathrm{icom}\left({\mathbf{licom}}\right)$int64int32nag_int array
2:     $\mathrm{ifail}$int64int32nag_int scalar
${\mathbf{ifail}}={\mathbf{0}}$ unless the function detects an error (see Error Indicators and Warnings).

## Error Indicators and Warnings

Errors or warnings detected by the function:
${\mathbf{ifail}}=1$
Constraint: ${\mathbf{neq}}\ge 1$.
${\mathbf{ifail}}=2$
Constraint: ${\mathbf{ml}}\le {\mathbf{neq}}-1$.
Constraint: ${\mathbf{ml}}\ge 0$.
${\mathbf{ifail}}=3$
Constraint: ${\mathbf{mu}}\le {\mathbf{neq}}-1$.
Constraint: ${\mathbf{mu}}\ge 0$.
${\mathbf{ifail}}=4$
Either the initialization function has not been called prior to the first call of this function or the communication array has become corrupted.
${\mathbf{ifail}}=5$
On entry, licom is too small.
${\mathbf{ifail}}=-99$
${\mathbf{ifail}}=-399$
Your licence key may have expired or may not have been installed correctly.
${\mathbf{ifail}}=-999$
Dynamic memory allocation failed.

Not applicable.

None.

## Example

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

% 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

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