deSolve (version 1.10-2)

events: Implementing Events 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.

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. The function should return the modified y-values, 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.

Also, if events are specified by an R-function, this requires either input of the time of the events, a vector 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 (note that from version 1.9.1 the following routines have root-finding capability: lsoda, lsode, lsodes, radau). 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 roots are found, then the output will have attribute troot which will contain the times at which a root was found. There will be at most events$maxroot such values. The default is 100. See two last examples; also see example of ccl4model. If the evenst are specified by a **data.frame** (argument events$data), this should contain the following columns (and in that order):

1. var

{the state variable name or number that is affected by the event} time{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) }

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:

1. 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 } 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. }

See Also

forcings, for how to implement forcing functions.

Examples

## =============================================================================
## 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 = 1:9) )
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)))