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_errest (d02za)

Purpose

nag_ode_ivp_stiff_errest (d02za) calculates the weighted norm of the local error estimate from inside a monitr called from an integrator in sub-chapter D02M–N (e.g., see nag_ode_ivp_stiff_exp_fulljac (d02nb)).

Syntax

[result, ifail] = d02za(v, w, 'neq', neq)
[result, ifail] = nag_ode_ivp_stiff_errest(v, w, 'neq', neq)

Description

nag_ode_ivp_stiff_errest (d02za) is for use with the forward communication integrators nag_ode_ivp_stiff_exp_fulljac (d02nb), nag_ode_ivp_stiff_exp_bandjac (d02nc), nag_ode_ivp_stiff_exp_sparjac (d02nd), nag_ode_ivp_stiff_imp_fulljac (d02ng), nag_ode_ivp_stiff_imp_bandjac (d02nh) and nag_ode_ivp_stiff_imp_sparjac (d02nj) and the reverse communication integrators nag_ode_ivp_stiff_exp_revcom (d02nm) and nag_ode_ivp_stiff_imp_revcom (d02nn). It must be used only inside monitr (if this option is selected) for the forward communication functions or on the equivalent return for the reverse communication functions. It may be used to evaluate the norm of the scaled local error estimate, vv, where the weights used are contained in ww and the norm used is as defined by an earlier call to the integrator setup function (nag_ode_ivp_stiff_dassl (d02mv), nag_ode_ivp_stiff_bdf (d02nv) or nag_ode_ivp_stiff_blend (d02nw)). Its use is described under the description of monitr in the specifications for the forward communication integrators mentioned above.

References

None.

Parameters

Compulsory Input Parameters

1:     v(neq) – double array
The vector, the weighted norm of which is to be evaluated by nag_ode_ivp_stiff_errest (d02za). v is calculated internally by the integrator being used.
2:     w(neq) – double array
The weights, calculated internally by the integrator, to be used in the norm evaluation.

Optional Input Parameters

1:     neq – int64int32nag_int scalar
Default: The dimension of the arrays v, w. (An error is raised if these dimensions are not equal.)
The number of differential equations, as defined for the integrator being used.

Input Parameters Omitted from the MATLAB Interface

None.

Output Parameters

1:     result – double scalar
The result of the function.
2:     ifail – int64int32nag_int scalar
ifail = 0ifail=0 unless the function detects an error (see [Error Indicators and Warnings]).

Error Indicators and Warnings

Note: nag_ode_ivp_stiff_errest (d02za) may return useful information for one or more of the following detected errors or 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.

W ifail = 1ifail=1
The value of the norm would either overflow or is close to overflowing. A value close to the square root of the largest number on the computer is returned.

Accuracy

The result is calculated close to machine precision except in the case when the function exits with ifail = 1ifail=1.

Further Comments

nag_ode_ivp_stiff_errest (d02za) should only be used within monitr associated with the integrators in sub-chapter D02M–N (e.g., see nag_ode_ivp_stiff_exp_fulljac (d02nb)). Its use and only valid calling sequence are fully documented in the description of monitr in the function documents for the integrators.

Example

function nag_ode_ivp_stiff_errest_example
t = 0;
tout = 10;
y = [1; 0; 0];
rwork = zeros(62,1);
rtol = [0.0001];
atol = [1e-07];
itol = int64(1);
inform = zeros(23, 1, 'int64');
ysave = zeros(3, 6);
wkjac = zeros(12, 1);
itask = int64(4);
itrace = int64(0);
[const, rwork, ifail] = ...
  nag_ode_ivp_stiff_bdf(int64(3), int64(6), int64(5), 'Newton', false, zeros(6), ...
     10, 1e-10, 10, 0, int64(200), int64(5), 'Average-L2', rwork);

[rwork, ifail] = ...
    nag_ode_ivp_stiff_fulljac_setup(int64(3), int64(3), 'Analytical', int64(12), rwork);

fprintf('\n      X        Y(1)         Y(2)         Y(3)\n');
fprintf(' %10.6f%13.7f%13.7f%13.7f\n', t, y(1), y(2), y(3));

[tOut, yOut, ydot, rworkOut, informOut, ysaveOut, wkjacOut, ifail] = ...
   nag_ode_ivp_stiff_exp_fulljac(t, tout, y, rwork, rtol, atol, itol, inform, ...
   @fcn, ysave, @jac, wkjac, @monitr, itask, itrace);
 tOut, yOut, ydot, informOut, ifail


function [f, ires] = fcn(neq, t, y, ires)
% Evaluate derivative vector.
f = zeros(3,1);
f(1) = -0.04d0*y(1) + 1.0d4*y(2)*y(3);
f(2) = 0.04d0*y(1) - 1.0d4*y(2)*y(3) - 3.0d7*y(2)*y(2);
f(3) = 3.0d7*y(2)*y(2);
function p = jac(neq, t, y, h, d, p)
% Evaluate the Jacobian.
p = zeros(neq, neq);
hxd = h*d;
p(1,1) = 1.0d0 - hxd*(-0.04d0);
p(1,2) = -hxd*(1.0d4*y(3));
p(1,3) = -hxd*(1.0d4*y(2));
p(2,1) = -hxd*(0.04d0);
p(2,2) = 1.0d0 - hxd*(-1.0d4*y(3)-6.0d7*y(2));
p(2,3) = -hxd*(-1.0d4*y(2));
p(3,2) = -hxd*(6.0d7*y(2));
p(3,3) = 1.0d0 - hxd*(0.0d0);
function [hnext, y, imon, inln, hmin, hmax] = ...
    monitr(neq, neqmax, t, hlast, hnext, y, ydot, ysave, r, acor, imon, hmin, hmax, nqu)

  inln=int64(0);
  if (imon == int64(1))
    [errloc, ifail] = nag_ode_ivp_stiff_errest(acor(1, 2), acor(1, 1));
    if (ifail ~= int64(0))
      imon = int64(-2);
    elseif (errloc > 5)
      fprintf(' %10.6f%13.7f%13.7f%13.7f ** WARNING scaled local error = %13.5f\n', t, y(1), y(2), y(3), errloc);
    else
      fprintf(' %10.6f%13.7f%13.7f%13.7f\n', t, y(1), y(2), y(3));
    end
  end
 

      X        Y(1)         Y(2)         Y(3)
   0.000000    1.0000000    0.0000000    0.0000000
   0.000108    0.9999957    0.0000042    0.0000001
   0.000215    0.9999914    0.0000083    0.0000003
   0.000305    0.9999878    0.0000115    0.0000007
   0.000394    0.9999842    0.0000145    0.0000013
   0.000550    0.9999780    0.0000193    0.0000027
   0.000707    0.9999717    0.0000234    0.0000049
   0.000863    0.9999655    0.0000267    0.0000078
   0.001134    0.9999546    0.0000308    0.0000145
   0.001338    0.9999465    0.0000328    0.0000207
   0.001541    0.9999384    0.0000342    0.0000274
   0.001744    0.9999302    0.0000351    0.0000347
   0.002024    0.9999190    0.0000357    0.0000453
   0.002304    0.9999078    0.0000360    0.0000561
   0.002584    0.9998966    0.0000362    0.0000671
   0.002865    0.9998855    0.0000363    0.0000782
   0.003252    0.9998700    0.0000364    0.0000936
   0.003639    0.9998545    0.0000365    0.0001090
   0.004026    0.9998390    0.0000365    0.0001245
   0.005346    0.9997864    0.0000365    0.0001772
   0.006665    0.9997337    0.0000365    0.0002298
   0.011496    0.9995413    0.0000364    0.0004223
   0.016328    0.9993492    0.0000364    0.0006144
   0.027384    0.9989107    0.0000363    0.0010529
   0.038440    0.9984742    0.0000362    0.0014896
   0.049496    0.9980395    0.0000362    0.0019243
   0.090362    0.9964493    0.0000359    0.0035149
   0.131228    0.9948843    0.0000356    0.0050802
   0.172093    0.9933438    0.0000353    0.0066209
   0.256544    0.9902354    0.0000348    0.0097299
   0.340995    0.9872231    0.0000342    0.0127427
   0.425446    0.9843009    0.0000337    0.0156653
   0.509897    0.9814639    0.0000332    0.0185029
   0.647901    0.9769992    0.0000325    0.0229684
   0.785904    0.9727312    0.0000318    0.0272371
   0.923908    0.9686436    0.0000311    0.0313253
   1.061912    0.9647221    0.0000305    0.0352474
   1.315223    0.9579148    0.0000294    0.0420558
   1.568533    0.9515559    0.0000285    0.0484156
   1.821844    0.9455915    0.0000276    0.0543809
   2.075154    0.9399766    0.0000268    0.0599966
   2.328465    0.9346727    0.0000261    0.0653013
   2.707880    0.9272386    0.0000251    0.0727363
   3.087296    0.9203375    0.0000242    0.0796383
   3.466712    0.9138972    0.0000234    0.0860794
   3.846128    0.9078595    0.0000227    0.0921178
   4.225543    0.9021762    0.0000220    0.0978018
   4.812256    0.8939921    0.0000211    0.1059867
   5.398969    0.8864397    0.0000203    0.1135400
   5.985682    0.8794262    0.0000196    0.1205542
   6.572395    0.8728779    0.0000190    0.1271031
   7.159109    0.8667346    0.0000184    0.1332470
   8.060884    0.8579697    0.0000176    0.1420126
   8.962660    0.8499047    0.0000169    0.1500784
   9.481330    0.8455416    0.0000166    0.1544418
  10.000000    0.8413577    0.0000162    0.1586261

tOut =

    10


yOut =

    0.8414
    0.0000
    0.1586


ydot =

   -0.0079
   -0.0000
    0.0079


informOut =

                   55
                   81
                   17
                    3
                    3
                   79
                    0
                    3
                    0
                    0
                    0
                    0
                    0
                    0
                    0
                    0
                    0
                    0
                    0
                    0
                    0
                    0
                    0


ifail =

                    0


function d02za_example
t = 0;
tout = 10;
y = [1; 0; 0];
rwork = zeros(62,1);
rtol = [0.0001];
atol = [1e-07];
itol = int64(1);
inform = zeros(23, 1, 'int64');
ysave = zeros(3, 6);
wkjac = zeros(12, 1);
itask = int64(4);
itrace = int64(0);
[const, rwork, ifail] = ...
  d02nv(int64(3), int64(6), int64(5), 'Newton', false, zeros(6), ...
     10, 1e-10, 10, 0, int64(200), int64(5), 'Average-L2', rwork);

[rwork, ifail] = d02ns(int64(3), int64(3), 'Analytical', int64(12), rwork);

fprintf('\n      X        Y(1)         Y(2)         Y(3)\n');
fprintf(' %10.6f%13.7f%13.7f%13.7f\n', t, y(1), y(2), y(3));

[tOut, yOut, ydot, rworkOut, informOut, ysaveOut, wkjacOut, ifail] = ...
   d02nb(t, tout, y, rwork, rtol, atol, itol, inform, ...
   @fcn, ysave, @jac, wkjac, @monitr, itask, itrace);
 tOut, yOut, ydot, informOut, ifail


function [f, ires] = fcn(neq, t, y, ires)
% Evaluate derivative vector.
f = zeros(3,1);
f(1) = -0.04d0*y(1) + 1.0d4*y(2)*y(3);
f(2) = 0.04d0*y(1) - 1.0d4*y(2)*y(3) - 3.0d7*y(2)*y(2);
f(3) = 3.0d7*y(2)*y(2);
function p = jac(neq, t, y, h, d, p)
% Evaluate the Jacobian.
p = zeros(neq, neq);
hxd = h*d;
p(1,1) = 1.0d0 - hxd*(-0.04d0);
p(1,2) = -hxd*(1.0d4*y(3));
p(1,3) = -hxd*(1.0d4*y(2));
p(2,1) = -hxd*(0.04d0);
p(2,2) = 1.0d0 - hxd*(-1.0d4*y(3)-6.0d7*y(2));
p(2,3) = -hxd*(-1.0d4*y(2));
p(3,2) = -hxd*(6.0d7*y(2));
p(3,3) = 1.0d0 - hxd*(0.0d0);
function [hnext, y, imon, inln, hmin, hmax] = ...
    monitr(neq, neqmax, t, hlast, hnext, y, ydot, ysave, r, acor, imon, hmin, hmax, nqu)

  inln=int64(0);
  if (imon == int64(1))
    [errloc, ifail] = d02za(acor(1, 2), acor(1, 1));
    if (ifail ~= int64(0))
      imon = int64(-2);
    elseif (errloc > 5)
      fprintf(' %10.6f%13.7f%13.7f%13.7f ** WARNING scaled local error = %13.5f\n',...
      t, y(1), y(2), y(3), errloc);
    else
      fprintf(' %10.6f%13.7f%13.7f%13.7f\n', t, y(1), y(2), y(3));
    end
  end
 

      X        Y(1)         Y(2)         Y(3)
   0.000000    1.0000000    0.0000000    0.0000000
   0.000108    0.9999957    0.0000042    0.0000001
   0.000215    0.9999914    0.0000083    0.0000003
   0.000305    0.9999878    0.0000115    0.0000007
   0.000394    0.9999842    0.0000145    0.0000013
   0.000550    0.9999780    0.0000193    0.0000027
   0.000707    0.9999717    0.0000234    0.0000049
   0.000863    0.9999655    0.0000267    0.0000078
   0.001134    0.9999546    0.0000308    0.0000145
   0.001338    0.9999465    0.0000328    0.0000207
   0.001541    0.9999384    0.0000342    0.0000274
   0.001744    0.9999302    0.0000351    0.0000347
   0.002024    0.9999190    0.0000357    0.0000453
   0.002304    0.9999078    0.0000360    0.0000561
   0.002584    0.9998966    0.0000362    0.0000671
   0.002865    0.9998855    0.0000363    0.0000782
   0.003252    0.9998700    0.0000364    0.0000936
   0.003639    0.9998545    0.0000365    0.0001090
   0.004026    0.9998390    0.0000365    0.0001245
   0.005346    0.9997864    0.0000365    0.0001772
   0.006665    0.9997337    0.0000365    0.0002298
   0.011496    0.9995413    0.0000364    0.0004223
   0.016328    0.9993492    0.0000364    0.0006144
   0.027384    0.9989107    0.0000363    0.0010529
   0.038440    0.9984742    0.0000362    0.0014896
   0.049496    0.9980395    0.0000362    0.0019243
   0.090362    0.9964493    0.0000359    0.0035149
   0.131228    0.9948843    0.0000356    0.0050802
   0.172093    0.9933438    0.0000353    0.0066209
   0.256544    0.9902354    0.0000348    0.0097299
   0.340995    0.9872231    0.0000342    0.0127427
   0.425446    0.9843009    0.0000337    0.0156653
   0.509897    0.9814639    0.0000332    0.0185029
   0.647901    0.9769992    0.0000325    0.0229684
   0.785904    0.9727312    0.0000318    0.0272371
   0.923908    0.9686436    0.0000311    0.0313253
   1.061912    0.9647221    0.0000305    0.0352474
   1.315223    0.9579148    0.0000294    0.0420558
   1.568533    0.9515559    0.0000285    0.0484156
   1.821844    0.9455915    0.0000276    0.0543809
   2.075154    0.9399766    0.0000268    0.0599966
   2.328465    0.9346727    0.0000261    0.0653013
   2.707880    0.9272386    0.0000251    0.0727363
   3.087296    0.9203375    0.0000242    0.0796383
   3.466712    0.9138972    0.0000234    0.0860794
   3.846128    0.9078595    0.0000227    0.0921178
   4.225543    0.9021762    0.0000220    0.0978018
   4.812256    0.8939921    0.0000211    0.1059867
   5.398969    0.8864397    0.0000203    0.1135400
   5.985682    0.8794262    0.0000196    0.1205542
   6.572395    0.8728779    0.0000190    0.1271031
   7.159109    0.8667346    0.0000184    0.1332470
   8.060884    0.8579697    0.0000176    0.1420126
   8.962660    0.8499047    0.0000169    0.1500784
   9.481330    0.8455416    0.0000166    0.1544418
  10.000000    0.8413577    0.0000162    0.1586261

tOut =

    10


yOut =

    0.8414
    0.0000
    0.1586


ydot =

   -0.0079
   -0.0000
    0.0079


informOut =

                   55
                   81
                   17
                    3
                    3
                   79
                    0
                    3
                    0
                    0
                    0
                    0
                    0
                    0
                    0
                    0
                    0
                    0
                    0
                    0
                    0
                    0
                    0


ifail =

                    0



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