[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