ddeabm_module Module

the number of (first order) differential equations to be integrated (>=0)

the expense of solving the problem is monitored by counting the number of steps attempted. when this exceeds maxnum, the counter is reset to zero and the user is informed about possible excessive work.

to define the system of first order differential equations which is to be solved. for the given values of x and the vector u(*)=(u(1),u(2),...,u(neq)), the subroutine must evaluate the neq components of the system of differential equations du/dx=df(x,u) and store the derivatives in the array uprime(*), that is, uprime(i) = du(i)/dx for equations i=1,...,neq.

a function for reporting the integration steps (optional). this can be associated in ddeabm_initialize or ddeabm_with_event_initialize, and is enabled using the mode argument in ddeabm_wrapper or ddeabm_with_event_wrapper

the user input rel tol

the user input abs tol

rel tol used internally

abs tol used internally

info array

how to choose the initial step h:

when initial_step_mode=3, the h value to use. (for the other modes, this will also contain the value used)

maximum step size magnitude. If <=0, there is no limit.

maximum order of the method (1 to ddeabm_max_order_limit)

replaced iwork(liw) in original code

replaced rwork(itstar) in original code

$2^i$ array

$| \gamma^{*}_i |$ array

(neq) contains the approximate derivative of the solution component y(i). in ddeabm, it is obtained by calling subroutine df to evaluate the differential equation using t and y(*) when idid=1 or 2, and by interpolation when idid=3.

(neq)

(neq)

(neq)

(neq)

(neq,max(max_order+4,6)). Columns 1:max_order+2 are the modified divided differences, and the last two are used for roundoff control.

(max_order)

(max_order)

(max_order)

(max_order)

(max_order)

(max_order+1)

(max_order+1)

(max_order)

the step size h to be attempted on the next step.

if the tolerances have been increased by the code (idid=-2), they were multiplied by eps.

contains the current value of the independent variable, i.e. the farthest point integration has reached. this will be different from t only when interpolation has been performed (idid=3).

(see info(4))

machine dependent parameter

machine dependent parameter

indicator for dsteps code

indicator for dsteps code

indicator for dsteps code

indicator for stiffness detection

indicator for intermediate-output

initialization indicator

counter for attempted steps

step counter for stiffness detection

termination flag

(max_order)

will be set in stop_integration. indicates an error and to abort the integration

state interpolation function

user can call this in df routine to stop the integration

call this after an interrupted task (idid = -1, -2, -3, -4) to continue the integration (equivalent to setting info(1)=1)

tolerance for root finding

time of last successful step (used to continue after a root finding)

state of last successful step. size (neq). (used to continue after a root finding)

time from which to resume event checking within the last step (the location of the last root found)

if the interval [t_resume, t_saved] of the last step has not yet been checked for events

event function: g(t,x)=0 at event

a function for determining if the root is bracketed. the default just checks for a sign change in the event function. the user can specify this to use other conditions.

state interpolation function

user can call this in df routine to stop the integration

call this after an interrupted task (idid = -1, -2, -3, -4) to continue the integration (equivalent to setting info(1)=1)

main routine for integration to an event (note that integrate can also still be used from ddeabm_class)

number of gfunc event equations

tolerance for root finding. this is a vector (size n_g_eqns)

time of last successful step (used to continue after a root finding)

state of last successful step. size (neq). (used to continue after a root finding)

time from which to resume event checking within the last step (the location of the last root found)

if the interval [t_resume, t_saved] of the last step has not yet been checked for events

event function: g(t,x)=0 at event

a function for determining if the root is bracketed. the default just checks for a sign change in the event function. the user can specify this to use other conditions.

state interpolation function

user can call this in df routine to stop the integration

call this after an interrupted task (idid = -1, -2, -3, -4) to continue the integration (equivalent to setting info(1)=1)

main routine for integration to an event (note that integrate can also still be used from ddeabm_class)

number of equations to be integrated

maximum number of integration steps allowed

derivative function dx/dt

relative tolerance for integration (see ddeabm)

absolution tolerance for integration (see ddeabm)

reporting function

how to choose the initial step h:

for initial_step_mode=3

maximum step size magnitude (>0). If not present, the step size is not limited.

maximum order of the method (1 to ddeabm_max_order_limit). If not present, ddeabm_default_max_order (12) is used. A lower value can improve stability on some problems. A higher value can reduce the number of steps for very tight tolerances when using higher precision (REAL128).

number of equations to be integrated

maximum number of integration steps allowed

derivative function $dx/dt$

relative tolerance for integration (see ddeabm)

absolution tolerance for integration (see ddeabm)

Event function $g(t,x)$. This should be a smooth function than can have values $<0$ and $\ge 0$. When $g=0$ (within the tolerance), then a root has been located and the integration will stop.

tolerance for the root finding

reporting function

root bracketing function. if not present, the default is used.

how to choose the initial step h:

for initial_step_mode=3

maximum step size magnitude (>0). If not present, the step size is not limited.

maximum order of the method (1 to ddeabm_max_order_limit). If not present, ddeabm_default_max_order (12) is used. A lower value can improve stability on some problems. A higher value can reduce the number of steps for very tight tolerances when using higher precision (REAL128).

number of equations to be integrated

maximum number of integration steps allowed

derivative function $dx/dt$

relative tolerance for integration (see ddeabm)

absolution tolerance for integration (see ddeabm)

number of event functions in g

Event function $g(t,x)$. This should be a smooth function than can have values $<0$ and $\ge 0$. When $g=0$ (within the tolerance), then a root has been located and the integration will stop.

tolerance for the root finding. This should be sized ng or 1 (in which case the value is used for all elements)

reporting function

root bracketing function. if not present, the default is used.

how to choose the initial step h:

for initial_step_mode=3

maximum step size magnitude (>0). If not present, the step size is not limited.

maximum order of the method (1 to ddeabm_max_order_limit). If not present, ddeabm_default_max_order (12) is used. A lower value can improve stability on some problems. A higher value can reduce the number of steps for very tight tolerances when using higher precision (REAL128).

time at which a solution is desired.

for some problems it may not be permissible to integrate past a point tstop because a discontinuity occurs there or the solution or its derivative is not defined beyond tstop. when you have declared a tstop point (see info(4)), you have told the code not to integrate past tstop. in this case any tout beyond tstop is invalid input. [not used if not present]

indicator reporting what the code did. you must monitor this integer variable to decide what action to take next.

Step mode: 1 - normal integration from t to tout, no reporting [default]. 2 - normal integration from t to tout, report each step.

Fixed time step to use for integration_mode=2. If not present, then default integrator steps are used. If integration_mode=1, then this is ignored.

max time at which a solution is desired. (if root not located, it will integrate to tmax)

for some problems it may not be permissible to integrate past a point tstop because a discontinuity occurs there or the solution or its derivative is not defined beyond tstop. when you have declared a tstop point (see info(4)), you have told the code not to integrate past tstop. in this case any tmax beyond tstop is invalid input. [not used if not present]

indicator reporting what the code did. you must monitor this integer variable to decide what action to take next. idid=1000 means a root was found. See ddeabm for other values.

value of the event function g(t,x) at the final time t (not set if the inputs are invalid: idid=-33)

Step mode: 1 - normal integration from t to tout, no reporting [default]. 2 - normal integration from t to tout, report each step.

Fixed time step to use for reporting and evaluation of the event function (the points are interpolated within each integrator step). If not present, then default integrator steps are used. The event function is also evaluated at the end of each integrator step.

to continue after a previous event location. This option can be used after a previous call to this routine has found an event, to continue integration to the next event without having to restart the integration. It will not report the initial point (which would have been reported as the last point of the previous call). The rest of the integration step in which the previous event was found is checked for events (using interpolation) before integration continues. Limitation: on that remaining interval, the event function is ignored (it is zero at the start), so a second root of it within the same integration step will not be found.

point at which solution is desired

interpolated state at tc

status flag:

max time at which a solution is desired. (if root not located, it will integrate to tmax)

for some problems it may not be permissible to integrate past a point tstop because a discontinuity occurs there or the solution or its derivative is not defined beyond tstop. when you have declared a tstop point (see info(4)), you have told the code not to integrate past tstop. in this case any tmax beyond tstop is invalid input. [not used if not present]

indicator reporting what the code did. you must monitor this integer variable to decide what action to take next. idid>1000 means a root was found for the (idid-1000)th g function. See ddeabm for other values.

value of the event functions g(t,x) at the final time t. Must be size ng (not set if the inputs are invalid: idid=-33)

Step mode: 1 - normal integration from t to tout, no reporting [default]. 2 - normal integration from t to tout, report each step.

Fixed time step to use for reporting and evaluation of the event functions (the points are interpolated within each integrator step). If not present, then default integrator steps are used. The event functions are also evaluated at the end of each integrator step.

to continue after a previous event location. This option can be used after a previous call to this routine has found an event, to continue integration to the next event without having to restart the integration. It will not report the initial point (which would have been reported as the last point of the previous call). The rest of the integration step in which the previous event was found is checked for events (using interpolation) before integration continues. Limitation: on that remaining interval, the event function(s) that triggered the previous event are ignored (they are zero at the start), so a second root of the same function within the same integration step will not be found. Roots of the other functions are found.

the number of (first order) differential equations to be integrated. (neq >= 1)

value of the independent variable. on input, set it to the initial point of the integration. the code changes its value. on output, the solution was successfully advanced to the output value of t.

the number of (first order) differential equations to be integrated.

the initial point of integration.

a value of the independent variable used to define the direction of integration. a reasonable choice is to set b to the first point at which a solution is desired. you can also use b, if necessary, to restrict the length of the first integration step because the algorithm will not compute a starting step length which is bigger than abs(b-a), unless b has been chosen too close to a. (it is presumed that dhstrt has been called with b different from a on the machine being used. also see the discussion about the parameter small.)

the vector of initial values of the neq solution components at the initial point a.

the vector of derivatives of the neq solution components at the initial point a. (defined by the differential equations in subroutine me%df)

the vector of error tolerances corresponding to the neq solution components. it is assumed that all elements are positive. following the first integration step, the tolerances are expected to be used by the integrator in an error test which roughly requires that abs(local error) <= etol for each vector component.

the order of the formula which will be used by the initial value method for taking the first integration step.

a small positive machine dependent constant which is used for protecting against computations with numbers which are too small relative to the precision of floating point arithmetic. small should be set to (approximately) the smallest positive double precision number such that (1.0_wp+small) > 1.0_wp. on the machine being used. the quantity small**(3.0_wp/8.0_wp) is used in computing increments of variables for approximating derivatives by differences. also the algorithm will not compute a starting step length which is smaller than 100.0_wp*small*abs(a).

a large positive machine dependent constant which is used for preventing machine overflows. a reasonable choice is to set big to (approximately) the square root of the largest double precision number which can be held in the machine.

work array of length neq which provide the routine with needed storage space.

work array of length neq which provide the routine with needed storage space.

work array of length neq which provide the routine with needed storage space.

work array of length neq which provide the routine with needed storage space.

appropriate starting step size to be attempted by the differential equation method.

point at which solution is desired

solution at xout

derivative of solution at xout

(max_order+1)

number of equations to be integrated

solution vector at x

independent variable

appropriate step size for next step. normally determined by code

local error tolerance

vector of weights for error criterion

logical variable set .true. for first step, .false. otherwise

step size used for last successful step

appropriate order for next step (determined by code)

order used for last successful step

logical variable set .true. when no step can be taken, .false. otherwise.

required for the interpolation subroutine dintp

required for the interpolation subroutine dintp

derivative of solution vector at x after successful step

required for the interpolation subroutine dintp

required for the interpolation subroutine dintp

required for the interpolation subroutine dintp

number of dsteps taken with size h

counter on attempted steps

2.*u where u is machine unit roundoff quantity

4.*u where u is machine unit roundoff quantity

required for the interpolation subroutine dintp

required for the interpolation subroutine dintp

required for the interpolation subroutine dintp

required for the interpolation subroutine dintp

required for the interpolation subroutine dintp

subroutine where error occurred

error message

error code

[Not used here] * -1: warning message (once), * 0: warning message, * 1: recoverable error, * 2: fatal error