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:
Mark 22: lrwork, liwork have been 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:     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 parameter neqf in a prior call to the setup function nag_ode_ivp_adams_setup (d02qw).
2:     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
• toldtwantT$\mathit{told}\le {\mathbf{twant}}\le \mathrm{T}$,
or if integration is proceeding in the negative direction
• toldtwantT$\mathit{told}\ge {\mathbf{twant}}\ge \mathrm{T}$,
where told$\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:     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: 1nwantneqf$1\le {\mathbf{nwant}}\le {\mathbf{neqf}}$.
4:     rwork(lrwork) – double array
This must be the same parameter 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:     iwork(liwork) – int64int32nag_int array
This must be the same parameter 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.

lrwork liwork

Output Parameters

1:     ywant(nwant) – double array
The calculated value of the i$\mathit{i}$th component of the solution at twant, for i = 1,2,,nwant$\mathit{i}=1,2,\dots ,{\mathbf{nwant}}$.
2:     ypwant(nwant) – double array
The calculated value of the i$\mathit{i}$th component of the derivative at twant, for i = 1,2,,nwant$\mathit{i}=1,2,\dots ,{\mathbf{nwant}}$.
3:     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$
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 statef = 'S'${\mathbf{statef}}=\text{'S'}$, one or more of the parameters lrwork, liwork and neqf does not match the same parameter supplied to nag_ode_ivp_adams_setup (d02qw), or nwant does not satisfy 1nwantneqf$1\le {\mathbf{nwant}}\le {\mathbf{neqf}}$.
W ifail = 2${\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 Section [Accuracy]).
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,  y(0) = 0,  y′(0) = 1 $y′′=-y, y(0)=0, y′(0)=1$
reposed as
 y1 ′ = y2 y2 ′ = − y1
$y1′=y2 y2′=-y1$
over the range [0,π / 2] $\left[0,\pi /2\right]$ with initial conditions y1 = 0${y}_{1}=0$ and y2 = 1${y}_{2}=1$ using vector error control (vectol = true${\mathbf{vectol}}=\mathbf{true}$) and nag_ode_ivp_adams_roots (d02qf) in one-step mode (onestp = true${\mathbf{onestp}}=\mathbf{true}$). nag_ode_ivp_adams_interp (d02qz) is used to provide solution values at intervals of π / 16$\pi /16$.
```function nag_ode_ivp_adams_interp_example
% Initialize variables and arrays.
neqf = 2;
neqg = 0;
latol = neqf;
lrtol = neqf;
lrwork = 23 + 23*neqf + 14*neqg;
liwork = 21 + 4*neqg;
rwork = zeros(lrwork);
iwork = zeros(liwork);
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 = 500;
t = tstart;
twant = tstart + tinc;
nwant = neqf;
y = [0.0; 1.0];

% nag_ode_ivp_adams_setup is a setup routine to be called prior to nag_ode_ivp_adams_roots.
[statef, alterg, rwork, iwork, ifail] = nag_ode_ivp_adams_setup(statef, int64(neqf), ...
vectol, atol, rtol, onestp, crit, tcrit, hmax, int64(maxstp), ...
int64(neqg), alterg, sophst, rwork, int64(iwork));
if ifail ~= 0
% Unsuccessful call.  Print message and exit.
error('Warning: nag_ode_ivp_adams_setup returned with ifail = %1d ',ifail);
end

% 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), y(2));
xarray(1) = t;
for ieqf = 1:neqf
yarray(1, ieqf) = y(ieqf);
end

j = 1;
while t < tout

% Integrate ODEs from t to tout.
[t, y, root, rwork, iwork, ifail] = nag_ode_ivp_adams_roots(@f, t, y, tout, ...
if ifail ~= 0
% Unsuccessful call.  Print message and exit.
error('Warning: nag_ode_ivp_adams_roots returned with ifail = %1d ',ifail);
end

% Now interpolate the solution at twant; keep looping till
% twant is bigger than t.
while twant <= t
[ywant, ypwant, ifail] = nag_ode_ivp_adams_interp(int64(neqf), twant, ...
int64(nwant), rwork, iwork);
if ifail ~= 0
% Unsuccessful call.  Print message and exit.
error('Warning: nag_ode_ivp_adams_interp returned with ifail = %1d ',ifail);
end
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));
for ieqf = 1:neqf
yarray(j, ieqf) = ywant(ieqf);
end
twant = tstart + j*tinc;
end
end
% Plot results.
fig = figure('Number', 'off');
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)
% Formatting for title and axis labels.
titleFmt = {'FontName', 'Helvetica', 'FontWeight', 'Bold', 'FontSize', 14};
labFmt = {'FontName', 'Helvetica', 'FontWeight', 'Bold', 'FontSize', 13};
set(gca, 'FontSize', 13); % for legend, axis tick labels, etc.
% 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]);
set(haxes(iaxis), 'FontSize', 13);
end
title('A Simple Problem with Sine Solution', titleFmt{:});
% Label the axes.
xlabel('x', labFmt{:});
ylabel(haxes(1), 'Solution', labFmt{:});
ylabel(haxes(2), 'Error', labFmt{:});
legend('y','y''','Error in y','Location','Best')
% Set some features of the three lines.
set(hline1, 'Linewidth', 0.25, 'Marker', '+', 'Line', '-');
set(hline2, 'Linewidth', 0.25, 'Marker', 'x', 'Line', '--');
set(hline3, 'Linewidth', 0.25, 'Marker', '*', 'Line', ':');
```
```

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

```
```function d02qz_example
% Initialize variables and arrays.
neqf = 2;
neqg = 0;
latol = neqf;
lrtol = neqf;
lrwork = 23 + 23*neqf + 14*neqg;
liwork = 21 + 4*neqg;
rwork = zeros(lrwork);
iwork = zeros(liwork);
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 = 500;
t = tstart;
twant = tstart + tinc;
nwant = neqf;
y = [0.0; 1.0];

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

% d02qw is a setup routine to be called prior to d02qf.
[statef, alterg, rwork, iwork, ifail] = d02qw(statef, int64(neqf), ...
vectol, atol, rtol, onestp, crit, tcrit, hmax, int64(maxstp), ...
int64(neqg), alterg, sophst, rwork, int64(iwork));
if ifail ~= 0
% Unsuccessful call.  Print message and exit.
error('Warning: d02qw returned with ifail = %1d ',ifail);
end

% 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), y(2));
xarray(1) = t;
for ieqf = 1:neqf
yarray(1, ieqf) = y(ieqf);
end

j = 1;
while t < tout

% Integrate ODEs from t to tout.
[t, y, root, rwork, iwork, ifail] = d02qf(@f, t, y, tout, ...
'd02qfz', int64(neqg), rwork, iwork);
if ifail ~= 0
% Unsuccessful call.  Print message and exit.
error('Warning: d02qf returned with ifail = %1d ',ifail);
end

% Now interpolate the solution at twant; keep looping till
% twant is bigger than t.
while twant <= t
[ywant, ypwant, ifail] = d02qz(int64(neqf), twant, ...
int64(nwant), rwork, iwork);
if ifail ~= 0
% Unsuccessful call.  Print message and exit.
error('Warning: d02qz returned with ifail = %1d ',ifail);
end
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));
for ieqf = 1:neqf
yarray(j, ieqf) = ywant(ieqf);
end
twant = tstart + j*tinc;
end
end
% Plot results.
fig = figure('Number', 'off');
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)
% Formatting for title and axis labels.
titleFmt = {'FontName', 'Helvetica', 'FontWeight', 'Bold', 'FontSize', 14};
labFmt = {'FontName', 'Helvetica', 'FontWeight', 'Bold', 'FontSize', 13};
set(gca, 'FontSize', 13); % for legend, axis tick labels, etc.
% 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]);
set(haxes(iaxis), 'FontSize', 13);
end
title('A Simple Problem with Sine Solution', titleFmt{:});
% Label the axes.
xlabel('x', labFmt{:});
ylabel(haxes(1), 'Solution', labFmt{:});
ylabel(haxes(2), 'Error', labFmt{:});
legend('y','y''','Error in y','Location','Best')
% Set some features of the three lines.
set(hline1, 'Linewidth', 0.25, 'Marker', '+', 'Line', '-');
set(hline2, 'Linewidth', 0.25, 'Marker', 'x', 'Line', '--');
set(hline3, 'Linewidth', 0.25, 'Marker', '*', 'Line', ':');
```
```
d02qz example program 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

```