[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