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_exp_revcom (d02nm)

Purpose

nag_ode_ivp_stiff_exp_revcom (d02nm) is a reverse communication function for integrating stiff systems of explicit ordinary differential equations.

Syntax

[t, y, ydot, rwork, inform, ysav, wkjac, jacpvt, imon, inln, ires, irevcm, ifail] = d02nm(t, tout, y, ydot, rwork, rtol, atol, itol, inform, ysav, wkjac, jacpvt, imon, inln, ires, irevcm, itask, itrace, 'neq', neq, 'sdysav', sdysav, 'nwkjac', nwkjac)
[t, y, ydot, rwork, inform, ysav, wkjac, jacpvt, imon, inln, ires, irevcm, ifail] = nag_ode_ivp_stiff_exp_revcom(t, tout, y, ydot, rwork, rtol, atol, itol, inform, ysav, wkjac, jacpvt, imon, inln, ires, irevcm, itask, itrace, 'neq', neq, 'sdysav', sdysav, 'nwkjac', nwkjac)
Note: the interface to this routine has changed since earlier releases of the toolbox:
Mark 22: njcpvt has been removed from the interface; neq has been made optional
.

Description

nag_ode_ivp_stiff_exp_revcom (d02nm) is a general purpose function for integrating the initial value problem for a stiff system of explicit ordinary differential equations,
y = g(t,y).
y=g(t,y).
An outline of a typical calling program is given below:
% 
% data initialization
% 
call linear algebra setup routine 
call integrator setup routine 
irevcm = int64(0);
[..., irevcm, ...] = d02nm();
while (irevcm > 0)
  if (irevcm == 8)
    supply the jacobian matrix                                  (i)
  elseif (irevcm == 9)
    perform monitoring tasks requested by the user             (ii)
  elseif (irecvm == 1 or irevcm >= 3 and irevcm <= 5)
    evaluate the derivative                                   (iii)
  elseif (irevcm == 10)
    indicates an unsuccessful step 
  end 
  [..., irevcm, ...] = d02nm();
end
% 
% post processing (optional linear algebra diagnostic call 
% (sparse case only), optional integrator diagnostic call) 
%
There are three major operations that may be required of the calling (sub)program on an intermeditate return (irevcm0irevcm0) from nag_ode_ivp_stiff_exp_revcom (d02nm); these are denoted (i), (ii) and (iii) above.
The following sections describe in greater detail exactly what is required of each of these operations.
(i) Supply the Jacobian Matrix
You need only provide this facility if the parameter jceval = 'A'jceval='A' (or jceval = 'F'jceval='F' if using sparse matrix linear algebra) in a call to the linear algebra setup function (see jceval in nag_ode_ivp_stiff_fulljac_setup (d02ns)). If the Jacobian matrix is to be evaluated numerically by the integrator, then the remainder of section (i) can be ignored.
We must define the system of nonlinear equations which is solved internally by the integrator. The time derivative, yy, has the form
y = (yz) / (hd),
y=(y-z)/(hd),
where hh is the current step size and dd is a parameter that depends on the integration method in use. The vector yy is the current solution and the vector zz depends on information from previous time steps. This means that d/(dy) (​ ​) = (hd) d/(dy) (​ ​) d dy (​ ​) = (hd) ddy (​ ​) .
The system of nonlinear equations that is solved has the form
yg(t,y) = 0
y-g(t,y)=0
but is solved in the form
r(t,y) = 0,
r(t,y)=0,
where the function rr is defined by
r(t,y) = (hd)((yz) / (hd)g(t,y)).
r(t,y)=(hd)((y-z)/(hd)-g(t,y)).
It is the Jacobian matrix (r)/(y) r y  that you must supply as follows:
(ri)/(yj) = 1(hd)(gi)/(yj) if ​i = j,
(ri)/(yj) = 0(hd)(gi)/(yj) otherwise,
ri yj =1-(hd) gi yj if ​i=j, ri yj =0-(hd) gi yj otherwise,
where tt, hh and dd are located in rwork(19)rwork19, rwork(16)rwork16 and rwork(20)rwork20 respectively and the array y contains the current values of the dependent variables. Only the nonzero elements of the Jacobian need be set, since the locations where it is to be stored are preset to zero.
Hereafter in this document this operation will be referred to as JAC.
(ii) Perform Tasks Requested by You
This operation is essentially a monitoring function and additionally provides the opportunity of changing the current values of y, HNEXT (the step size that the integrator proposes to take on the next step), HMIN (the minimum step size to be taken on the next step), and HMAX (the maximum step size to be taken on the next step). The scaled local error at the end of a timestep may be obtained by calling double function nag_ode_ivp_stiff_errest (d02za) as follows:
[errloc, ifail] = d02za(rwork(51+neq:51+neq+neq-1), rwork(51:51+neq-1));
% Check ifail before proceeding
The following gives details of the location within the array rwork of variables that may be of interest to you:
Variable Specification Location
tcurr the current value of the independent variable rwork(19)rwork19
hlast last step size successfully used by the integrator rwork(15)rwork15
hnext step size that the integrator proposes to take on the next step rwork(16)rwork16
hmin minimum step size to be taken on the next step rwork(17)rwork17
hmax maximum step size to be taken on the next step rwork(18)rwork18
nqu the order of the integrator used on the last step rwork(10)rwork10
You are advised to consult the description of monitr in nag_ode_ivp_stiff_exp_fulljac (d02nb) for details on what optional input can be made.
If y is changed, then imon must be set to 22 before return to nag_ode_ivp_stiff_exp_revcom (d02nm). If either of the values of HMIN or HMAX are changed, then imon must be set 33 before return to nag_ode_ivp_stiff_exp_revcom (d02nm). If HNEXT is changed, then imon must be set to 44 before return to nag_ode_ivp_stiff_exp_revcom (d02nm).
In addition you can force nag_ode_ivp_stiff_exp_revcom (d02nm) to evaluate the residual vector
yg(t,y)
y-g(t,y)
by setting imon = 0imon=0 and inln = 3inln=3 and then returning to nag_ode_ivp_stiff_exp_revcom (d02nm); on return to this monitoring operation the residual vector will be stored in rwork(50 + 2 × neq + i)rwork50+2×neq+i, for i = 1,2,,neqi=1,2,,neq.
Hereafter in this document this operation will be referred to as MONITR.
(iii) Evaluate the Derivative
This operation must evaluate the derivative vector for the explicit ordinary differential equation system defined by
y = g(t,y),
y=g(t,y),
where tt is located in rwork(19)rwork19.
Hereafter in this document this operation will be referred to as FCN.

References

See the D02M–N sub-chapter Introduction.

Parameters

Note: this function uses reverse communication. Its use involves an initial entry, intermediate exits and re-entries, and a final exit, as indicated by the parameter irevcm. Between intermediate exits and re-entries, all parameters other than ydot, rwork, wkjac, imon, inln and ires must remain unchanged.

Compulsory Input Parameters

1:     t – double scalar
On initial entry: tt, the value of the independent variable. The input value of t is used only on the first call as the initial point of the integration.
2:     tout – double scalar
On initial entry: the next value of tt at which a computed solution is desired. For the initial tt, the input value of tout is used to determine the direction of integration. Integration is permitted in either direction (see also itask).
Constraint: toutttoutt.
3:     y(neq) – double array
neq, the dimension of the array, must satisfy the constraint neq1neq1.
On initial entry: the values of the dependent variables (solution). On the first call the first neq elements of yy must contain the vector of initial values.
4:     ydot(neq) – double array
neq, the dimension of the array, must satisfy the constraint neq1neq1.
On intermediate re-entry: must be set to the derivatives as defined under the description of irevcm.
5:     rwork(50 + 4 × neq50+4×neq) – double array
On initial entry: must be the same array as used by one of the method setup functions nag_ode_ivp_stiff_dassl (d02mv), nag_ode_ivp_stiff_bdf (d02nv) or nag_ode_ivp_stiff_blend (d02nw), and by one of the storage setup functions nag_ode_ivp_stiff_bandjac_setup (d02nt), nag_ode_ivp_stiff_sparjac_setup (d02nu) or nag_ode_ivp_stiff_bdf (d02nv). The contents of rwork must not be changed between any call to a setup function and the first call to nag_ode_ivp_stiff_exp_revcom (d02nm).
On intermediate re-entry: elements of rwork must be set to quantities as defined under the description of irevcm.
6:     rtol( : :) – double array
Note: the dimension of the array rtol must be at least 11 if itol = 1itol=1 or 22, and at least neqneq otherwise.
On initial entry: the relative local error tolerance.
Constraint: rtol(i)0.0rtoli0.0 for all relevant ii (see itol).
7:     atol( : :) – double array
Note: the dimension of the array atol must be at least 11 if itol = 1itol=1 or 33, and at least neqneq otherwise.
On initial entry: the absolute local error tolerance.
Constraint: atol(i)0.0atoli0.0 for all relevant ii (see itol).
8:     itol – int64int32nag_int scalar
On initial entry: a value to indicate the form of the local error test. itol indicates to nag_ode_ivp_stiff_exp_revcom (d02nm) whether to interpret either or both of rtol or atol as a vector or a scalar. The error test to be satisfied is ei / wi < 1.0ei/wi<1.0, where wiwi is defined as follows:
itol rtol atol wiwi
1 scalar scalar rtol(1) × |yi| + atol(1)rtol1×|yi|+atol1
2 scalar vector rtol(1) × |yi| + atol(i)rtol1×|yi|+atoli
3 vector scalar rtol(i) × |yi| + atol(1)rtoli×|yi|+atol1
4 vector vector rtol(i) × |yi| + atol(i)rtoli×|yi|+atoli
eiei is an estimate of the local error in yiyi, computed internally, and the choice of norm to be used is defined by a previous call to an integrator setup function.
Constraint: itol = 1itol=1, 22, 33 or 44.
9:     inform(2323) – int64int32nag_int array
10:   ysav(ldysav,sdysav) – double array
ldysav, the first dimension of the array, must satisfy the constraint ldysavneqldysavneq.
On initial entry: the second dimension of the array ysav as declared in the (sub)program from which nag_ode_ivp_stiff_exp_revcom (d02nm) is called. An appropriate value for sdysav is described in the specifications of the integrator setup functions nag_ode_ivp_stiff_bdf (d02nv) and nag_ode_ivp_stiff_blend (d02nw). This value must be the same as that supplied to the integrator setup function.
11:   wkjac(nwkjac) – double array
On intermediate re-entry: elements of the Jacobian as defined under the description of irevcm. If a numerical Jacobian was requested then wkjac is used for workspace.
12:   jacpvt(njcpvt) – int64int32nag_int array
On initial entry: the dimension of the array jacpvt as declared in the (sub)program from which nag_ode_ivp_stiff_exp_revcom (d02nm) is called. The actual size depends on the linear algebra method used. An appropriate value for njcpvt is described in the specifications of the linear algebra setup functions nag_ode_ivp_stiff_bandjac_setup (d02nt) and nag_ode_ivp_stiff_sparjac_setup (d02nu) for banded and sparse matrix linear algebra respectively. This value must be the same as that supplied to the linear algebra setup function. When full matrix linear algebra is chosen, the array jacpvt is not used and hence njcpvt should be set to 11.
13:   imon – int64int32nag_int scalar
On intermediate re-entry: may be reset to determine subsequent action in nag_ode_ivp_stiff_exp_revcom (d02nm).
imon = 2imon=-2
Integration is to be halted. A return will be made from nag_ode_ivp_stiff_exp_revcom (d02nm) to the calling (sub)program with ifail = 12ifail=12.
imon = 1imon=-1
Allow nag_ode_ivp_stiff_exp_revcom (d02nm) to continue with its own internal strategy. The integrator will try up to three restarts unless imon1imon-1.
imon = 0imon=0
Return to the internal nonlinear equation solver, where the action taken is determined by the value of inln.
imon = 1imon=1
Normal exit to nag_ode_ivp_stiff_exp_revcom (d02nm) to continue integration.
imon = 2imon=2
Restart the integration at the current time point. The integrator will restart from order 11 when this option is used. The solution y, provided by the monitr operation (see Section [Description]), will be used for the initial conditions.
imon = 3imon=3
Try to continue with the same step size and order as was to be used before entering the monitr operation (see Section [Description]). hmin and hmax may be altered if desired.
imon = 4imon=4
Continue the integration but using a new value of hnext and possibly new values of hmin and hmax.
14:   inln – int64int32nag_int scalar
On intermediate re-entry: with imon = 0imon=0 and irevcm = 9irevcm=9, inln specifies the action to be taken by the internal nonlinear equation solver. By setting inln = 3inln=3 and returning to nag_ode_ivp_stiff_exp_revcom (d02nm), the residual vector is evaluated and placed in rwork(50 + 2 × neq + i)rwork50+2×neq+i, for i = 1,2,,neqi=1,2,,neq and then the monitr operation (see Section [Description]) is invoked again. At present this is the only option available: inln must not be set to any other value.
15:   ires – int64int32nag_int scalar
On intermediate re-entry: should be unchanged unless one of the following actions is required of nag_ode_ivp_stiff_exp_revcom (d02nm) in which case ires should be set accordingly.
ires = 2ires=2
Indicates to nag_ode_ivp_stiff_exp_revcom (d02nm) that control should be passed back immediately to the calling (sub)program with the error indicator set to ifail = 11ifail=11.
ires = 3ires=3
Indicates to nag_ode_ivp_stiff_exp_revcom (d02nm) that an error condition has occurred in the solution vector, its time derivative or in the value of tt. The integrator will use a smaller time step to try to avoid this condition. If this is not possible nag_ode_ivp_stiff_exp_revcom (d02nm) returns to the calling (sub)program with the error indicator set to ifail = 7ifail=7.
ires = 4ires=4
Indicates to nag_ode_ivp_stiff_exp_revcom (d02nm) to stop its current operation and to enter the monitr operation (see Section [Description]) immediately.
16:   irevcm – int64int32nag_int scalar
On initial entry: must contain 00.
On intermediate re-entry: should remain unchanged.
Constraint: irevcm = 0irevcm=0, 11, 33, 44, 55, 88, 99 or 1010.
17:   itask – int64int32nag_int scalar
On initial entry: the task to be performed by the integrator.
itask = 1itask=1
Normal computation of output values of y(t)y(t) at t = toutt=tout (by overshooting and interpolating).
itask = 2itask=2
Take one step only and return.
itask = 3itask=3
Stop at the first internal integration point at or beyond t = toutt=tout and return.
itask = 4itask=4
Normal computation of output values of y(t)y(t) at t = toutt=tout but without overshooting t = tcritt=tcrit. tcrit must be specified as an option in one of the integrator setup functions before the first call to the integrator, or specified in the optional input function before a continuation call. tcrit (e.g., see nag_ode_ivp_stiff_bdf (d02nv)) may be equal to or beyond tout, but not before it in the direction of integration.
itask = 5itask=5
Take one step only and return, without passing tcrit (e.g., see nag_ode_ivp_stiff_bdf (d02nv)). tcrit must be specified under itask = 4itask=4.
Constraint: 1itask51itask5.
18:   itrace – int64int32nag_int scalar
On initial entry: the level of output that is printed by the integrator. itrace may take the value 1-1, 00, 11, 22 or 33.
itrace < 1itrace<-1
1-1 is assumed and similarly if itrace > 3itrace>3, then 33 is assumed.
itrace = 1itrace=-1
No output is generated.
itrace = 0itrace=0
Only warning messages are printed on the current error message unit (see nag_file_set_unit_error (x04aa)).
itrace > 0itrace>0
Warning messages are printed as above, and on the current advisory message unit (see nag_file_set_unit_advisory (x04ab)) output is generated which details Jacobian entries, the nonlinear iteration and the time integration. The advisory messages are given in greater detail the larger the value of itrace.

Optional Input Parameters

1:     neq – int64int32nag_int scalar
Default: The dimension of the arrays y, ydot and the first dimension of the array ysav. (An error is raised if these dimensions are not equal.)
On initial entry: the number of differential equations to be solved.
Constraint: neq1neq1.
2:     sdysav – int64int32nag_int scalar
Default: The second dimension of the array ysav.
On initial entry: the second dimension of the array ysav as declared in the (sub)program from which nag_ode_ivp_stiff_exp_revcom (d02nm) is called. An appropriate value for sdysav is described in the specifications of the integrator setup functions nag_ode_ivp_stiff_bdf (d02nv) and nag_ode_ivp_stiff_blend (d02nw). This value must be the same as that supplied to the integrator setup function.
3:     nwkjac – int64int32nag_int scalar
Default: The dimension of the array wkjac.
On initial entry: the dimension of the array wkjac as declared in the (sub)program from which nag_ode_ivp_stiff_exp_revcom (d02nm) is called. The actual size depends on the linear algebra method used. An appropriate value for nwkjac is described in the specifications of the linear algebra setup functions nag_ode_ivp_stiff_fulljac_setup (d02ns), nag_ode_ivp_stiff_bandjac_setup (d02nt) and nag_ode_ivp_stiff_sparjac_setup (d02nu) for full, banded and sparse matrix linear algebra respectively. This value must be the same as that supplied to the linear algebra setup function.

Input Parameters Omitted from the MATLAB Interface

ldysav njcpvt

Output Parameters

1:     t – double scalar
On final exit: the value at which the computed solution yy is returned (usually at tout).
2:     y(neq) – double array
On final exit: the computed solution vector evaluated at t (usually t = toutt=tout).
3:     ydot(neq) – double array
On final exit: the time derivatives yy of the vector yy at the last integration point.
4:     rwork(50 + 4 × neq50+4×neq) – double array
On intermediate exit: contains information for JAC, FCN and MONITR operations as described in Section [Description] and the parameter irevcm.
5:     inform(2323) – int64int32nag_int array
6:     ysav(ldysav,sdysav) – double array
7:     wkjac(nwkjac) – double array
On intermediate exit: the Jacobian is overwritten.
8:     jacpvt(njcpvt) – int64int32nag_int array
9:     imon – int64int32nag_int scalar
On intermediate exit: used to pass information between nag_ode_ivp_stiff_exp_revcom (d02nm) and the MONITR operation (see Section [Description]). With irevcm = 9irevcm=9, imon contains a flag indicating under what circumstances the return from nag_ode_ivp_stiff_exp_revcom (d02nm) occurred:
imon = 2imon=-2
Exit from nag_ode_ivp_stiff_exp_revcom (d02nm) after ires = 4ires=4 caused an early termination (this facility could be used to locate discontinuities).
imon = 1imon=-1
The current step failed repeatedly.
imon = 0imon=0
Exit from nag_ode_ivp_stiff_exp_revcom (d02nm) after a call to the internal nonlinear equation solver.
imon = 1imon=1
The current step was successful.
10:   inln – int64int32nag_int scalar
On intermediate exit: contains a flag indicating the action to be taken, if any, by the internal nonlinear equation solver.
11:   ires – int64int32nag_int scalar
On intermediate exit: with irevcm = 1irevcm=1, 22, 33, 44 or 55, ires contains the value 11.
12:   irevcm – int64int32nag_int scalar
On intermediate exit: indicates what action you must take before re-entering. The possible exit values of irevcm are 11, 33, 44, 55, 88, 99 or 1010, which should be interpreted as follows:
irevcm = 1irevcm=1, 33, 44 and 55
Indicates that an FCN operation (see Section [Description]) is required: y = g(t,y)y=g(t,y) must be supplied, where y(i)yi is located in yiyi, for i = 1,2,,neqi=1,2,,neq.
For irevcm = 1irevcm=1 or 33, yiyi should be placed in location rwork(50 + 2 × neq + i)rwork50+2×neq+i, for i = 1,2,,neqi=1,2,,neq.
For irevcm = 4irevcm=4, yiyi should be placed in location rwork(50 + neq + i)rwork50+neq+i, for i = 1,2,,neqi=1,2,,neq.
For irevcm = 5irevcm=5, yiyi should be placed in location ydot(i)ydoti, for i = 1,2,,neqi=1,2,,neq.
irevcm = 8irevcm=8
Indicates that a JAC operation (see Section [Description]) is required: the Jacobian matrix must be supplied.
If full matrix linear algebra is being used, then the (i,j)(i,j)th element of the Jacobian must be stored in wkjac((j1) × neq + i)wkjac((j-1)×neq+i).
If banded matrix linear algebra is being used then the (i,j)(i,j)th element of the Jacobian
must be stored in wkjac((i1) × mB + k)wkjac((i-1)×mB+k), where mB = mL + mU + 1mB=mL+mU+1 and
k = min (mLi + 1,0) + jk=min(mL-i+1,0)+j; here mLmL and mUmU are the number of subdiagonals and superdiagonals, respectively, in the band.
If sparse matrix linear algebra is being used then nag_ode_ivp_stiff_sparjac_enq (d02nr) must be called to determine which column of the Jacobian is required and where it should be stored.
 [j, iplace] = d02nr(inform); 
will return in j the number of the column of the Jacobian that is required and will set iplace = 1iplace=1 or 22. If iplace = 1iplace=1, then the (i,j)(i,j)th element of the Jacobian must be stored in rwork(50 + 2 × neq + i)rwork50+2×neq+i; otherwise it must be stored in rwork(50 + neq + i)rwork50+neq+i.
irevcm = 9irevcm=9
Indicates that a MONITR operation (see Section [Description]) can be performed.
irevcm = 10irevcm=10
Indicates that the current step was not successful, due to error test failure or convergence test failure. The only information supplied to you on this return is the current value of the independent variable tt, located in rwork(19)rwork19. No values must be changed before re-entering nag_ode_ivp_stiff_exp_revcom (d02nm); this facility enables you to determine the number of unsuccessful steps.
On final exit: irevcm = 0irevcm=0 indicated the user-specified task has been completed or an error has been encountered (see the descriptions for itask and ifail).
13:   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, the integrator detected an illegal input, or that a linear algebra and/or integrator setup function has not been called prior to the call to the integrator. If itrace0itrace0, the form of the error will be detailed on the current error message unit (see nag_file_set_unit_error (x04aa)).
  ifail = 2ifail=2
The maximum number of steps specified has been taken (see the description of optional inputs in the integrator setup functions and the optional input continuation function, nag_ode_ivp_stiff_contin (d02nz)).
W ifail = 3ifail=3
With the given values of rtol and atol no further progress can be made across the integration range from the current point t. The components y(1),y(2),,y(neq)y1,y2,,yneq contain the computed values of the solution at the current point t.
W ifail = 4ifail=4
There were repeated error test failures on an attempted step, before completing the requested task, but the integration was successful as far as t. The problem may have a singularity, or the local error requirements may be inappropriate.
W ifail = 5ifail=5
There were repeated convergence test failures on an attempted step, before completing the requested task, but the integration was successful as far as t. This may be caused by an inaccurate Jacobian matrix or one which is incorrectly computed.
W ifail = 6ifail=6
Some error weight wiwi became zero during the integration (see the description of itol). Pure relative error control (atol(i) = 0.0)(atoli=0.0) was requested on a variable (the iith) which has now vanished. The integration was successful as far as t.
  ifail = 7ifail=7
The FCN operation (see Section [Description]) set the error flag ires = 3ires=3 continually despite repeated attempts by the integrator to avoid this.
  ifail = 8ifail=8
Not used for this integrator.
  ifail = 9ifail=9
A singular Jacobian (r)/(y) r y  has been encountered. This error exit is unlikely to be taken when solving explicit ordinary differential equations. You should check the problem formulation and Jacobian calculation.
  ifail = 10ifail=10
An error occurred during Jacobian formulation or back-substitution (a more detailed error description may be directed to the current error message unit, see nag_file_set_unit_error (x04aa)).
W ifail = 11ifail=11
The FCN operation (see Section [Description]) signalled the integrator to halt the integration and return by setting ires = 2ires=2. Integration was successful as far as t.
W ifail = 12ifail=12
The MONITR operation (see Section [Description]) set imon = 2imon=-2 and so forced a return but the integration was successful as far as t.
W ifail = 13ifail=13
The requested task has been completed, but it is estimated that a small change in rtol and atol is unlikely to produce any change in the computed solution. (Only applies when you are not operating in one step mode, that is when itask2itask2 or 55.)
  ifail = 14ifail=14
The values of rtol and atol are so small that nag_ode_ivp_stiff_exp_revcom (d02nm) is unable to start the integration.

Accuracy

The accuracy of the numerical solution may be controlled by a careful choice of the parameters rtol and atol, and to a much lesser extent by the choice of norm. You are advised to use scalar error control unless the components of the solution are expected to be poorly scaled. For the type of decaying solution typical of many stiff problems, relative error control with a small absolute error threshold will be most appropriate (that is, you are advised to choose itol = 1itol=1 with atol(1)atol1 small but positive).

Further Comments

The cost of computing a solution depends critically on the size of the differential system and to a lesser extent on the degree of stiffness of the problem; also on the type of linear algebra being used. For further details see Section [Further Comments] of the documents for nag_ode_ivp_stiff_exp_fulljac (d02nb) (full matrix), nag_ode_ivp_stiff_exp_bandjac (d02nc) (banded matrix) or nag_ode_ivp_stiff_exp_sparjac (d02nd) (sparse matrix).
In general, you are advised to choose the backward differentiation formula option (setup function nag_ode_ivp_stiff_bdf (d02nv)) but if efficiency is of great importance and especially if it is suspected that (g)/(y) g y  has complex eigenvalues near the imaginary axis for some part of the integration, you should try the BLEND option (setup function nag_ode_ivp_stiff_blend (d02nw)).

Example

This example solves the well-known stiff Robertson problem
a = 0.04a + 1.0e4bc3.0e7b2
b = 0.04a1.0e4bc3.0e7b2
c = 0.04a1.0e4bc3.0e7b2
a = -0.04a+1.0e4bc-3.0e7b2 b = -0.04a-1.0e4bc-3.0e7b2 c = -0.04a-1.0e4bc-3.0e7b2
over the range [0,10][0,10] with initial conditions a = 1.0a=1.0 and b = c = 0.0b=c=0.0 and with scalar error control (itol = 1itol=1). The integration proceeds until tout = 10.0tout=10.0 is passed, providing C1C1 interpolation at intervals of 2.02.0 through a MONITR operation. The integration method used is the BDF method (setup function nag_ode_ivp_stiff_bdf (d02nv)) with a modified Newton method. The Jacobian is a full matrix, which is specified using the setup function nag_ode_ivp_stiff_fulljac_setup (d02ns); this Jacobian is to be calculated numerically.
function nag_ode_ivp_stiff_exp_revcom_example
% Initialize variables and arrays.
neq = 3;
neqmax = neq;
nrw = 50+4*neqmax;
ninf = 23;
nwkjac = neqmax*(neqmax + 1);
njcpvt = 1;
maxord = 5;
sdysav = maxord+1;
maxstp = 200;
mxhnil = 5;
ldysav = neq;

h0 = 0.0;
hmax = 10.0;
hmin = 1.0e-10;
tcrit = 0.0;
petzld = false;

const = zeros(6, 1);
rwork = zeros(nrw, 1);
wkjac = zeros(nwkjac, 1);
ydot = zeros(neq, 1);
ysave = zeros(ldysav, sdysav);
inform = zeros(ninf, 1);
jacpvt = zeros(njcpvt, 1);
algequ = zeros(neq, 1);

fprintf('nag_ode_ivp_stiff_exp_revcom example program results\n');

% Set initial parameter values.  We integrate to tout by overshooting
% (itask = 1) using B.D.F. formulae with a Newton method.  Default values
% for the const array are used.  We employ scalar tolerances and the
% Jacobian is evaluated numerically.  On the reverse communication call
% (equivalent to calling the monitr routine in the forward communication
% version), we carry out interpolation using nag_ode_ivp_stiff_c1_interp.
t = 0;
tout = 10.0;
itask = 1;
iout = 1;
xout = 2.0;
y = [1; 0; 0];
itol = 1;
rtol = [1e-04];
atol = [1e-07];

% nag_ode_ivp_stiff_bdf is a setup routine (specifically for the backward
% differentiation formula (BDF) integrator) to be called prior to
% nag_ode_ivp_stiff_exp_revcom.
[const, rwork, ifail] = nag_ode_ivp_stiff_bdf(int64(neqmax), int64(sdysav), ...
    int64(maxord), 'Newton', petzld, const, tcrit, hmin, ...
    hmax, h0, int64(maxstp), int64(mxhnil), 'Average-L2', rwork);
if ifail ~= 0
    % Unsuccessful call.  Print message and exit.
    error('nag_ode_ivp_stiff_bdf returned with ifail = %1d ',ifail);
end

% nag_ode_ivp_stiff_fulljac_setup is a setup routine (specifically for
% full matrix linear algebra) to be called prior to a routine from
% sub-chapter d02n-m (e.g. nag_ode_ivp_stiff_exp_revcom, as here).  We
% specify (jceval = 'Numerical') that the Jacobian is to be evaluated
% numerically.
[rwork, ifail] = nag_ode_ivp_stiff_fulljac_setup(int64(neq), int64(neqmax), 'Numerical', ...
    int64(nwkjac), rwork);
lacorb = 50 + neq;
lacor1 = lacorb + 1;
lacor2 = lacorb + 2;
lacor3 = lacorb + 3;
lsavrb = lacorb + neq;
lsavr1 = lsavrb + 1;
lsavr2 = lsavrb + 2;
lsavr3 = lsavrb + 3;

% Output header and initial results.
fprintf('\n      x           y(1)          y(2)          y(3)\n');
fprintf(' %8.3f %13.5f %13.5f %13.5f\n', t, y(1), y(2), y(3));

% Prepare to store results for plotting.
ncall = 1;
ykeep = y';
tkeep = [t];

% nag_ode_ivp_stiff_exp_revcom uses reverse communication, involving initial entry, intermediate
% exits and re-entries, and a final exit.  These variables are only used
% on intermediate re-entry, but they need to be initialized (here, to a
% bogus value) before calling the routine for the first time.
irevcm = -999;
imon = -999;
inln = -999;
ires = -999;

% Ask for warning messages only.
itrace = 0;

% The reverse communication process is controlled by the value of irevcm
% (which is 0 for the final exit).
while (irevcm ~= 0)
    [t, y, ydot, rwork, inform, ysave, wkjac, jacpvt, imon, inln, ires, ...
        irevcm, ifail] = nag_ode_ivp_stiff_exp_revcom(t, tout, y, ydot, rwork, rtol, atol, ...
        int64(itol), int64(inform), ysave, wkjac, int64(jacpvt), ...
        int64(imon), int64(inln), int64(ires), int64(irevcm), ...
        int64(itask), int64(itrace), 'neq', int64(neq), 'sdysav', ...
        int64(sdysav), 'nwkjac', int64(nwkjac));

    % Evaluate the derivative, then use the value of irevcm to decide
    % where it is to be placed.
    f1 = -0.04*y(1) + 1.0e4*y(2)*y(3);
    f2 =  0.04*y(1) - 1.0e4*y(2)*y(3) - 3.0e7*y(2)*y(2);
    f3 =  3.0e7*y(2)*y(2);
    if (irevcm == 1 || irevcm == 3 )
        rwork(lsavr1) = f1;
        rwork(lsavr2) = f2;
        rwork(lsavr3) = f3;
    elseif (irevcm == 4)
        rwork(lacor1) = f1;
        rwork(lacor2) = f2;
        rwork(lacor3) = f3;
    elseif (irevcm == 5)
        ydot(1) = f1;
        ydot(2) = f2;
        ydot(3) = f3;
    elseif (irevcm == 9)
        % If the current step was successful, perform monitoring tasks.
        if (imon == 1)

            % Extract useful information about the current step.
            tc = rwork(19);
            hlast = rwork(15);
            hnext = rwork(16);
            nqu = int64(rwork(10));

            % For low values of the independent variable, store the results
            % for plotting.
            if tc < 2.0
                ncall = ncall + 1;
                tkeep(ncall,1) = tc;
                ykeep(ncall,:) = y;
            end

            % xout is the point at which we'll do interpolation on the
            % solution for output.  Check to see if we're near it.
            while (iout < 6 && tc-hlast < xout && xout <= tc)
                [sol, ifail] = nag_ode_ivp_stiff_c1_interp(xout, int64(neq), ysave, ...
                    rwork(lacor1:lacor1+neq-1), tc, nqu, hlast, hnext);
                if (ifail ~= 0)
                    imon = -2;
                else

                    % Output interpolated result, and save it for plotting.
                    fprintf(' %8.3f %13.5f %13.5f %13.5f\n', xout, ...
                        sol(1), sol(2), sol(3));
                    ncall = ncall + 1;
                    tkeep(ncall,1) = xout;
                    ykeep(ncall,:) = sol;

                    iout = iout + 1;
                    xout = double(iout)*2.0;
                end
            end

            % Other values are illegal - e.g. irevcm = 8 asks for the
            % Jacobian to be calculated (analogously to the derivative
            % above), but we've already specified (in the call to nag_ode_ivp_stiff_fulljac_setup)
            % that this is to be evaluated numerically, so this value of
            % irevcm should never occur here.
        elseif irevcm == 2 || irevcm == 6 || irevcm == 7 || irevcm == 8
            fprintf('Warning: illegal value of irevcm = %d in nag_ode_ivp_stiff_exp_revcom\n', ...
                irevcm);
            return
        end
    end
end

% nag_ode_ivp_stiff_integ_diag is an integrator diagnostic routine which can be called
% after nag_ode_ivp_stiff_exp_revcom.  Output results.
[hu, h, tcur, tolsf, nst, nre, nje, nqu, nq, niter, imxer, algequ, ...
    ifail] = nag_ode_ivp_stiff_integ_diag(int64(neq), int64(neqmax), rwork, inform);
fprintf('\n hused = %1.5e   hnext = %1.5e   tcur = %1.5e\n', hu, h, tcur);
fprintf(' nst = %1.0f   nre = %1.0f   nje = %1.0f\n',nst, nre, nje);
fprintf(' nqu = %1.0f   nq = %1.0f   niter = %1.0f\n',nqu, nq ,niter);
fprintf(' Max Err Comp = %1.0f  \n\n', imxer);
% Plot results.
fig = figure('Number', 'off');
display_plot(tkeep, ykeep)

function display_plot(x, y)
% Formatting for title and axis labels.
titleFmt = {'FontName', 'Helvetica', 'FontWeight', 'Bold', 'FontSize', 14};
labFmt = {'FontName', 'Helvetica', 'FontWeight', 'Bold', 'FontSize', 13};
set(gca, 'FontSize', 13); % for legend, axis tick labels, etc.
% Plot one of the curves and then add the other two.
hline1 = semilogx(x, y(:,1));
hold on
[haxes, hline2, hline3] = plotyy(x, y(:,3),x, y(:,2), ...
    'semilogx', 'semilogx');
% Set the axis limits and the tick specifications to beautify the plot.
set(haxes(1), 'YLim', [-0.05 1.2]);
set(haxes(1), 'XMinorTick', 'on', 'YMinorTick', 'on');
set(haxes(1), 'YTick', [0.0 0.2 0.4 0.6 0.8 1.0 1.2]);
set(haxes(2), 'YLim', [0.0 4e-5]);
set(haxes(2), 'YMinorTick', 'on');
set(haxes(2), 'YTick', [5.0e-6 1.0e-5 1.5e-5 2.0e-5 2.5e-5 3.0e-5 3.5e-5]);
for iaxis = 1:2
    % These properties must be the same for both sets of axes.
    set(haxes(iaxis), 'XLim', [0 12]);
    set(haxes(iaxis), 'XTick', [0.0001 0.001 0.01 0.1 1 10]);
    set(haxes(iaxis), 'FontSize', 13);
end
set(gca, 'box', 'off'); % so ticks aren't shown on opposite axes.
% Add title.
title(['Stiff Robertson Problem (Explicit ODE): BDF with Newton ', ...
    'Iterations'], titleFmt{:});
% Label the x axis, and both y axes.
xlabel('x', labFmt{:});
ylabel(haxes(1),'Solution (a,c)', labFmt{:});
ylabel(haxes(2),'Solution (b)', labFmt{:});
% Add a legend.
legend('a','c','b','Location','Best');
% Set some features of the three lines.
set(hline1, 'Marker', '+','Linestyle','-');
set(hline2, 'Marker', '*','Linestyle',':');
set(hline3, 'Marker', 'x','Linestyle','--');
 
nag_ode_ivp_stiff_exp_revcom example program results

      x           y(1)          y(2)          y(3)
    0.000       1.00000       0.00000       0.00000
    2.000       0.94161       0.00003       0.05836
    4.000       0.90551       0.00002       0.09446
    6.000       0.87926       0.00002       0.12072
    8.000       0.85854       0.00002       0.14144
   10.000       0.84135       0.00002       0.15863

 hused = 9.01776e-01   hnext = 9.01776e-01   tcur = 1.07662e+01
 nst = 55   nre = 128   nje = 16
 nqu = 4   nq = 4   niter = 78
 Max Err Comp = 3  


function d02nm_example
% Initialize variables and arrays.
neq = 3;
neqmax = neq;
nrw = 50+4*neqmax;
ninf = 23;
nwkjac = neqmax*(neqmax + 1);
njcpvt = 1;
maxord = 5;
sdysav = maxord+1;
maxstp = 200;
mxhnil = 5;
ldysav = neq;

h0 = 0.0;
hmax = 10.0;
hmin = 1.0e-10;
tcrit = 0.0;
petzld = false;

const = zeros(6, 1);
rwork = zeros(nrw, 1);
wkjac = zeros(nwkjac, 1);
ydot = zeros(neq, 1);
ysave = zeros(ldysav, sdysav);
inform = zeros(ninf, 1);
jacpvt = zeros(njcpvt, 1);
algequ = zeros(neq, 1);

fprintf('d02nm example program results\n');

% Set initial parameter values.  We integrate to tout by overshooting
% (itask = 1) using B.D.F. formulae with a Newton method.  Default values
% for the const array are used.  We employ scalar tolerances and the
% Jacobian is evaluated numerically.  On the reverse communication call
% (equivalent to calling the monitr routine in the forward communication
% version), we carry out interpolation using d02xf.
t = 0;
tout = 10.0;
itask = 1;
iout = 1;
xout = 2.0;
y = [1; 0; 0];
itol = 1;
rtol = [1e-04];
atol = [1e-07];

% d02nv is a setup routine (specifically for the backward
% differentiation formula (BDF) integrator) to be called prior to d02nm.
[const, rwork, ifail] = d02nv(int64(neqmax), int64(sdysav), ...
    int64(maxord), 'Newton', petzld, const, tcrit, hmin, ...
    hmax, h0, int64(maxstp), int64(mxhnil), 'Average-L2', rwork);
if ifail ~= 0
    % Unsuccessful call.  Print message and exit.
    error('d02nv returned with ifail = %1d ',ifail);
end

% d02ns is a setup routine (specifically for full matrix linear algebra)
% to be called prior to a routine from sub-chapter d02n-m (e.g. d02nm,
% as here).  We specify (jceval = 'Numerical') that the Jacobian is to be
% evaluated numerically.
[rwork, ifail] = d02ns(int64(neq), int64(neqmax), 'Numerical', ...
    int64(nwkjac), rwork);
lacorb = 50 + neq;
lacor1 = lacorb + 1;
lacor2 = lacorb + 2;
lacor3 = lacorb + 3;
lsavrb = lacorb + neq;
lsavr1 = lsavrb + 1;
lsavr2 = lsavrb + 2;
lsavr3 = lsavrb + 3;

% Output header and initial results.
fprintf('\n      x           y(1)          y(2)          y(3)\n');
fprintf(' %8.3f %13.5f %13.5f %13.5f\n', t, y(1), y(2), y(3));

% Prepare to store results for plotting.
ncall = 1;
ykeep = y';
tkeep = [t];

% d02nm uses reverse communication, involving initial entry, intermediate
% exits and re-entries, and a final exit.  These variables are only used
% on intermediate re-entry, but they need to be initialized (here, to a
% bogus value) before calling the routine for the first time.
irevcm = -999;
imon = -999;
inln = -999;
ires = -999;

% Ask for warning messages only.
itrace = 0;

% The reverse communication process is controlled by the value of irevcm
% (which is 0 for the final exit).
while (irevcm ~= 0)
    [t, y, ydot, rwork, inform, ysave, wkjac, jacpvt, imon, inln, ires, ...
        irevcm, ifail] = d02nm(t, tout, y, ydot, rwork, rtol, atol, ...
        int64(itol), int64(inform), ysave, wkjac, int64(jacpvt), ...
        int64(imon), int64(inln), int64(ires), int64(irevcm), ...
        int64(itask), int64(itrace), 'neq', int64(neq), 'sdysav', ...
        int64(sdysav), 'nwkjac', int64(nwkjac));

    % Evaluate the derivative, then use the value of irevcm to decide
    % where it is to be placed.
    f1 = -0.04*y(1) + 1.0e4*y(2)*y(3);
    f2 =  0.04*y(1) - 1.0e4*y(2)*y(3) - 3.0e7*y(2)*y(2);
    f3 =  3.0e7*y(2)*y(2);
    if (irevcm == 1 || irevcm == 3 )
        rwork(lsavr1) = f1;
        rwork(lsavr2) = f2;
        rwork(lsavr3) = f3;
    elseif (irevcm == 4)
        rwork(lacor1) = f1;
        rwork(lacor2) = f2;
        rwork(lacor3) = f3;
    elseif (irevcm == 5)
        ydot(1) = f1;
        ydot(2) = f2;
        ydot(3) = f3;
    elseif (irevcm == 9)
        % If the current step was successful, perform monitoring tasks.
        if (imon == 1)

            % Extract useful information about the current step.
            tc = rwork(19);
            hlast = rwork(15);
            hnext = rwork(16);
            nqu = int64(rwork(10));

            % For low values of the independent variable, store the results
            % for plotting.
            if tc < 2.0
                ncall = ncall + 1;
                tkeep(ncall,1) = tc;
                ykeep(ncall,:) = y;
            end

            % xout is the point at which we'll do interpolation on the
            % solution for output.  Check to see if we're near it.
            while (iout < 6 && tc-hlast < xout && xout <= tc)
                [sol, ifail] = d02xk(xout, int64(neq), ysave, ...
                    rwork(lacor1:lacor1+neq-1), tc, nqu, hlast, hnext);
                if (ifail ~= 0)
                    imon = -2;
                else

                    % Output interpolated result, and save it for plotting.
                    fprintf(' %8.3f %13.5f %13.5f %13.5f\n', xout, ...
                        sol(1), sol(2), sol(3));
                    ncall = ncall + 1;
                    tkeep(ncall,1) = xout;
                    ykeep(ncall,:) = sol;

                    iout = iout + 1;
                    xout = double(iout)*2.0;
                end
            end

            % Other values are illegal - e.g. irevcm = 8 asks for the
            % Jacobian to be calculated (analogously to the derivative
            % above), but we've already specified (in the call to d02ns)
            % that this is to be evaluated numerically, so this value of
            % irevcm should never occur here.
        elseif irevcm == 2 || irevcm == 6 || irevcm == 7 || irevcm == 8
            fprintf('Warning: illegal value of irevcm = %d in d02nm\n', ...
                irevcm);
            return
        end
    end
end

% d02ny is an integrator diagnostic routine which can be called
% after d02nm.  Output results.
[hu, h, tcur, tolsf, nst, nre, nje, nqu, nq, niter, imxer, algequ, ...
    ifail] = d02ny(int64(neq), int64(neqmax), rwork, inform);
fprintf('\n hused = %1.5e   hnext = %1.5e   tcur = %1.5e\n', hu, h, tcur);
fprintf(' nst = %1.0f   nre = %1.0f   nje = %1.0f\n',nst, nre, nje);
fprintf(' nqu = %1.0f   nq = %1.0f   niter = %1.0f\n',nqu, nq ,niter);
fprintf(' Max Err Comp = %1.0f  \n\n', imxer);
% Plot results.
fig = figure('Number', 'off');
display_plot(tkeep, ykeep)

function display_plot(x, y)
% Formatting for title and axis labels.
titleFmt = {'FontName', 'Helvetica', 'FontWeight', 'Bold', 'FontSize', 14};
labFmt = {'FontName', 'Helvetica', 'FontWeight', 'Bold', 'FontSize', 13};
set(gca, 'FontSize', 13); % for legend, axis tick labels, etc.
% Plot one of the curves and then add the other two.
hline1 = semilogx(x, y(:,1));
hold on
[haxes, hline2, hline3] = plotyy(x, y(:,3),x, y(:,2), ...
    'semilogx', 'semilogx');
% Set the axis limits and the tick specifications to beautify the plot.
set(haxes(1), 'YLim', [-0.05 1.2]);
set(haxes(1), 'XMinorTick', 'on', 'YMinorTick', 'on');
set(haxes(1), 'YTick', [0.0 0.2 0.4 0.6 0.8 1.0 1.2]);
set(haxes(2), 'YLim', [0.0 4e-5]);
set(haxes(2), 'YMinorTick', 'on');
set(haxes(2), 'YTick', [5.0e-6 1.0e-5 1.5e-5 2.0e-5 2.5e-5 3.0e-5 3.5e-5]);
for iaxis = 1:2
    % These properties must be the same for both sets of axes.
    set(haxes(iaxis), 'XLim', [0 12]);
    set(haxes(iaxis), 'XTick', [0.0001 0.001 0.01 0.1 1 10]);
    set(haxes(iaxis), 'FontSize', 13);
end
set(gca, 'box', 'off'); % so ticks aren't shown on opposite axes.
% Add title.
title(['Stiff Robertson Problem (Explicit ODE): BDF with Newton ', ...
    'Iterations'], titleFmt{:});
% Label the x axis, and both y axes.
xlabel('x', labFmt{:});
ylabel(haxes(1),'Solution (a,c)', labFmt{:});
ylabel(haxes(2),'Solution (b)', labFmt{:});
% Add a legend.
legend('a','c','b','Location','Best');
% Set some features of the three lines.
set(hline1, 'Marker', '+','Linestyle','-');
set(hline2, 'Marker', '*','Linestyle',':');
set(hline3, 'Marker', 'x','Linestyle','--');
 
d02nm example program results

      x           y(1)          y(2)          y(3)
    0.000       1.00000       0.00000       0.00000
    2.000       0.94161       0.00003       0.05836
    4.000       0.90551       0.00002       0.09446
    6.000       0.87926       0.00002       0.12072
    8.000       0.85854       0.00002       0.14144
   10.000       0.84135       0.00002       0.15863

 hused = 9.01776e-01   hnext = 9.01776e-01   tcur = 1.07662e+01
 nst = 55   nre = 128   nje = 16
 nqu = 4   nq = 4   niter = 78
 Max Err Comp = 3  



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