[R] deSolve question
Chris Campbell
ccampbell at mango-solutions.com
Thu Jun 20 11:47:38 CEST 2013
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}}
More information about the R-help
mailing list