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_rkts_interp (d02ps)


    1  Purpose
    2  Syntax
    7  Accuracy
    9  Example


nag_ode_ivp_rkts_interp (d02ps) computes the solution of a system of ordinary differential equations using interpolation anywhere on an integration step taken by nag_ode_ivp_rkts_onestep (d02pf).


[ywant, ypwant, wcomm, user, iwsav, rwsav, ifail] = d02ps(n, twant, ideriv, nwant, f, wcomm, iwsav, rwsav, 'lwcomm', lwcomm, 'user', user)
[ywant, ypwant, wcomm, user, iwsav, rwsav, ifail] = nag_ode_ivp_rkts_interp(n, twant, ideriv, nwant, f, wcomm, iwsav, rwsav, 'lwcomm', lwcomm, 'user', user)


nag_ode_ivp_rkts_interp (d02ps) and its associated functions (nag_ode_ivp_rkts_onestep (d02pf), nag_ode_ivp_rkts_setup (d02pq), nag_ode_ivp_rkts_reset_tend (d02pr), nag_ode_ivp_rkts_diag (d02pt) and nag_ode_ivp_rkts_errass (d02pu)) 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_rkts_onestep (d02pf) computes the solution at the end of an integration step. Using the information computed on that step nag_ode_ivp_rkts_interp (d02ps) computes the solution by interpolation at any point on that step. It cannot be used if method=3 or -3 was specified in the call to setup function nag_ode_ivp_rkts_setup (d02pq).


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


Compulsory Input Parameters

1:     n int64int32nag_int scalar
n, the number of ordinary differential equations in the system to be solved by the integration function.
Constraint: n1.
2:     twant – double scalar
t, the value of the independent variable where a solution is desired.
3:     ideriv int64int32nag_int scalar
Determines whether the solution and/or its first derivative are to be computed
compute approximate solution.
compute approximate first derivative.
compute approximate solution and first derivative.
Constraint: ideriv=0, 1 or 2.
4:     nwant int64int32nag_int scalar
The number of components of the solution to be computed. The first nwant components are evaluated.
Constraint: 1nwantn.
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_rkts_onestep (d02pf).
[yp, user] = f(t, n, y, user)

Input Parameters

1:     t – double scalar
t, the current value of the independent variable.
2:     n int64int32nag_int scalar
n, the number of ordinary differential equations in the system to be solved.
3:     yn – double array
The current values of the dependent variables, yi, for i=1,2,,n.
4:     user – Any MATLAB object
f is called from nag_ode_ivp_rkts_interp (d02ps) with the object supplied to nag_ode_ivp_rkts_interp (d02ps).

Output Parameters

1:     ypn – double array
The values of fi, for i=1,2,,n.
2:     user – Any MATLAB object
6:     wcommlwcomm – double array
This array stores information that can be utilized on subsequent calls to nag_ode_ivp_rkts_interp (d02ps).
7:     iwsav130 int64int32nag_int array
8:     rwsav32×n+350 – double array
These must be the same arrays supplied in a previous call nag_ode_ivp_rkts_onestep (d02pf). They must remain unchanged between calls.

Optional Input Parameters

1:     lwcomm int64int32nag_int scalar
Default: the dimension of the array wcomm.
Length of wcomm.
If in a previous call to nag_ode_ivp_rkts_setup (d02pq):
  • method=1 or -1 then lwcomm must be at least 1.
  • method=2 or -2 then lwcomm must be at least n+maxn,5×nwant.
  • method=3 or -3 then wcomm and lwcomm are not referenced.
2:     user – Any MATLAB object
user is not used by nag_ode_ivp_rkts_interp (d02ps), but is passed to f. Note that for large objects it may be more efficient to use a global variable which is accessible from the m-files than to use user.

Output Parameters

1:     ywantnwant – double array
An approximation to the first nwant components of the solution at twant if ideriv=0 or 2. Otherwise ywant is not defined.
2:     ypwantnwant – double array
An approximation to the first nwant components of the first derivative at twant if ideriv=1 or 2. Otherwise ypwant is not defined.
3:     wcommlwcomm – double array
Communication array, used to store information between calls to nag_ode_ivp_rkts_interp (d02ps).
4:     user – Any MATLAB object
5:     iwsav130 int64int32nag_int array
6:     rwsav32×n+350 – double array
Information about the integration for use on subsequent calls to nag_ode_ivp_rkts_onestep (d02pf), nag_ode_ivp_rkts_interp (d02ps) or other associated functions.
7:     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:
Constraint: 1nwantn.
Constraint: for method=-1 or 1, lwcomm1.
Constraint: for method=-2 or 2, lwcommn+maxn,5×nwant.
Constraint: ideriv=0, 1 or 2.
method=-3 or 3 in setup, but interpolation is not available for this method. Either use method=-2 or 2 in setup or use reset function to force the integrator to step to particular points.
On entry, a previous call to the setup function has not been made or the communication arrays have become corrupted, or a catastrophic error has already been detected elsewhere.
You cannot continue integrating the problem.
On entry, n=_, but the value passed to the setup function was .
You cannot call this function after the integrator has returned an error.
You cannot call this function before you have called the step integrator.
You cannot call this function when you have specified, in the setup function, that the range integrator will be used.
An unexpected error has been triggered by this routine. Please contact NAG.
Your licence key may have expired or may not have been installed correctly.
Dynamic memory allocation failed.


The computed values will be of a similar accuracy to that computed by nag_ode_ivp_rkts_onestep (d02pf).

Further Comments



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_rkts_onestep (d02pf) is used to integrate the problem one step at a time and nag_ode_ivp_rkts_interp (d02ps) 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 low order Runge–Kutta method (method=-1) is also used with tolerances tol=1.0e−4 and tol=1.0e−5 in turn so that solutions may be compared.
function d02ps_example

fprintf('d02ps 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);
nwant  = int64(1);
tol    =  1.0E-2;
npts   = 32;
wcomm  = zeros(n+5*nwant, 1);
twant  = zeros(npts+1, 1);
ywant  = zeros(npts, nwant);
ypwant = zeros(npts, nwant);

% 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
  [iwsav, rwsav, ifail] = d02pq(tstart, tend, yinit, tol, thresh, method);

  % 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;

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

  % Integrate one step from t, updating t.
    [t, y, yp, user, iwsav, rwsav, ifail] = d02pf(@f, n, iwsav, rwsav);

    % Interpolate at required additional points up to t
    while twant(j+1) < t
      % Interpolate and print solution at t = twant.
      ideriv = int64(2);
      [ywant(j, :), ypwant(j, :), wcomm, user, iwsav, rwsav, ifail] = ...
        d02ps(n, twant(j), ideriv, nwant, @f, wcomm, iwsav, rwsav);

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

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


% Plot results
fig1 = figure;
title('Simple Sine Solution, tol=1e-4');
hold on;
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, user] = f(t, n, y, user)
  yp = [y(2); -y(1)];
d02ps example results

Calculation with TOL =  1.0e-03

     Number of evaluations of f = 152

Calculation with TOL =  1.0e-04

     Number of evaluations of f = 231

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