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_rk_interp (d02px)

Purpose

nag_ode_ivp_rk_interp (d02px) computes the solution of a system of ordinary differential equations using interpolation anywhere on an integration step taken by nag_ode_ivp_rk_onestep (d02pd).
Note: this function is scheduled to be withdrawn, please see d02px in Advice on Replacement Calls for Withdrawn/Superseded Routines..

Syntax

[ywant, ypwant, work, wrkint, ifail] = d02px(neq, twant, reqest, nwant, f, work, wrkint, 'lenint', lenint)
[ywant, ypwant, work, wrkint, ifail] = nag_ode_ivp_rk_interp(neq, twant, reqest, nwant, f, work, wrkint, 'lenint', lenint)

Description

nag_ode_ivp_rk_interp (d02px) and its associated functions (nag_ode_ivp_rk_onestep (d02pd), nag_ode_ivp_rk_setup (d02pv), nag_ode_ivp_rk_reset_tend (d02pw), nag_ode_ivp_rk_diag (d02py) and nag_ode_ivp_rk_errass (d02pz)) solve the initial value problem for a first-order system of ordinary differential equations. The functions, based on Runge–Kutta methods and derived from RKSUITE (see Brankin et al. (1991)), integrate
y = f(t,y)  given  y(t0) = y0
y=f(t,y)  given  y(t0)=y0
where yy is the vector of nn solution components and tt is the independent variable.
nag_ode_ivp_rk_onestep (d02pd) computes the solution at the end of an integration step. Using the information computed on that step nag_ode_ivp_rk_interp (d02px) computes the solution by interpolation at any point on that step. It cannot be used if method = 3method=3 was specified in the call to setup function nag_ode_ivp_rk_setup (d02pv).

References

Brankin R W, Gladwell I and Shampine L F (1991) RKSUITE: A suite of Runge–Kutta codes for the initial value problems for ODEs SoftReport 91-S1 Southern Methodist University

Parameters

Compulsory Input Parameters

1:     neq – int64int32nag_int scalar
nn, the number of ordinary differential equations in the system to be solved by the integration function.
Constraint: neq1neq1.
2:     twant – double scalar
tt, the value of the independent variable where a solution is desired.
3:     reqest – string (length ≥ 1)
Determines whether the solution and/or its first derivative are to be computed.
reqest = 'S'reqest='S'
Compute the approximate solution only.
reqest = 'D'reqest='D'
Compute the approximate first derivative of the solution only.
reqest = 'B'reqest='B'
Compute both the approximate solution and its first derivative.
Constraint: reqest = 'S'reqest='S', 'D''D' or 'B''B'.
4:     nwant – int64int32nag_int scalar
The number of components of the solution to be computed. The first nwant components are evaluated.
Constraint: 1nwantn1nwantn, where nn is specified by neq in the prior call to nag_ode_ivp_rk_setup (d02pv).
5:     f – function handle or string containing name of m-file
f must evaluate the functions fifi (that is the first derivatives yiyi) for given values of the arguments t,yit,yi. It must be the same procedure as supplied to nag_ode_ivp_rk_onestep (d02pd).
[yp] = f(t, y)

Input Parameters

1:     t – double scalar
tt, the current value of the independent variable.
2:     y( : :) – double array
The current values of the dependent variables, yiyi, for i = 1,2,,ni=1,2,,n.

Output Parameters

1:     yp( : :) – double array
The values of fifi, for i = 1,2,,ni=1,2,,n.
6:     work( : :) – double array
Note: the dimension of the array work must be at least lenwrklenwrk (see nag_ode_ivp_rk_setup (d02pv)).
This must be the same array as supplied to nag_ode_ivp_rk_onestep (d02pd) and must remain unchanged between calls.
7:     wrkint(lenint) – double array
lenint, the dimension of the array, must satisfy the constraint
  • lenint1lenint1 if method = 1method=1 in the prior call to nag_ode_ivp_rk_setup (d02pv)
  • lenintn + 5 × nwantlenintn+5×nwant if method = 2method=2 and nn is specified by neq in the prior call to nag_ode_ivp_rk_setup (d02pv)
  • .
    Must be the same array as supplied in previous calls, if any, and must remain unchanged between calls to nag_ode_ivp_rk_interp (d02px).

    Optional Input Parameters

    1:     lenint – int64int32nag_int scalar
    Default: The dimension of the array wrkint.
    The dimension of the array wrkint as declared in the (sub)program from which nag_ode_ivp_rk_interp (d02px) is called.
    Constraints:

    Input Parameters Omitted from the MATLAB Interface

    None.

    Output Parameters

    1:     ywant( : :) – double array
    Note: the dimension of the array ywant must be at least nwantnwant if reqest = 'S'reqest='S' or 'B''B', and at least 11 otherwise.
    An approximation to the first nwant components of the solution at twant if reqest = 'S'reqest='S' or 'B''B'. Otherwise ywant is not defined.
    2:     ypwant( : :) – double array
    Note: the dimension of the array ypwant must be at least nwantnwant if reqest = 'D'reqest='D' or 'B''B', and at least 11 otherwise.
    An approximation to the first nwant components of the first derivative at twant if reqest = 'D'reqest='D' or 'B''B'. Otherwise ypwant is not defined.
    3:     work( : :) – double array
    Note: the dimension of the array work must be at least lenwrklenwrk (see nag_ode_ivp_rk_setup (d02pv)).
    Contains information about the integration for use on subsequent calls to nag_ode_ivp_rk_onestep (d02pd) or other associated functions.
    4:     wrkint(lenint) – double array
    The contents are modified.
    5:     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:
      ifail = 1ifail=1
    On entry, an invalid input value for nwant or lenint was detected or an invalid call to nag_ode_ivp_rk_interp (d02px) was made, for example without a previous call to the integration function nag_ode_ivp_rk_onestep (d02pd), or after an error return from nag_ode_ivp_rk_onestep (d02pd), or if nag_ode_ivp_rk_onestep (d02pd) was being used with method = 3method=3. You cannot continue integrating the problem.

    Accuracy

    The computed values will be of a similar accuracy to that computed by nag_ode_ivp_rk_onestep (d02pd).

    Further Comments

    None.

    Example

    This example solves the equation
    y = y ,   y(0) = 0,   y(0) = 1
    y = -y ,   y(0)=0,   y(0)=1
    reposed as
    y1 = y2
    y1 = y2
    y2 = y1
    y2 = -y1
    over the range [0,2π][0,2π] with initial conditions y1 = 0.0y1=0.0 and y2 = 1.0y2=1.0. Relative error control is used with threshold values of 1.0e−81.0e−8 for each solution component. nag_ode_ivp_rk_onestep (d02pd) is used to integrate the problem one step at a time and nag_ode_ivp_rk_interp (d02px) is used to compute the first component of the solution and its derivative at intervals of length π / 8π/8 across the range whenever these points lie in one of those integration steps. A moderate order Runge–Kutta method (method = 2method=2) is also used with tolerances tol = 1.0e−3tol=1.0e−3 and tol = 1.0e−4tol=1.0e−4 in turn so that solutions may be compared. The value of ππ is obtained by using nag_math_pi (x01aa).
    Note that the length of work is large enough for any valid combination of input arguments to nag_ode_ivp_rk_setup (d02pv) and the length of wrkint is large enough for any valid value of the parameter nwant.
    function nag_ode_ivp_rk_interp_example
    % Initialize variables and arrays.
    tstart = 0;
    hstart = 0;
    ystart = [0; 1];
    tend = 2*pi;
    thres = [1e-08; 1e-08];
    method = 2;
    task = 'Complex Task';
    errass = false;
    neq = 2;
    lenwrk = 32*neq;
    twant = pi/8;
    reqest = 'Both';
    nwant = 1;
    wrkint = zeros(7, 1);
    npts = 16;
    lenint = neq+5*nwant;
    tinc = tend/npts;
    
    fprintf('nag_ode_ivp_rk_interp example program results\n\n');
    
    % We run through the calculation twice with two tolerance values; we
    % output and plot the results for both of them.
    tol = [1e-03; 1e-04];
    for itol = 1:2
    
        % nag_ode_ivp_rk_setup is a setup routine to be called prior to nag_ode_ivp_rk_onestep.
        [work, ifail] = nag_ode_ivp_rk_setup(tstart, ystart, tend, tol(itol), thres, ...
            int64(method), task, errass, int64(lenwrk), 'neq', ...
            int64(neq));
        if ifail ~= 0
            % Unsuccessful call.  Print message and exit.
            error('Warning: nag_ode_ivp_rk_setup returned with ifail = %1d ',ifail);
        end
    
        fprintf('Calculation with tol = %1.3e\n\n',tol(itol));
        disp('  t           y1           y1''');
        fprintf('%1.3f     %8.4f     %8.4f\n', tstart, ystart(1), ystart(2));
    
        % Initialize the storage counter, store the current values.
        ncall = 1;
        y1keep(ncall) = ystart(1);
        y2keep(ncall) = ystart(2);
        tkeep(ncall) = tstart;
    
        % Initialize the current t values, and start looping over t. tnow
        % nmarks the end of the integration step for nag_ode_ivp_rk_onestep; twant is the
        % point on that step at which nag_ode_ivp_rk_interp calculates the solution
        % by interpolation.
        tnow = tstart;
        twant = tinc;
        while (tnow < tend)
    
            [tnow, ynow, ypnow, work, ifail] = nag_ode_ivp_rk_onestep(@f, ...
                int64(neq), work);
            if ifail ~= 0
                % Unsuccessful call.  Print message and exit.
                error('Warning: nag_ode_ivp_rk_onestep returned with ifail = %1d ',...
                    ifail);
            end
    
            while (twant < tnow)
    
                % Call nag_ode_ivp_rk_interp to obtain solution by interpolation on results
                % from nag_ode_ivp_rk_onestep.  Note that the function file must be the same
                % for both routines.
                [ywant, ypwant, work, wrkint, ifail] = nag_ode_ivp_rk_interp(int64(neq), ...
                    twant, reqest, int64(nwant), @f, work, wrkint,...
                    'lenint', int64(lenint));
    
                % Output results, then store them for plotting.
                fprintf('%1.3f     %8.4f     %8.4f\n', twant, ywant, ypwant);
                ncall = ncall+1;
                y1keep(ncall) = ywant;
                y2keep(ncall) = ypwant;
                tkeep(ncall) = twant;
    
                twant = twant + tinc;
            end
        end
    
        % nag_ode_ivp_rk_diag is a diagnostic routine.
        [totfcn, stpcst, waste, stpsok, hnext, ifail] = nag_ode_ivp_rk_diag;
    
        fprintf(['\nCost of the integration in evaluations of F is ' ...
            '%1.0f\n\n'], totfcn);
    
        % Plot results.
        fig = figure('Number', 'off');
        display_plot(tkeep, y1keep, y2keep, tol(itol))
    end
    
    function [yp] = f(t, y)
    % Evaluate derivative vector.
    yp = zeros(2, 1);
    yp(1) =  y(2);
    yp(2) = -y(1);
    function display_plot(tkeep, y1keep, y2keep, tol)
    % 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, y1keep, '-+', tkeep, y2keep, ':x')
    % Add title.
    title(['Simple Sine Solution, TOL = ',num2str(tol)], titleFmt{:});
    % Label the axes.
    xlabel('t', labFmt{:});
    ylabel('Solution', labFmt{:});
    % Add a legend.
    legend('solution','derivative','Location','Best');
    
     
    nag_ode_ivp_rk_interp example program results
    
    Calculation with tol = 1.000e-03
    
      t           y1           y1'
    0.000       0.0000       1.0000
    0.393       0.3827       0.9239
    0.785       0.7071       0.7071
    1.178       0.9239       0.3826
    1.571       1.0000      -0.0001
    1.963       0.9238      -0.3828
    2.356       0.7070      -0.7073
    2.749       0.3825      -0.9240
    3.142      -0.0002      -0.9999
    3.534      -0.3829      -0.9238
    3.927      -0.7072      -0.7069
    4.320      -0.9239      -0.3823
    4.712      -0.9999       0.0004
    5.105      -0.9236       0.3830
    5.498      -0.7068       0.7073
    5.890      -0.3823       0.9239
    
    Cost of the integration in evaluations of F is 68
    
    Calculation with tol = 1.000e-04
    
      t           y1           y1'
    0.000       0.0000       1.0000
    0.393       0.3827       0.9239
    0.785       0.7071       0.7071
    1.178       0.9239       0.3827
    1.571       1.0000      -0.0000
    1.963       0.9239      -0.3827
    2.356       0.7071      -0.7071
    2.749       0.3827      -0.9239
    3.142      -0.0000      -1.0000
    3.534      -0.3827      -0.9239
    3.927      -0.7071      -0.7071
    4.320      -0.9239      -0.3827
    4.712      -1.0000       0.0000
    5.105      -0.9238       0.3827
    5.498      -0.7071       0.7071
    5.890      -0.3826       0.9239
    
    Cost of the integration in evaluations of F is 105
    
    
    
    function d02px_example
    % Initialize variables and arrays.
    tstart = 0;
    hstart = 0;
    ystart = [0; 1];
    tend = 2*pi;
    thres = [1e-08; 1e-08];
    method = 2;
    task = 'Complex Task';
    errass = false;
    neq = 2;
    lenwrk = 32*neq;
    twant = pi/8;
    reqest = 'Both';
    nwant = 1;
    wrkint = zeros(7, 1);
    npts = 16;
    lenint = neq+5*nwant;
    tinc = tend/npts;
    
    fprintf('d02px example program results\n\n');
    
    % We run through the calculation twice with two tolerance values; we
    % output and plot the results for both of them.
    tol = [1e-03; 1e-04];
    for itol = 1:2
    
        % d02pv is a setup routine to be called prior to d02pd.
        [work, ifail] = d02pv(tstart, ystart, tend, tol(itol), thres, ...
            int64(method), task, errass, int64(lenwrk), 'neq', ...
            int64(neq));
        if ifail ~= 0
            % Unsuccessful call.  Print message and exit.
            error('Warning: d02pv returned with ifail = %1d ',ifail);
        end
    
        fprintf('Calculation with tol = %1.3e\n\n',tol(itol));
        disp('  t           y1           y1''');
        fprintf('%1.3f     %8.4f     %8.4f\n', tstart, ystart(1), ystart(2));
    
        % Initialize the storage counter, store the current values.
        ncall = 1;
        y1keep(ncall) = ystart(1);
        y2keep(ncall) = ystart(2);
        tkeep(ncall) = tstart;
    
        % Initialize the current t values, and start looping over t. tnow
        % nmarks the end of the integration step for d02pd; twant is the
        % point on that step at which d02px calculates the solution
        % by interpolation.
        tnow = tstart;
        twant = tinc;
        while (tnow < tend)
    
            [tnow, ynow, ypnow, work, ifail] = d02pd(@f, ...
                int64(neq), work);
            if ifail ~= 0
                % Unsuccessful call.  Print message and exit.
                error('Warning: d02pd returned with ifail = %1d ',...
                    ifail);
            end
    
            while (twant < tnow)
    
                % Call d02px to obtain solution by interpolation on results
                % from d02pd.  Note that the function file must be the same
                % for both routines.
                [ywant, ypwant, work, wrkint, ifail] = d02px(int64(neq), ...
                    twant, reqest, int64(nwant), @f, work, wrkint,...
                    'lenint', int64(lenint));
    
                % Output results, then store them for plotting.
                fprintf('%1.3f     %8.4f     %8.4f\n', twant, ywant, ypwant);
                ncall = ncall+1;
                y1keep(ncall) = ywant;
                y2keep(ncall) = ypwant;
                tkeep(ncall) = twant;
    
                twant = twant + tinc;
            end
        end
    
        % d02py is a diagnostic routine.
        [totfcn, stpcst, waste, stpsok, hnext, ifail] = d02py;
    
        fprintf(['\nCost of the integration in evaluations of F is ' ...
            '%1.0f\n\n'], totfcn);
    
        % Plot results.
        fig = figure('Number', 'off');
        display_plot(tkeep, y1keep, y2keep, tol(itol))
    end
    
    function [yp] = f(t, y)
    % Evaluate derivative vector.
    yp = zeros(2, 1);
    yp(1) =  y(2);
    yp(2) = -y(1);
    function display_plot(tkeep, y1keep, y2keep, tol)
    % 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, y1keep, '-+', tkeep, y2keep, ':x')
    % Add title.
    title(['Simple Sine Solution, TOL = ',num2str(tol)], titleFmt{:});
    % Label the axes.
    xlabel('t', labFmt{:});
    ylabel('Solution', labFmt{:});
    % Add a legend.
    legend('solution','derivative','Location','Best');
    
     
    d02px example program results
    
    Calculation with tol = 1.000e-03
    
      t           y1           y1'
    0.000       0.0000       1.0000
    0.393       0.3827       0.9239
    0.785       0.7071       0.7071
    1.178       0.9239       0.3826
    1.571       1.0000      -0.0001
    1.963       0.9238      -0.3828
    2.356       0.7070      -0.7073
    2.749       0.3825      -0.9240
    3.142      -0.0002      -0.9999
    3.534      -0.3829      -0.9238
    3.927      -0.7072      -0.7069
    4.320      -0.9239      -0.3823
    4.712      -0.9999       0.0004
    5.105      -0.9236       0.3830
    5.498      -0.7068       0.7073
    5.890      -0.3823       0.9239
    
    Cost of the integration in evaluations of F is 68
    
    Calculation with tol = 1.000e-04
    
      t           y1           y1'
    0.000       0.0000       1.0000
    0.393       0.3827       0.9239
    0.785       0.7071       0.7071
    1.178       0.9239       0.3827
    1.571       1.0000      -0.0000
    1.963       0.9239      -0.3827
    2.356       0.7071      -0.7071
    2.749       0.3827      -0.9239
    3.142      -0.0000      -1.0000
    3.534      -0.3827      -0.9239
    3.927      -0.7071      -0.7071
    4.320      -0.9239      -0.3827
    4.712      -1.0000       0.0000
    5.105      -0.9238       0.3827
    5.498      -0.7071       0.7071
    5.890      -0.3826       0.9239
    
    Cost of the integration in evaluations of F is 105
    
    
    

    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