These are just warnings from the package telling you what is and is not included in your plot:
## =============================================================================
## 1. EVENTS in a data.frame
## =============================================================================
library(deSolve)
## 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
#> var time value method
#> 1 v1 1 1 add
#> 2 v2 1 2 mult
#> 3 v2 5 3 rep
#> 4 v1 9 4 add
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
#> var time value method
#> 1 v1 1 0.60130622 add
#> 2 v1 2 0.26361379 add
#> 3 v1 3 0.96842125 add
#> 4 v1 4 0.02222490 add
#> 5 v1 5 0.03510969 add
#> 6 v1 6 0.30194760 add
#> 7 v1 7 0.55619074 add
#> 8 v1 8 0.49451620 add
#> 9 v1 9 0.42082526 add
#> 10 v1 10 0.24796604 add
#> 11 v2 1 0.93000955 add
#> 12 v2 2 0.46625803 add
#> 13 v2 3 0.50769950 add
#> 14 v2 4 0.11666074 add
#> 15 v2 5 0.64572750 add
#> 16 v2 6 0.74749996 add
#> 17 v2 7 0.53995838 add
#> 18 v2 8 0.46295596 add
#> 19 v2 9 0.69547626 add
#> 20 v2 10 0.62424576 add
out <- ode(func = derivs, y = yini, times = times, parms = NULL,
events = list(data = eventdat))
#> Warning in checkevents(events, times, Ynames, dllname): Time of events
#> ('time' column of 'events') was not ordered.
plot(out)
## =============================================================================
## 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")]
#> $troot
#> [1] 0.7853982 3.9269908 7.0685835 10.2101761 13.3517688
#>
#> $nroot
#> [1] 5
#>
#> $indroot
#> [1] 1 1 1 1 1
## =============================================================================
## 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)))
#> Warning in checkevents(events, times, Ynames, dllname): Not all event times
#> 'events$time' are in output 'times' so they are automatically included.
#> Warning in checkevents(events, times, Ynames, dllname): Some time steps
#> were very close to events - only the event times are used in these cases.
plot(out, out2)
Created on 2018-10-16 by the reprex package (v0.2.1.9000)