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_interp (d02mz)

## Purpose

nag_ode_ivp_stiff_interp (d02mz) interpolates components of the solution of a system of first-order differential equations from information provided by those integrators in Sub-chapter D02M–N using methods set up by calls to nag_ode_ivp_stiff_dassl (d02mv), nag_ode_ivp_stiff_bdf (d02nv) or nag_ode_ivp_stiff_blend (d02nw).

## Syntax

[sol, ifail] = d02mz(tsol, m, neq, ysav, rwork, 'sdysav', sdysav)
[sol, ifail] = nag_ode_ivp_stiff_interp(tsol, m, neq, ysav, rwork, 'sdysav', sdysav)

## Description

nag_ode_ivp_stiff_interp (d02mz) evaluates the first m components of the solution of a system of ordinary differential equations at any point using natural polynomial interpolation based on information generated by the integrator. This information must be passed unchanged to nag_ode_ivp_stiff_interp (d02mz). nag_ode_ivp_stiff_interp (d02mz) should not normally be used to extrapolate outside the range of values obtained from the above function.

## References

See the D02M–N Sub-chapter Introduction.

## Parameters

### Compulsory Input Parameters

1:     $\mathrm{tsol}$ – double scalar
The point at which the first m components of the solution are to be evaluated. tsol should not normally be an extrapolation point. Extrapolation is permitted but not recommended.
2:     $\mathrm{m}$int64int32nag_int scalar
The number of components of the solution whose values are required.
Constraint: $1\le {\mathbf{m}}\le {\mathbf{neq}}$.
3:     $\mathrm{neq}$int64int32nag_int scalar
The value used for the argument neq when calling the integrator.
Constraint: $1\le {\mathbf{neq}}\le \mathit{ldysav}$.
4:     $\mathrm{ysav}\left(\mathit{ldysav},{\mathbf{sdysav}}\right)$ – double array
ldysav, the first dimension of the array, must satisfy the constraint $\mathit{ldysav}\ge 1$.
The values provided in the array ysav on return from the integrator.
5:     $\mathrm{rwork}\left(50+4×{\mathbf{neq}}\right)$ – double array
The values provided in the array rwork on return from the integrator.

### Optional Input Parameters

1:     $\mathrm{sdysav}$int64int32nag_int scalar
Default: the second dimension of the array ysav.
The value used for the argument sdysav when calling the integrator.

### Output Parameters

1:     $\mathrm{sol}\left({\mathbf{m}}\right)$ – double array
The calculated value of the solution at tsol.
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:

Cases prefixed with W are classified as warnings and do not generate an error of type NAG:error_n. See nag_issue_warnings.

${\mathbf{ifail}}=1$
 On entry, ${\mathbf{m}}<1$, or $\mathit{ldysav}<1$, or ${\mathbf{neq}}<1$, or ${\mathbf{m}}>{\mathbf{neq}}$, or ${\mathbf{neq}}>\mathit{ldysav}$.
${\mathbf{ifail}}=2$
On entry, when accessing an element of the array rwork an unexpected quantity was found. You have not passed the correct array to nag_ode_ivp_stiff_interp (d02mz) or has overwritten elements of this array.
W  ${\mathbf{ifail}}=3$
On entry, nag_ode_ivp_stiff_interp (d02mz) has been called for extrapolation. Before returning with this error exit, the value of the solution at tsol is calculated and placed in sol.
${\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.

## Accuracy

The solution values returned will be of a similar accuracy to those computed by the integrator.

None.

## Example

This example solves the well-known stiff Robertson problem written in implicit form
 $r1 = -0.04a + 1.0E4bc - a′ r2 = 0.04a - 1.0E4bc - 3.0E7⁢b2 - b′ r3 = 3.0E7⁢b2 - c′$
with initial conditions $a=1.0$ and $b=c=0.0$ over the range $\left[0,0.1\right]$ with vector error control (${\mathbf{itol}}=4$), the BDF method (setup function nag_ode_ivp_stiff_bdf (d02nv)) and functional iteration. The Jacobian is calculated numerically if the functional iteration encounters difficulty and the integration is in one-step mode (${\mathbf{itask}}=2$), with natural interpolation to calculate the solution at intervals of $0.02$ using nag_ode_ivp_stiff_interp (d02mz) externally. nag_ode_ivp_stiff_exp_fulljac_dummy_monit (d02nby) is used for monitr.
```function d02mz_example

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

% Integrate to using B.D.F formulae with a functional iteration method
neq    = int64(3);
maxord = int64(5);
sdysav = int64(maxord+1);
petzld = false;
const  = zeros(6,1);
tcrit  = 0;
hmin   = 1e-10;
hmax =   10;
h0     = 0;
maxstp = int64(200);
mxhnil = int64(5);
rwork  = zeros(16+20*neq,1);
[const, rwork, ifail] = d02nv(...
neq, sdysav, maxord, 'functional-iteration', ...
petzld, const, tcrit, hmin, hmax, h0, maxstp, ...
mxhnil, 'average-l2', rwork);

% Use numerical Jacobian if functional iteration encounters any difficulty.
nwkjac = int64(neq * (neq + 1));
[rwork, ifail] = d02ns(neq, neq, 'numerical', nwkjac, rwork);

algequ = zeros(1, neq);
lderiv = zeros(1, 2);

% Variable and array initializations for integrator
% Initial conditions
t      = 0;
tout   = 0.1;
y      = [1; 0; 0];
ydot   = zeros(1, neq);
% vector tolerances
itol   = int64(4);
rtol   = [1.0e-4; 1.0e-3; 1.0e-4];
atol   = [1.0e-7; 1.0e-8; 1.0e-7];
inform = zeros(23, 1, 'int64');
ysav   = zeros(neq, sdysav);
wkjac  = zeros(1, nwkjac);
lderiv = [false; false];
itrace = int64(0);

% Intialize for solution output
xout = 0.02e0;
fprintf('\n     x           y(1)          y(2)          y(3)\n');
fprintf('%8.3f %13.5f %13.5f %13.5f\n', t, y);

while (t < tout)
[t, tout, y, ydot, rwork, inform, ysav, wkjac, lderiv, ifail] = ...
d02ng(...
t, tout, y, ydot, rwork, rtol, atol, itol, inform, ...
@resid, ysav, 'd02ngz', wkjac, 'd02nby', lderiv, itask, ...
itrace);
while (xout <= t)
[sol, ifail] = d02mz(xout, neq, neq, ysav, rwork);
fprintf('%8.3f %13.5f %13.5f %13.5f\n', xout, sol);
xout = xout + 0.02e0;
end
end

function [r, ires] = resid(neq, t, y, ydot, ires)
r = -ydot;
if ires == int64(1)
r(1) = -0.04e0*y(1) + 1.0e4*y(2)*y(3)                   + r(1);
r(2) =  0.04e0*y(1) - 1.0e4*y(2)*y(3) - 3.0e7*y(2)*y(2) + r(2);
r(3) =                                  3.0e7*y(2)*y(2) + r(3);
end
```
```d02mz example results

x           y(1)          y(2)          y(3)
0.000       1.00000       0.00000       0.00000
0.020       0.99920       0.00004       0.00076
0.040       0.99841       0.00004       0.00155
0.060       0.99763       0.00004       0.00234
0.080       0.99685       0.00004       0.00311
0.100       0.99608       0.00004       0.00389
```

Chapter Contents
Chapter Introduction
NAG Toolbox

© The Numerical Algorithms Group Ltd, Oxford, UK. 2009–2015