[R-sig-dyn-mod] NEGATIVE VALUES SIMULATED WITH DESOLVE

[Karline Soetaert]

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