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_roots_revcom (d02qg)

Purpose

nag_ode_ivp_adams_roots_revcom (d02qg) is a reverse communication function for integrating a non-stiff system of first-order ordinary differential equations using a variable-order variable-step Adams method. A root-finding facility is provided.

Syntax

[t, y, root, irevcm, trvcm, yrvcm, yprvcm, kgrvcm, rwork, iwork, ifail] = d02qg(t, y, tout, neqg, irevcm, grvcm, kgrvcm, rwork, iwork, 'neqf', neqf)
[t, y, root, irevcm, trvcm, yrvcm, yprvcm, kgrvcm, rwork, iwork, ifail] = nag_ode_ivp_adams_roots_revcom(t, y, tout, neqg, irevcm, grvcm, kgrvcm, rwork, iwork, 'neqf', neqf)
Note: the interface to this routine has changed since earlier releases of the toolbox:
Mark 22: lrwork, liwork have been removed from the interface
.

Description

Given the initial values x,y1,y2,,yneqfx,y1,y2,,yneqf nag_ode_ivp_adams_roots_revcom (d02qg) integrates a non-stiff system of first-order differential equations of the type
yi = fi(x,y1,y2,,yneqf),  i = 1,2,,neqf,
yi=fi(x,y1,y2,,yneqf),  i=1,2,,neqf,
from x = tx=t to x = toutx=tout using a variable-order variable-step Adams method. You define the system by reverse communication, evaluating fifi in terms of xx and y1,y2,,yneqfy1,y2,,yneqf, and y1,y2,,yneqfy1,y2,,yneqf are supplied at x = tx=t by nag_ode_ivp_adams_roots_revcom (d02qg). The function is capable of finding roots (values of xx) of prescribed event functions of the form
gj (x,y,y) = 0 ,   j = 1,2,,neqg.
gj (x,y,y) = 0 ,   j=1,2,,neqg.
Each gjgj is considered to be independent of the others so that roots are sought of each gjgj individually. The root reported by the function will be the first root encountered by any gjgj. 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 gjgj are looked for at the ends of each integration step. You also define each gjgj 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).
For a description of the practical implementation of an Adams formula see Shampine and Gordon (1975).

References

Shampine L F and Gordon M K (1975) Computer Solution of Ordinary Differential Equations – The Initial Value Problem W H Freeman & Co., San Francisco
Shampine L F and Watts H A (1979) DEPAC – design of a user oriented package of ODE solvers Report SAND79-2374 Sandia National Laboratory
Watts H A (1985) RDEAM – An Adams ODE code with root solving capability Report SAND85-1595 Sandia National Laboratory

Parameters

Note: this function uses reverse communication. Its use involves an initial entry, intermediate exits and re-entries, and a final exit, as indicated by the parameter irevcm. Between intermediate exits and re-entries, all parameters other than grvcm and rwork must remain unchanged.

Compulsory Input Parameters

1:     t – double scalar
On initial entry: that is after a call to nag_ode_ivp_adams_setup (d02qw) with statef = 'S'statef='S', t must be set to the initial value of the independent variable xx.
2:     y(neqf) – double array
neqf, the dimension of the array, must satisfy the constraint neqf1neqf1.
On initial entry: the initial values of the solution y1,y2,,yneqfy1,y2,,yneqf.
3:     tout – double scalar
On initial entry: the next value of xx at which a computed solution is required. For the initial t, the input value of tout is used to determine the direction of integration. Integration is permitted in either direction. If tout = ttout=t on exit, tout must be reset beyond t in the direction of integration, before any continuation call.
4:     neqg – int64int32nag_int scalar
On initial entry: the number of event functions which you are defining for root-finding. If root-finding is not required the value for neqg must be 00. Otherwise it must be the same value as the parameter neqg used in the prior call to nag_ode_ivp_adams_setup (d02qw).
5:     irevcm – int64int32nag_int scalar
On initial entry: must have the value 00.
6:     grvcm – double scalar
On initial entry: need not be set.
On intermediate re-entry: with irevcm = 9irevcm=9, 1010, 1111 or 1212, grvcm must contain the value of gk(x,y,y)gk(x,y,y), where kk is given by kgrvcm.
7:     kgrvcm – int64int32nag_int scalar
On intermediate re-entry: with irevcm = 9irevcm=9, 1010, 1111 or 1212, kgrvcm must remain unchanged from a previous call to nag_ode_ivp_adams_roots_revcom (d02qg).
8:     rwork(lrwork) – double array
This must be the same parameter rwork as supplied to nag_ode_ivp_adams_setup (d02qw). It is used to pass information from nag_ode_ivp_adams_setup (d02qw) to nag_ode_ivp_adams_roots_revcom (d02qg), and from nag_ode_ivp_adams_roots_revcom (d02qg) to nag_ode_ivp_adams_diag (d02qx), nag_ode_ivp_adams_rootdiag (d02qy) and nag_ode_ivp_adams_interp (d02qz). Therefore the contents of this array must not be changed before the call to nag_ode_ivp_adams_roots_revcom (d02qg) or calling any of the functions nag_ode_ivp_adams_diag (d02qx), nag_ode_ivp_adams_rootdiag (d02qy) and nag_ode_ivp_adams_interp (d02qz).
9:     iwork(liwork) – int64int32nag_int array
This must be the same parameter iwork as supplied to nag_ode_ivp_adams_setup (d02qw). It is used to pass information from nag_ode_ivp_adams_setup (d02qw) to nag_ode_ivp_adams_roots_revcom (d02qg), and from nag_ode_ivp_adams_roots_revcom (d02qg) to nag_ode_ivp_adams_diag (d02qx), nag_ode_ivp_adams_rootdiag (d02qy) and nag_ode_ivp_adams_interp (d02qz). Therefore the contents of this array must not be changed before the call to nag_ode_ivp_adams_roots_revcom (d02qg) or calling any of the functions nag_ode_ivp_adams_diag (d02qx), nag_ode_ivp_adams_rootdiag (d02qy) and nag_ode_ivp_adams_interp (d02qz).

Optional Input Parameters

1:     neqf – int64int32nag_int scalar
Default: The dimension of the array y.
On initial entry: the number of first-order ordinary differential equations to be solved by nag_ode_ivp_adams_roots_revcom (d02qg). It must contain the same value as the parameter neqf used in the prior call to nag_ode_ivp_adams_setup (d02qw).
Constraint: neqf1neqf1.

Input Parameters Omitted from the MATLAB Interface

lrwork liwork

Output Parameters

1:     t – double scalar
On final exit: the value of xx at which yy has been computed. This may be an intermediate output point, a root, tout or a point at which an error has occurred. If the integration is to be continued, possibly with a new value for tout, t must not be changed.
2:     y(neqf) – double array
On final exit: 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 tout, these values must not be changed.
3:     root – logical scalar
On final exit: if root-finding was required (neqg > 0neqg>0 on entry), then root specifies whether or not the output value of the parameter t is a root of one of the event functions. If root = falseroot=false, then no root was detected, whereas root = trueroot=true indicates a root and you should make a call to nag_ode_ivp_adams_rootdiag (d02qy) for further information.
If root-finding was not required (neqg = 0neqg=0 on entry), then root = falseroot=false.
4:     irevcm – int64int32nag_int scalar
On intermediate exit: specifies what action you must take before re-entering nag_ode_ivp_adams_roots_revcom (d02qg) with irevcm unchanged.
irevcm = 1irevcm=1, 22, 33, 44, 55, 66 or 77
Indicates that you must supply y = f(x,y)y=f(x,y), where xx is given by trvcm and yiyi is returned in y(i)yi, for i = 1,2,,neqfi=1,2,,neqf when yrvcm = 0yrvcm=0 and rwork(yrvcm + i1)rworkyrvcm+i-1, for i = 1,2,,neqfi=1,2,,neqf when yrvcm0yrvcm0. yiyi should be placed in location rwork(yprvcm + i1)rworkyprvcm+i-1, for i = 1,2,,neqfi=1,2,,neqf.
irevcm = 8irevcm=8
Indicates that the current step was not successful due to error test failure. The only information supplied to you on this return is the current value of the independent variable t, as given by trvcm. No values must be changed before re-entering nag_ode_ivp_adams_roots_revcom (d02qg). This facility enables you to determine the number of unsuccessful steps.
irevcm = 9irevcm=9, 1010, 1111 or 1212
Indicates that you must supply gk(x,y,y)gk(x,y,y), where kk is given by kgrvcm, xx is given by trvcm, yiyi is given by y(i)yi and yiyi is given by rwork(yprvcm1 + i)rworkyprvcm-1+i. The result gkgk should be placed in the variable grvcm.
On final exit: has the value 00, which indicates that an output point or root has been reached or an error has occurred (see ifail).
5:     trvcm – double scalar
On intermediate exit: the current value of the independent variable.
6:     yrvcm – int64int32nag_int scalar
On intermediate exit: with irevcm = 1irevcm=1, 22, 33, 44, 55, 66, 77, 99, 1010, 1111 or 1212, yrvcm specifies the locations of the dependent variables yy for use in evaluating the differential system or the event functions.
yrvcm = 0yrvcm=0
yiyi is given by y(i)yi, for i = 1,2,,neqfi=1,2,,neqf.
yrvcm0yrvcm0
yiyi is given by rwork(yrvcm + i1)rworkyrvcm+i-1, for i = 1,2,,neqfi=1,2,,neqf.
7:     yprvcm – int64int32nag_int scalar
On intermediate exit: with irevcm = 1irevcm=1, 22, 33, 44, 55, 66 or 77, yprvcm specifies the positions in rwork at which you should place the derivatives yy. yiyi should be placed in location rwork(yprvcm + i1)rworkyprvcm+i-1, for i = 1,2,,neqfi=1,2,,neqf.
With irevcm = 9irevcm=9, 1010, 1111 or 1212, yprvcm specifies the locations of the derivatives yy for use in evaluating the event functions. yiyi is given by rwork(yprvcm + i1)rworkyprvcm+i-1, for i = 1,2,,neqfi=1,2,,neqf.
8:     kgrvcm – int64int32nag_int scalar
On intermediate exit: with irevcm = 9irevcm=9, 1010, 1111 or 1212, kgrvcm specifies which event function gk(x,y,y)gk(x,y,y) you must evaluate.
9:     rwork(lrwork) – double array
10:   iwork(liwork) – int64int32nag_int array
11:   ifail – int64int32nag_int scalar
ifail = 0ifail=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 = 1ifail=1
On entry, the integrator detected an illegal input, or nag_ode_ivp_adams_setup (d02qw) has not been called before the call to the integrator.
This error may be caused by overwriting elements of rwork and iwork.
W ifail = 2ifail=2
The maximum number of steps has been attempted (at a cost of about 22 derivative evaluations per step). (See parameter maxstp in nag_ode_ivp_adams_setup (d02qw).) 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.
  ifail = 3ifail=3
The step size needed to satisfy the error requirements is too small for the machine precision being used. (See parameter tolfac in nag_ode_ivp_adams_diag (d02qx).)
W ifail = 4ifail=4
Some error weight wiwi became zero during the integration (see parameters vectol, rtol and atol in nag_ode_ivp_adams_setup (d02qw).) Pure relative error control (atol = 0.0atol=0.0) was requested on a variable (the iith) which has now become zero. (See parameter badcmp in nag_ode_ivp_adams_diag (d02qx).) The integration was successful as far as t.
W ifail = 5ifail=5
The problem appears to be stiff (see the D02 Chapter Introduction for a discussion of the term ‘stiff’). Although it is inefficient to use this integrator to solve stiff problems, integration may be continued by resetting ifail and calling the function again.
W ifail = 6ifail=6
A change in sign of an event function has been detected but the root-finding process appears to have converged to a singular point t rather than a root. Integration may be continued by resetting ifail and calling the function again.
  ifail = 7ifail=7
The code has detected two successive error exits at the current value of t and cannot proceed. Check all input variables.

Accuracy

The accuracy of integration is determined by the parameters 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 parameters. 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 parameters vectol, rtol and atol in the function document for nag_ode_ivp_adams_setup (d02qw).
The accuracy of any roots located will depend on the accuracy of integration and may also be restricted by the numerical properties of g(x,y,y)g(x,y,y). When evaluating gg you should try to write the code so that unnecessary cancellation errors will be avoided.

Further Comments

If nag_ode_ivp_adams_roots_revcom (d02qg) fails with ifail = 3ifail=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 = 3ifail=3, but usually an error exit with ifail = 2ifail=2 or 55. 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 parameter 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 = truecrit=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

This example solves the following system (for a projectile)
y = tanφ
v = ( 0.032 tanφ )/v(0.02v)/( cosφ )
φ = (0.032)/(v2)
y=tanϕ v= -0.032 tanϕ v - 0.02v cosϕ ϕ= -0.032 v2
over an interval [0.0,10.0][0.0,10.0] starting with values y = 0.5y=0.5, v = 0.5v=0.5 and φ = π / 5ϕ=π/5 using scalar error control (vectol = falsevectol=false) until the first point where y = 0.0y=0.0 is encountered.
Also, nag_ode_ivp_adams_roots_revcom (d02qg) is used to produce output at intervals of 2.02.0.
function nag_ode_ivp_adams_roots_revcom_example
% Initialize variables and arrays.
t = 0;
y = [0.5; 0.5; 0.2*pi];
neqg = 1;
grvcm = 0;
kgrvcm = 0;
rwork = zeros(106, 1);
iwork = zeros(25, 1);

statef = 'S';
neqf = 3;
vectol = false;
atol = [1e-4];
rtol = [1e-7];
onestp = false;
crit = true;
tcrit = 10;
hmax = 2;
maxstp = 500;
alterg = false;
sophst = true;
tinc = 2;

% Initialize storage counter, and define storage arrays.
ncall = 1;
tkeep = zeros(1,1);
ykeep = zeros(1,length(y));

fprintf('nag_ode_ivp_adams_roots_revcom example program results\n\n');

% nag_ode_ivp_adams_setup is a setup routine to be called prior to nag_ode_ivp_adams_roots_revcom.
[statef, alterg, rwork, iwork, ifail] = nag_ode_ivp_adams_setup(statef, int64(neqf), ...
    vectol, atol, rtol, onestp, crit, tcrit, hmax, int64(maxstp), ...
    int64(neqg), alterg, sophst, rwork, int64(iwork));
if ifail ~= 0
    % Illegal input.  Print message and exit.
    error('Warning: nag_ode_ivp_adams_setup returned with ifail = %1d ',ifail);
end

% 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 nag_ode_ivp_adams_roots_revcom for this t value.
    irevcm = 0;
    [t, y, root, irevcm, trvcm, yrvcm, yprvcm, kgrvcm, rwork, iwork, ...
        ifail] = nag_ode_ivp_adams_roots_revcom(t, y, tout, int64(neqg), int64(irevcm), grvcm, ...
        int64(kgrvcm), rwork, int64(iwork));

    % nag_ode_ivp_adams_roots_revcom uses reverse communication.  Keep re-entering the function
    % until a final exit (irevcm = 0) occurs.
    while (irevcm ~= 0)
        if (irevcm < 8)
            % Set current values for derivatives.
            if (yrvcm == 0)
                rwork(yprvcm) = tan(y(3));
                rwork(yprvcm+1) = -0.032*tan(y(3))/y(2) - 0.02*y(2)/cos(y(3));
                rwork(yprvcm+2) = -0.032/y(2)^2;
            else
                rwork(yprvcm) = tan(rwork(yrvcm+2));
                rwork(yprvcm+1) = -0.032*tan(rwork(yrvcm+2))/rwork(yrvcm+1) - ...
                    0.02*rwork(yrvcm+1)/cos(rwork(yrvcm+2));
                rwork(yprvcm+2) = -0.032/rwork(yrvcm+1)^2;
            end
        elseif (irevcm > 8)
            grvcm = y(1);
        end

        [t, y, root, irevcm, trvcm, yrvcm, yprvcm, kgrvcm, rwork, ...
            iwork, ifail] = nag_ode_ivp_adams_roots_revcom(t, y, tout, int64(neqg), ...
            int64(irevcm), grvcm, int64(kgrvcm), rwork, int64(iwork));
    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.
fig = figure('Number', 'off');
display_plot(tkeep, ykeep)

function display_plot(tkeep, ykeep)
% Formatting for title and axis labels.
titleFmt = {'FontName', 'Helvetica', 'FontWeight', 'Bold', 'FontSize', 14};
labFmt = {'FontName', 'Helvetica', 'FontWeight', 'Bold', 'FontSize', 13};
set(gca, 'FontSize', 13); % for legend, axis tick labels, etc.
% 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', titleFmt{:});
% Label the axes.
xlabel('x', labFmt{:});
ylabel('Solution', labFmt{:});
% Add a legend.
legend('height','velocity','angle','Location','Best');
 
nag_ode_ivp_adams_roots_revcom example program results

   t        y(1)      y(2)      y(3)
 0.0000    0.5000    0.5000    0.6283
 2.0000    1.5493    0.4055    0.3066
 4.0000    1.7424    0.3743   -0.1289
 6.0000    1.0058    0.4173   -0.5507
 7.2879   -0.0000    0.4749   -0.7601
10.0000   -3.6287    0.6333   -1.0515

function d02qg_example
% Initialize variables and arrays.
t = 0;
y = [0.5; 0.5; 0.2*pi];
neqg = 1;
grvcm = 0;
kgrvcm = 0;
rwork = zeros(106, 1);
iwork = zeros(25, 1);

statef = 'S';
neqf = 3;
vectol = false;
atol = [1e-4];
rtol = [1e-7];
onestp = false;
crit = true;
tcrit = 10;
hmax = 2;
maxstp = 500;
alterg = false;
sophst = true;
tinc = 2;

% Initialize storage counter, and define storage arrays.
ncall = 1;
tkeep = zeros(1,1);
ykeep = zeros(1,length(y));

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

% d02qw is a setup routine to be called prior to d02qg.
[statef, alterg, rwork, iwork, ifail] = d02qw(statef, int64(neqf), ...
    vectol, atol, rtol, onestp, crit, tcrit, hmax, int64(maxstp), ...
    int64(neqg), alterg, sophst, rwork, int64(iwork));
if ifail ~= 0
    % Illegal input.  Print message and exit.
    error('Warning: d02qw returned with ifail = %1d ',ifail);
end

% 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 = 0;
    [t, y, root, irevcm, trvcm, yrvcm, yprvcm, kgrvcm, rwork, iwork, ...
        ifail] = d02qg(t, y, tout, int64(neqg), int64(irevcm), grvcm, ...
        int64(kgrvcm), rwork, int64(iwork), 'ifail', int64(-1));

    % d02qg uses reverse communication.  Keep re-entering the function
    % until a final exit (irevcm = 0) occurs.
    while (irevcm ~= 0)
        if (irevcm < 8)
            % Set current values for derivatives.
            if (yrvcm == 0)
                rwork(yprvcm) = tan(y(3));
                rwork(yprvcm+1) = -0.032*tan(y(3))/y(2) - 0.02*y(2)/cos(y(3));
                rwork(yprvcm+2) = -0.032/y(2)^2;
            else
                rwork(yprvcm) = tan(rwork(yrvcm+2));
                rwork(yprvcm+1) = -0.032*tan(rwork(yrvcm+2))/rwork(yrvcm+1) - ...
                    0.02*rwork(yrvcm+1)/cos(rwork(yrvcm+2));
                rwork(yprvcm+2) = -0.032/rwork(yrvcm+1)^2;
            end
        elseif (irevcm > 8)
            grvcm = y(1);
        end

        [t, y, root, irevcm, trvcm, yrvcm, yprvcm, kgrvcm, rwork, ...
            iwork, ifail] = d02qg(t, y, tout, int64(neqg), ...
            int64(irevcm), grvcm, int64(kgrvcm), rwork, int64(iwork));
    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.
fig = figure('Number', 'off');
display_plot(tkeep, ykeep)

function display_plot(tkeep, ykeep)
% Formatting for title and axis labels.
titleFmt = {'FontName', 'Helvetica', 'FontWeight', 'Bold', 'FontSize', 14};
labFmt = {'FontName', 'Helvetica', 'FontWeight', 'Bold', 'FontSize', 13};
set(gca, 'FontSize', 13); % for legend, axis tick labels, etc.
% 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', titleFmt{:});
% Label the axes.
xlabel('x', labFmt{:});
ylabel('Solution', labFmt{:});
% Add a legend.
legend('height','velocity','angle','Location','Best');
 
d02qg example program results

   t        y(1)      y(2)      y(3)
 0.0000    0.5000    0.5000    0.6283
 2.0000    1.5493    0.4055    0.3066
 4.0000    1.7424    0.3743   -0.1289
 6.0000    1.0058    0.4173   -0.5507
 7.2879   -0.0000    0.4749   -0.7601
10.0000   -3.6287    0.6333   -1.0515


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–2013