R help - Apr 2013 - lsoda question from deSolve package

Dear List,
Wonder if you have some thoughts on the following question using lsoda in
desolve:
I have the following data and function:
require(deSolve)
times <- c(0:24)
tin  <- 0.5
D <- 400
V    <- 26.3
k <-0.056
k12  <- 0.197118
k21  <- 0.022665
yini <- c(dy1 = 0,dy2 = 0)
 events <- data.frame(var = "dy1",time = c(10,15),value =
c(200,100),method = "add")
pkmod <- function(t, y, p) {
  if (t < tin) R <- (D/tin) else R <- 0
  dy1 <- R - (p["k12"]+p["k"])* y[1] +
p["k21"]*y[2]
  dy2 <- p["k12"]* y[1] - p["k21"] *y[2]
  list(c(dy1, dy2))
}
 
p <- c(k=k, k12=k12, k21=k21, V=V)
result <- lsoda( yini, times, pkmod, p,rtol = 1e-08, atol = 1e-08,events =
list(data = events))
result <- data.frame(result)
names(result) <- c("time","Amount1", "Amount2")
plot(result$Amount1 ~ result$time, type="b", main="Central",
xlab="time", ylab="Amount")
plot(result$Amount2 ~ result$time, type="b",
main="Peripheral", xlab="time", ylab="Amount")
What I would like to do is to set up the events for the differential equations
based on the following logic:
The way this model is written currently is so that it adds a D of 500 over 0.5
hours to the system at a rate R starting at time zero. The input is into state
“dy1”. Then, based on the “events” data frame created, we will add another 200
and 100 (value=c(200,100)) at times 10 hours and 15 hours (time = c(10,15)),
respectively, and instantaneously to the system via “dy1”. This is how far I
got, but this is not exactly what I want:-) What I would like to achieve is to
add the 200 to the system over 0.6 hours, ie: an input at a rate of R1 =200/0.6
that starts at hour 10 and lasts for 0.6 hours, and to add the 100 to the system
over 0.75 hours, ie: an input at a rate of R2 =100/0.75 that starts at hour 15
and lasts for 0.75 hours, as opposed to adding it instantaneously, which is how
it is happening now to the best of my understanding. The trick is probably in
the way the events data frame is set up, but I cannot get my head around it.
Your help
 is greatly appreciated,
 
Sincerely
 Andras 
	[[alternative HTML version deleted]]

Hi Andreas,

adding a varying amount of something over time is even easier than 
implementing an "immediate" event. It can be implemented as a
"forcing",
see ?forcings and example below.

Hope it helps,

Thomas


require(deSolve)

pkmod <- function(t, y, p) {
   inp <- forc(t)

   if (t < tin) R <- (D/tin) else R <- 0
   dy1 <- R - (p["k12"]+p["k"])* y[1] +
p["k21"]*y[2]   + inp
   dy2 <- p["k12"]* y[1] - p["k21"] *y[2]
   list(c(dy1, dy2), inp=inp)
}

## an example, not your data
intimes <- c(0, 2, 5, 10, 24)
input   <- c(0, 200, 100, 0, 0)
forc <- approxfun(intimes, input, method="constant")

times <- seq(0, 24, 0.1)
tin  <- 0.5
D    <- 400
V    <- 26.3
yini <- c(Central = 0, Peripheral = 0)
p    <- c(k=0.05, k12=0.197118, k21=0.022665, V=26.3)

result <- lsoda(yini, times, pkmod, p, rtol = 1e-08, atol = 1e-08)
plot(result)


-- 
Thomas Petzoldt
Technische Universitaet Dresden
Faculty of Environmental Sciences
Institute of Hydrobiology
01062 Dresden, Germany
http://tu-dresden.de/Members/thomas.petzoldt