[R] lsoda question from deSolve package

[Andras Farkas]

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]]


______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.