cOde (version 1.1.1)

odeC: Interface to ode()

Description

Interface to ode()

Usage

odeC(y, times, func, parms, ...)

Arguments

y

named vector of type numeric. Initial values for the integration

times

vector of type numeric. Integration times. If includeTimeZero is TRUE (see funC), the times vector is augmented by t = 0. If nGridpoints (see funC) was set >= 0, uniformly distributed time points between the first and last time point are introduced and the solution is returned for these time points, too. Any additional time points that are introduced during integration (e.g. event time points) are returned unless nGridpoints = -1 (the default).

func

return value from funC()

parms

named vector of type numeric.

...

further arguments going to ode()

Value

matrix with times and states

Details

See deSolve-package for a full description of possible arguments

Examples

# NOT RUN {
######################################################################
## Ozone formation and decay, modified by external forcings
######################################################################

library(deSolve)
data(forcData)
forcData$value <- forcData$value + 1

# O2 + O <-> O3
f <- c(
  O3 = " (build_O3 + u_build) * O2 * O - (decay_O3 + u_degrade) * O3",
  O2 = "-(build_O3 + u_build) * O2 * O + (decay_O3 + u_degrade) * O3",
  O  = "-(build_O3 + u_build) * O2 * O + (decay_O3 + u_degrade) * O3"
)

# Generate ODE function
forcings <- c("u_build", "u_degrade")
func <- funC(f, forcings = forcings, modelname = "test", 
             fcontrol = "nospline", nGridpoints = 10)

# Initialize times, states, parameters and forcings
times <- seq(0, 8, by = .1)
yini <- c(O3 = 0, O2 = 3, O = 2)
pars <- c(build_O3 = 1/6, decay_O3 = 1)

forc <- setForcings(func, forcData)

# Solve ODE
out <- odeC(y = yini, times = times, func = func, parms = pars, 
            forcings = forc)

# Plot solution

par(mfcol=c(1,2))
t1 <- unique(forcData[,2])
M1 <- matrix(forcData[,3], ncol=2)
t2 <- out[,1]
M2 <- out[,2:4]
M3 <- out[,5:6]

matplot(t1, M1, type="l", lty=1, col=1:2, xlab="time", ylab="value", 
	main="forcings", ylim=c(0, 4))
matplot(t2, M3, type="l", lty=2, col=1:2, xlab="time", ylab="value", 
	main="forcings", add=TRUE)

legend("topleft", legend = c("u_build", "u_degrade"), lty=1, col=1:2)
matplot(t2, M2, type="l", lty=1, col=1:3, xlab="time", ylab="value", 
	main="response")
legend("topright", legend = c("O3", "O2", "O"), lty=1, col=1:3)



######################################################################
## Ozone formation and decay, modified by events
######################################################################


f <- c(
  O3 = " (build_O3 + u_build) * O2 * O - 
         (decay_O3 + u_degrade) * O3",
  O2 = "-(build_O3 + u_build) * O2 * O + 
         (decay_O3 + u_degrade) * O3",
  O  = "-(build_O3 + u_build) * O2 * O + 
         (decay_O3 + u_degrade) * O3",
  u_build = "0",    # piecewise constant
  u_degrade = "0"   # piecewise constant
)

# Define parametric events
events.pars <- data.frame(
  var = c("u_degrade", "u_degrade", "u_build"),
  time = c("t_on", "t_off", "2"),
  value = c("plus", "minus", "2"),
  method = "replace"
)

# Declar parameteric events when generating funC object
func <- funC(f, forcings = NULL, events = events.pars, modelname = "test", 
             fcontrol = "nospline", nGridpoints = -1)

# Set Parameters
yini <- c(O3 = 0, O2 = 3, O = 2, u_build = 1, u_degrade = 1)
times <- seq(0, 8, by = .1)
pars <- c(build_O3 = 1/6, decay_O3 = 1, t_on = exp(rnorm(1, 0)), t_off = 6, plus = 3, minus = 1)

# Solve ODE with additional fixed-value events
out <- odeC(y = yini, times = times, func = func, parms = pars)


# Plot solution

par(mfcol=c(1,2))
t2 <- out[,1]
M2 <- out[,2:4]
M3 <- out[,5:6]


matplot(t2, M3, type="l", lty=2, col=1:2, xlab="time", ylab="value", 
        main="events")
legend("topleft", legend = c("u_build", "u_degrade"), lty=1, col=1:2)
matplot(t2, M2, type="l", lty=1, col=1:3, xlab="time", ylab="value", 
        main="response")
legend("topright", legend = c("O3", "O2", "O"), lty=1, col=1:3)


######################################################################
## Ozone formation and decay, modified by events triggered by root
######################################################################


f <- c(
  O3 = " (build_O3 + u_build) * O2 * O - 
         (decay_O3 + u_degrade) * O3",
  O2 = "-(build_O3 + u_build) * O2 * O + 
         (decay_O3 + u_degrade) * O3",
  O  = "-(build_O3 + u_build) * O2 * O + 
         (decay_O3 + u_degrade) * O3",
  u_build = "0",    # piecewise constant
  u_degrade = "0"   # piecewise constant
)

# Define parametric events
events.pars <- data.frame(
  var = c("u_degrade", "u_degrade", "u_build", "O3"),
  time = c("t_on", "t_off", "2", "t_thres_O3"),
  value = c("plus", "minus", "2", "0"),
  root = c(NA, NA, NA, "O3 - thres_O3"),
  method = "replace"
)

# Declar parameteric events when generating funC object
func <- funC(f, forcings = NULL, events = events.pars, modelname = "test", 
             fcontrol = "nospline", nGridpoints = -1)

# Set Parameters
yini <- c(O3 = 0, O2 = 3, O = 2, u_build = 1, u_degrade = 1)
times <- seq(0, 8, by = .01)
pars <- c(build_O3 = 1/6, decay_O3 = 1, 
          t_on = exp(rnorm(1, 0)), t_off = 6, plus = 3, minus = 1, 
          thres_O3 = 0.5, t_thres_O3 = 1)

# Solve ODE with additional fixed-value events
out <- odeC(y = yini, times = times, func = func, parms = pars, method = "lsode")

class(out) <- c("deSolve")
plot(out)
# Plot solution

par(mfcol=c(1,2))
t2 <- out[,1]
M2 <- out[,2:4]
M3 <- out[,5:6]


matplot(t2, M3, type="l", lty=2, col=1:2, xlab="time", ylab="value", 
        main="events")
legend("topleft", legend = c("u_build", "u_degrade"), lty=1, col=1:2)
matplot(t2, M2, type="l", lty=1, col=1:3, xlab="time", ylab="value", 
        main="response")
legend("topright", legend = c("O3", "O2", "O"), lty=1, col=1:3)




# }