Hello, I am using odesolve to simulate a group of people moving through time and transmitting infections to one another.
In Matlab, there is a NonNegative option which tells the Matlab solver to keep the vector elements of the ODE solution non-negative at all times. What is the right way to do this in R? Thanks, Jeremy P.S., Below is a simplified version of the code I use to try to do this, but I am not sure that it is theoretically right dynmodel <- function(t,y,p) { ## Initialize parameter values birth <- p$mybirth(t) death <- p$mydeath(t) recover <- p$myrecover beta <- p$mybeta vaxeff <- p$myvaxeff vaccinated <- p$myvax(t) vax <- vaxeff*vaccinated/100 ## If the state currently has negative quantities (shouldn't have), then reset to reasonable values for computing meaningful derivatives for (i in 1:length(y)) { if (y[i]<0) { y[i] <- 0 } } S <- y[1] I <- y[2] R <- y[3] N <- y[4] shat <- (birth*(1-vax)) - (death*S) - (beta*S*I/N) ihat <- (beta*S*I/N) - (death*I) - (recover*I) rhat <- (birth*(vax)) + (recover*I) - (death*R) ## Do we overshoot into negative space, if so shrink derivative to bring state to 0 ## then rescale the components that take the derivative negative if (shat+S<0) { shat_old <- shat shat <- -1*S scaled_transmission <- (shat/shat_old)*(beta*S*I/N) ihat <- scaled_transmission - (death*I) - (recover*I) } if (ihat+I<0) { ihat_old <- ihat ihat <- -1*I scaled_recovery <- (ihat/ihat_old)*(recover*I) rhat <- scaled_recovery +(birth*(vax)) - (death*R) } if (rhat+R<0) { rhat <- -1*R } nhat <- shat + ihat + rhat if (nhat+N<0) { nhat <- -1*N } ## return derivatives list(c(shat,ihat,rhat,nhat),c(0)) } ______________________________________________ R-help@stat.math.ethz.ch 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.