[R-sig-dyn-mod] metapopulation SIR model

[Thomas Petzoldt]

Dear Vincent,

R is a vectorized language. Loops are often not necessary, and named 
vectors work differently. Try the example below. In addition, it would 
be also possible to write the model completely in matrix formulation.

Thomas



n <- 3
parameters <- c(beta=1000,gamma=365/13)
S <- rep(0.99, times=n)
I <- rep(0.001, times=n)
R <- rep(0.009, times=n)
state <- c(S, I, R)

sirmodel <- function(t, state, parameters) {
   S <- state[1:n]
   I <- state[(n+1):(2*n)]
   R <- state[(2*n+1):(3*n)]

   with (as.list(parameters),{
     dS <- -beta*S*I
     dI <- beta*S*I-gamma*I
     dR <- gamma*I
   # return the rate of change
   list(c(dS, dI, dR))
})
}
times <-seq(from=0,to=60/365,by=1/365/4)
print(system.time(
  out <- ode(y = state, times = times, func = sirmodel, parms = parameters)
  ))
head(out)
par(oma = c(0, 0, 3, 0))
plot(out, xlab = "time", ylab = "-", mfrow=c(3,3))