# Implementing Events and Roots in Differential Equation Models.

### Description

An `event`

occurs when the value of a state variable is suddenly
changed, e.g. because a value is added, subtracted, or multiplied. The
integration routines cannot deal easily with such state variable
changes. Typically these events occur only at specific times. In
`deSolve`

, events can be imposed by means of an input data.frame,
that specifies at which time and how a certain state variable is altered,
or via an event function.

Roots occur when a root function becomes zero. By default when a root is found, the simulation either stops (no event), or triggers an event.

### Details

The `events`

are specified by means of argument `events`

passed to the integration routines.

`events`

should be a list that contains one of the following:

func: an R-function or the name of a function in compiled code that specifies the event,

data: a data.frame that specifies the state variables, times, values and types of the events. Note that the event times must also be part of the integration output times, else the event will not take place. As from version 1.9.1, this is checked by the solver, and a warning message is produced if event times are missing in times; see also

`cleanEventTimes`

for utility functions to check and solve such issues.time: when events are specified by an event function: the times at which the events take place. Note that these event times must also be part of the integration output times exactly, else the event would not take place. As from version 1.9.1 this is checked by the solver, and an error message produced if event times are missing in times; see also

`cleanEventTimes`

for utility functions to check and solve such issues.root: when events are specified by a function and triggered by a root, this logical should be set equal to

`TRUE`

terminalroot when events are triggered by a root, the default is that the simulation continues after the event is executed. In

`terminalroot`

, we can specify which roots should terminate the simulation.maxroot: when

`root = TRUE`

, the maximal number of times at with a root is found and that are kept; defaults to 100. If the number of roots >`maxroot`

, then only the first`maxroot`

will be outputted.ties: if events, as specified by a data.frame are "ordered", set to "ordered", the default is "notordered". This will save some computational time.

In case the events are specified by means of an **R** **function**
(argument `events$func`

),
it must be defined as: `function(t, y, parms, ...)`

.
`t`

is the current time point in the integration,
`y`

is the current estimate of the variables in the ODE system.
If the initial values `y`

has a `names`

attribute, the
names will be available inside `events$func`

. `parms`

is a
vector or list of parameters; `...`

(optional) are any other
arguments passed to the function via the call to the integration method.
The event function should return the y-values (some of which modified),
as a *vector*.

If `events$func`

is a string, this indicates that the events are
specified by a `function`

in compiled code. This function has as
arguments, the number of state variables, the time, and the state
variable vector. See package vignette "compiledCode" for more details.

In case events are specified by an **R**-function,
this requires either: input of the *time* of the events, a vector as
defined in `events$time`

OR the specification of a *root* function. In the
latter case, the model must be solved with an integration routine
with root-finding capability

The root function itself should be specified with argument `rootfunc`

.
In this case, the integrator is informed that the simulation it to be
continued after a root is found by
setting `events$root`

equal to `TRUE`

.

If the events are specified by a **data frame**
(argument `events$data`

), this should
contain the following columns (and in that order):

var the state variable

*name*or*number*that is affected by the eventtime the time at which the event is to take place; the solvers will check if the time is embraced by the simulation time

value the value, magnitude of the event

method which event is to take place; should be one of ("replace", "add", "multiply"); also allowed is to specify the number (1 = replace, 2 = add, 3 = multiply)

For instance, the following line

`"v1" 10 2 "add"`

will cause the value 2 to be added to a state variable, called `"v1"`

at
`time = 10`

.

From deSolve version 1.9.1 the following routines have **root-finding** capability:
lsoda, lsode, lsodes, and radau.
For the first 3 integration methods, the
root finding algorithm is based on the algorithm in solver LSODAR, and is
implemented in FORTRAN. For radau, the root solving algorithm is written
in C-code, and it works slightly different. Thus, some problems involving roots
may be more efficiently solved with either lsoda, lsode, or lsodes, while other
problems are more efficiently solved with radau.

If a root function is defined, but not an event function, then by default the solver will stop at a root. If this is not desirable, e.g. because we want to record the position of many roots, then a dummy "event" function can be defined which returns the values of the state variables - unaltered.

If roots and events are combined, and roots are found, then the output will have attribute
`troot`

which will contain the `times`

at which a root was found (and
the event trigerred).
There will be at most `events$maxroot`

such values. The default is 100.

See two last examples; also see example of `ccl4model`

.

### Author(s)

Karline Soetaert,

### See Also

forcings, for how to implement forcing functions.

lsodar, for more examples of roots

### Examples

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 | ```
## =============================================================================
## 1. EVENTS in a data.frame
## =============================================================================
## derivative function: derivatives set to 0
derivs <- function(t, var, parms) {
list(dvar = rep(0, 2))
}
yini <- c(v1 = 1, v2 = 2)
times <- seq(0, 10, by = 0.1)
eventdat <- data.frame(var = c("v1", "v2", "v2", "v1"),
time = c(1, 1, 5, 9) ,
value = c(1, 2, 3, 4),
method = c("add", "mult", "rep", "add"))
eventdat
out <- vode(func = derivs, y = yini, times = times, parms = NULL,
events = list(data = eventdat))
plot(out)
##
eventdat <- data.frame(var = c(rep("v1", 10), rep("v2", 10)),
time = c(1:10, 1:10),
value = runif(20),
method = rep("add", 20))
eventdat
out <- ode(func = derivs, y = yini, times = times, parms = NULL,
events = list(data = eventdat))
plot(out)
## =============================================================================
## 2. EVENTS in a function
## =============================================================================
## derivative function: rate of change v1 = 0, v2 reduced at first-order rate
derivs <- function(t, var, parms) {
list(c(0, -0.5 * var[2]))
}
# events: add 1 to v1, multiply v2 with random number
eventfun <- function(t, y, parms){
with (as.list(y),{
v1 <- v1 + 1
v2 <- 5 * runif(1)
return(c(v1, v2))
})
}
yini <- c(v1 = 1, v2 = 2)
times <- seq(0, 10, by = 0.1)
out <- ode(func = derivs, y = yini, times = times, parms = NULL,
events = list(func = eventfun, time = c(1:9, 2.2, 2.4)) )
plot(out, type = "l")
## =============================================================================
## 3. EVENTS triggered by a root function
## =============================================================================
## derivative: simple first-order decay
derivs <- function(t, y, pars) {
return(list(-0.1 * y))
}
## event triggered if state variable = 0.5
rootfun <- function (t, y, pars) {
return(y - 0.5)
}
## sets state variable = 1
eventfun <- function(t, y, pars) {
return(y = 1)
}
yini <- 2
times <- seq(0, 100, 0.1)
## uses ode to solve; root = TRUE specifies that the event is
## triggered by a root.
out <- ode(times = times, y = yini, func = derivs, parms = NULL,
events = list(func = eventfun, root = TRUE),
rootfun = rootfun)
plot(out, type = "l")
## time of the root:
troot <- attributes(out)$troot
points(troot, rep(0.5, length(troot)))
## =============================================================================
## 4. More ROOT examples: Rotation function
## =============================================================================
Rotate <- function(t, x, p )
list(c( x[2],
-x[1] ))
## Root = when second state variable = 0
rootfun <- function(t, x, p) x[2]
## "event" returns state variables unchanged
eventfun <- function(t, x, p) x
times <- seq(from = 0, to = 15, by = 0.1)
## 1. No event: stops at first root
out1 <- ode(func = Rotate, y = c(5, 5), parms = 0,
times = times, rootfun = rootfun)
tail(out1)
## 2. Continues till end of times and records the roots
out <- ode(func = Rotate, y = c(5, 5), parms = 0,
times = times, rootfun = rootfun,
events = list(func = eventfun, root = TRUE) )
plot(out)
troot <- attributes(out)$troot # time of roots
points(troot,rep(0, length (troot)))
## Multiple roots: either one of the state variables = 0
root2 <- function(t, x, p) x
out2 <- ode(func = Rotate, y = c(5, 5), parms = 0,
times = times, rootfun = root2,
events = list(func = eventfun, root = TRUE) )
plot(out2, which = 2)
troot <- attributes(out2)$troot
indroot <- attributes(out2)$indroot # which root was found
points(troot, rep(0, length (troot)), col = indroot, pch = 18, cex = 2)
## Multiple roots and stop at first time root 1.
out3 <- ode(func = Rotate, y = c(5, 5), parms = 0,
times = times, rootfun = root2,
events = list(func = eventfun, root = TRUE, terminalroot = 1))
## =============================================================================
## 5. Stop at 5th root - only works with radau.
## =============================================================================
Rotate <- function(t, x, p )
list(c( x[2],
-x[1],
0 ))
## Root = when second state variable = 0
root3 <- function(t, x, p) c(x[2], x[3] - 5)
event3 <- function (t, x, p) c(x[1:2], x[3]+1)
times <- seq(0, 15, 0.1)
out3 <- ode(func = Rotate, y = c(x1 = 5, x2 = 5, nroot = 0),
parms = 0, method = "radau",
times = times, rootfun = root3,
events = list(func = event3, root = TRUE, terminalroot = 2))
plot(out3)
attributes(out3)[c("troot", "nroot", "indroot")]
## =============================================================================
## 6 Event in R-code, model function in compiled code - based on vode example
## =============================================================================
times <- 1:365
Flux <- cbind(times, sin(pi*times/365)^2) # forcing function
# run without events
out <- ode(y = c(C = 1), times, func = "scocder", parms = c(k=0.01),
dllname = "deSolve", initforc = "scocforc", forcings = Flux,
initfunc = "scocpar", nout = 2, outnames = c("Mineralisation", "Depo"))
# Event halves the concentration
EventMin <- function(t, y , p) y/2
out2 <- ode(y = c(C = 1), times, func = "scocder", parms = c(k=0.01),
dllname = "deSolve", initforc = "scocforc", forcings = Flux,
initfunc = "scocpar", nout = 2, outnames = c("Mineralisation", "Depo"),
events = list (func = EventMin, time = c(50.1, 200, 210.5)))
plot(out, out2)
``` |