PDF version (NAG web site
, 64bit version, 64bit version)
NAG Toolbox: nag_ode_ivp_2nd_rkn_interp (d02lz)
Purpose
nag_ode_ivp_2nd_rkn_interp (d02lz) interpolates components of the solution of a nonstiff system of secondorder differential equations from information provided by the integrator
nag_ode_ivp_2nd_rkn (d02la), when the loworder method has been used.
Syntax
Note: the interface to this routine has changed since earlier releases of the toolbox:
At Mark 22: 
lrwork was removed from the interface 
Description
nag_ode_ivp_2nd_rkn_interp (d02lz) evaluates the first
${\mathbf{nwant}}$ (
$\text{}\le {\mathbf{neq}}$) components of the solution of a nonstiff system of secondorder ordinary differential equations at any point using a special Runge–Kutta–Nystrom formula (see
Dormand and Prince (1986)) and information generated by
nag_ode_ivp_2nd_rkn (d02la) when the loworder method has been used. This information must be presented unchanged to
nag_ode_ivp_2nd_rkn_interp (d02lz).
nag_ode_ivp_2nd_rkn_interp (d02lz) should not normally be used to extrapolate outside the range of the values from
nag_ode_ivp_2nd_rkn (d02la).
References
Dormand J R and Prince P J (1986) Runge–Kutta–Nystrom triples Mathematical Report TPCS8605 Teesside Polytechnic
Parameters
Compulsory Input Parameters
 1:
$\mathrm{t}$ – double scalar

$t$, the current value at which the solution and its derivative have been computed (as returned in argument
t on output from
nag_ode_ivp_2nd_rkn (d02la)).
 2:
$\mathrm{y}\left({\mathbf{neq}}\right)$ – double array

The
$\mathit{i}$th component of the solution at
$t$, for
$\mathit{i}=1,2,\dots ,{\mathbf{neq}}$, as returned from
nag_ode_ivp_2nd_rkn (d02la).
 3:
$\mathrm{yp}\left({\mathbf{neq}}\right)$ – double array

The
$\mathit{i}$th component of the derivative at
$t$, for
$\mathit{i}=1,2,\dots ,{\mathbf{neq}}$, as returned from
nag_ode_ivp_2nd_rkn (d02la).
 4:
$\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{neq}}$.
 5:
$\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
or if integration is proceeding in the negative direction
where
told is the previous integration point which is held in an element of the array
rwork and is, to within rounding,
${\mathbf{t}}{\mathbf{hused}}$. (
hused is given by
nag_ode_ivp_2nd_rkn_diag (d02ly).) Extrapolation is permitted but not recommended, and
${\mathbf{ifail}}={\mathbf{2}}$ is returned whenever extrapolation is attempted.
 6:
$\mathrm{rwork}\left(\mathit{lrwork}\right)$ – double array

This
must be the same argument
rwork as supplied to
nag_ode_ivp_2nd_rkn (d02la). It is used to pass information from
nag_ode_ivp_2nd_rkn (d02la) to
nag_ode_ivp_2nd_rkn_interp (d02lz) and therefore the contents of this array
must not be changed before calling
nag_ode_ivp_2nd_rkn_interp (d02lz).
Optional Input Parameters
 1:
$\mathrm{neq}$ – int64int32nag_int scalar

Default:
the dimension of the arrays
y,
yp. (An error is raised if these dimensions are not equal.)
The number of secondorder ordinary differential equations being solved by
nag_ode_ivp_2nd_rkn (d02la). It must contain the same value as the argument
neq in a prior call to
nag_ode_ivp_2nd_rkn (d02la).
Output Parameters
 1:
$\mathrm{ywant}\left({\mathbf{nwant}}\right)$ – double array

The calculated value of the
$\mathit{i}$th component of the solution at $t={\mathbf{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 $t={\mathbf{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).
If
nag_ode_ivp_2nd_rkn_interp (d02lz) is to be used for extrapolation at
twant,
ifail should be set to
$1$ before entry. It is then essential to test the value of
ifail on exit.
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$

Illegal input detected, i.e., one of the following conditions:
This error exit can be caused if elements of
rwork have been overwritten.
 W ${\mathbf{ifail}}=2$

nag_ode_ivp_2nd_rkn_interp (d02lz) 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 error number (see
Accuracy).
 ${\mathbf{ifail}}=3$

nag_ode_ivp_2nd_rkn (d02la) last used the high order method to integrate the system of differential equations. Interpolation is not permitted with this method.
 ${\mathbf{ifail}}=99$
An unexpected error has been triggered by this routine. Please
contact
NAG.
 ${\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 error in interpolation is of a similar order to the error arising from the integration using
nag_ode_ivp_2nd_rkn (d02la) with the lower order method.
The same order of accuracy can be expected when extrapolating using nag_ode_ivp_2nd_rkn_interp (d02lz). 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
See
Example in
nag_ode_ivp_2nd_rkn (d02la).
Open in the MATLAB editor:
d02lz_example
function d02lz_example
fprintf('d02lz example results\n\n');
neq = 2;
h = 0;
tol = 1e4;
thres = zeros(neq, 1);
thresp = zeros(neq, 1);
maxstp = int64(1000);
start = true;
onestp = true;
high = false;
rwork = zeros(16+20*neq,1);
[startOut, rwork, ifail] = d02lx...
(h, tol, thres, thresp, maxstp, start, onestp, high, rwork);
t = 0;
tend = 20;
y = [0.5; 0];
yp = [0; sqrt(3)];
ydp = zeros(neq, 1);
tnext = 2;
nwant = int64(2);
fprintf('\n T Y(1) Y(2)\n');
fprintf('%4.1f %10.5f %10.5f\n', t, y(1), y(2));
while (t < tend && ifail == 0)
[t, y, yp, ydp, rwork, ifail] = d02la(@fcn, t, tend, y, yp, ydp, rwork);
while (tnext <= t && ifail == 0)
[ywant, ypwant, ifail] = d02lz(t, y, yp, nwant, tnext, rwork);
fprintf('%4.1f %10.5f %10.5f\n', tnext, ywant(1:2));
tnext = tnext + 2;
end
end
if (ifail == 0)
[hnext, hused, hstart, nsucc, nfail, natt, thres, thresp, ifail] = ...
d02ly(nwant, rwork);
fprintf('\n Number of successful steps = %d\n', nsucc);
fprintf(' Number of failed steps = %d\n', nfail);
end
function ydp = fcn(neq, t, y)
r = sqrt(y(1)^2+y(2)^2)^3;
ydp(1) = y(1)/r;
ydp(2) = y(2)/r;
d02lz example results
T Y(1) Y(2)
0.0 0.50000 0.00000
2.0 1.20573 0.61357
4.0 1.33476 0.47685
6.0 0.35748 0.44558
8.0 1.03762 0.73022
10.0 1.42617 0.32658
12.0 0.05515 0.72032
14.0 0.82880 0.81788
16.0 1.48103 0.16788
18.0 0.26719 0.84223
20.0 0.57803 0.86339
Number of successful steps = 108
Number of failed steps = 16
PDF version (NAG web site
, 64bit version, 64bit version)
© The Numerical Algorithms Group Ltd, Oxford, UK. 2009–2015