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

 Contents

    1  Purpose
    2  Syntax
    7  Accuracy
    9  Example

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_withdraw_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=ft,y  given  yt0=y0  
where y is the vector of n solution components and t 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=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
n, the number of ordinary differential equations in the system to be solved by the integration function.
Constraint: neq1.
2:     twant – double scalar
t, 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'
Compute the approximate solution only.
reqest='D'
Compute the approximate first derivative of the solution only.
reqest='B'
Compute both the approximate solution and its first derivative.
Constraint: reqest='S', 'D' or 'B'.
4:     nwant int64int32nag_int scalar
The number of components of the solution to be computed. The first nwant components are evaluated.
Constraint: 1nwantn, where n 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 fi (that is the first derivatives yi) for given values of the arguments t,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
t, the current value of the independent variable.
2:     y: – double array
The current values of the dependent variables, yi, for i=1,2,,n.

Output Parameters

1:     yp: – double array
The values of fi, for i=1,2,,n.
6:     work: – double array
The dimension of the array work must be at least lenwrk (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:     wrkintlenint – double array
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.
Constraints:

Output Parameters

1:     ywant: – double array
The dimension of the array ywant will be nwant if reqest='S' or 'B' and 1 otherwise
An approximation to the first nwant components of the solution at twant if reqest='S' or 'B'. Otherwise ywant is not defined.
2:     ypwant: – double array
The dimension of the array ypwant will be nwant if reqest='D' or 'B' and 1 otherwise
An approximation to the first nwant components of the first derivative at twant if reqest='D' or 'B'. Otherwise ypwant is not defined.
3:     work: – double array
The dimension of the array work will be lenwrk (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:     wrkintlenint – double array
The contents are modified.
5:     ifail int64int32nag_int scalar
ifail=0 unless the function detects an error (see Error Indicators and Warnings).

Error Indicators and Warnings

Errors or warnings detected by the function:
   ifail=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=3. You cannot continue integrating the problem.
   ifail=-99
An unexpected error has been triggered by this routine. Please contact NAG.
   ifail=-399
Your licence key may have expired or may not have been installed correctly.
   ifail=-999
Dynamic memory allocation failed.

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 ,   y0=0,   y0=1  
reposed as
y1 = y2  
y2 = -y1  
over the range 0,2π with initial conditions y1=0.0 and y2=1.0. Relative error control is used with threshold values of 1.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 across the range whenever these points lie in one of those integration steps. A moderate order Runge–Kutta method (method=2) is also used with tolerances tol=1.0e−3 and tol=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 argument nwant.
function d02px_example


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

% Set initial conditions and input
method = int64(2);
tstart = 0;
tend   = 2*pi;
yinit  = [0;1];
hstart = 0;
thresh = [1e-08; 1e-08];
n      = int64(2);
task   = 'Complex Task';
errass = false;
lenwrk = int64(32*n);
reqest = 'Both';
nwant  = int64(1);
tol    =  1.0E-2;
npts   = 32;
work   = zeros(n+5*nwant, 1);
twant  = zeros(npts+1, 1);
ywant  = zeros(npts, nwant);
ypwant = zeros(npts, nwant);
lenint = n+5*nwant;
wrkint = zeros(lenint,1);

% Set control for printing solution
tinc   = (tend-tstart)/npts;

% Run through the calculation twice with two tolerance values
for i = 1:2

  tol = tol*0.1;

  % Call setup function
  [work, ifail] = d02pv(tstart, yinit, tend, tol, thresh, method, task, ...
                        errass, lenwrk);

  % Set up first point at which solution is required
  twant(1) = tstart;
  twant(2) = tstart + tinc;
  ywant(1,1) = yinit(1);
  ypwant(1,1) = yinit(2);
  t = tstart;
  j=1;

  % Integrate by steps until tend is reached or error is encountered.
  while t < tend

  % Integrate one step from t, updating t.
    [t, y, yp, work, ifail] = d02pd(@f, n, work);

    % Interpolate at required additional points up to t
    while twant(j+1) < t
      j=j+1;
      % Interpolate and print solution at t = twant.
      [ywant(j,:), ypwant(j,:), work, wrkint, ifail] = d02px(n, twant(j), ...
                                             reqest, nwant, @f, work, wrkint);

      % Set next required solution point
      twant(j+1) = twant(j) + tinc;
    end
  end

  % Get integration statistics.
  [fevals, stepcost, waste, stepsok, hnext, fail] = d02py;
  fprintf('\nCalculation with TOL = %8.1e\n\n', tol);
  fprintf('     Number of evaluations of f = %d\n', fevals);

end

% Plot results
fig1 = figure;
title('Simple Sine Solution, tol=1e-4');
hold on;
xlabel('t');
ylabel('Solution');
plot(twant(1:npts), ywant(:, 1), '-xr');
text(3, 0.25, 'Solution', 'Color', 'r');
plot(twant(1:npts), ypwant(:, 1), '-xg');
text(1.45, 0.25, 'Derivative', 'Color', 'g');
hold off



function [yp] = f(t, y)
% Evaluate derivative vector.
yp = zeros(2, 1);
yp(1) =  y(2);
yp(2) = -y(1);
d02px example results


Calculation with TOL =  1.0e-03

     Number of evaluations of f = 68

Calculation with TOL =  1.0e-04

     Number of evaluations of f = 105
d02px_fig1.png

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