## general solver for simObj models and analytically solved odeModel's
##
## NOTE 1: iteration returns the new state
## (as opposed to ODE solvers which return the first derivatives)
## NOTE 2: the special parameter DELTAT is set here to ensure
## consistency with the times vector.
setGeneric("iteration",
function(y, times=FALSE, func=FALSE, parms=FALSE, animate=FALSE, ...)
standardGeneric("iteration")
)
setMethod("iteration", "simObj",
function(y, times=NULL, func=NULL, parms=NULL, animate=FALSE, ...) {
init <- y@init
times <- fromtoby(y@times)
func <- y@main
parms <- y@parms
inputs <- y@inputs
equations <- y@equations
environment(func) <- environment()
attach(equations)
on.exit(detach(equations))
parms$DELTAT <- 0
out <- list(func(times[1], init, parms))
for (i in 2:length(times)) {
time <- times[i]
parms$DELTAT <- times[i] - times[i-1]
init <- func(time, init, parms)
out <- c(out, list(init))
if (animate) {
y@out <- out
plot(y, index=i, ...)
}
}
out
}
)
setMethod("iteration", "odeModel",
function(y, times=NULL, func=NULL, parms=NULL, animate=FALSE, ...) {
init <- y@init
times <- fromtoby(y@times)
func <- y@main
parms <- y@parms
inputs <- y@inputs
equations <- y@equations
environment(func) <- environment()
attach(equations)
on.exit(detach(equations))
n <- length(init)
parms <- c(parms, DELTAT = 0)
nm <- c("time", if (!is.null(attr(init, "names")))
names(init) else as.character(1:n))
out <- unlist(func(times[1], init, parms))
for (i in 2:length(times)) {
time <- times[i]
parms["DELTAT"] <- times[i] - times[i - (i>1)] # is zero if i=1
init <- unlist(func(time, init, parms))
out <- rbind(out, init)
if (animate) {
y@out <- out
plot(y, index=i, ...)
}
}
row.names(out) <- NULL
out <- as.data.frame(cbind(times, out))
names(out) <- nm
out
}
)