Events for dsolve[numeric]
Round-off and simple triggers
Discrete variables associated to triggers
Event triggering (details)
dsolve(..., numeric, events=eventlist)
dsolve(..., numeric, events=eventlist, discrete_variables=discvars)
dsolve(..., numeric, events=eventlist, event_maxiter=miter)
dsolve(..., numeric, events=eventlist, event_pre=mode)
dsolve(..., numeric, events=eventlist, event_iterate=form)
dsolve(..., numeric, events=eventlist, event_initial=truefalse)
list of events
set of discrete variables
maximum number of event iterations for a single trigger
operating mode of pre, one of iteration (default), initial, assignment, or change
iteration control criteria, one of trigger or prechange
Events, in combination with discrete variables, are a mechanism that can be used in dsolve/numeric to handle stopping criteria, reset conditions, zero order holds, and most other events that occur in ODE and DAE system simulation. Note that events are restricted to operations that act on variable values, and cannot be used to change the essential form of the ODE system or exchange equations.
These are available in dsolve,numeric for the default non-stiff (rkf45, rkf45_dae, ck45, ck45_dae) and stiff (rosenbrock and rosenbrock_dae) solvers. Note that events replace the stop condition mechanism previously available for these solvers.
The sections below describe the behavior and features of the event implementation in dsolve[numeric].
Events are specified by the events=[event1,event2,...] option in the call to dsolve[numeric], where each event is a list of [trigger,action] pairs, where the trigger describes the trigger of the event, and the action describes the action to perform when the event is triggered.
The following are the options for the form of the trigger part of the description:
y(t): simple root finding trigger
f(t,y(t)): root finding trigger
[f(t,y(t)),c(t,y(t))<0]: root finding with a conditional trigger
[f(t,y(t)), And(c1(t,y(t)) > 0, c2(t,y(t)) > 0)]: root finding with a compound conditional trigger
[0, c(t,y(t)) < 0]: discrete event with a conditional trigger
[0, And(c1(t,y(t)) > 0, c2(t,y(t)) > 0)]: discrete event with a compound conditional trigger
f(t,y(t)) = lo..hi: range trigger, that fires once the expression goes outside the range
t=start: timed trigger that fires when t=start
t=[start, δ]: timed trigger that fires when t=..., start-δ, start, start+δ, start+2*δ, ...
Any events where the trigger depends on continuous variables fire when the trigger becomes satisfied during the integration process (on the so-called leading edge). This means that if a continuous condition is satisfied at the initial point, it will not fire (not by default, though this can be changed via fireinitial described below), but will fire if the condition is satisfied at any time past the initial point.
This is not the case for timed triggers, which do fire at the initial point, so, for example, the trigger t-t0 will not fire at the initial point, while the trigger t=t0 will.
Discrete events (where either the rootfind part is zero, or the event contains only discrete variables) can fire only during an event iteration that has been initiated by a continuous trigger (see the following sections), or if it has been specified that events iteration should occur each time the integrator resumes (the event_initial option).
The following are the options for the form of the action part of the description:
none: do nothing, but this forces a point to be included at the trigger point for an adaptive plot mode (for example, the range option of adaptive odeplot output), and induces an event iteration which can trigger discrete events.
halt: halt integration with a stop message when triggered
tobegin: (iteration) jump to the start of the event iteration, incrementing the iteration counter and updating pre.
toend: (iteration) jump to the end of the event iteration, performing a termination check before iterating.
breakiteration: (iteration) halt event iteration, resuming continuous integration.
delayhalt: (iteration) schedule a halt of the integration on completion of event iteration.
nomark: (iteration) do not tag event as having executed in an event iteration (only useful if event_iterate=trigger and described in more detail in the discussion of event iteration).
u(t) = -u(t): execute statement when triggered
[u(t) = -u(t), y(t) = -y(t)]: execute statements when triggered
[If(y(t) < 0, y(t) = -y(t), halt)]: execute statements conditionally when triggered
The conditional part of a trigger can use the inert (capitalized) form of the Not, And, Or, Xor, and Implies logical operations. The action can use these, the If conditional, the pre operator (see below), halt, tobegin, toend, breakiteration, delayhalt, nomark and none; within the action, equations are interpreted as assignments. Note that delayhalt and nomark must be the last statement in any embedded code, as when executed they force an immediate return to the controlling routine (which is obviously also true for halt, tobegin, toend, and breakiteration). Note that only none cannot be embedded in code, and only nochange must be embedded in code. Note that the action can also use temporary variables in the associated computation, so long as the names are distinct from the dependent and independent variable names for the problem. It can also use derivative values up to one order lower than the order of the derivative in the problem (that is, for a second order ODE in y(t), up to diff(y(t), t) can be used). As an example, a valid action for a second order ODE in y(t) could be [t1 = sqrt(y(x)^2 + diff(y(t),t)^2), x(t) = x(t)*t1, y(t) = y(t)*t1].
When a root-finding trigger fires, the value of the expression in the root finder is only an approximate zero due to numerical round-off. What this means is that f(t,y(t)) may evaluate, with respect to the current solution, to a value slightly above or slightly below zero. This may not be desirable for some computations where it is critical that the value of the trigger expression is either exactly zero, as approaching from the left, or as just past the trigger point. This is primarily desired to be able to apply events for an ODE system that has been separated into disjoint cases dependent on the values of particular triggers (in which case you always want to use a form that provides the values just past the trigger point).
A mechanism has been provided for this, but only for simple root finding triggers (i.e. the expression for the root find must be a dependent variable of the problem). In the event that an expression of the dependent and independent variables is required, an index-1 relation of the form newvar(t) = f(t,y(t)) can be added, and the trigger can then be based on newvar(t).
The side function, in place of the root finding dependent variable y(t) can be used in the forms side(y(t),pre) or side(y(t),post) to obtain the value just prior, at, or just after the root.
For example, for a root finding trigger side(y(t),post) where y(t) is initially negative, at the time at which the event is triggered the value of y(t) will be a small positive quantity.
Discrete variables are a way of maintaining discrete data as part of a numerical ODE/DAE computation in dsolve[numeric], and can be as simple as a boolean flag that indicates that a switch is either in the on or off position, an integer value that counts how many times a particular criterion is met, or a sampling of continuous data, as would be required for a zero-order hold.
The values of discrete variables can only ever be changed in events, since they are not continuous (and hence will not have an associated differential equation). As a result, derivatives of discrete variables cannot appear in the input to dsolve[numeric].
Discrete variables depend upon the independent variable (since their value can change within event handling), so must be declared as functions of the independent variable. These can be specified to dsolve[numeric] by using the option discrete_variables=[v1,v2,...], where for any i, vi can be a dependent variable (for example, u(t)) or a dependent variable with a declared type (u(t)::integer). Only the types boolean, integer and float are currently available. Note: it is strongly recommended that any discrete variables that are not of type float (the default) be explicitly declared.
Discrete variables must also be given an initial value (even if it is not needed) in the input system.
As an example of the use of a discrete variable for a zero order hold on the continuous dependent variable y(t) sampled at 12 second intervals, specify discrete_variables=[u(t)::float], events=[[t = [0,1/2], u(t) = y(t)]] and an initial condition on u(t).
Note that discrete variables that are declared as boolean can be used directly in the conditionals of any event construct (for example, if u(t) was declared as boolean, then a conditional can be of the form If(u(t),...)).
Another mechanism for handling events that is most useful in cases where a trigger expression asymptotically approaches (i.e. touches) zero is the ability to associate a discrete (boolean) variable to a trigger.
This is done using the requires option, and is only available for simple triggers (i.e. basic root finding triggers, supporting increasing, but with no additional options, and no trigger condition). Also note that increasing triggers must occur in pairs, where one trigger depends on the requires variable, and the other depends on the Not of the requires variable, and the sum of the arguments to the corresponding increasing (or the difference of the increasing and decreasing arguments) must be a constant.
The value of the requires option specifies the variable associated to the trigger, and the value that the variable should have when the trigger fires in an increasing direction (which is the opposite to the value the variable will have after the firing of the trigger). The requires variable will be automatically updated to the opposite value as specified for the trigger firing prior to execution of the corresponding action, and the action cannot update the value of the requires variable (though it can use it).
As a simple example, consider a problem with continuous dependent variable y⁡t, where the discrete variable py⁡t is true (1) if 0<y⁡t and false (0) if y⁡t<0 (and note that the value when y⁡t=0 can be either). The event specification that can be used to tie py⁡t into y⁡t is simply y⁡t,none,requires=¬py⁡t, and the problem must start with an initial condition that is consistent (i.e. if 0<y⁡t then py⁡t=1). This trigger will then fire when y⁡t passes through zero in the positive direction with py⁡t=0, and will also fire when y⁡t passes through zero in the negative direction with py⁡t=1.
The same behavior can be accomplished with a somewhat less efficient pair of events: increasing⁡y⁡t,requires=¬py⁡t and increasing⁡−y⁡t,requires=py⁡t. The first event sets py⁡t=1 when y⁡t passes through zero in the positive direction with py⁡t=0, and the second event sets py⁡t=0 when y⁡t passes through zero in the negative direction with py⁡t=1.
The latter mechanism is most useful when providing a band to avoid spurious events triggering as a result of round off (so one would change the trigger expressions to increasing⁡y⁡t−ε and increasing⁡−y⁡t−ε respectively to provide an ε width band about the trigger expression).
Events that are based on continuous criteria (such as root finding or events that fire for a fixed time or sequence of times) can fire only once at any given time. This is not the case for events that trigger on discrete variables.
Events that are based on discrete variables can be utilized to code an algorithm (of sorts) into the event process, but note that for the events to be initiated, the condition must be satisfied at a point after the trigger of a continuous event (for example, a sign change or a certain time being reached).
If event_initial is set to true (the default is false) event iteration also occurs at the initial point, and at any point when resuming integration after a halt action has been performed.
The way in which sequences of events run is very relevant, and quite simply, they run in the order that they appear in the input events specification. If any event has triggered during an iteration, or the value of variables has changed (which criteria is in use is controlled by the event_iterate option), then the events are iterated through again, from the beginning, until event_maxiter is reached. The event_maxiter option is provided to prevent infinite event loops (which are surprisingly easy to create), and can be set to any positive integer value, but cannot be disabled.
There are two main ways in which event iteration can terminate. One of these is when a pass is made over all events, and none of the events trigger, and this is the default termination type (specified by event_iterate=trigger). The other termination type (event_iterate=prechange) compares the values of all variables (discrete and continuous) to their pre values (see below) after each complete pass over the events, and if no differences occur, the iteration terminates.
There are also several control actions that only make sense from within an event iteration:
The tobegin action exits the current event, and goes back to the start of the iteration (updating pre if relevant) while the toend action exits the current event, and jumps to the end of the iteration. Note that the only difference between these two behaviors is that when event_iterate=prechange, the toend action performs a termination check, while the tobegin action does not.
The nomark action is only useful for event_iterate=trigger, and exits the current event telling the control code to act as though the event did not execute, so as to not affect the trigger termination condition (e.g. if event_iterate=trigger, and the only action to run in an iteration ended with nomark, then the iteration would terminate).
The delayhalt action exits the current action and schedules a halt (a stop condition) for when the iteration terminates. This is in contrast to the halt action which exits the current action and halts immediately, exiting both the iteration and the continuous integration.
The breakiteration action exits the current event and immediately halts the event iteration, resuming continuous integration (unless a delayhalt was scheduled).
Another common mechanism required in event processing for simulation is the ability to examine a prior value of a variable. This can be used to check if a variable has changed, or a boolean variable has changed from true to false or false to true.
The mechanism available for this in dsolve[numeric] event processing is the pre operator. Note that pre has no meaning outside of event processing, as the value of pre is keyed by other events. The condition to detect if the value of, for example, u(t) has changed during event processing can be written as u(t) <> pre(u(t)).
There are some subtleties in dealing with pre, and these have to do with the point in the computation that the pre value refers to. There are four options for this, controlled by the event_pre optional argument:
iteration (default): pre(u(t)) refers to the value of u(t) at the start of the current iteration.
initial: pre(u(t)) refers to the value of u(t) at the start of the events at the current point (this value will never change during event iteration).
assignment: pre(u(t)) refers to the value of u(t) prior to the last assignment to u(t) since the start of the event iteration at the current point.
change: pre(u(t)) refers to the value of u(t) prior to the last nontrivial assignment to u(t) since the start of the events at the current point.
There are several simple aliases that are related to pre variables for use in the conditional part of the trigger of an event, and these are as follows:
change(u(t)): returns true if u(t) <> pre(u(t)), false otherwise.
risingedge(u(t)): returns true if pre(u(t)) is false and u(t) is true, or if pre(u(t)) < u(t) for non-boolean u(t).
fallingedge(u(t)): returns true if pre(u(t)) is true and u(t) is false, or if pre(u(t)) > u(t) for non-boolean u(t).
edge(u(t)): same as risingedge(u(t)).
In addition to the typical description of an event [trigger, action], additional options can be specified to control the behavior of events. Two options are currently available:
retrigger=truefalse: This option specifies whether an event should be allowed to trigger more than once at a specific time. By default this is false for continuous events, and true for discrete events. Warning: setting this to true for a continuous event will result in an event that can never be exited, which means the integration must halt at that point. This cannot be set to true for a time event.
fireinitial=truefalse: This option specifies whether an event is allowed to trigger at the initial point. By default this is false for root-finding events, true for time events, and irrelevant for discrete events. This option may be needed for problems where the event is very close to, but not at, the initial point of the problem.
requires=discrete_variable: This option allows association of a discrete variable to a trigger expression, and is described in detail earlier in the Discrete variables associated to triggers section.
In addition to these options, a root-finding trigger (given as an expression) can be restricted to roots having a positive rate of change with respect to the independent variable (increasing) or a negative rate of change with respect to the independent variable (decreasing) by wrapping the root-finding trigger expression in a function with the name increasing or decreasing respectively. For example, the trigger increasing(y(t)) will only fire when y(t) passes through zero at a point t0 with D(y)(t0)>0.
The general dsolve option event_doublecross=tol provides a tolerance such that any root-finding events that are satisfied twice within an interval relative to tol will be ignored, and also that the solver is allowed to step an amount relative to tol past the actual trigger point so that the overall state of the solution in a halt is just past the trigger. This is most useful for cases where the events are being used to transition from one configuration to another, and use of the values just within the next configuration are needed to detect the correct configuration.
There are many features available for dealing with events for procedure-form numerical ODE/DAE solutions that are an extension of the features previously available for stop conditions. This includes the abilities to query which event triggered a stop, continue numerical integration by clearing a continuous stop, and enable/disable events. The first example below demonstrates some of these features.
The features previously available through stop conditions include:
eventstop: query the last event that fired (0 if no event has fired)
eventclear: clear the last event allowing integration to continue
eventcount: the total number of events that have fired since integration started at the initial point
eventdisable=set: set of events to be disabled for further solution
eventenable=set: set of events to be enabled for further solution
eventstatus: query enabled/disabled events
In addition to the existing stop condition features a new event control and query feature, eventfired, has been added.
A call to a procedure-form dsolve numeric procedure with the argument eventfired=[n1,n2,n3], where n1,n2,n3 are numerical references to the location of a root-finding or time-interval event in the initial events specification list, will return a list [v1,v2,v3] of the times at which the event most recently fired. If the event has never fired, a time on the opposite side of the initial point (as the current integration direction) will be returned instead.
The calling sequence eventfired=[n1=v1,n2=v2,n3=v3] can be used to specify when each event has last fired (and as a result the events will not fire until one of the triggers fire pas the relevant specified value). When this form is used, the prior values of eventfired are returned in a list.
Any attempt to set or query eventfired for discrete or range conditions, or for an index outside of the bounds of the initial events list, will result in an error. Any attempt to set or query eventfired for a root-finding event with retrigger set to true will result in a warning as eventfired is not used for retriggerable events.
Until an integration direction is chosen (i.e. by requesting a solution value), eventfired cannot be used (as a point prior to the initial condition cannot be determined). It is important to note that the integration direction is an absolute direction relative to the initial value for the problem, not relative to the last computed point (dsolve procedures only ever integrate away from the initial condition for a problem). In order to accomplish this without requiring an initial computation, one can call the procedure with direction=1 to set a positive integration direction, or direction=-1 to set a negative integration direction. Note that setting the integration direction for a procedure that has already been used to compute a solution clears the prior computation, forcing recomputation from the initial point for any subsequent solution requests. If a direction cannot be determined (i.e. this is not done, and no computation has been performed) an error will result when calling with eventfired.
Consider a simple harmonic oscillator with an event that halts the computation whenever the solution passes through zero:
dsn ≔ dsolve⁡ⅆ2ⅆt2⁢y⁡t+y⁡t=0,y⁡0=0,D⁡y⁡0=1,numeric,events=y⁡t,halt
dsn ≔ procx_rkf45...end proc
The exact solution to this equation is y⁡t=sin⁡t, so we expect the computation to halt at approximately t=Pi:
This halt can be cleared, and we can resume computation until we either reach 10 or the next event:
In order to reset the integrator to repeat the computation again, it must be called with a time value that is slightly larger than the initial time, and the computation can then be repeated:
Now consider the same simple harmonic oscillator with two events that correspond to the solution passing through zero in an upwards direction and a downwards direction:
dsn ≔ dsolve⁡ⅆ2ⅆt2⁢y⁡t+y⁡t=0,y⁡0=0,D⁡y⁡0=1,numeric,events=y⁡t,0<ⅆⅆt⁢y⁡t,halt,y⁡t,ⅆⅆt⁢y⁡t<0,halt
This is effectively the same as the first example, except now the halt for the solution passing through zero from below and from above are controlled by different events.
You can query which events are enabled/disabled by using eventstatus.
Now if you integrate in the positive direction, you get a halt from event #2 at approximately Pi.
You can then clear this event, and continue integration, which causes a halt from event #1 at approximately 2*Pi.
Now, if you clear this event and continue integration, you would expect a halt from event #2 at approximately 3*Pi, but if you disable this event, this halt will not occur.
Further, you can disable the remaining event, and integrate to any time.
When events are re-enabled, the problem is reset to initial values, as though the procedure has never been called, and events have never been cleared.
Note that when multiple events are present, and a halt is cleared, the event processing will resume immediately after the last executed event (so if event #2 triggered a halt, and was cleared, event processing will resume at event #3).
Next, consider a bouncing ball problem where the collision with the ground (at y=0) is assumed to be lossless, and upon collision, the velocity in the y direction is simply reversed. This can be posed (for an initial position of (0,1) and an initial velocity of (1,0)) as follows:
dsn ≔ dsolve⁡ⅆ2ⅆt2⁢x⁡t=0,x⁡0=0,D⁡x⁡0=1,ⅆ2ⅆt2⁢y⁡t=−9.8,y⁡0=1,D⁡y⁡0=0,numeric,events=y⁡t,ⅆⅆt⁢y⁡t=−ⅆⅆt⁢y⁡t,range=0..3
If you then plot the position of the ball using refinement, you can see the behavior of the ball as a function of time.
As a second example, consider a problem in 2 dependent variables where you want to swap the values of variables once a condition is met. This is a little tricky, since the values in the statements take effect immediately.
dsys ≔ ⅆⅆt⁢x1⁡t=1,x1⁡0=−1,ⅆⅆt⁢x2⁡t=−1,x2⁡0=1
Now one solution would be to specify the event action as [x1(t) = x2(t), x2(t) = x1(t)], but this would not give the desired result, since the value of x1(t) will be changed by the first statement before the value of x2(t) can be set in the second statement, and both variables would then have the value that x2(t) had at the beginning of the event.
To solve this issue, you need to add a temporary variable to allow the values to be swapped as [t1 = x1(t), x1(t) = x2(t), x2(t) = t1], which is done below.
dsn ≔ dsolve⁡dsys,numeric,events=x1⁡t=1,t1=x1⁡t,x1⁡t=x2⁡t,x2⁡t=t1,range=0..6
Download Help Document