naginterfaces.library.ode.ivp_​stiff_​exp_​sparjac

naginterfaces.library.ode.ivp_stiff_exp_sparjac(t, tout, y, comm, rtol, atol, itol, fcn, itask, itrace, jac=None, monitr=None, data=None, io_manager=None)

ivp_stiff_exp_sparjac is a direct communication function for integrating stiff systems of explicit ordinary differential equations when the Jacobian is a sparse matrix.

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

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

Parameters

tfloat

$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

The next value of $t$ at which a computed solution is desired. For the initial $t$, the input value of $tout$ is used to determine the direction of integration. Integration is permitted in either direction (see also $itask$).

yfloat, array-like, shape $\left(\text{neq}\right)$

The values of the dependent variables (solution). On the first call the first $\text{neq}$ elements of $y$ must contain the vector of initial values.

commdict, communication object, modified in place

Communication structure.

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

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

Note: the required length for this argument is determined as follows: if $itol\text{in}(1,2)$: $1$; otherwise: $\text{neq}$.

The relative local error tolerance.

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

Note: the required length for this argument is determined as follows: if $itol\text{in}(1,3)$: $1$; otherwise: $\text{neq}$.

The absolute local error tolerance.

itolint

A value to indicate the form of the local error test. $itol$ indicates to ivp_stiff_exp_sparjac whether to interpret either or both of $rtol$ or $atol$ as a vector or a scalar. The error test to be satisfied is $\parallel e_i/w_i\parallel <1.0$, where $w_i$ is defined as follows:

$itol$	$rtol$	$atol$	$w_i$
1	scalar	scalar	$rtol\left[0\right]\times \left|y_i\right|+atol\left[0\right]$
2	scalar	vector	$rtol\left[0\right]\times \left|y_i\right|+atol[i-1]$
3	vector	scalar	$rtol[i-1]\times \left|y_i\right|+atol\left[0\right]$
4	vector	vector	$rtol[i-1]\times \left|y_i\right|+atol[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.

fcncallable (f, ires) = fcn(t, y, ires, data=None)

$fcn$ must evaluate the derivative vector for the explicit ordinary differential equation system, defined by $y^{\prime }=g(t,y)$.

Parameters

tfloat

$t$, the current value of the independent variable.

yfloat, ndarray, shape $\left(\text{neq}\right)$

The value of $y_{\text{i}}$, for $\text{i}=1,2,\dots ,\text{neq}$.

iresint

$ires=1$.

dataarbitrary, optional, modifiable in place

User-communication data for callback functions.

Returns

ffloat, array-like, shape $\left(\text{neq}\right)$

The value $y_{\text{i}}^{\prime }$, given by $y_{\text{i}}^{\prime }=g_{\text{i}}(t,y)$, for $\text{i}=1,2,\dots ,\text{neq}$.

iresint

You may set $ires$ as follows to indicate certain conditions in $fcn$ to the integrator:

$ires=1$

Indicates a normal return from $fcn$, that is $ires$ has not been altered by you and integration continues.

$ires=2$

Indicates to the integrator that control should be passed back immediately to the calling (sub)program with the error indicator set to $errno$ = 11.

$ires=3$

Indicates to the integrator 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 the integrator returns to the calling (sub)program with the error indicator set to $errno$ = 7.

$ires=4$

Indicates to the integrator to stop its current operation and to enter $monitr$ immediately with argument $imon=-2$.

itaskint

The task to be performed by the integrator.

$itask=1$

Normal computation of output values of $y\left(t\right)$ at $t=tout$ (by overshooting and interpolating).

$itask=2$

Take one step only and return.

$itask=3$

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

$itask=4$

Normal computation of output values of $y\left(t\right)$ at $t=tout$ but without overshooting $t=\text{tcrit}$ (e.g., see ivp_stiff_dassl()). $\text{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. $\text{tcrit}$ may be equal to or beyond $tout$, but not before it, in the direction of integration.

$itask=5$

Take one step only and return, without passing $\text{tcrit}$ (e.g., see ivp_stiff_dassl()). $\text{tcrit}$ must be specified as under $itask=4$.

itraceint

The level of output that is printed by the integrator. $itrace$ may take the value $-1$, $0$, $1$, $2$ or $3$.

$itrace<-1$

$-1$ is assumed and similarly if $itrace>3$, $3$ is assumed.

$itrace=-1$

No output is generated.

$itrace=0$

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

$itrace>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 $itrace$.

jacNone or callable pdj = jac(t, y, h, d, j, pdj, data=None), optional

Note: if this argument is None then a NAG-supplied facility will be used.

$jac$ must evaluate the Jacobian of the system.

If this option is not required, the actual argument for $jac$ must be None. You must indicate to the integrator whether this option is to be used by setting the argument $\text{jceval}$ appropriately in a call to the sparse linear algebra setup function ivp_stiff_sparjac_setup().

First we must define the system of nonlinear equations which is solved internally by the integrator.

The time derivative, $y^{\prime }$, generated internally, has the form

$$y^{\prime }=(y-z)/\left(hd\right)\text{,}$$

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}^{\prime }}\left(\right)=\left(hd\right)\frac{d}{dy}\left(\right)$.

The system of nonlinear equations that is solved has the form

$$y^{\prime }-g(t,y)=0$$

but it is solved in the form

$$r(t,y)=0\text{,}$$

where $r$ is the function defined by

$$r(t,y)=\left(hd\right)\left((y-z)/\left(hd\right)-g(t,y)\right)\text{.}$$

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

$$\begin{array}{c}\begin{array}{cc}\frac{\partial r_i}{\partial y_j}=1-\left(hd\right)\frac{\partial g_i}{\partial y_j}\text{,}& \text{if}i=j\text{,}\\ & \\ \frac{\partial r_i}{\partial y_j}=-\left(hd\right)\frac{\partial g_i}{\partial y_j}\text{,}& \text{otherwise.}\end{array}\end{array}$$

Parameters

tfloat

$t$, the current value of the independent variable.

yfloat, ndarray, shape $\left(\text{neq}\right)$

$y_{\text{i}}$, for $\text{i}=1,2,\dots ,\text{neq}$, the current solution component.

hfloat

The current step size.

dfloat

The argument $d$ which depends on the integration method.

jint

The column of the Jacobian that $jac$ must return in the array $pdj$.

pdjfloat, ndarray, shape $\left(\text{neq}\right)$

Is set to zero.

dataarbitrary, optional, modifiable in place

User-communication data for callback functions.

Returns

pdjfloat, array-like, shape $\left(\text{neq}\right)$

$pdj[i-1]$ should be set to the $(i,j)$th element of the Jacobian, where $j$ is given by $j$. Only nonzero elements of this array need be set, since it is preset to zero before the call to $jac$.

monitrNone or callable (hnext, y, imon, inln, hmin, hmax) = monitr(t, hlast, hnext, y, ydot, comm, r, acor, imon, hmin, hmax, nqu, data=None), optional

Note: if this argument is None then a NAG-supplied facility will be used.

$monitr$ performs tasks requested by you.

If this option is not required, the actual argument for $monitr$ must be None.

Parameters

tfloat

The current value of the independent variable.

hlastfloat

The last step size successfully used by the integrator.

hnextfloat

The step size that the integrator proposes to take on the next step.

yfloat, ndarray, shape $\left(\text{neq}\right)$

$y$, the values of the dependent variables evaluated at $t$.

ydotfloat, ndarray, shape $\left(\text{neq}\right)$

The time derivatives $y^{\prime }$ of the vector $y$.

commdict, communication object, modifiable in place

Communication structure.

This argument may be passed as argument $\text{comm}$ in a call to ivp_stiff_nat_interp() or $\text{comm}$ in a call to ivp_stiff_c1_interp() to enable you to carry out interpolation using either of those functions.

rfloat, ndarray, shape $\left(\text{neq}\right)$

If $imon=0$ and $inln=3$, the first $\text{neq}$ elements contain the residual vector, $y^{\prime }-g(t,y)$.

acorfloat, ndarray, shape $(\text{neq},2)$

With $imon=1$, $acor[i-1,0]$ contains the weight used for the $i$th equation when the norm is evaluated, and $acor[i-1,1]$ contains the estimated local error for the $i$th equation. The scaled local error at the end of a timestep may be obtained by calling ivp_stiff_errest().

imonint

A flag indicating under what circumstances $monitr$ was called:

$imon=-2$

Entry from the integrator after $ires=4$ (set in $fcn$) caused an early termination (this facility could be used to locate discontinuities).

$imon=-1$

The current step failed repeatedly.

$imon=0$

Entry after a call to the internal nonlinear equation solver (see $inln$).

$imon=1$

The current step was successful.

hminfloat

The minimum step size to be taken on the next step.

hmaxfloat

The maximum step size to be taken on the next step.

nquint

The order of the integrator used on the last step. This is supplied to enable you to carry out interpolation using either of the functions ivp_stiff_nat_interp() or ivp_stiff_c1_interp().

dataarbitrary, optional, modifiable in place

User-communication data for callback functions.

Returns

hnextfloat

The next step size to be used. If this is different from the input value, $imon$ must be set to $4$.

yfloat, array-like, shape $\left(\text{neq}\right)$

These values must not be changed unless $imon$ is set to $2$.

imonint

May be reset to determine subsequent action in ivp_stiff_exp_sparjac.

$imon=-2$

Integration is to be halted. A return will be made from the integrator to the calling (sub)program with $errno$ = 12.

$imon=-1$

Allow the integrator to continue with its own internal strategy. The integrator will try up to three restarts unless $imon$ is set $\ne -1$ on exit.

$imon=0$

Return to the internal nonlinear equation solver, where the action taken is determined by the value of $inln$ (see $inln$).

$imon=1$

Normal exit to the integrator to continue integration.

$imon=2$

Restart the integration at the current time point. The integrator will restart from order $1$ when this option is used. The solution $y$, provided by $monitr$, will be used for the initial conditions.

$imon=3$

Try to continue with the same step size and order as was to be used before the call to $monitr$. $hmin$ and $hmax$ may be altered if desired.

$imon=4$

Continue the integration but using a new value of $hnext$ and possibly new values of $hmin$ and $hmax$.

inlnint

The action to be taken by the internal nonlinear equation solver when $monitr$ is exited with $imon=0$. By setting $inln=3$ and returning to the integrator, the residual vector is evaluated and placed in the array $r$, and then $monitr$ is called again. At present this is the only option available: $inln$ must not be set to any other value.

hminfloat

The minimum step size to be used. If this is different from the input value, $imon$ must be set to $3$ or $4$.

hmaxfloat

The maximum step size to be used. If this is different from the input value, $imon$ must be set to $3$ or $4$. If $hmax$ is set to zero, no limit is assumed.

dataarbitrary, optional

User-communication data for callback functions.

io_managerFileObjManager, optional

Manager for I/O in this routine.

Returns

tfloat

The value at which the computed solution $y$ is returned (usually at $tout$).

yfloat, ndarray, shape $\left(\text{neq}\right)$

The computed solution vector, evaluated at $t$ (usually $t=tout$).

ydotfloat, ndarray, shape $\left(\text{neq}\right)$

The time derivatives $y^{\prime }$ of the vector $y$ at the last integration point.

Raises

NagValueError

(errno $1$)

$monitr$ set $imon=⟨value⟩$.

Constraint: $-2\le imon\le 4$.

(errno $1$)

Failure during internal time interpolation. $tout$ and the current time are too close.

$itask=⟨value⟩$ and $tout=⟨value⟩$ and the current time is $⟨value⟩$.

(errno $1$)

On entry, $itask=4$ or $5$ and $\text{tcrit}$ is before the current time in the direction of integration.

$itask=⟨value⟩$, $\text{tcrit}=⟨value⟩$ and the current time is $⟨value⟩$.

(errno $1$)

$itask=3$ and $tout$ is more than an integration step behind the current time.

$tout=⟨value⟩$, current time minus step size: $⟨value⟩$.

(errno $1$)

The initial stepsize, $⟨value⟩$, is too small.

(errno $1$)

On entry, $itask=4$ or $5$ and $\text{tcrit}$ is before $tout$ in the direction of integration.

$itask=⟨value⟩$, $\text{tcrit}=⟨value⟩$ and $tout=⟨value⟩$.

(errno $1$)

On entry, $atol=⟨value⟩$.

Constraint: $atol\ge 0.0$.

(errno $1$)

On entry, $rtol=⟨value⟩$.

Constraint: $rtol\ge 0.0$.

(errno $1$)

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

$\text{sdysav}=⟨value⟩$, $\text{sdysav}$ (setup) $=⟨value⟩$.

(errno $1$)

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

(errno $1$)

On entry, $tout$ is less than $t$ with respect to the direction of integration given by the sign of $\text{h0}$ in a prior call to a setup function.

$tout=⟨value⟩$, $t=⟨value⟩$ and $\text{h0}=⟨value⟩$.

(errno $1$)

On entry, an illegal (negative) maximum number of steps was provided in a prior call to a setup function. $\text{maxstp}=⟨value⟩$.

(errno $1$)

On entry, $itol=⟨value⟩$.

Constraint: $1\le itol\le 4$.

(errno $1$)

On entry, $\text{neq}=⟨value⟩$ and $\text{ldysav}=⟨value⟩$.

Constraint: $\text{neq}\le \text{ldysav}$.

(errno $1$)

On entry, $\text{neq}=⟨value⟩$.

Constraint: $\text{neq}\ge 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, $tout$ is too close to $t$ to start integration.

$tout=⟨value⟩$ and $t=⟨value⟩$.

(errno $1$)

On entry, $itask=⟨value⟩$.

Constraint: $1\le itask\le 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=⟨value⟩$.

(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$)

$monitr$ set $inln=⟨value⟩$.

Constraint: $inln=3$.

(errno $1$)

$monitr$ appears to have overwritten the solution vector.

Further integration will not be attempted.

(errno $2$)

At time $⟨value⟩$ the maximum number of allowed steps on this call was taken before reaching the next output point $tout=⟨value⟩$.

Maximum number of steps $=⟨value⟩$.

(errno $3$)

Too much accuracy requested for precision of the machine at time $⟨value⟩$.

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

(errno $7$)

The derivative evaluation procedure set the error flag $ires=3$ repeatedly despite attempts by the integrator to avoid this.

(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: $⟨value⟩$. Amount provided: $⟨value⟩$.

(errno $10$)

Larger integer workspace required.

Provided: $⟨value⟩$; required: $⟨value⟩$.

(errno $10$)

Not enough real store provided for sparse matrix solver.

Units of store needed: $⟨value⟩$. Amount provided: $⟨value⟩$.

(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 $⟨value⟩$.

(errno $15$)

The linear algebra setup function for sparse Jacobian was not called first.

Warns

NagAlgorithmicWarning

(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$)

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.

(errno $6$)

At time $⟨value⟩$, error weight $⟨value⟩$ became zero. Check the values of $atol$, $rtol$ and $itol$ supplied.

(errno $6$)

Weight number $i=⟨value⟩$ used in the local error test is too small. Check the values of $rtol$ and $atol$.

$atol[i-1]$ and $y[i-1]$ may both be zero.

Weight $i=⟨value⟩$.

(errno $11$)

$fcn$ set $ires=2$, which signals that the integration should terminate. $ires=⟨value⟩$ at time $⟨value⟩$.

(errno $12$)

A return was forced by setting $imon=-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 $rtol$ and $atol$ 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 $itask\ne 2$ or $5$.)

Notes

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

ivp_stiff_exp_sparjac is a general purpose function for integrating the initial value problem for a stiff system of explicit ordinary differential equations,

$$y^{\prime }=g(t,y)\text{.}$$

It is designed specifically for the case where the Jacobian $\frac{\partial g}{\partial y}$ is a sparse matrix.

Both interval and step oriented modes of operation are available and also modes designed to permit intermediate output within an interval oriented mode.

A typical calling program for ivp_stiff_exp_sparjac would involve calling the sparse matrix linear algebra setup function ivp_stiff_sparjac_setup(), the Backward Differentiation Formula (BDF) integrator setup function ivp_stiff_bdf(), its diagnostic counterpart ivp_stiff_integ_diag(), and the sparse linear algebra diagnostic function ivp_stiff_sparjac_diag() as well as the call to ivp_stiff_exp_sparjac. The linear algebra setup function ivp_stiff_sparjac_setup() and one of the integrator setup functions, ivp_stiff_bdf() or ivp_stiff_blend(), must be called prior to the call of ivp_stiff_exp_sparjac. Either or both of the integrator diagnostic function ivp_stiff_integ_diag(), or the sparse matrix linear algebra diagnostic function ivp_stiff_sparjac_diag(), may be called after the call to ivp_stiff_exp_sparjac. There is also a function, ivp_stiff_contin(), designed to permit you to change step size on a continuation call to ivp_stiff_exp_sparjac without restarting the integration process.