R help - Apr 2014 - Difficulty coding time-forced functions in deSolve

Hi all,

I'm trying to use deSolve to solve a series of differential equations
with rk4 mimicking a SEIR model,  while including an event/function
that is not solely time-dependent.

Explicitly:
I want to introduce vaccination 7 days after the proportion of I2/N2
reaches 0.01.


Here is the code I am using:


require(deSolve)

require(sfsmisc)


SEIR <- function(t, x, p) {

with(as.list(c(x,p)),{





dS<-b*N-d*S-beta*S*I/N-V(I, N, t)*S


dE<- -d*E+beta*S*I/N -epsilon*E-V(I, N, t)*E


dI<- -d*I+epsilon*E-gamma*I-mu*I-V(I, N, t)*I


dR<--d*R+gamma*I+V(I, N, t)*S+V(I, N, t)*E+V(I, N, t)*I


dN<-dS+dE+dI+dR


list(c(dS, dE, dI, dR, dN))




})

}

V <-function(I, N, t) {ifelse(t >=8 & I[t-7]/N[t-7]>0.01, 0.25, 0)}

num_years <- 10.0

time_limit <-num_years*365.00


Ni <-1.0E3

b <-1/(10.0*365)

d <-b

beta <-0.48

epsilon <-1/4

gamma <-1/4

mu <--log(1-0.25)*gamma


parms <-c(Ni=Ni, b=b, d=d, beta=beta, epsilon=epsilon, gamma=gamma, mu=mu)

xstart <-c(S=999, E=0, I=1, R=0, N=1000)


 times   <- seq(0.0, time_limit, 1.0)

        tol <- 1e-16

        my.atol <- rep(tol,5)

        my.rtol <- 1e-12

 out_rk4 <-  as.data.frame(rk4(xstart, times, SEIR, parms))



  outfilename <- paste("Basic SEIR.csv")



   write.csv(out_rk4,file=outfilename,row.names=FALSE, col.names=FALSE)



If I remove function V and the associated parts within the
differential equations, the model runs just fine. If I define V as
V<-function(I, N) {ifelse(I/N >0.01, 0.25, 0)
the model functions just fine.

Any pointers as to how I can code a function that relies on solutions
from previous time steps?

thank you in advance,
Aimee.

On 4/2/2014 5:51 AM, Aimee Kopolow wrote:
> Any pointers as to how I can code a function that relies on solutions
> from previous time steps?
Such a system would be called a delay differential equation (DDE). It 
can be solved with the dede function, see ?dede for details.

However if you want to model something like this:

 > Explicitly:
 > I want to introduce vaccination 7 days after the proportion of I2/N2
 > reaches 0.01.


Than this is called "root finding", that can be combined with events, 
see example "EVENTS triggered by a root function" in ?events.

More can be found in the papers listed at:

http://desolve.r-forge.r-project.org

... or you may consider to ask the R-sig-dynamic models mailing list:

https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models

Hope it helps

Thomas Petzoldt

Dr. Thomas Petzoldt
Technische Universitaet Dresden
Faculty of Environmental Sciences
Institute of Hydrobiology
01062 Dresden, Germany

http://tu-dresden.de/Members/thomas.petzoldt

All,

I'm getting this:
> sprintf("%.17f", 0.8)
[1] "0.80000000000000004"

Where does the `4` at the end come from? Shouldn't it be zero at the
end? Maybe I'm missing something.
> sessionInfo()
R version 3.0.2 (2013-09-25)
Platform: x86_64-redhat-linux-gnu (64-bit)

locale:
 [1] LC_CTYPE=en_US.utf8       LC_NUMERIC=C
LC_TIME=en_US.utf8
 [4] LC_COLLATE=en_US.utf8     LC_MONETARY=en_US.utf8
LC_MESSAGES=en_US.utf8
 [7] LC_PAPER=en_US.utf8       LC_NAME=C                 LC_ADDRESS=C

[10] LC_TELEPHONE=C            LC_MEASUREMENT=en_US.utf8
LC_IDENTIFICATION=C

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base


Thanks,

M

Hi,

does the following help you?

Thomas


require(deSolve)

SEIR <- function(t, x, p) {

   if (t < 7)
     xlag <- x
   else
     xlag <- lagvalue(t - 7)

   V <- ifelse(xlag[3]/sum(xlag) > 0.01, 0.25, 0)

   ## uncomment the following for printing some internal information
   #cat("t=", t, " -- ")
   #cat(xlag, " -- ", V, "\n")

   N <- sum(x)
   with(as.list(c(x,p)),{
     dS<-  b*N - d*S - beta*S*I/N - V*S
     dE<- -d*E+beta*S*I/N - epsilon*E - V*E
     dI<- -d*I + epsilon*E - gamma*I - mu*I - V*I
     dR<- -d*R + gamma*I + V*S + V*E + V*I
     list(c(dS, dE, dI, dR), N = unname(N), V = unname(V))
   })
}

num_years  <- 1
time_limit <- num_years*365.00

b       <- 1/(10.0*365)
d       <- b
beta    <- 0.48
epsilon <- 1/4
gamma   <- 1/4
mu      <- -log(1-0.25)*gamma
parms   <- c(b=b, d=d, beta=beta, epsilon=epsilon, gamma=gamma, mu=mu)

xstart <-c(S=999, E=0, I=1, R=0)
times   <- seq(0.0, time_limit, 1.0)

out <-  dede(xstart, times, SEIR, parms, rtol=1e-8, atol=1e-8)

plot(out)