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

## Purpose

nag_ode_ivp_adams_interp (d02qz) interpolates components of the solution of a non-stiff system of first-order differential equations from information provided by the integrator functions nag_ode_ivp_adams_roots (d02qf) or nag_ode_ivp_adams_roots_revcom (d02qg).

## Syntax

[ywant, ypwant, ifail] = d02qz(neqf, twant, nwant, rwork, iwork)
[ywant, ypwant, ifail] = nag_ode_ivp_adams_interp(neqf, twant, nwant, rwork, iwork)
Note: the interface to this routine has changed since earlier releases of the toolbox:
 At Mark 22: lrwork and liwork were removed from the interface

## Description

nag_ode_ivp_adams_interp (d02qz) evaluates the first nwant components of the solution of a non-stiff system of first-order ordinary differential equations at any point using the method of Watts and Shampine (1986) and information generated by nag_ode_ivp_adams_roots (d02qf) or nag_ode_ivp_adams_roots_revcom (d02qg). nag_ode_ivp_adams_interp (d02qz) should not normally be used to extrapolate outside the current range of the values produced by the integration function.

## References

Watts H A and Shampine L F (1986) Smoother interpolants for Adams codes SIAM J. Sci. Statist. Comput. 7 334–345

## Parameters

### Compulsory Input Parameters

1:     $\mathrm{neqf}$int64int32nag_int scalar
The number of first-order ordinary differential equations being solved by the integration function. It must contain the same value as the argument neqf in a prior call to the setup function nag_ode_ivp_adams_setup (d02qw).
2:     $\mathrm{twant}$ – double scalar
The point at which components of the solution and derivative are to be evaluated. twant should not normally be an extrapolation point, that is twant should satisfy
• $\mathit{told}\le {\mathbf{twant}}\le \mathrm{T}$,
or if integration is proceeding in the negative direction
• $\mathit{told}\ge {\mathbf{twant}}\ge \mathrm{T}$,
where $\mathit{told}$ is the previous integration point and is, to within rounding, tcurrhlast (see nag_ode_ivp_adams_diag (d02qx)). Extrapolation is permitted but not recommended and ${\mathbf{ifail}}={\mathbf{2}}$ is returned whenever extrapolation is attempted.
3:     $\mathrm{nwant}$int64int32nag_int scalar
The number of components of the solution and derivative whose values at twant are required. The first nwant components are evaluated.
Constraint: $1\le {\mathbf{nwant}}\le {\mathbf{neqf}}$.
4:     $\mathrm{rwork}\left(\mathit{lrwork}\right)$ – double array
This must be the same argument rwork as supplied to nag_ode_ivp_adams_setup (d02qw) and to nag_ode_ivp_adams_roots (d02qf) or nag_ode_ivp_adams_roots_revcom (d02qg). It is used to pass information from these functions to nag_ode_ivp_adams_interp (d02qz). Therefore its contents must not be changed before a call to nag_ode_ivp_adams_interp (d02qz).
5:     $\mathrm{iwork}\left(\mathit{liwork}\right)$int64int32nag_int array
This must be the same argument iwork as supplied to nag_ode_ivp_adams_setup (d02qw) and to nag_ode_ivp_adams_roots (d02qf) or nag_ode_ivp_adams_roots_revcom (d02qg). It is used to pass information from these functions to nag_ode_ivp_adams_interp (d02qz). Therefore its contents must not be changed before a call to nag_ode_ivp_adams_interp (d02qz).

None.

### Output Parameters

1:     $\mathrm{ywant}\left({\mathbf{nwant}}\right)$ – double array
The calculated value of the $\mathit{i}$th component of the solution at twant, for $\mathit{i}=1,2,\dots ,{\mathbf{nwant}}$.
2:     $\mathrm{ypwant}\left({\mathbf{nwant}}\right)$ – double array
The calculated value of the $\mathit{i}$th component of the derivative at twant, for $\mathit{i}=1,2,\dots ,{\mathbf{nwant}}$.
3:     $\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$
An integration function (nag_ode_ivp_adams_roots (d02qf) or nag_ode_ivp_adams_roots_revcom (d02qg)) has not been called, no integration steps have been taken since the last call to nag_ode_ivp_adams_setup (d02qw) with ${\mathbf{statef}}=\text{'S'}$, one or more of the arguments lrwork, liwork and neqf does not match the same argument supplied to nag_ode_ivp_adams_setup (d02qw), or nwant does not satisfy $1\le {\mathbf{nwant}}\le {\mathbf{neqf}}$.
W  ${\mathbf{ifail}}=2$
nag_ode_ivp_adams_interp (d02qz) has been called for extrapolation. The values of the solution and its derivative at twant have been calculated and placed in ywant and ypwant before returning with this warning (see Accuracy).
${\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.
These error exits may be caused by overwriting elements of rwork and iwork.

## Accuracy

The error in interpolation is of a similar order to the error arising from the integration. The same order of accuracy can be expected when extrapolating using nag_ode_ivp_adams_interp (d02qz). However, the actual error in extrapolation will, in general, be much larger than for interpolation.

When interpolation for only a few components is required then it is more efficient to order the components of interest so that they are numbered first.

## Example

This example solves the equation
 $y′′=-y, y0=0, y′0=1$
reposed as
 $y1′=y2 y2′=-y1$
over the range $\left[0,\pi /2\right]$ with initial conditions ${y}_{1}=0$ and ${y}_{2}=1$ using vector error control (${\mathbf{vectol}}=\mathit{true}$) and nag_ode_ivp_adams_roots (d02qf) in one-step mode (${\mathbf{onestp}}=\mathit{true}$). nag_ode_ivp_adams_interp (d02qz) is used to provide solution values at intervals of $\pi /16$.
```function d02qz_example

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

% Initialize variables and arrays.
neqf = int64(2);
neqg = int64(0);
latol = neqf;
lrtol = neqf;
lrwork = 23 + 23*neqf + 14*neqg;
liwork = 21 + 4*neqg;
rwork = zeros(lrwork);
iwork = zeros(liwork,'int64');
tstart = 0;
hmax = 2.0;
statef = 's';
vectol = true;
alterg = false;
sophst = false;
atol = [1.0e-08; 1.0e-08];
rtol = [1.0e-04; 1.0e-04];
onestp = true;
crit = true;
tinc = 0.0625*pi;
tcrit = 8.0*tinc;
tout = tcrit;
maxstp = int64(500);
t = tstart;
twant = tstart + tinc;
nwant = neqf;
y = [0 1];

% d02qw is a setup routine to be called prior to d02qf.
[statef, alterg, rwork, iwork, ifail] = ...
d02qw( ...
statef, neqf, vectol, atol, rtol, onestp, crit, tcrit, hmax, ...
maxstp, neqg, alterg, sophst, rwork, iwork);

% Prepare to store results for plotting.
xarray = zeros(1);
yarray = zeros(1, neqf+1);

% Output initial results.
fprintf('  t         y_1      y_2\n');
fprintf('%6.4f   %7.4f  %7.4f\n', t, y(1:2));
xarray(1) = t;
yarray(1, 1:neqf) = y(1:neqf);

j = 1;
while t < tout

% Integrate ODEs from t to tout.
[t, y, root, rwork, iwork, ifail] = ...
d02qf( ...
@f, t, y, tout, 'd02qfz', neqg, rwork, iwork);

% Now interpolate the solution at twant; keep looping till
% twant is bigger than t.
while twant <= t
[ywant, ypwant, ifail] = ...
d02qz( ...
neqf, twant, nwant, rwork, iwork);
fprintf('%6.4f   %7.4f  %7.4f\n', twant, ywant(1), ywant(2));
j = j+1;
xarray(j) = twant;
yarray(j, neqf+1) = abs(ywant(1) - sin(twant));
yarray(j, 1:neqf) = ywant(1:neqf);
twant = tstart + j*tinc;
end
end

% Plot results.
fig1 = figure;
display_plot(xarray, yarray);

function f = f(neqf, x, y)
% Evaluate derivative vector.

f = zeros(neqf, 1);
f(1) =  y(2);
f(2) = -y(1);

function display_plot(x, y)
% Plot two of the curves, then add the other one.
[haxes, hline1, hline2] = plotyy(x, y(:,1), x, y(:,3), 'plot', 'semilogy');
% We want the third curve to be plotted on the left-hand y-axis.
hold(haxes(1), 'on');
hline3 = plot(x, y(:,2));
% Set the axis limits and the tick specifications to beautify the plot.
set(haxes(1), 'YLim', [0 1.1]);
set(haxes(1), 'YMinorTick', 'on');
set(haxes(1), 'YTick', [0.0 0.25 0.5 0.75 1.0]);
set(haxes(2), 'YMinorTick', 'on');
set(haxes(2), 'YTick', [1.0e-10 1.0e-8 1.0e-6 1.0e-4 1.0e-2 1.0]);
for iaxis = 1:2
% These properties must be the same for both sets of axes.
set(haxes(iaxis), 'XLim', [0 1.6]);
end
title('A Simple Problem with Sine Solution');
% Label the axes.
xlabel('x');
ylabel(haxes(1), 'Solution');
yh2 = ylabel(haxes(2), 'Error');
set(yh2,'position',[1.52,1e-7]);
legend('y','y''','Error in y','Location','West')
% Set some features of the three lines.
set(hline1, 'Linewidth', 0.25, 'Marker', '+', 'LineStyle', '-');
set(hline2, 'Linewidth', 0.25, 'Marker', 'x', 'LineStyle', '--');
set(hline3, 'Linewidth', 0.25, 'Marker', '*', 'LineStyle', ':', ...
'Color', 'Magenta');
```
```d02qz example results

t         y_1      y_2
0.0000    0.0000   1.0000
0.1963    0.1951   0.9808
0.3927    0.3827   0.9239
0.5890    0.5556   0.8315
0.7854    0.7071   0.7071
0.9817    0.8315   0.5556
1.1781    0.9239   0.3827
1.3744    0.9808   0.1951
1.5708    1.0000  -0.0000
```