naginterfaces.library.ode.ivp_stiff_imp_revcom¶

naginterfaces.library.ode.ivp_stiff_imp_revcom(t, tout, y, ydot, comm, rtol, atol, itol, wkjac, imon, inln, ires, irevcm, lderiv, itask, itrace, io_manager=None)[source]¶

ivp_stiff_imp_revcom is a reverse communication function for integrating stiff systems of implicit ordinary differential equations coupled with algebraic equations.

For full information please refer to the NAG Library document for d02nn

https://support.nag.com/numeric/nl/nagdoc_30.3/flhtml/d02/d02nnf.html

Parameters

tfloat

On initial entry: $t$ , 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.

toutfloat

On initial entry: the next value of $t$ at which a computed solution is desired. For the initial $t$ , the input value of $t o u t$ is used to determine the direction of integration. Integration is permitted in either direction (see also $i t a s k$ ).

yfloat, ndarray, shape $(neq)$ , modified in place

On initial entry: the values of the dependent variables (solution). On the first call the first $neq$ elements of $y$ must contain the vector of initial values.

On final exit: the computed solution vector evaluated at $t$ (usually $t = t o u t$ ).

ydotfloat, ndarray, shape $(neq)$ , modified in place

On initial entry: if $l d e r i v [0] = T r u e$ , $y d o t$ must contain approximations to the time derivatives $y^{'}$ of the vector $y$ . If $l d e r i v [0] = F a l s e$ , $y d o t$ need not be set on entry.

On final exit: contains the time derivatives $y^{'}$ of the vector $y$ at the last integration point.

commdict, communication object, modified in place

Communication structure.

This argument must have been initialized by prior calls to one of ivp_stiff_fulljac_setup(), ivp_stiff_bandjac_setup() or ivp_stiff_sparjac_setup() and one of ivp_stiff_dassl(), ivp_stiff_bdf() or ivp_stiff_blend().

rtolfloat, array-like, shape $(:)$

Note: the required length for this argument is determined as follows: if $i t o l in (1, 2)$ : $1$ ; otherwise: $neq$ .

On initial entry: the relative local error tolerance.

atolfloat, array-like, shape $(:)$

Note: the required length for this argument is determined as follows: if $i t o l in (1, 3)$ : $1$ ; otherwise: $neq$ .

On initial entry: the absolute local error tolerance.

itolint

On initial entry: a value to indicate the form of the local error test. $i t o l$ indicates to ivp_stiff_imp_revcom whether to interpret either or both of $r t o l$ or $a t o l$ as a vector or a scalar. The error test to be satisfied is $∥ e_{i} / w_{i} ∥ < 1.0$ , where $w_{i}$ is defined as follows:

$i t o l$	$r t o l$	$a t o l$	$w_{i}$
1	scalar	scalar	$r t o l [0] \times \| y_{i} \| + a t o l [0]$
2	scalar	vector	$r t o l [0] \times \| y_{i} \| + a t o l [i - 1]$
3	vector	scalar	$r t o l [i - 1] \times \| y_{i} \| + a t o l [0]$
4	vector	vector	$r t o l [i - 1] \times \| y_{i} \| + a t o l [i - 1]$

$e_{i}$ is an estimate of the local error in $y_{i}$ , computed internally, and the choice of norm to be used is defined by a previous call to an integrator setup function.

wkjacfloat, ndarray, shape $(nwkjac)$ , modified in place

On intermediate entry: elements of the Jacobian as defined under the description of $i r e v c m$ . If a numerical Jacobian was requested then $w k j a c$ is used for workspace.

On intermediate exit: the Jacobian is overwritten.

imonint

On intermediate entry: may be reset to determine subsequent action in ivp_stiff_imp_revcom.

$i m o n = - 2$

Integration is to be halted. A return will be made from ivp_stiff_imp_revcom to the calling (sub)program with $e r r n o$ = 12.

$i m o n = - 1$

Allow ivp_stiff_imp_revcom to continue with its own internal strategy. The integrator will try up to three restarts unless $i m o n \neq - 1$ .

$i m o n = 0$

Return to the internal nonlinear equation solver, where the action taken is determined by the value of $i n l n$ .

$i m o n = 1$

Normal exit to ivp_stiff_imp_revcom to continue integration.

$i m o n = 2$

Restart the integration at the current time point. The integrator will restart from order $1$ when this option is used. The internal initialization module solves for new values of $y$ and $y^{'}$ by using the values supplied in $y$ and $y d o t$ by the $m o n i t r$ operation (see Notes) as initial estimates.

$i m o n = 3$

Try to continue with the same step size and order as was to be used before entering the $m o n i t r$ operation (see Notes). $h m i n$ and $h m a x$ may be altered if desired.

$i m o n = 4$

Continue the integration but using a new value of $h n e x t$ and possibly new values of $h m i n$ and $h m a x$ .

inlnint

On intermediate entry: with $i m o n = 0$ and $i r e v c m = 9$ , $i n l n$ specifies the action to be taken by the internal nonlinear equation solver. By setting $i n l n = 3$ and returning to ivp_stiff_imp_revcom, the residual vector is evaluated and placed in $c o m m [^{'} r w o r k^{'}] [50 + 2 \times neq + i - 1]$ , for $i = 1, 2, \dots, neq$ and then the $m o n i t r$ operation (see Notes) is invoked again. At present this is the only option available: $i n l n$ must not be set to any other value.

iresint

On intermediate entry: should be unchanged unless one of the following actions is required of ivp_stiff_imp_revcom in which case $i r e s$ should be set accordingly.

$i r e s = 2$

Indicates to ivp_stiff_imp_revcom that control should be passed back immediately to the calling (sub)program with the error indicator set to $e r r n o$ = 11.

$i r e s = 3$

Indicates to ivp_stiff_imp_revcom that an error condition has occurred in the solution vector, its time derivative or in the value of $t$ . The integrator will use a smaller time step to try to avoid this condition. If this is not possible ivp_stiff_imp_revcom returns to the calling (sub)program with the error indicator set to $e r r n o$ = 7.

$i r e s = 4$

Indicates to ivp_stiff_imp_revcom to stop its current operation and to enter the $m o n i t r$ operation (see Notes) immediately.

irevcmint

On initial entry: must contain $0$ .

On intermediate entry: should remain unchanged.

lderivbool, ndarray, shape $(2)$ , modified in place

On initial entry: $l d e r i v [0]$ must be set to $T r u e$ if you have supplied both an initial $y$ and an initial $y^{'}$ . $l d e r i v [0]$ must be set to $F a l s e$ if only the initial $y$ has been supplied.

$l d e r i v [1]$ must be set to $T r u e$ if the integrator is to use a modified Newton method to evaluate the initial $y$ and $y^{'}$ .

Note that $y$ and $y^{'}$ , if supplied, are used as initial estimates.

This method involves taking a small step at the start of the integration, and if $i t a s k = 6$ on entry, $t$ and $t o u t$ will be set to the result of taking this small step. $l d e r i v [1]$ must be set to $F a l s e$ if the integrator is to use functional iteration to evaluate the initial $y$ and $y^{'}$ , and if this fails a modified Newton method will then be attempted. $l d e r i v [1] = T r u e$ is recommended if there are implicit equations or the initial $y$ and $y^{'}$ are zero.

On final exit: $l d e r i v [0]$ is normally unchanged. However if $i t a s k = 6$ and internal initialization was successful then $l d e r i v [0] = T r u e$ .

$l d e r i v [1] = T r u e$ , if implicit equations were detected.

Otherwise $l d e r i v [1] = F a l s e$ .

itaskint

On initial entry: the task to be performed by the integrator.

$i t a s k = 1$

Normal computation of output values of $y (t)$ at $t = t o u t$ (by overshooting and interpolating).

$i t a s k = 2$

Take one step only and return.

$i t a s k = 3$

Stop at the first internal integration point at or beyond $t = t o u t$ and return.

$i t a s k = 4$

Normal computation of output values of $y (t)$ at $t = t o u t$ but without overshooting $t = 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 ivp_stiff_bdf()) may be equal to or beyond $t o u t$ , but not before it in the direction of integration.

$i t a s k = 5$

Take one step only and return, without passing $tcrit$ (e.g., see ivp_stiff_bdf()). $tcrit$ must be specified under $i t a s k = 4$ .

$i t a s k = 6$

The integrator will solve for the initial values of $y$ and $y^{'}$ only and then return to the calling (sub)program without doing the integration. This option can be used to check the initial values of $y$ and $y^{'}$ . Functional iteration or a ‘small’ backward Euler method used in conjunction with a damped Newton iteration is used to calculate these values (see $l d e r i v$ ). Note that if a backward Euler step is used then the value of $t$ will have been advanced a short distance from the initial point.

Note: if ivp_stiff_imp_revcom is recalled with a different value of $i t a s k$ (and $t o u t$ altered) then the initialization procedure is repeated, possibly leading to different initial conditions.

itraceint

On initial entry: the level of output that is printed by the integrator. $i t r a c e$ may take the value $- 1$ , $0$ , $1$ , $2$ or $3$ .

$i t r a c e < - 1$

$- 1$ is assumed and similarly if $i t r a c e > 3$ , $3$ is assumed.

$i t r a c e = - 1$

No output is generated.

$i t r a c e = 0$

Only warning messages are printed on the current error message unit (see FileObjManager).

$i t r a c e > 0$

Warning messages are printed as above, and on the file object associated with the advisory I/O unit (see FileObjManager) 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 $i t r a c e$ .

io_managerFileObjManager, optional

Manager for I/O in this routine.

Returns

tfloat

On final exit: the value at which the computed solution $y$ is returned (usually at $t o u t$ ).

toutfloat

Is unaltered unless $i t a s k = 6$ and $l d e r i v [1] = T r u e$ on entry (see also $i t a s k$ and $l d e r i v$ ) in which case $t o u t$ will be set to the result of taking a small step at the start of the integration.

imonint

On intermediate exit: used to pass information between ivp_stiff_imp_revcom and the $m o n i t r$ operation (see Notes). With $i r e v c m = 9$ , $i m o n$ contains a flag indicating under what circumstances the return from ivp_stiff_imp_revcom occurred:

$i m o n = - 2$

Exit from ivp_stiff_imp_revcom after $i r e s = 4$ (set in the $r e s i d$ operation (see Notes) caused an early termination (this facility could be used to locate discontinuities).

$i m o n = - 1$

The current step failed repeatedly.

$i m o n = 0$

Exit from ivp_stiff_imp_revcom after a call to the internal nonlinear equation solver.

$i m o n = 1$

The current step was successful.

inlnint

On intermediate exit: contains a flag indicating the action to be taken, if any, by the internal nonlinear equation solver.

iresint

On intermediate exit: with $i r e v c m = 1$ , $2$ , $3$ , $4$ , $5$ , $6$ , $7$ or $11$ , $i r e s$ specifies the form of the residual to be returned by the $r e s i d$ operation (see Notes).

If $i r e s = 1$ , $- r = g (t, y) - A (t, y) y^{'}$ must be returned.

If $i r e s = - 1$ , $-^r = - A (t, y) y^{'}$ must be returned.

irevcmint

On intermediate exit: indicates what action you must take before re-entering ivp_stiff_imp_revcom. The possible exit values of $i r e v c m$ are $1$ , $2$ , $3$ , $4$ , $5$ , $6$ , $7$ , $8$ , $9$ , $10$ or $11$ which should be interpreted as follows:

$i r e v c m = 1$ , $2$ , $3$ , $4$ , $5$ , $6$ , $7$ or $11$

Indicates that a $r e s i d$ operation (see Notes) is required: you must supply the residual of the system. For each of these values of $i r e v c m$ , $y_{i}$ is located in $y [i - 1]$ , for $i = 1, 2, \dots, neq$ .

For $i r e v c m = 1$ , $3$ , $6$ or $11$ , $y_{i}^{'}$ is located in $y d o t [i - 1]$ and $r_{i}$ should be stored in $c o m m [^{'} r w o r k^{'}] [50 + 2 \times neq + i - 1]$ , for $i = 1, 2, \dots, neq$ .

For $i r e v c m = 2$ , $y_{i}^{'}$ is located in $c o m m [^{'} r w o r k^{'}] [50 + neq + i - 1]$ and $r_{i}$ should be stored in $c o m m [^{'} r w o r k^{'}] [50 + 2 \times neq + i - 1]$ , for $i = 1, 2, \dots, neq$ .

For $i r e v c m = 4$ or $7$ , $y_{i}^{'}$ is located in $y d o t [i - 1]$ and $r_{i}$ should be stored in $c o m m [^{'} r w o r k^{'}] [50 + neq + i - 1]$ , for $i = 1, 2, \dots, neq$ .

For $i r e v c m = 5$ , $y_{i}^{'}$ is located in $c o m m [^{'} r w o r k^{'}] [50 + 2 \times neq + i - 1]$ and $r_{i}$ should be stored in $y d o t [i - 1]$ , for $i = 1, 2, \dots, neq$ .

$i r e v c m = 8$

Indicates that a $j a c$ operation (see Notes) is required: you must supply the Jacobian matrix.

If full matrix linear algebra is being used, the $(i, j)$ th element of the Jacobian must be stored in $w k j a c [(j - 1) \times neq + i - 1]$ .

If banded matrix linear algebra is being used, the $(i, j)$ th element of the Jacobian must be stored in $w k j a c [(i - 1) \times m_{B} + k - 1]$ , where $m_{B} = m_{L} + m_{U} + 1$ and $k = m i n (m_{L} - i + 1, 0) + j$ ; here $m_{L}$ and $m_{U}$ are the number of subdiagonals and superdiagonals, respectively, in the band.

If sparse matrix linear algebra is being used, ivp_stiff_sparjac_enq() must be called to determine which column of the Jacobian is required and where it should be stored.

$i r e v c m = 9$

Indicates that a $m o n i t r$ operation (see Notes) can be performed.

$i r e v c m = 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 variable $t$ , located in $c o m m [^{'} r w o r k^{'}] [18]$ . No values must be changed before re-entering ivp_stiff_imp_revcom; this facility enables you to determine the number of unsuccessful steps.

On final exit: $i r e v c m = 0$ indicating that the user-specified task has been completed or an error has been encountered (see the descriptions for $i t a s k$ and $errno$ ).

Raises

NagValueError

(errno $1$ )

On entry, $i r e s = ⟨ v a l u e ⟩$ and $\frac{d y}{d t} = 0.0$ for all elements.

Check the evaluation of the residual for this value of $i r e s$ .

(errno $1$ )

On re-entry, $i r e s$ has been set to an illegal value during initialization.

$i r e s = ⟨ v a l u e ⟩$ at time $⟨ v a l u e ⟩$ .

(errno $1$ )

On re-entry, $i m o n = ⟨ v a l u e ⟩$ .

Constraint: $- 2 \leq i m o n \leq 4$ .

(errno $1$ )

Failure during internal time interpolation. $t o u t$ and the current time are too close.

$i t a s k = ⟨ v a l u e ⟩$ and $t o u t = ⟨ v a l u e ⟩$ and the current time is $⟨ v a l u e ⟩$ .

(errno $1$ )

On entry, $i t a s k = 4$ or $5$ and $tcrit$ is before the current time in the direction of integration.

$i t a s k = ⟨ v a l u e ⟩$ , $tcrit = ⟨ v a l u e ⟩$ and the current time is $⟨ v a l u e ⟩$ .

(errno $1$ )

$i t a s k = 3$ and $t o u t$ is more than an integration step behind the current time.

$t o u t = ⟨ v a l u e ⟩$ , current time minus step size: $⟨ v a l u e ⟩$ .

(errno $1$ )

The initial stepsize, $h = ⟨ v a l u e ⟩$ , is too small.

(errno $1$ )

Weight number $i = ⟨ v a l u e ⟩$ used in the local error test is too small. Check the values of $r t o l$ and $a t o l$ .

$a t o l [i - 1]$ and $y [i - 1]$ may both be zero.

Weight $i = ⟨ v a l u e ⟩$ .

(errno $1$ )

On entry, $i t a s k = 4$ or $5$ and $tcrit$ is before $t o u t$ in the direction of integration.

$i t a s k = ⟨ v a l u e ⟩$ , $tcrit = ⟨ v a l u e ⟩$ and $t o u t = ⟨ v a l u e ⟩$ .

(errno $1$ )

On entry, $a t o l = ⟨ v a l u e ⟩$ .

Constraint: $a t o l \geq 0.0$ .

(errno $1$ )

On entry, $r t o l = ⟨ v a l u e ⟩$ .

Constraint: $r t o l \geq 0.0$ .

(errno $1$ )

Either the value of $nwkjac$ is not the same as the value supplied to the setup function or a communication array has become corrupted.

$nwkjac = ⟨ v a l u e ⟩$ , $nwkjac$ (setup) $= ⟨ v a l u e ⟩$ .

(errno $1$ )

Either the value of $sdysav$ is not the same as the value supplied to the setup function or a communication array has become corrupted.

$sdysav = ⟨ v a l u e ⟩$ , $sdysav$ (setup) $= ⟨ v a l u e ⟩$ .

(errno $1$ )

On entry, an illegal (negative) minimum stepsize was provided in a prior call to a setup function. $hmin = ⟨ v a l u e ⟩$ .

(errno $1$ )

On entry, $t o u t$ is less than $t$ with respect to the direction of integration given by the sign of $h0$ in a prior call to a setup function.

$t o u t = ⟨ v a l u e ⟩$ , $t = ⟨ v a l u e ⟩$ and $h0 = ⟨ v a l u e ⟩$ .

(errno $1$ )

On entry, an illegal (negative) maximum number of steps was provided in a prior call to a setup function. $maxstp = ⟨ v a l u e ⟩$ .

(errno $1$ )

On entry, $i t o l = ⟨ v a l u e ⟩$ .

Constraint: $1 \leq i t o l \leq 4$ .

(errno $1$ )

On entry, $neq = ⟨ v a l u e ⟩$ and $ldysav = ⟨ v a l u e ⟩$ .

Constraint: $neq \leq ldysav$ .

(errno $1$ )

On entry, $neq = ⟨ v a l u e ⟩$ .

Constraint: $neq \geq 1$ .

(errno $1$ )

Either the function was entered on a continuation call without a prior call of this function, or a communication array has become corrupted.

(errno $1$ )

On entry, $t o u t$ is too close to $t$ to start integration.

$t o u t = ⟨ v a l u e ⟩$ and $t = ⟨ v a l u e ⟩$ .

(errno $1$ )

On entry, $i t a s k = ⟨ v a l u e ⟩$ .

Constraint: $1 \leq i t a s k \leq 5$ .

(errno $1$ )

Either the integrator setup function has not been called prior to the first call of this function, or a communication array has become corrupted.

(errno $1$ )

On entry, an illegal (negative) maximum stepsize was provided in a prior call to a setup function. $hmax = ⟨ v a l u e ⟩$ .

(errno $1$ )

Either the linear algebra setup function has not been called prior to the first call of this function, or a communication array has become corrupted.

(errno $1$ )

On re-entry, $i n l n = ⟨ v a l u e ⟩$ for the case $i r e v c m = 9$ and $i m o n = 0$ .

Constraint: $i n l n = 3$ .

(errno $1$ )

On entry, $i r e v c m = ⟨ v a l u e ⟩$ .

Constraint: $0 \leq i r e v c m \leq 11$ .

(errno $1$ )

On re-entry, the solution vector appears to have been overwritten.

Further integration will not be attempted.

(errno $2$ )

At time $⟨ v a l u e ⟩$ the maximum number of allowed steps on this call was taken before reaching the next output point $t o u t = ⟨ v a l u e ⟩$ .

Maximum number of steps $= ⟨ v a l u e ⟩$ .

(errno $7$ )

On re-entry, $i r e s = 3$ which signals that an error condition has occurred in the solution vector, its time derivative or in the value of $t$ . It was not possible to remove this condition.

$i r e s = ⟨ v a l u e ⟩$ at $t = ⟨ v a l u e ⟩$ .

(errno $8$ )

Attempt was made to reduce the step size to a value less than the minimum step size during the calculation of initial values.

Minimum stepsize: $⟨ v a l u e ⟩$ .

(errno $8$ )

The user problem has one or more inconsistencies between the $i r e s = 1$ and $i r e s = - 1$ parts. Integration will not be attempted.

(errno $8$ )

The residual function returned an error when calculating the initial values of the solution and its time derivative.

(errno $8$ )

Workspace error occurred when trying to form the Jacobian matrix in calculating the initial values of the solution and its time derivative.

(errno $9$ )

A singular Jacobian has been encountered. You should check the problem formulation and Jacobian calculation.

(errno $10$ )

Workspace error occurred when trying to form the Jacobian matrix in calculating the initial values of the solution and its time derivative.

(errno $10$ )

Not enough integer store provided for sparse matrix solver.

Units of store needed: $⟨ v a l u e ⟩$ . Amount provided: $⟨ v a l u e ⟩$ .

(errno $10$ )

Larger integer workspace required.

Provided: $⟨ v a l u e ⟩$ ; required: $⟨ v a l u e ⟩$ .

(errno $10$ )

Not enough real store provided for sparse matrix solver.

Units of store needed: $⟨ v a l u e ⟩$ . Amount provided: $⟨ v a l u e ⟩$ .

(errno $14$ )

On entry, too much accuracy requested for precision of the machine at the start of problem. The tolerances should be checked; the requested accuracy should be reduced by a factor of at least $⟨ v a l u e ⟩$ .

Warns

NagAlgorithmicWarning

(errno $3$ )

Too much accuracy requested for precision of the machine at time $⟨ v a l u e ⟩$ .

The tolerances should be checked; the requested accuracy should be reduced by a factor of at least $⟨ v a l u e ⟩$ .

(errno $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.

(errno $5$ )

Nonlinear solver failed to converge using a damped Newton method to solve for initial values.

Damping factor: $⟨ v a l u e ⟩$ ; convergence rate: $⟨ v a l u e ⟩$ .

(errno $6$ )

At time $⟨ v a l u e ⟩$ , error weight $⟨ v a l u e ⟩$ became zero. Check the values of $a t o l$ , $r t o l$ and $i t o l$ supplied.

(errno $11$ )

On re-entry, $i r e s = 2$ , which signals that the integration should terminate. $i r e s = ⟨ v a l u e ⟩$ at time $⟨ v a l u e ⟩$ .

(errno $12$ )

A return was forced by setting $i m o n = - 2$ , but the integration was successful as far as $t$ .

(errno $13$ )

The requested task has been completed, but it is estimated that a small change in $r t o l$ and $a t o l$ is unlikely to produce any change in the computed solution. (This ONLY applies when you are NOT operating in one step mode; that is, when $i t a s k \neq 2$ or $5$ .)

Notes

No equivalent traditional C interface for this routine exists in the NAG Library.

ivp_stiff_imp_revcom is a general purpose function for integrating the initial value problem for a stiff system of implicit ordinary differential equations coupled with algebraic equations, written in the form

A (t, y) y^{'} = g (t, y) .

An outline of a typical calling program is

call the linear algebra setup routine

call the integrator setup routine

set $i r e v c m = 0$

start looping

call ivp_stiff_imp_revcom

if $i r e v c m = 8$ then

supply the Jacobian matrix

else if $i r e v c m = 9$ then

perform monitoring tasks requested by the user

else if $i r e v c m = 10$ then

indicates an unsuccessful step

else if $i r e v c m > 0$ then

evaluate the residual

else

stop looping

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 function on an intermediate return ( $i r e v c m \neq 0$ ) from ivp_stiff_imp_revcom; these are denoted (i), (ii) and (iii).

The following sections describe in greater detail exactly what is required of each of these operations.

Supply the Jacobian matrix

You need only provide this facility if the argument $jceval ='A'$ (or $jceval ='F'$ if using sparse matrix linear algebra) in a call to the linear algebra setup function (see $jceval$ in ivp_stiff_sparjac_setup()). 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, $y^{'}$ , has the form

$y^{'} = (y - z) / (h d),$

where $h$ is the current step size and $d$ is an argument that depends on the integration method in use. The vector $y$ is the current solution and the vector $z$ depends on information from previous time steps. This means that $\frac{d}{d y^{'}} () = (h d) \frac{d}{d y} ()$ .

The system of nonlinear equations that is solved has the form

$A (t, y) y^{'} - g (t, y) = 0$

but is solved in the form

$f (t, y) = 0,$

where $f$ is the function defined by

$f (t, y) = (h d) (A (t, y) (y - z) / (h d) - g (t, y)) .$

It is the Jacobian matrix $\frac{\partial r}{\partial y}$ that you must supply as follows:

$\frac{\partial f_{i}}{\partial y_{j}} = a_{i j} (t, y) + h d \frac{\partial}{\partial y_{j}} (neq \sum k = 1 a_{i k} (t, y) {y^{'}}_{k} - g_{i} (t, y)),$

where $t$ , $h$ and $d$ are located in $c o m m [^{'} r w o r k^{'}] [18]$ , $c o m m [^{'} r w o r k^{'}] [15]$ and $c o m m [^{'} r w o r k^{'}] [19]$ respectively and the arrays $y$ and $y d o t$ contain the current solution and time derivatives respectively. Only the nonzero elements of the Jacobian need be set, since the locations where it is to be stored are preset to zero.

In this document this operation is referred to as JAC.

Perform tasks requested by you

This operation is essentially a monitoring function and additionally provides the opportunity of changing the current values of $y$ , $y d o t$ , $h n e x t$ (the step size that the integrator proposes to take on the next step), $h m i n$ (the minimum step size to be taken on the next step), and $h m a x$ (the maximum step size to be taken on the next step). The scaled local error at the end of a time step may be obtained by calling ivp_stiff_errest().

The following gives details of the location within the array $c o m m$ [‘rwork’] of variables that may be of interest to you:

Variable	Specification	Location
$t c u r r$	the current value of the independent variable	$c o m m [^{'} r w o r k^{'}] [18]$
$h l a s t$	last step size successfully used by the integrator	$c o m m [^{'} r w o r k^{'}] [14]$
$h n e x t$	step size that the integrator proposes to take on the next step	$c o m m [^{'} r w o r k^{'}] [15]$
$h m i n$	minimum step size to be taken on the next step	$c o m m [^{'} r w o r k^{'}] [16]$
$h m a x$	maximum step size to be taken on the next step	$c o m m [^{'} r w o r k^{'}] [17]$
$n q u$	the order of the integrator used on the last step	$c o m m [^{'} r w o r k^{'}] [9]$

You are advised to consult the description of $monitr$ in ivp_stiff_imp_fulljac() for details on what optional input can be made.

If either $y$ or $y d o t$ are changed, then $i m o n$ must be set to $2$ before return to ivp_stiff_imp_revcom. If either of the values $h m i n$ or $h m a x$ are changed, then $i m o n$ must be set $\geq 3$ before return to ivp_stiff_imp_revcom. If $h n e x t$ is changed, then $i m o n$ must be set to $4$ before return to ivp_stiff_imp_revcom.

In addition you can force ivp_stiff_imp_revcom to evaluate the residual vector

A (t, y) y^{'} - g (t, y)

by setting $i m o n = 0$ and $i n l n = 3$ and then returning to ivp_stiff_imp_revcom; on return to this monitoring operation the residual vector will be stored in $c o m m [^{'} r w o r k^{'}] [50 + 2 \times neq + i - 1]$ , for $i = 1, 2, \dots, neq$ .

In this document this operation is referred to as MONITR.

Evaluate the residual

This operation must evaluate the residual

$- r = g (t, y) - A (t, y) y^{'}$

in one case and the reduced residual

$-^r = - A (t, y) y^{'}$

in another, where $t$ is located in $c o m m [^{'} r w o r k^{'}] [18]$ . The form of the residual that is returned is determined by the value of $i r e s$ returned by ivp_stiff_imp_revcom. If $i r e s = - 1$ , then the residual defined by equation (2) above must be returned; if $i r e s = 1$ , then the residual returned by equation (1) above must be returned.

In this document this operation is referred to as RESID.

NAG and Python

Return to Front

naginterfaces.library.ode.ivp_stiff_imp_revcom¶

naginterfaces.library.ode.ivp_​stiff_​imp_​revcom¶

naginterfaces.library.ode.ivp_stiff_imp_revcom¶