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:     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:     m – int64int32nag_int scalar
The number of components of the solution whose values are required.
Constraint: 1mneq$1\le {\mathbf{m}}\le {\mathbf{neq}}$.
3:     neq – int64int32nag_int scalar
The value used for the parameter neq when calling the integrator.
Constraint: 1neqldysav$1\le {\mathbf{neq}}\le \mathit{ldysav}$.
4:     ysav(ldysav,sdysav) – double array
ldysav, the first dimension of the array, must satisfy the constraint ldysav1$\mathit{ldysav}\ge 1$.
The values provided in the array ysav on return from the integrator.
5:     rwork(50 + 4 × neq$50+4×{\mathbf{neq}}$) – double array
The values provided in the array rwork on return from the integrator.

### Optional Input Parameters

1:     sdysav – int64int32nag_int scalar
Default: The second dimension of the array ysav.
The value used for the parameter sdysav when calling the integrator.

ldysav

### Output Parameters

1:     sol(m) – double array
The calculated value of the solution at tsol.
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:

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

ifail = 1${\mathbf{ifail}}=1$
 On entry, m < 1${\mathbf{m}}<1$, or ldysav < 1$\mathit{ldysav}<1$, or neq < 1${\mathbf{neq}}<1$, or m > neq${\mathbf{m}}>{\mathbf{neq}}$, or neq > ldysav${\mathbf{neq}}>\mathit{ldysav}$.
ifail = 2${\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 ifail = 3${\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.

## Accuracy

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

None.

## Example

```function nag_ode_ivp_stiff_interp_example
%
% Integrate to tout by overshooting tout in one step mode (itask=2)
% using B.D.F formulae with a functional iteration method.
% Default values for the array const are used. Employ vector
% tolerances and the Jacobian is evaluated internally, if necessary.
% monitr subroutine replaced by NAG dummy routine nag_ode_ivp_stiff_exp_fulljac_dummy_monit.
% Interpolation outside nag_ode_ivp_stiff_imp_fulljac using nag_ode_ivp_stiff_interp.
%
neq = 3;
atol = zeros(1, neq);
rtol = zeros(1, neq);
sol = zeros(1, neq);
nwkjac = int64(neq * (neq + 1));
neq = int64(neq);
wkjac = zeros(1, nwkjac);
y = zeros(1, neq);
ydot = zeros(1, neq);
algequ = zeros(1, neq);
lderiv = zeros(1, 2);
inform = zeros(23, 1, 'int64');
maxord = int64(5);
sdysav = int64(6);
ysav = zeros(neq, sdysav);
petzld = false;
const = zeros(6,1);
h0 = 0.0e0;
hmax = 10.0e0;
hmin = 1.0e-10;
tcrit = 0.0e0;
maxstp = int64(200);
mxhnil = int64(5);
rwork = zeros(62,1);
t = 0.0e0;
tout = 0.1e0;
y = [1; 0; 0];
lderiv = [false; false];
itol = int64(4);
rtol = [1.0e-4; 1.0e-3; 1.0e-4];
atol = [1.0e-7; 1.0e-8; 1.0e-7];
sol = zeros(6, neq);
sol(1,:) = y;
fprintf('\nnag_ode_ivp_stiff_interp example program results\n');
[const, rwork, ifail] = ...
nag_ode_ivp_stiff_bdf(neq, sdysav, maxord, 'functional-iteration', ...
petzld, const, tcrit, hmin, hmax, h0, maxstp, ...
mxhnil, 'average-l2', rwork);
if ifail == int64(0)
% Linear algebra setup required (in case functional iteration
% encounters any difficulty).
[rwork, ifail] = ...
nag_ode_ivp_stiff_fulljac_setup(neq, neq, 'numerical', nwkjac, rwork);
xout = 0.02e0;
iout = 1;
fprintf('\n');
fprintf('     x           y(1)          y(2)          y(3)\n');
fprintf('%8.3f %13.5f %13.5f %13.5f\n', t, y(1:double(neq)));
itrace = int64(0);
tcur = 1;
hu = 0;
while (tcur-hu >= xout) || (xout > tcur) && (ifail == 0)
[t, tout, y, ydot, rwork, inform, ysav, wkjac, lderiv, ifail] = ...
nag_ode_ivp_stiff_imp_fulljac(t, tout, y, ydot, rwork, rtol, atol, itol, inform, ...
@resid, ysav, 'nag_ode_ivp_stiff_imp_fulljac_dummy_jac', wkjac, 'nag_ode_ivp_stiff_exp_fulljac_dummy_monit', lderiv, itask, ...
itrace);
[hu, h, tcur, tolsf, nst, nre, nje, nqu, nq, niter, imxer, algequ, ifail] = ...
nag_ode_ivp_stiff_integ_diag(neq, neq, rwork, inform);
while (iout < 6) && (tcur-hu < xout) && (xout <= tcur) && (ifail == 0)
[sol(iout+1,:), ifail] = nag_ode_ivp_stiff_interp(xout, neq, neq, ysav, rwork);
fprintf('%8.3f %13.5f %13.5f %13.5f\n', xout, sol(iout+1,1:double(neq)));
iout = iout + 1;
xout = iout * 0.02e0;
end
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
```
```

nag_ode_ivp_stiff_interp example program 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

```
```function d02mz_example
%
% Integrate to tout by overshooting tout in one step mode (itask=2)
% using B.D.F formulae with a functional iteration method.
% Default values for the array const are used. Employ vector
% tolerances and the Jacobian is evaluated internally, if necessary.
% monitr subroutine replaced by NAG dummy routine d02nby.
% Interpolation outside d02ng using d02mz.
%
neq = 3;
atol = zeros(1, neq);
rtol = zeros(1, neq);
sol = zeros(1, neq);
nwkjac = int64(neq * (neq + 1));
neq = int64(neq);
wkjac = zeros(1, nwkjac);
y = zeros(1, neq);
ydot = zeros(1, neq);
algequ = zeros(1, neq);
lderiv = zeros(1, 2);
inform = zeros(23, 1, 'int64');
maxord = int64(5);
sdysav = int64(6);
ysav = zeros(neq, sdysav);
petzld = false;
const = zeros(6,1);
h0 = 0.0e0;
hmax = 10.0e0;
hmin = 1.0e-10;
tcrit = 0.0e0;
maxstp = int64(200);
mxhnil = int64(5);
rwork = zeros(62,1);
t = 0.0e0;
tout = 0.1e0;
y = [1; 0; 0];
lderiv = [false; false];
itol = int64(4);
rtol = [1.0e-4; 1.0e-3; 1.0e-4];
atol = [1.0e-7; 1.0e-8; 1.0e-7];
sol = zeros(6, neq);
sol(1,:) = y;
fprintf('\nd02mz example program results\n');
[const, rwork, ifail] = d02nv(neq, sdysav, maxord, 'functional-iteration', ...
petzld, const, tcrit, hmin, hmax, h0, maxstp, ...
mxhnil, 'average-l2', rwork);
if ifail == int64(0)
% Linear algebra setup required (in case functional iteration
% encounters any difficulty).
[rwork, ifail] = d02ns(neq, neq, 'numerical', nwkjac, rwork);
xout = 0.02e0;
iout = 1;
fprintf('\n');
fprintf('     x           y(1)          y(2)          y(3)\n');
fprintf('%8.3f %13.5f %13.5f %13.5f\n', t, y(1:double(neq)));
itrace = int64(0);
tcur = 1;
hu = 0;
while (tcur-hu >= xout) || (xout > tcur) && (ifail == 0)
[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);
[hu, h, tcur, tolsf, nst, nre, nje, nqu, nq, niter, imxer, algequ, ifail] = ...
d02ny(neq, neq, rwork, inform);
while (iout < 6) && (tcur-hu < xout) && (xout <= tcur) && (ifail == 0)
[sol(iout+1,:), ifail] = d02mz(xout, neq, neq, ysav, rwork);
fprintf('%8.3f %13.5f %13.5f %13.5f\n', xout, sol(iout+1,1:double(neq)));
iout = iout + 1;
xout = iout * 0.02e0;
end
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 program 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

```