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.

Reply via email to