[R-sig-dyn-mod] NEGATIVE VALUES SIMULATED WITH DESOLVE
Andras Farkas
motyocska at yahoo.com
Sun Jun 30 22:33:25 CEST 2013
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
More information about the R-sig-dynamic-models
mailing list