hide long namesshow long names
hide short namesshow short names
Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int
PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

NAG Toolbox: nag_ode_ivp_adams_interp (d02qz)

▸▿ Contents

    1  Purpose
    2  Syntax
    3  Description
    4  References
▸▿ 5  Parameters
    6  Error Indicators and Warnings
    7  Accuracy
    8  Further Comments
    9  Example

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:     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:     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
  • told≤twant≤T,
or if integration is proceeding in the negative direction
  • told≥twant≥T,
where told is the previous integration point and is, to within rounding, tcurr – hlast (see nag_ode_ivp_adams_diag (d02qx)). Extrapolation is permitted but not recommended and ifail=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: 1≤nwant≤neqf.
4:     rworklrwork – 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:     iworkliwork – 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).

Optional Input Parameters

None.

Output Parameters

1:     ywantnwant – double array
The calculated value of the ith component of the solution at twant, for i=1,2,…,nwant.
2:     ypwantnwant – double array
The calculated value of the ith component of the derivative at twant, for i=1,2,…,nwant.
3:     ifail – int64int32nag_int scalar
ifail=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
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', 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≤nwant≤neqf.
W  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).
   ifail=-99
An unexpected error has been triggered by this routine. Please contact NAG.
   ifail=-399
Your licence key may have expired or may not have been installed correctly.
   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.

Further Comments

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 0,π/2  with initial conditions y1=0 and y2=1 using vector error control (vectol=true) and nag_ode_ivp_adams_roots (d02qf) in one-step mode (onestp=true). nag_ode_ivp_adams_interp (d02qz) is used to provide solution values at intervals of π/16.

Open in the MATLAB editor: d02qz_example

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
  % Add title.
  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]);
  % Add a legend.
  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
d02qz_fig1.png
PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox
© The Numerical Algorithms Group Ltd, Oxford, UK. 2009–2015