NAG Toolbox: nag_ode_ivp_stiff_fulljac_setup (d02ns)

Purpose

Syntax

Description

References

Parameters

Compulsory Input Parameters

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

The number of differential equations.

Constraint: $1\le \mathbf{neq}\le \mathbf{neqmax}$.

2:     $\mathrm{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:     $\mathrm{jceval}$ – string (length ≥ 1)

Specifies the technique to be used to compute the Jacobian.

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

$\mathbf{jceval}=\text{'A'}$

You must supply a (sub)program to evaluate the Jacobian on a call to the integrator.

$\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 argument 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: $\mathbf{jceval}=\text{'N'}$, $\text{'A'}$ or $\text{'D'}$.

4:     $\mathrm{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: $\mathbf{nwkjac}\ge \mathbf{neqmax}\times \left(\mathbf{neqmax}+1\right)$.

5:     $\mathrm{rwork}\left(50+4\times \mathbf{neqmax}\right)$ – 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.

Optional Input Parameters

None.

Output Parameters

1:     $\mathrm{rwork}\left(50+4\times \mathbf{neqmax}\right)$ – double 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

   $\mathbf{ifail}=1$

On entry,	$\mathbf{neq}<1$,
or	$\mathbf{neq}>\mathbf{neqmax}$,
or	$\mathbf{nwkjac}<\mathbf{neqmax}\times \left(\mathbf{neqmax}+1\right)$,
or	$\mathbf{jceval}\ne \text{'N'}$, $\text{'A'}$ or $\text{'D'}$.

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

function d02ns_example


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

% Initialize variables and arrays for setup routine

n      = int64(3);
ord    = int64(5);
sdy    = int64(ord + 1);
petzld = false;
con    = zeros(6);
tcrit  = 10;
hmin   = 1.0e-10;
hmax   = 10;
h0     = 0;
maxstp = int64(200);
mxhnil = int64(5);
lrwork = 50 + 4*n;
rwork  = zeros(lrwork);

% Setup integration method using d02nv.
[const, rwork, ifail] = d02nv(n, sdy, ord, 'Newton', petzld, con, tcrit, ...
                              hmin, hmax, h0, maxstp, mxhnil, 'Average-L2', ...
                              rwork);
% Setup for analytical Jacbian using d02ns
nwkjac = int64(n*(n+1));
[rwork, ifail] = d02ns(n, n, 'Analytical', nwkjac, rwork);

% Initialize variables and arrays for integration
t      = 0;
tout   = 10;
y      = [1; 0; 0];
rtol   = [0.0001];
atol   = [1e-07];
itol   = int64(1);
inform = zeros(23, 1, 'int64');
ysave  = zeros(n, sdy);
wkjac  = zeros(nwkjac, 1);
itask  = int64(4);
itrace = int64(0);

% Integrate ODE from t=0 to t=tout, no monitoring, using d02nb
[t, y, ydot, rwork, inform, ysave, wkjac, ifail] = d02nb(t, tout, y, rwork, ...
          rtol, atol, itol, inform, @fcn, ysave, @jac, wkjac, 'd02nby', ...
          itask, itrace);

fprintf('Solution y and derivative y'' at t = %7.4f is:\n',t);
fprintf('\n %10s %10s\n','y','y''');
for i=1:n
  fprintf(' %10.4f %10.4f\n',y(i),ydot(i));
end



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) = -hxd*(-0.04d0    ) + 1;
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) = -hxd*(-1.0d4*y(3)-6.0d7*y(2)) + 1;
p(2,3) = -hxd*(-1.0d4*y(2));
p(3,2) = -hxd*( 6.0d7*y(2));
p(3,3) = -hxd*( 0.0d0     ) + 1;

d02ns example results

Solution y and derivative y' at t = 10.0000 is:

          y         y'
     0.8414    -0.0079
     0.0000    -0.0000
     0.1586     0.0079