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_stiff_interp (d02mz)


    1  Purpose
    2  Syntax
    7  Accuracy
    9  Example


nag_ode_ivp_stiff_interp (d02mz) interpolates components of the solution of a system of first-order differential equations from information provided by those integrators in Sub-chapter D02M–N using methods set up by calls to nag_ode_ivp_stiff_dassl (d02mv), nag_ode_ivp_stiff_bdf (d02nv) or nag_ode_ivp_stiff_blend (d02nw).


[sol, ifail] = d02mz(tsol, m, neq, ysav, rwork, 'sdysav', sdysav)
[sol, ifail] = nag_ode_ivp_stiff_interp(tsol, m, neq, ysav, rwork, 'sdysav', sdysav)


nag_ode_ivp_stiff_interp (d02mz) evaluates the first m components of the solution of a system of ordinary differential equations at any point using natural polynomial interpolation based on information generated by the integrator. This information must be passed unchanged to nag_ode_ivp_stiff_interp (d02mz). nag_ode_ivp_stiff_interp (d02mz) should not normally be used to extrapolate outside the range of values obtained from the above function.


See the D02M–N Sub-chapter Introduction.


Compulsory Input Parameters

1:     tsol – double scalar
The point at which the first m components of the solution are to be evaluated. tsol should not normally be an extrapolation point. Extrapolation is permitted but not recommended.
2:     m int64int32nag_int scalar
The number of components of the solution whose values are required.
Constraint: 1mneq.
3:     neq int64int32nag_int scalar
The value used for the argument neq when calling the integrator.
Constraint: 1neqldysav.
4:     ysavldysavsdysav – double array
ldysav, the first dimension of the array, must satisfy the constraint ldysav1.
The values provided in the array ysav on return from the integrator.
5:     rwork50+4×neq – double array
The values provided in the array rwork on return from the integrator.

Optional Input Parameters

1:     sdysav int64int32nag_int scalar
Default: the second dimension of the array ysav.
The value used for the argument sdysav when calling the integrator.

Output Parameters

1:     solm – double array
The calculated value of the solution at tsol.
2:     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:

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

On entry,m<1,
On entry, when accessing an element of the array rwork an unexpected quantity was found. You have not passed the correct array to nag_ode_ivp_stiff_interp (d02mz) or has overwritten elements of this array.
W  ifail=3
On entry, nag_ode_ivp_stiff_interp (d02mz) has been called for extrapolation. Before returning with this error exit, the value of the solution at tsol is calculated and placed in sol.
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 solution values returned will be of a similar accuracy to those computed by the integrator.

Further Comments



This example solves the well-known stiff Robertson problem written in implicit form
r1 = -0.04a + 1.0E4bc - a r2 = 0.04a - 1.0E4bc - 3.0E7b2 - b r3 = 3.0E7b2 - c  
with initial conditions a=1.0 and b=c=0.0 over the range 0,0.1 with vector error control (itol=4), the BDF method (setup function nag_ode_ivp_stiff_bdf (d02nv)) and functional iteration. The Jacobian is calculated numerically if the functional iteration encounters difficulty and the integration is in one-step mode (itask=2), with natural interpolation to calculate the solution at intervals of 0.02 using nag_ode_ivp_stiff_interp (d02mz) externally. nag_ode_ivp_stiff_exp_fulljac_dummy_monit (d02nby) is used for monitr.
function d02mz_example

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

% Integrate to using B.D.F formulae with a functional iteration method
neq    = int64(3);
maxord = int64(5);
sdysav = int64(maxord+1);
petzld = false;
const  = zeros(6,1);
tcrit  = 0;
hmin   = 1e-10;
hmax =   10;
h0     = 0;
maxstp = int64(200);
mxhnil = int64(5);
rwork  = zeros(16+20*neq,1);
[const, rwork, ifail] = d02nv(...
                              neq, sdysav, maxord, 'functional-iteration', ...
                              petzld, const, tcrit, hmin, hmax, h0, maxstp, ...
                              mxhnil, 'average-l2', rwork);

% Use numerical Jacobian if functional iteration encounters any difficulty.
nwkjac = int64(neq * (neq + 1));
[rwork, ifail] = d02ns(neq, neq, 'numerical', nwkjac, rwork);

algequ = zeros(1, neq);
lderiv = zeros(1, 2);

% Variable and array initializations for integrator
% Initial conditions
t      = 0;
tout   = 0.1;
y      = [1; 0; 0];
ydot   = zeros(1, neq);
% vector tolerances
itol   = int64(4);
rtol   = [1.0e-4; 1.0e-3; 1.0e-4];
atol   = [1.0e-7; 1.0e-8; 1.0e-7];
inform = zeros(23, 1, 'int64');
ysav   = zeros(neq, sdysav);
wkjac  = zeros(1, nwkjac);
lderiv = [false; false];
itask  = int64(2);
itrace = int64(0);

% Intialize for solution output
xout = 0.02e0;
fprintf('\n     x           y(1)          y(2)          y(3)\n');
fprintf('%8.3f %13.5f %13.5f %13.5f\n', t, y);

while (t < tout)
   [t, tout, y, ydot, rwork, inform, ysav, wkjac, lderiv, ifail] = ...
             t, tout, y, ydot, rwork, rtol, atol, itol, inform, ...
             @resid, ysav, 'd02ngz', wkjac, 'd02nby', lderiv, itask, ...
   while (xout <= t)
      [sol, ifail] = d02mz(xout, neq, neq, ysav, rwork);
      fprintf('%8.3f %13.5f %13.5f %13.5f\n', xout, sol);
      xout = xout + 0.02e0;

function [r, ires] = resid(neq, t, y, ydot, ires)
  r = -ydot;
  if ires == int64(1)
    r(1) = -0.04e0*y(1) + 1.0e4*y(2)*y(3)                   + r(1);
    r(2) =  0.04e0*y(1) - 1.0e4*y(2)*y(3) - 3.0e7*y(2)*y(2) + r(2);
    r(3) =                                  3.0e7*y(2)*y(2) + r(3);
d02mz example results

     x           y(1)          y(2)          y(3)
   0.000       1.00000       0.00000       0.00000
   0.020       0.99920       0.00004       0.00076
   0.040       0.99841       0.00004       0.00155
   0.060       0.99763       0.00004       0.00234
   0.080       0.99685       0.00004       0.00311
   0.100       0.99608       0.00004       0.00389

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