[R] deSolve question

[Chris Campbell]

Dear Andras

In algebra (and computing) terms are operated upon in order of priority. Terms in brackets are evaluated first. In your original code, - pars$k * state[1] is not divided by pars$v since division happens before addition. When you use the brackets that expression becomes part of the term that is divided.

Therefore the outputs of these functions should not be the same.

Best wishes

Chris

Chris Campbell, PhD    
Tel. +44 (0) 1249 705 450 | Mobile. +44 (0) 7929 628 349    
mailto:ccampbell at mango-solutions.com | http://www.mango-solutions.com    
Mango Solutions, 2 Methuen Park, Chippenham, Wiltshire , SN14 OGB UK    

-----Original Message-----
From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On Behalf Of Andras Farkas
Sent: 19 June 2013 00:46
To: r-help at r-project.org; r-sig-dynamic-models at r-project.org
Subject: [R] deSolve question



Dear All

wonder if you could provide some insights on the following: currently I have this code which produces the expected results:

require(deSolve)
pars <- list(k = 0.08,v=15)
intimes <- c(0,0.5,12)
input   <- c(800,0,0)
forc <- approxfun(intimes, input, method="constant", rule=2)

derivs <- function(t, state, pars) {
  inp <- forc(t)
  dy1 <- - pars$k * state[1] + inp/pars$v
  return(list(c(dy1)))
}

model <- function(pars, times=seq(0, 13, by = 0.01)) {
  state <- c(y = 0)
  return(ode(y = state, times = times, func = derivs, parms = pars,
             method="lsoda"))
}

plot(model(pars))

intuitively it would make sense to me that if I change the derivs as:

derivs <- function(t, state, pars) {
  inp <- forc(t)
  #note the added ()
  dy1 <- (- pars$k * state[1] + inp)/pars$v
  return(list(c(dy1)))
}

than I would get the same results, but that is not the case. I need to "relocate" pars$v to another position within derivs so that I could implement it in a forcing function to be able to change its value during derivation. Appreciate your help,

Andras

______________________________________________
R-help at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

--

LEGAL NOTICE\ \ This message is intended for the use of ...{{dropped:18}}