PDF version (NAG web site
, 64bit version, 64bit version)
NAG Toolbox: nag_ode_ivp_2nd_rkn (d02la)
Purpose
nag_ode_ivp_2nd_rkn (d02la) is a function for integrating a nonstiff system of secondorder ordinary differential equations using Runge–Kutta–Nystrom techniques.
Syntax
[
t,
y,
yp,
ydp,
rwork,
ifail] = d02la(
fcn,
t,
tend,
y,
yp,
ydp,
rwork, 'neq',
neq)
[
t,
y,
yp,
ydp,
rwork,
ifail] = nag_ode_ivp_2nd_rkn(
fcn,
t,
tend,
y,
yp,
ydp,
rwork, 'neq',
neq)
Note: the interface to this routine has changed since earlier releases of the toolbox:
At Mark 22: 
lrwork was removed from the interface 
Description
Given the initial values
$x,{y}_{1},{y}_{2},\dots ,{y}_{{\mathbf{neq}}},{y}_{1}^{\prime},{y}_{2}^{\prime},\dots ,{y}_{{\mathbf{neq}}}^{\prime}$ nag_ode_ivp_2nd_rkn (d02la) integrates a nonstiff system of secondorder differential equations of the type
from
$x={\mathbf{t}}$ to
$x={\mathbf{tend}}$ using a Runge–Kutta–Nystrom formula pair. The system is defined by
fcn, which evaluates
${f}_{i}$ in terms of
$x$ and
${y}_{1},{y}_{2},\dots ,{y}_{{\mathbf{neq}}}$, where
${y}_{1},{y}_{2},\dots ,{y}_{{\mathbf{neq}}}$ are supplied at
$x$.
There are two Runge–Kutta–Nystrom formula pairs implemented in this function. The lower order method is intended if you have moderate accuracy requirements and may be used in conjunction with the interpolation function
nag_ode_ivp_2nd_rkn_interp (d02lz) to produce solutions and derivatives at userspecified points. The higher order method is intended if you have high accuracy requirements.
In onestep mode the function returns approximations to the solution, derivative and
${f}_{i}$ at each integration point. In interval mode these values are returned at the end of the integration range. You select the order of the method, the mode of operation, the error control and various optional inputs by a prior call to
nag_ode_ivp_2nd_rkn_setup (d02lx).
For a description of the Runge–Kutta–Nystrom formula pairs see
Dormand et al. (1986a) and
Dormand et al. (1986b) and for a description of their practical implementation see
Brankin et al. (1989).
References
Brankin R W, Dormand J R, Gladwell I, Prince P J and Seward W L (1989) Algorithm 670: A Runge–Kutta–Nystrom Code ACM Trans. Math. Software 15 31–40
Dormand J R, El–Mikkawy M E A and Prince P J (1986a) Families of Runge–Kutta–Nystrom formulae Mathematical Report TPMR 861 Teesside Polytechnic
Dormand J R, El–Mikkawy M E A and Prince P J (1986b) High order embedded Runge–Kutta–Nystrom formulae Mathematical Report TPMR 862 Teesside Polytechnic
Parameters
Compulsory Input Parameters
 1:
$\mathrm{fcn}$ – function handle or string containing name of mfile

fcn must evaluate the functions
${f}_{i}$ (that is the second derivatives
${y}_{i}^{\prime \prime}$) for given values of its arguments
$x$,
${y}_{1},{y}_{2},\dots ,{y}_{{\mathbf{neq}}}$.
[f] = fcn(neq, t, y)
Input Parameters
 1:
$\mathrm{neq}$ – int64int32nag_int scalar

The number of differential equations.
 2:
$\mathrm{t}$ – double scalar

$x$, the value of the argument.
 3:
$\mathrm{y}\left({\mathbf{neq}}\right)$ – double array

${y}_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,{\mathbf{neq}}$, the value of the argument.
Output Parameters
 1:
$\mathrm{f}\left({\mathbf{neq}}\right)$ – double array

The value of
${f}_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,{\mathbf{neq}}$.
 2:
$\mathrm{t}$ – double scalar

The initial value of the independent variable $x$.
Constraint:
${\mathbf{t}}\ne {\mathbf{tend}}$.
 3:
$\mathrm{tend}$ – double scalar

The end point of the range of integration. If
${\mathbf{tend}}<{\mathbf{t}}$ on initial entry, integration will proceed in the negative direction.
tend may be reset, in the direction of integration, before any continuation call.
 4:
$\mathrm{y}\left({\mathbf{neq}}\right)$ – double array

The initial values of the solution ${y}_{1},{y}_{2},\dots ,{y}_{{\mathbf{neq}}}$.
 5:
$\mathrm{yp}\left({\mathbf{neq}}\right)$ – double array

The initial values of the derivatives ${y}_{1}^{\prime},{y}_{2}^{\prime},\dots ,{y}_{{\mathbf{neq}}}^{\prime}$.
 6:
$\mathrm{ydp}\left({\mathbf{neq}}\right)$ – double array

Must be unchanged from a previous call to nag_ode_ivp_2nd_rkn (d02la).
 7:
$\mathrm{rwork}\left(\mathit{lrwork}\right)$ – double array

This
must be the same argument
rwork as supplied to
nag_ode_ivp_2nd_rkn_setup (d02lx). It is used to pass information from
nag_ode_ivp_2nd_rkn_setup (d02lx) to
nag_ode_ivp_2nd_rkn (d02la), and from
nag_ode_ivp_2nd_rkn (d02la) to both
nag_ode_ivp_2nd_rkn_diag (d02ly) and
nag_ode_ivp_2nd_rkn_interp (d02lz). Therefore the contents of this array
must not be changed before the call to
nag_ode_ivp_2nd_rkn (d02la) or calling either of the functions
nag_ode_ivp_2nd_rkn_diag (d02ly) and
nag_ode_ivp_2nd_rkn_interp (d02lz).
Optional Input Parameters
 1:
$\mathrm{neq}$ – int64int32nag_int scalar

Default:
the dimension of the arrays
y,
yp,
ydp. (An error is raised if these dimensions are not equal.)
The number of secondorder ordinary differential equations to be solved by
nag_ode_ivp_2nd_rkn (d02la). It must contain the same value as the argument
neq used in a prior call to
nag_ode_ivp_2nd_rkn_setup (d02lx).
Constraint:
${\mathbf{neq}}\ge 1$.
Output Parameters
 1:
$\mathrm{t}$ – double scalar

The value of the independent variable, which is usually
tend, unless an error has occurred or the code is operating in onestep mode. If the integration is to be continued, possibly with a new value for
tend,
t must not be changed.
 2:
$\mathrm{y}\left({\mathbf{neq}}\right)$ – double array

The computed values of the solution at the exit value of
t. If the integration is to be continued, possibly with a new value for
tend, these values must not be changed.
 3:
$\mathrm{yp}\left({\mathbf{neq}}\right)$ – double array

The computed values of the derivatives at the exit value of
t. If the integration is to be continued, possibly with a new value for
tend, these values must not be changed.
 4:
$\mathrm{ydp}\left({\mathbf{neq}}\right)$ – double array

The computed values of the second derivative at the exit value of
t, unless illegal input is detected, in which case the elements of
ydp may not have been initialized. If the integration is to be continued, possibly with a new value for
tend, these values must not be changed.
 5:
$\mathrm{rwork}\left(\mathit{lrwork}\right)$ – double array

 6:
$\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:
 ${\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.
 ${\mathbf{ifail}}=2$

The maximum number of steps has been attempted. (See argument
maxstp in
nag_ode_ivp_2nd_rkn_setup (d02lx).) If integration is to be continued then you need only reset
ifail and call the function again and a further
maxstp steps will be attempted.
 ${\mathbf{ifail}}=3$

In order to satisfy the error requirements, the step size needed is too small for the
machine precision being used.
 ${\mathbf{ifail}}=4$

The code has detected two successive error exits at the current value of $x$ and cannot proceed. Check all input variables.
 ${\mathbf{ifail}}=5$

The code has detected inefficient use of the integration method. The step size has been reduced by a significant amount too often in order to hit the output points specified by
tend. (Of the last
$100$ or more successful steps more than
$10\%$ are steps with sizes that have had to be reduced by a factor of greater than a half.)
 ${\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 accuracy of integration is determined by the arguments
tol,
thres and
thresp in a prior call to
nag_ode_ivp_2nd_rkn_setup (d02lx). Note that only the local error at each step is controlled by these arguments. The error estimates obtained are not strict bounds but are usually reliable over one step. Over a number of steps the overall error may accumulate in various ways, depending on the system. The code is designed so that a reduction in
tol should lead to an approximately proportional reduction in the error. You are strongly recommended to call
nag_ode_ivp_2nd_rkn (d02la) with more than one value for
tol and to compare the results obtained to estimate their accuracy.
The accuracy obtained depends on the type of error test used. If the solution oscillates around zero a relative error test should be avoided, whereas if the solution is exponentially increasing an absolute error test should not be used. For a description of the error test see the specifications of the arguments
tol,
thres and
thresp in function document
nag_ode_ivp_2nd_rkn_setup (d02lx).
Further Comments
If
nag_ode_ivp_2nd_rkn (d02la) fails with
${\mathbf{ifail}}={\mathbf{3}}$ then the value of
tol may be so small that a solution cannot be obtained, in which case the function should be called again with a larger value for
tol. If the accuracy requested is really needed then you should consider whether there is a more fundamental difficulty. For example:
(a) 
in the region of a singularity the solution components will usually be of a large magnitude. nag_ode_ivp_2nd_rkn (d02la) could be used in onestep mode to monitor the size of the solution with the aim of trapping the solution before the singularity. In any case numerical integration cannot be continued through a singularity, and analytical treatment may be necessary; 
(b) 
if the solution contains fast oscillatory components, the function will require a very small step size to preserve stability. This will usually be exhibited by excessive computing time and sometimes an error exit with ${\mathbf{ifail}}={\mathbf{3}}$. The Runge–Kutta–Nystrom methods are not efficient in such cases and you should consider reposing your problem as a system of firstorder ordinary differential equations and then using a function from Subchapter D02M–N with the Blend formulae (see nag_ode_ivp_stiff_dassl (d02mv)). 
nag_ode_ivp_2nd_rkn (d02la) can be used for producing results at short intervals (for example, for tabulation), in two ways. By far the less efficient is to call
nag_ode_ivp_2nd_rkn (d02la) successively over short intervals,
$t+\left(i1\right)\times h$ to
$t+i\times h$, although this is the only way if the higher order method has been selected and precisely
not what it is intended for. A more efficient way,
only for use when the lower order method has been selected, is to use
nag_ode_ivp_2nd_rkn (d02la) in onestep mode. The output values of arguments
y,
yp,
ydp,
t and
rwork are set correctly for a call to
nag_ode_ivp_2nd_rkn_interp (d02lz) to compute the solution and derivative at the required points.
Example
This example solves the following system (the two body problem)
over the range
$\left[0,20\right]$ with initial conditions
${y}_{1}=1.0\epsilon $,
${y}_{2}=0.0$,
${y}_{1}^{\prime}=0.0$ and
${y}_{2}^{\prime}=\sqrt{\left(\frac{1+\epsilon}{1\epsilon}\right)}$ where
$\epsilon $, the eccentricity, is
$0.5$. The system is solved using the lower order method with relative local error tolerances
$\text{1.0e\u22124}$ and
$\text{1.0e\u22125}$ and default threshold tolerances.
nag_ode_ivp_2nd_rkn (d02la) is used in onestep mode (
${\mathbf{onestp}}=\mathit{true}$) and
nag_ode_ivp_2nd_rkn_interp (d02lz) provides solution values at intervals of
$2.0$.
Open in the MATLAB editor:
d02la_example
function d02la_example
fprintf('d02la example results\n\n');
ecc = 0.5;
tend = 20;
tstart = 0;
tinc = 0.1;
neq = 2;
lrwork = 16+20*neq;
nwant = neq;
thres = zeros(neq, 1);
thresp = zeros(neq, 1);
fprintf('d02la example program results \n\n');
for itol = 4:5
tol = 10^(itol);
thres(1) = 0.0;
thresp(1) = 0.0;
h = 0.0;
maxstp = 0;
start = true;
high = false;
onestp = true;
rwork = zeros(lrwork,1);
t = 0;
y = [1ecc; 0];
yp = [0; sqrt((1+ecc)/(1ecc))];
ydp = zeros(neq, 1);
[start, rwork, ifail] = d02lx(h, tol, thres, thresp, int64(maxstp),...
start, onestp, high, rwork, 'neq', int64(neq));
fprintf('Calculation with tol = %1.1e \n\n',tol);
t = tstart;
tnext = t + tinc;
if itol == 4
ncall = 1;
tkeep = tnext;
ykeep = y;
nkeep = 1;
end
while (t < tend)
[tout, y, yp, ydp, rwork, ifail] = d02la(@fcn, t, tend, y,...
yp, ydp, rwork);
t = tout;
while (tnext <= t)
[ywant, ypwant, ifail] = d02lz(t, y, yp, int64(neq), tnext,...
rwork, 'neq', int64(neq));
ncall = ncall+1;
if itol == 4
tkeep(ncall) = tnext;
ykeep(:,ncall) = ywant;
end
tnext = tnext + tinc;
end
end
[hnext, hused, hstart, nsucc, nfail, natt, thres, thresp, ifail] =...
d02ly(int64(neq), rwork);
fprintf('\nNumber of successful steps = %1.0f\n',nsucc);
fprintf('Number of failed steps = %1.0f\n',nfail);
fprintf('\n\n');
end
nhalf = 1;
while ykeep(2,nhalf) >= 0
nhalf = nhalf+1;
end
fig1 = figure;
display_plot(ykeep, nhalf);
function ydp = fcn(neq, t, y)
r = sqrt(y(1)^2+y(2)^2)^3;
f = zeros(2,1);
ydp(1) = y(1)/r;
ydp(2) = y(2)/r;
function display_plot(ykeep, nhalf)
norb1 = 2*(nhalf1); norb2 = 2*norb1; norb3 = 3*norb1;
plot(ykeep(1,1:norb1), ykeep(2,1:norb1), '+', ...
ykeep(1,norb1+1:norb2), ykeep(2,norb1+1:norb2)+0.1, 'x', ...
ykeep(1,norb2+1:norb3), ykeep(2,norb2+1:norb3)+0.2, ':*');
title({'Secondorder ODE Solution using RungeKuttaNystrom.', ...
'The Twobody Problem (using shifts to distinguish orbits)'});
xlabel('x'); ylabel('y');
legend('1st orbit', '2nd orbit+(0,0.1)', '3rd orbit+(0,0.2)', ...
'Location', 'Best');
d02la example results
d02la example program results
Calculation with tol = 1.0e04
Number of successful steps = 108
Number of failed steps = 16
Calculation with tol = 1.0e05
Number of successful steps = 169
Number of failed steps = 15
PDF version (NAG web site
, 64bit version, 64bit version)
© The Numerical Algorithms Group Ltd, Oxford, UK. 2009–2015