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_ivp_stiff_fulljac_setup (d02ns)

## Purpose

nag_ode_ivp_stiff_fulljac_setup (d02ns) is a setup function which must be called prior to an integrator in sub-chapter D02M–N, if full matrix linear algebra is required.

## Syntax

[rwork, ifail] = d02ns(neq, neqmax, jceval, nwkjac, rwork)
[rwork, ifail] = nag_ode_ivp_stiff_fulljac_setup(neq, neqmax, jceval, nwkjac, rwork)

## Description

nag_ode_ivp_stiff_fulljac_setup (d02ns) defines the linear algebra to be used as full matrix linear algebra, permits you to specify the method for calculating the Jacobian and checks the validity of certain input values.

## References

See the D02M–N sub-chapter Introduction.

## Parameters

### Compulsory Input Parameters

1:     neq – int64int32nag_int scalar
The number of differential equations.
Constraint: 1neqneqmax$1\le {\mathbf{neq}}\le {\mathbf{neqmax}}$.
2:     neqmax – int64int32nag_int scalar
A bound on the maximum number of differential equations to be solved during the integration.
Constraint: ${\mathbf{neqmax}}\ge {\mathbf{neq}}$.
3:     jceval – string (length ≥ 1)
Specifies the technique to be used to compute the Jacobian.
jceval = 'N'${\mathbf{jceval}}=\text{'N'}$
The Jacobian is to be evaluated numerically by the integrator. If this option is used, then the actual argument corresponding to jac in the call to nag_ode_ivp_stiff_exp_fulljac (d02nb) or nag_ode_ivp_stiff_imp_fulljac (d02ng) must be either nag_ode_ivp_stiff_exp_fulljac_dummy_jac (d02nbz) or nag_ode_ivp_stiff_imp_fulljac_dummy_jac (d02ngz) respectively.
jceval = 'A'${\mathbf{jceval}}=\text{'A'}$
You must supply a (sub)program to evaluate the Jacobian on a call to the integrator.
jceval = 'D'${\mathbf{jceval}}=\text{'D'}$
The default choice is to be made. In this case 'D' is interpreted as 'N'.
Only the first character of the actual parameter jceval is passed to nag_ode_ivp_stiff_fulljac_setup (d02ns); hence it is permissible for the actual argument to be more descriptive ‘Numerical’, ‘Analytical’ or ‘Default’ on a call to nag_ode_ivp_stiff_fulljac_setup (d02ns).
Constraint: jceval = 'N'${\mathbf{jceval}}=\text{'N'}$, 'A'$\text{'A'}$ or 'D'$\text{'D'}$.
4:     nwkjac – int64int32nag_int scalar
The size of the workspace array wkjac, which you are supplying to the integrator, as declared in the (sub)program from which nag_ode_ivp_stiff_fulljac_setup (d02ns) is called.
Constraint: nwkjacneqmax × (neqmax + 1)${\mathbf{nwkjac}}\ge {\mathbf{neqmax}}×\left({\mathbf{neqmax}}+1\right)$.
5:     rwork(50 + 4 × neqmax$50+4×{\mathbf{neqmax}}$) – double array
This must be the same workspace array as the array rwork supplied to the integrator. It is used to pass information from the setup function to the integrator and therefore the contents of this array must not be changed before calling the integrator.

None.

None.

### Output Parameters

1:     rwork(50 + 4 × neqmax$50+4×{\mathbf{neqmax}}$) – double array
2:     ifail – int64int32nag_int scalar
${\mathrm{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:
ifail = 1${\mathbf{ifail}}=1$
 On entry, neq < 1${\mathbf{neq}}<1$, or ${\mathbf{neq}}>{\mathbf{neqmax}}$, or nwkjac < neqmax × (neqmax + 1)${\mathbf{nwkjac}}<{\mathbf{neqmax}}×\left({\mathbf{neqmax}}+1\right)$, or jceval ≠ 'N'${\mathbf{jceval}}\ne \text{'N'}$, 'A'$\text{'A'}$ or 'D'$\text{'D'}$.

## Accuracy

Not applicable.

nag_ode_ivp_stiff_fulljac_setup (d02ns) must be called as a setup function before a call to either nag_ode_ivp_stiff_exp_fulljac (d02nb) or nag_ode_ivp_stiff_imp_fulljac (d02ng) and may be called as the linear algebra setup function before a call to either nag_ode_ivp_stiff_exp_revcom (d02nm) or nag_ode_ivp_stiff_imp_revcom (d02nn).

## Example

```function nag_ode_ivp_stiff_fulljac_setup_example
t = 0;
tout = 10;
y = [1; 0; 0];
rwork = zeros(62, 1);
rtol = [0.0001];
atol = [1e-07];
itol = int64(1);
inform = zeros(23, 1, 'int64');
ysave = zeros(3, 6);
wkjac = zeros(12, 1);
itrace = int64(0);
[const, rwork, ifail] = ...
nag_ode_ivp_stiff_bdf(int64(3), int64(6), int64(5), 'Newton', false, zeros(6), ...
10, 1e-10, 10, 0, int64(200), int64(5), 'Average-L2', rwork);
[rwork, ifail] = ...
nag_ode_ivp_stiff_fulljac_setup(int64(3), int64(3), 'Numerical', int64(12), rwork);
[tOut, yOut, ydot, rworkOut, informOut, ysaveOut, wkjacOut, ifail] = ...
nag_ode_ivp_stiff_exp_fulljac(t, tout, y, rwork, rtol, atol, itol, inform, @fcn, ...
ysave, @jac, wkjac, 'nag_ode_ivp_stiff_exp_fulljac_dummy_monit', itask, itrace)

function [f, ires] = fcn(neq, t, y, ires)
% Evaluate derivative vector.
f = zeros(3,1);
f(1) = -0.04d0*y(1) + 1.0d4*y(2)*y(3);
f(2) = 0.04d0*y(1) - 1.0d4*y(2)*y(3) - 3.0d7*y(2)*y(2);
f(3) = 3.0d7*y(2)*y(2);
function p = jac(neq, t, y, h, d, p)
% Evaluate the Jacobian.
p = zeros(neq, neq);
hxd = h*d;
p(1,1) = 1.0d0 - hxd*(-0.04d0);
p(1,2) = -hxd*(1.0d4*y(3));
p(1,3) = -hxd*(1.0d4*y(2));
p(2,1) = -hxd*(0.04d0);
p(2,2) = 1.0d0 - hxd*(-1.0d4*y(3)-6.0d7*y(2));
p(2,3) = -hxd*(-1.0d4*y(2));
p(3,2) = -hxd*(6.0d7*y(2));
p(3,3) = 1.0d0 - hxd*(0.0d0);
```
```

tOut =

10

yOut =

0.8414
0.0000
0.1586

ydot =

-0.0079
-0.0000
0.0079

rworkOut =

1.0e+06 *

0.0000
0.0002
0.0000
0
0
0.0000
0.0000
0.0000
0
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0
0.0000
0.0000
0.0000
0.0000
0.0010
0.0000
0.0000
0.0000
0
0
0.0000
0.0000
0.0000
0.0000
0.0000
0
0
0
0
0.0000
0
0
0
0
0
0.0000
0.0000
0
0.0000
0
0.0118
9.8370
0.0643
-0.0000
0.0000
0.0000
0.0000
-0.0000
-0.0000
0.0000
0.0000
0.0000

informOut =

55
132
17
3
3
79
0
3
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0

ysaveOut =

0.8414   -0.0041    0.0001   -0.0000    0.0000   -0.0000
0.0000   -0.0000    0.0000   -0.0000    0.0000   -0.0000
0.1586    0.0041   -0.0001    0.0000   -0.0000    0.0000

wkjacOut =

1.0100
0
-0.0100
-380.6996
-247.4940
629.1936
-0.0408
-0.0040
2.5831
1.0000
3.0000
3.0000

ifail =

0

```
```function d02ns_example
t = 0;
tout = 10;
y = [1; 0; 0];
rwork = zeros(62, 1);
rtol = [0.0001];
atol = [1e-07];
itol = int64(1);
inform = zeros(23, 1, 'int64');
ysave = zeros(3, 6);
wkjac = zeros(12, 1);
itrace = int64(0);
[const, rwork, ifail] = ...
d02nv(int64(3), int64(6), int64(5), 'Newton', false, zeros(6), ...
10, 1e-10, 10, 0, int64(200), int64(5), 'Average-L2', rwork);
[rwork, ifail] = d02ns(int64(3), int64(3), 'Numerical', int64(12), rwork);
[tOut, yOut, ydot, rworkOut, informOut, ysaveOut, wkjacOut, ifail] = ...
d02nb(t, tout, y, rwork, rtol, atol, itol, inform, @fcn, ...
ysave, @jac, wkjac, 'd02nby', itask, itrace)

function [f, ires] = fcn(neq, t, y, ires)
% Evaluate derivative vector.
f = zeros(3,1);
f(1) = -0.04d0*y(1) + 1.0d4*y(2)*y(3);
f(2) = 0.04d0*y(1) - 1.0d4*y(2)*y(3) - 3.0d7*y(2)*y(2);
f(3) = 3.0d7*y(2)*y(2);
function p = jac(neq, t, y, h, d, p)
% Evaluate the Jacobian.
p = zeros(neq, neq);
hxd = h*d;
p(1,1) = 1.0d0 - hxd*(-0.04d0);
p(1,2) = -hxd*(1.0d4*y(3));
p(1,3) = -hxd*(1.0d4*y(2));
p(2,1) = -hxd*(0.04d0);
p(2,2) = 1.0d0 - hxd*(-1.0d4*y(3)-6.0d7*y(2));
p(2,3) = -hxd*(-1.0d4*y(2));
p(3,2) = -hxd*(6.0d7*y(2));
p(3,3) = 1.0d0 - hxd*(0.0d0);
```
```

tOut =

10

yOut =

0.8414
0.0000
0.1586

ydot =

-0.0079
-0.0000
0.0079

rworkOut =

1.0e+06 *

0.0000
0.0002
0.0000
0
0
0.0000
0.0000
0.0000
0
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0
0.0000
0.0000
0.0000
0.0000
0.0010
0.0000
0.0000
0.0000
0
0
0.0000
0.0000
0.0000
0.0000
0.0000
0
0
0
0
0.0000
0
0
0
0
0
0.0000
0.0000
0
0.0000
0
0.0118
9.8370
0.0643
-0.0000
0.0000
0.0000
0.0000
-0.0000
-0.0000
0.0000
0.0000
0.0000

informOut =

55
132
17
3
3
79
0
3
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0

ysaveOut =

0.8414   -0.0041    0.0001   -0.0000    0.0000   -0.0000
0.0000   -0.0000    0.0000   -0.0000    0.0000   -0.0000
0.1586    0.0041   -0.0001    0.0000   -0.0000    0.0000

wkjacOut =

1.0100
0
-0.0100
-380.6996
-247.4940
629.1936
-0.0408
-0.0040
2.5831
1.0000
3.0000
3.0000

ifail =

0

```