demo/sirode.R

sir.ode <- function(times,init,parms){
  with(as.list(c(parms,init)), {
  # ODEs
  dS <- -beta*S*I
  dI <- beta*S*I-gamma*I
  dR <- gamma*I
  list(c(dS,dI,dR))
  })
}
parms <- c(beta=0.1,gamma=0.05)
init <- c(S=0.99,I=0.01,R=0) 
times <- seq(0,200,by=0.01)
sir.out <- lsoda(init,times,sirode,parms)
# Turn output into 'long' format
sir.out.long <- melt(as.data.frame(sir.out),"time")
# Now plot
ggplot(sir.out.long,aes(x=time,y=value,colour=variable,group=variable))+
  # Add line
  geom_line(lwd=2)+
  #Add labels
  xlab("Time")+ylab("Number")
reconhub/epicookbook documentation built on May 27, 2019, 4:02 a.m.