[R] Difficulty coding time-forced functions in deSolve

Aimee Kopolow alj27 at georgetown.edu
Wed Apr 2 05:51:34 CEST 2014

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.

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

Here is the code I am using:



SEIR <- function(t, 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


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,

More information about the R-help mailing list