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

Daniel Reed reeddc at umich.edu
Fri Sep 13 21:57:11 CEST 2013


Hi Andras:

I default to using vode because I've found that it's usually faster than other numerical integrators. This may well vary between models or applications, but I think it's a good rule of thumb. That said, in this case the reduction in runtime is probably negligible since it's such a simple model. Not a bad habit to get into, though.

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 3:41 PM, Andras Farkas <motyocska at yahoo.com> wrote:

> 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
> 
> 
> _______________________________________________
> 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