d02qg:: Ordinary Differential EquationsIntegrators for Stiff Ordinary Differential Systems (NAG Toolbox)

Given the initial values

x, y_{1}, y_{2}, \dots, y_{neqf}

nag_ode_ivp_adams_roots_revcom (d02qg) integrates a non-stiff system of first-order differential equations of the type

y_{i}^{'} = f_{i} (x, y_{1}, y_{2}, \dots, y_{neqf}), i = 1, 2, \dots, neqf,

from

x = t

x = tout

using a variable-order variable-step Adams' method. You define the system by reverse communication, evaluating

f_{i}

in terms of

x

and

y_{1}, y_{2}, \dots, y_{neqf}

, and

y_{1}, y_{2}, \dots, y_{neqf}

are supplied at

x = t

by nag_ode_ivp_adams_roots_revcom (d02qg). The function is capable of finding roots (values of

x

) of prescribed event functions of the form

g_{j} (x, y, y^{'}) = 0, j = 1, 2, \dots, neqg .

Each

g_{j}

is considered to be independent of the others so that roots are sought of each

g_{j}

individually. The root reported by the function will be the first root encountered by any

g_{j}

. Two techniques for determining the presence of a root in an integration step are available: the sophisticated method described in Watts (1985) and a simplified method whereby sign changes in each

g_{j}

are looked for at the ends of each integration step. You also define each

g_{j}

by reverse communication. In one-step mode the function returns an approximation to the solution at each integration point. In interval mode this value is returned at the end of the integration range. If a root is detected this approximation is given at the root. You select the mode of operation, the error control, the root-finding technique and various optional inputs by a prior call to the setup function nag_ode_ivp_adams_setup (d02qw).

References

Parameters

Compulsory Input Parameters

Optional Input Parameters

Output Parameters

Error Indicators and Warnings

Cases prefixed with W are classified as warnings and do not generate an error of type NAG:error_n. See nag_issue_warnings.

Accuracy

The accuracy of integration is determined by the arguments vectol, rtol and atol in a prior call to nag_ode_ivp_adams_setup (d02qw). 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 properties of the differential equation system. The code is designed so that a reduction in the tolerances should lead to an approximately proportional reduction in the error. You are strongly recommended to call nag_ode_ivp_adams_roots_revcom (d02qg) with more than one set of tolerances 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. If different accuracies are required for different components of the solution then a component-wise error test should be used. For a description of the error test see the specifications of the arguments vectol, rtol and atol in the function document for nag_ode_ivp_adams_setup (d02qw).

Further Comments

If nag_ode_ivp_adams_roots_revcom (d02qg) fails with

ifail = 3

then the combination of atol and rtol may be so small that a solution cannot be obtained, in which case the function should be called again with larger values for rtol and/or atol (see nag_ode_ivp_adams_setup (d02qw)). 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_adams_roots_revcom (d02qg) could be used in one-step 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)

for ‘stiff’ equations, where the solution contains rapidly decaying 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

ifail = 3

, but usually an error exit with

ifail = 2

5

. The Adams' methods are not efficient in such cases and you should consider using a function from the Sub-chapter D02M–N. A high proportion of failed steps (see argument nfail in nag_ode_ivp_adams_diag (d02qx)) may indicate stiffness but there may be other reasons for this phenomenon.

nag_ode_ivp_adams_roots_revcom (d02qg) can be used for producing results at short intervals (for example, for graph plotting); you should set

crit = true

and tcrit to the last output point required in a prior call to nag_ode_ivp_adams_setup (d02qw) and then set tout appropriately for each output point in turn in the call to nag_ode_ivp_adams_roots_revcom (d02qg)

Example

function d02qg_example


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

% Initialize setup variables and arrays.
statef = 'S';
neqf   = int64(3);
vectol = false;
atol   = [1e-4];
rtol   = [1e-7];
onestp = false;
crit   = true;
tcrit  = 10;
hmax   = 2;
maxstp = int64(500);
neqg   = int64(1);
alterg = false;
sophst = true;
rwork  = zeros(106, 1);
iwork  = zeros(25, 1, 'int64');

% d02qw is a setup routine to be called prior to d02qg.
[statef, alterg, rwork, iwork, ifail] = ...
         d02qw(...
               statef, neqf, vectol, atol, rtol, onestp, crit, tcrit, ...
               hmax, maxstp, neqg, alterg, sophst, rwork, iwork);

%Initialize integrator input variables
t    = 0;
y    = [0.5; 0.5; 0.2*pi];
grv  = 0;
kgrv = int64(0);

% Initialize output variables and arrays.
tinc = 2;
ncall = 1;
tkeep = zeros(1,1);
ykeep = zeros(1,neqf);


% Output initial results, then store them for plotting.
fprintf('   t        y(1)      y(2)      y(3)\n');
fprintf('%7.4f   %7.4f   %7.4f   %7.4f\n', t, y(1), y(2), y(3));
tkeep(ncall)   = t;
ykeep(ncall,:) = y;

for j=1:5

    tout = double(j)*tinc;
    % Initial call of d02qg for this t value.
    irevcm     = int64(0);
    first_call = true;

    % d02qg uses reverse communication.  Keep re-entering the function
    % until a final exit (irevcm = 0) occurs.
    while (irevcm ~= 0 || first_call)
      first_call = false;
        [t, y, root, irevcm, trv, yrv, yprv, kgrv, rwork, iwork, ifail] = ...
            d02qg(...
                  t, y, tout, neqg, irevcm, grv, kgrv, rwork, iwork);
        if (irevcm > 0 && irevcm < 8)
            % Set current values for derivatives.
            if (yrv == 0)
                rwork(yprv)   = tan(y(3));
                rwork(yprv+1) = -0.032*tan(y(3))/y(2) - 0.02*y(2)/cos(y(3));
                rwork(yprv+2) = -0.032/y(2)^2;
            else
                y2 = rwork(yrv+1);
                y3 = rwork(yrv+2);
                rwork(yprv)   = tan(y3);
                rwork(yprv+1) = -0.032*tan(y3)/y2 - 0.02*y2/cos(y3);
                rwork(yprv+2) = -0.032/y2^2;
            end
        elseif (irevcm > 8)
            grv = y(1);
        end
    end

    % Output results, then store them for plotting.
    fprintf('%7.4f   %7.4f   %7.4f   %7.4f\n', t, y(1), y(2), y(3));
    ncall          = ncall + 1;
    tkeep(ncall)   = t;
    ykeep(ncall,:) = y;
end
% Plot results.
fig1 = figure;
display_plot(tkeep, ykeep)


function display_plot(tkeep, ykeep)
% Plot the results.
plot(tkeep, ykeep(:,1), '-+', tkeep, ykeep(:,2), '--x', tkeep, ...
     ykeep(:,3), ':*');
%Add title.
title('ODE Solution using Adams Method with Root-finding');
% Label the axes.
xlabel('x');
ylabel('Solution');
% Add a legend.
legend('height','velocity','angle','Location','Best');

NAG Toolbox: nag_ode_ivp_adams_roots_revcom (d02qg)

▸▿ Contents

Purpose

Syntax

Description