Hello Abdel,
I'm trying to model the spread of two viruses between different states, which
are i1 and i2, and i12 if you got both viruses, but you can go back to the
previous state with given probabilities (alpha, beta). The gamma probabilities
are just additional infection (without contact) and delta an interaction factor
between virus 1 and 2.
I tried lowering the time steps, and, as you said, i1 and i2 are going
negative, but it stops after a few steps. In effect, that's not what I'm
looking for. I just want to model the dynamic of this system, and do some king
of sensitivity analysis. I tried different set of parameters, but it still
gives me the same error message. Maybe this has to do with my equations ? But I
really doubt it.
Thank you for your time
Alex
From: abdel.hallo...@gmail.com
Date: Sat, 20 Feb 2016 11:01:00 -0600
Subject: Re: [R] Problems with the deSolve package
To: alexandresu...@hotmail.fr
CC: r-help@r-project.org
I think your parameters are off. If you look at the simul data frame, it gives
you a bunch of NaNs after the first initialization. If you put lower the
timesteps s.t.
>
times<-seq(0,200, by=0.01)
it begins to run but soon your values diverge, i1 & i2 going negative while i12
goes way high. Not sure what you are modeling but I assume those values aren't
to be like that. Try again with different parameters and see.
On Fri, Feb 19, 2016 at 2:42 AM, Alexandre Suire <alexandresu...@hotmail.fr>
wrote:
Hello R-users,
I'm trying to build a SIR-like model under R,
using the "deSolve" package. I'm trying to do simulations of its dynamic
over time, with three differential equations. I'm also looking to
calculate the equilibrium state.
So far, my code looks like this
library(deSolve)
#This is my system, with three differential equations
system<-function(times, init, parameters){
with(as.list(c(init, parameters)),{
di1_dt=(alpha1*(N-i1-i2-i12)*(i1+i12))+(beta2*i12+gamma1*(N-i1-i2-i12))-(beta1*i1)-(delta*alpha2*i1*(i2+i12))
di2_dt=(alpha2*(N-i1-i2-i12)*(i2+i12))+(beta1*i12+gamma2*(N-i1-i2-i12))-(beta2*i2)-(delta*alpha1*i2*(i1+i12))
di12_dt=(delta*alpha2*i1*(i12+i2))+(delta*alpha1*i2*(i12*i1))+(delta*gamma1*i1)+(delta*gamma2*i2)-((beta1+beta2)*i12)
return(list(c(di1_dt,di2_dt,di12_dt)))
})
}
# Initials values and parameters
init<-c(i1=10, i2=10, i12=0)
parameters<-c(alpha1=0.7, alpha2=0.5, beta1=0.5, beta2=0.3, gamma1=0.5,
gamma2=0.5, delta=0.5, N=100)
times<-seq(0,200, by=1)
simul <- as.data.frame(ode(y = init, times = times, func = system, parms =
parameters, method="ode45"))
simul$time <- NULL
head(simul,200)
#Plotting the results
matplot(times,
simul, type = "l", xlab = "Time", ylab = "i1 i2 i12", main =
"Simulation", lwd = 2, lty = 2, bty = "l", col=c("darkblue",
"darkred","mediumseagreen"))
legend("bottomright", c("i1", "i2","i12"), lty=2,lwd=2, col = c("darkblue",
"darkred", "mediumseagreen"))
At
first, I just tried studying with only the first two equations, and it
seems to work completely fine, but when I wanted to add the 3rd
equation, I sometimes get this message, even when I juggle the
parameters, when i launch the line:
#simul <- as.data.frame(ode(y = init, times = times, func = system, parms =
parameters))
Warning messages:
1: In lsode(y, times, func, parms, mf = 10, ...) :
an excessive amount of work (> maxsteps ) was done, but integration was not
successful - increase maxsteps
2: In lsode(y, times, func, parms, mf = 10, ...) :
Returning early. Results are accurate, as far as they go
Have I overlooked something ? I tried to use methods="ode45" and
methods="adams", without any sucess.
Thank you for your time
Alex
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.
