[R-sig-dyn-mod] NEGATIVE VALUES SIMULATED WITH DESOLVE
Karline Soetaert
Karline.Soetaert at nioz.nl
Mon Jul 1 09:22:06 CEST 2013
Andras,
It is a flaw in your model equation for y which says:
dy/dt = -k*z + inp
The decline of y is unrelated to y, and can proceed even if y = 0 or negative.
Karline
-----Original Message-----
From: r-sig-dynamic-models-bounces at r-project.org [mailto:r-sig-dynamic-models-bounces at r-project.org] On Behalf Of Andras Farkas
Sent: zondag 30 juni 2013 22:33
To: r-sig-dynamic-models at r-project.org
Subject: [R-sig-dyn-mod] NEGATIVE VALUES SIMULATED WITH DESOLVE
Dear List,
coul you please provide some insights for the following:
code is:
pars <-list(k1=0.06787964, k2= 0.02672658,k3= 0.03684916,v1=0.3296415,v2= 0.2737872,v3= 0.1936819)
input <-c(180 ,0, 180, 0, 180, 0, 180,0, 180,0, 180,0, 180,0, 180,0, 180, 0,180 , 0 ,180, 0 ,180, 0, 180, 0, 180, 0, 180, 0, 180 , 0 , 0)
intimes <-c( 0.0 , 0.5, 72.0, 72.5, 144.0, 144.5, 216.0, 216.5, 288.0, 288.5,
360.0, 360.5, 432.0, 432.5, 504.0, 504.5, 576.0, 576.5, 648.0 , 648.5 ,720.0, 720.5, 792.0, 792.5, 864.0, 864.5, 936.0, 936.5, 1008.0, 1008.5 , 1080.0, 1080.5, 1152.0) forc <- approxfun(intimes, input, method="constant", rule=2)
kintimes <- c(118.6 ,162.1, 382.6)
model <- function(pars, times=seq(0, 527, by = 0.05)) {
state <- c(z = 0,y = 0)
return(ode(y = state, times = times, func = ab, parms = pars,method="lsoda",rtol = 1e-11,atol =1e-11)) }
ab <-function(t, state, pars) {
inp <- forc(t)
if (t<=kintimes[1]) k =pars$k1 else if (t<=kintimes[2]) k =pars$k2 else k =pars$k3
if (t<=kintimes[1]) v =pars$v1 else if (t<=kintimes[2]) v =pars$v2 else v =pars$v3
dy1 <- - k * state[1] + inp
dc <- dy1/(v*55.99)
return(list(c(dy1,dc)))
}
summary(model(pars))
as you can see, y will include negative values, which is not desirabel. I played around a bit with rtol and atol, but no consitent success. I also tried to set the ab function as a forcing function, but still gettin the negative values. Any thoughts on what i may be doing wrong?
thanks,
Andras
_______________________________________________
R-sig-dynamic-models mailing list
R-sig-dynamic-models at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models
More information about the R-sig-dynamic-models
mailing list