[R] Problems with the deSolve package
Abdel Halloway
abdel.halloway at gmail.com
Wed Feb 24 16:21:55 CET 2016
I think your delta is too high. If you reduce your delta (0.1 for example),
you should be able to get good results.
On Mon, Feb 22, 2016 at 3:37 AM, Alexandre Suire <alexandresuire at hotmail.fr>
wrote:
> 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.halloway at gmail.com
> Date: Sat, 20 Feb 2016 11:01:00 -0600
> Subject: Re: [R] Problems with the deSolve package
> To: alexandresuire at hotmail.fr
> CC: r-help at 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 <
> alexandresuire at 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 at 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]]
More information about the R-help
mailing list