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)

Purpose

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).

Syntax

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

Description

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.

References

See the D02M–N sub-chapter Introduction.

Parameters

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: 1mneq1mneq.
3:     neq – int64int32nag_int scalar
The value used for the parameter neq when calling the integrator.
Constraint: 1neqldysav1neqldysav.
4:     ysav(ldysav,sdysav) – double array
ldysav, the first dimension of the array, must satisfy the constraint ldysav1ldysav1.
The values provided in the array ysav on return from the integrator.
5:     rwork(50 + 4 × neq50+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 parameter sdysav when calling the integrator.

Input Parameters Omitted from the MATLAB Interface

ldysav

Output Parameters

1:     sol(m) – double array
The calculated value of the solution at tsol.
2:     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,m < 1m<1,
orldysav < 1ldysav<1,
orneq < 1neq<1,
orm > neqm>neq,
orneq > ldysavneq>ldysav.
  ifail = 2ifail=2
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 = 3ifail=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.

Accuracy

The solution values returned will be of a similar accuracy to those computed by the integrator.

Further Comments

None.

Example

function nag_ode_ivp_stiff_interp_example
%
% Integrate to tout by overshooting tout in one step mode (itask=2)
% using B.D.F formulae with a functional iteration method.
% Default values for the array const are used. Employ vector
% tolerances and the Jacobian is evaluated internally, if necessary.
% monitr subroutine replaced by NAG dummy routine nag_ode_ivp_stiff_exp_fulljac_dummy_monit.
% Interpolation outside nag_ode_ivp_stiff_imp_fulljac using nag_ode_ivp_stiff_interp.
%
neq = 3;
atol = zeros(1, neq);
rtol = zeros(1, neq);
sol = zeros(1, neq);
nwkjac = int64(neq * (neq + 1));
neq = int64(neq);
wkjac = zeros(1, nwkjac);
y = zeros(1, neq);
ydot = zeros(1, neq);
algequ = zeros(1, neq);
lderiv = zeros(1, 2);
inform = zeros(23, 1, 'int64');
maxord = int64(5);
sdysav = int64(6);
ysav = zeros(neq, sdysav);
petzld = false;
const = zeros(6,1);
h0 = 0.0e0;
hmax = 10.0e0;
hmin = 1.0e-10;
tcrit = 0.0e0;
maxstp = int64(200);
mxhnil = int64(5);
rwork = zeros(62,1);
t = 0.0e0;
tout = 0.1e0;
itask = int64(2);
y = [1; 0; 0];
lderiv = [false; false];
itol = int64(4);
rtol = [1.0e-4; 1.0e-3; 1.0e-4];
atol = [1.0e-7; 1.0e-8; 1.0e-7];
sol = zeros(6, neq);
sol(1,:) = y;
fprintf('\nnag_ode_ivp_stiff_interp example program results\n');
[const, rwork, ifail] = ...
    nag_ode_ivp_stiff_bdf(neq, sdysav, maxord, 'functional-iteration', ...
                                              petzld, const, tcrit, hmin, hmax, h0, maxstp, ...
                                              mxhnil, 'average-l2', rwork);
if ifail == int64(0)
  % Linear algebra setup required (in case functional iteration
  % encounters any difficulty).
  [rwork, ifail] = ...
    nag_ode_ivp_stiff_fulljac_setup(neq, neq, 'numerical', nwkjac, rwork);
  xout = 0.02e0;
  iout = 1;
  fprintf('\n');
  fprintf('     x           y(1)          y(2)          y(3)\n');
  fprintf('%8.3f %13.5f %13.5f %13.5f\n', t, y(1:double(neq)));
  itrace = int64(0);
  tcur = 1;
  hu = 0;
  while (tcur-hu >= xout) || (xout > tcur) && (ifail == 0)
    [t, tout, y, ydot, rwork, inform, ysav, wkjac, lderiv, ifail] = ...
       nag_ode_ivp_stiff_imp_fulljac(t, tout, y, ydot, rwork, rtol, atol, itol, inform, ...
       @resid, ysav, 'nag_ode_ivp_stiff_imp_fulljac_dummy_jac', wkjac, 'nag_ode_ivp_stiff_exp_fulljac_dummy_monit', lderiv, itask, ...
        itrace);
    [hu, h, tcur, tolsf, nst, nre, nje, nqu, nq, niter, imxer, algequ, ifail] = ...
       nag_ode_ivp_stiff_integ_diag(neq, neq, rwork, inform);
    while (iout < 6) && (tcur-hu < xout) && (xout <= tcur) && (ifail == 0)
      [sol(iout+1,:), ifail] = nag_ode_ivp_stiff_interp(xout, neq, neq, ysav, rwork);
      fprintf('%8.3f %13.5f %13.5f %13.5f\n', xout, sol(iout+1,1:double(neq)));
      iout = iout + 1;
      xout = iout * 0.02e0;
    end
  end
end


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);
  end
 

nag_ode_ivp_stiff_interp example program 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

function d02mz_example
%
% Integrate to tout by overshooting tout in one step mode (itask=2)
% using B.D.F formulae with a functional iteration method.
% Default values for the array const are used. Employ vector
% tolerances and the Jacobian is evaluated internally, if necessary.
% monitr subroutine replaced by NAG dummy routine d02nby.
% Interpolation outside d02ng using d02mz.
%
neq = 3;
atol = zeros(1, neq);
rtol = zeros(1, neq);
sol = zeros(1, neq);
nwkjac = int64(neq * (neq + 1));
neq = int64(neq);
wkjac = zeros(1, nwkjac);
y = zeros(1, neq);
ydot = zeros(1, neq);
algequ = zeros(1, neq);
lderiv = zeros(1, 2);
inform = zeros(23, 1, 'int64');
maxord = int64(5);
sdysav = int64(6);
ysav = zeros(neq, sdysav);
petzld = false;
const = zeros(6,1);
h0 = 0.0e0;
hmax = 10.0e0;
hmin = 1.0e-10;
tcrit = 0.0e0;
maxstp = int64(200);
mxhnil = int64(5);
rwork = zeros(62,1);
t = 0.0e0;
tout = 0.1e0;
itask = int64(2);
y = [1; 0; 0];
lderiv = [false; false];
itol = int64(4);
rtol = [1.0e-4; 1.0e-3; 1.0e-4];
atol = [1.0e-7; 1.0e-8; 1.0e-7];
sol = zeros(6, neq);
sol(1,:) = y;
fprintf('\nd02mz example program results\n');
[const, rwork, ifail] = d02nv(neq, sdysav, maxord, 'functional-iteration', ...
                              petzld, const, tcrit, hmin, hmax, h0, maxstp, ...
                              mxhnil, 'average-l2', rwork);
if ifail == int64(0)
  % Linear algebra setup required (in case functional iteration
  % encounters any difficulty).
  [rwork, ifail] = d02ns(neq, neq, 'numerical', nwkjac, rwork);
  xout = 0.02e0;
  iout = 1;
  fprintf('\n');
  fprintf('     x           y(1)          y(2)          y(3)\n');
  fprintf('%8.3f %13.5f %13.5f %13.5f\n', t, y(1:double(neq)));
  itrace = int64(0);
  tcur = 1;
  hu = 0;
  while (tcur-hu >= xout) || (xout > tcur) && (ifail == 0)
    [t, tout, y, ydot, rwork, inform, ysav, wkjac, lderiv, ifail] = ...
       d02ng(t, tout, y, ydot, rwork, rtol, atol, itol, inform, ...
       @resid, ysav, 'd02ngz', wkjac, 'd02nby', lderiv, itask, ...
        itrace);
    [hu, h, tcur, tolsf, nst, nre, nje, nqu, nq, niter, imxer, algequ, ifail] = ...
       d02ny(neq, neq, rwork, inform);
    while (iout < 6) && (tcur-hu < xout) && (xout <= tcur) && (ifail == 0)
      [sol(iout+1,:), ifail] = d02mz(xout, neq, neq, ysav, rwork);
      fprintf('%8.3f %13.5f %13.5f %13.5f\n', xout, sol(iout+1,1:double(neq)));
      iout = iout + 1;
      xout = iout * 0.02e0;
    end
  end
end


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);
  end
 

d02mz example program 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–2013