[R-sig-dyn-mod] changing state variable in simulations

[Andras Farkas]

Thanks a lot Daniel for the detailed break down of the code and modifications.... Just would like to ask one question: I see you have changed the solver to vode, may I ask why would that be your choice in this case?

thanks again

Andras

 
----- Original Message -----
From: Daniel Reed <reeddc at umich.edu>
To: Special Interest Group for Dynamic Simulation Models in R <r-sig-dynamic-models at r-project.org>
Cc: 
Sent: Friday, September 13, 2013 3:28 PM
Subject: Re: [R-sig-dyn-mod] changing state variable in simulations

Hi Andras:

How about this?

#Load library
require(deSolve)

#Define    parameters and initial conditions
plist<-matrix(c(0.08,0.07,0.09,15,20,25,3,5,6),ncol=3 )
colnames(plist)=c("K","V","state")

#Define    "forc" function prior to "derivs", so it isn't redefined every time "derivs" is    called
intimes<-c(0,0.5,12)
input<-c(1200,0,0)
forc<-approxfun(intimes, input, method="constant",rule=2)

#Derivatives to be integrated
derivs<-function(t,state,pars){
  with(as.list(pars),{
    inp <- forc(t)
    dy1 <- inp/V - K * state[1]
    return(list(c(dy1)))
  })
}

#Define    output times
times=seq(0,13,by=0.05)

#Model wrapper
model<-function(pars) {
  state<-pars[3]
  return(ode(y=state,times=times,func=derivs,parms=pars[1:2],method="vode",rtol=1e-8,atol=1e-8,hmax=0.05))
}

#Empty matrix
modelout<-matrix(ncol=nrow(plist),nrow=length(times))

#Loop through parameters & initial conditions
for(i in 1:nrow(plist))
  modelout[,i]<-model(plist[i,])[,2]

#Plot output
plot(times,modelout[,1],type="l",lty="solid")
lines(times,modelout[,2],lty="dashed")
lines(times,modelout[,3],lty="dotted")

Cheers,
Daniel


_______________________________
Daniel C. Reed, PhD.
Postdoctoral Fellow,
Dept. of Earth & Environmental Sciences,
University of Michigan,
Ann Arbor, MI, USA.
email: reeddc at umich.edu
web: www.danielreed.org



On Sep 13, 2013, at 1:34 PM, Andras Farkas <motyocska at yahoo.com> wrote:

>  Dear All,
>  
> have this simple model:
>  
> require(deSolve)
> plist <-matrix(c(0.08,0.07,0.09,15,20,25),ncol=2 )
> colnames(plist)=c("K","V")
> derivs <- function(t, state, pars) {
>  with(as.list(pars), {
>    intimes <- c(0,0.5,12)
>    input  <- c(1200,0,0)
>    forc <- approxfun(intimes, input, method="constant",rule=2)
>    inp <- forc(t)
>    dy1 <- inp/V - K * state[1] 
>    return(list(c(dy1)))
>  })
> 
> }
> model <- function(pars, times=seq(0, 13, by = 0.05)) {
>  state <- c(a = 5 )
>  return(ode(y = state, times = times, func = derivs, parms = plist[i,],method="lsoda",rtol = 1e-8,atol =1e-8,hmax=0.05))
> }
> modelout <- list()
> for (i in 1:nrow(plist))
>  modelout[[i]] <- model(plist)
> 
> out <-data.frame(modelout)
> time <-out$time
> out <-out[, grepl("^a", colnames(out))]
> plot(time,out[,3],type="l")
> 
>  
> Please provide some insights for the following: I would like to change the initial state variable during simulations for each time we move down one row in plist, ie: just like the parameters K and V change from row to row, I would also like the initial state variable to change by "putting" the initial state variables into plist like this:
>  
>  plist <-matrix(c(0.08,0.07,0.09,15,20,25,3,5,6),ncol=3 )
> colnames(plist)=c("K","V","state")
> 
> I would greatly appreciate your thoughts on how to incorporate this into the above model,
>  
> 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

_______________________________________________
R-sig-dynamic-models mailing list
R-sig-dynamic-models at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models