NAG Toolbox: nag_ode_ivp_2nd_rkn (d02la)

Purpose

Syntax

Description

References

Parameters

Compulsory Input Parameters

1:     $\mathrm{fcn}$ – function handle or string containing name of m-file

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 second-order 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 one-step 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

   $\mathbf{ifail}=1$

Illegal input detected, i.e., one of the following conditions:

• on any call, $\mathbf{t}=\mathbf{tend}$, or the value of neq or lrwork has been altered;
• on a continuation call, the direction of integration has been changed;
• nag_ode_ivp_2nd_rkn_setup (d02lx) had not been called previously, or the previous call to nag_ode_ivp_2nd_rkn_setup (d02lx) resulted in an error exit.

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

Further Comments

Example

Open in the MATLAB editor: d02la_example

function d02la_example


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

% Initialize variables and arrays.
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 = [1-ecc; 0];
    yp = [0; sqrt((1+ecc)/(1-ecc))];
    ydp = zeros(neq, 1);

    % d02lx is a setup routine to be called prior to d02la.
    [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);
    % fprintf('   t          y_1           y_2 ');
    t = tstart;
    tnext = t + tinc;
    % fprintf('\n %4.1f    %10.5f    %10.5f\n', t, y);

    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)
            % d02lz interpolates components of solution provided by d02la.

             [ywant, ypwant, ifail] = d02lz(t, y, yp, int64(neq), tnext,...
                rwork, 'neq', int64(neq));

            ncall = ncall+1;
            % fprintf(' %4.1f    %10.5f    %10.5f\n', tnext, ywant);

            if itol == 4
                tkeep(ncall) = tnext;
                ykeep(:,ncall) = ywant;
            end

            tnext = tnext + tinc;
        end
    end

    % d02ly is a diagnostic routine for use after calling d02la.

    [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

% Find the index of the point halfway through the first orbit (i.e. the
% first point at which the y coord is negative).  This is used to calculate
% the end points of the orbits in display_plot.
nhalf = 1;
while ykeep(2,nhalf) >= 0
    nhalf = nhalf+1;
end
% Plot results.
fig1 = figure;
display_plot(ykeep, nhalf);


function ydp = fcn(neq, t, y)
% Evaluate second derivatives.
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)

% Find the end points of the three orbits.
norb1 = 2*(nhalf-1); norb2 = 2*norb1; norb3 = 3*norb1;

% Plot the orbits, including displacements.
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, ':*');
% Add title.
title({'Second-order ODE Solution using Runge-Kutta-Nystrom.', ...
    'The Two-body Problem (using shifts to distinguish orbits)'});
% Label the axes.
xlabel('x'); ylabel('y');
% Add a legend.
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.0e-04 


Number of successful steps = 108
Number of failed     steps = 16


Calculation with tol = 1.0e-05 


Number of successful steps = 169
Number of failed     steps = 15