CODEWITHSUNDEEP
rode
EBP
EBP

Reputation: 23

deSolve in R gives incorrect intermediate values

I am having problems with the ode solver in R from deSolve.

I have found that the solver gives non-zero intermediate values for variables that should be always equal to zero. It doesn't look like a problem until I start trying to debug with tests like if(R>0){browser()}, which gets triggered.

My code is below. Thanks in advance!

Ellen

library(deSolve)

simpleSIR <- function(t,states,par){
    with(as.list(c(states,par)),{

    S=states[1]
    I=states[2]
    R=states[3] 

    newinfections = beta*I*S

    dS <- -newinfections
    dI <- newinfections - gamma*I
    dR <- gamma*I

    if(R>0)
    {
        print(paste("why is this not zero?",R))
    }

    return(list(c(dS,dI,dR)))   
})}


par=list(beta=0.3,gamma=0.0)

init=c(0.99,0.01,0)

times <- seq(0,500,by = 1)


out <- as.data.frame(ode(y = init, times = times, func = simpleSIR, parms = par,maxsteps=2000))


plot(out[,2],type="l")
lines(out[,3],type="l",col="red")
lines(out[,4],type="l",col="blue")

Upvotes: 2

Views: 489

Answers (2)

David Heckmann
David Heckmann

Reputation: 2939

I think your problem lies in the fact that the default method for numerical integration in ode() is lsoda. This method can switch between solvers for stiff and non-stiff problems. In case you have switches to the stiff solvers, the jacobian is evaluted numerically, which might lead to the numerical errors you see. You can see that this might be the reason in the following code:

out <- deSolve::ode(y = init, times = times, func = simpleSIR, parms = par,maxsteps=2000)
deSolve::diagnostics.deSolve(out)

"[...]14 The number of Jacobian evaluations and LU decompositions so far: 23 [...]"

which corresponds to the number of print messages (23) your original code produces. You will get rid of the problem by using a non-stiff solver like RK4:

out.rk4 <- deSolve::ode(y = init, times = times, func = simpleSIR,method = "rk4", parms = par,maxsteps=2000)

if you insist on using lsoda, you might want to try supplying lsoda with the jacobian you calculated analytically.

Upvotes: 1

J_F
J_F

Reputation: 10372

First of all, I don't know if you understand the differences between state-variables and parameters.

R (in your case) is a state-variable. That means, that at the begin of the simulation (at t = 0) R = 0. But with time development this can change. See your state-variable S!

Your parameters beta and gamma, they will not change during the calculations.

Then you can simplify your code. Use this lines instead of your definition of par and init:

par <- c(beta = 0.3, gamma = 0.0)

init <- c(S = 0.99, I = 0.01, R = 0)

And now you can remove the following lines, because with with(as.list(c(states, par)) the names of the parameters and state-parameters are available with their names S, I and R.

S = states[1]
I = states[2]
R = states[3]

And you can simplify your plot, too:

matplot(out[,1], out[,2:4], type = "l", col = c("black", "red", "blue"))

Regards,

J_F

Upvotes: 1

Related Questions