## ODE forcings

``````## =============================================================================
## R-Code to solve example 3.3.2 from the book
## K. Soetaert, J. Cash and F. Mazzia, 2012.
## Solving differential equations in R. UseR, Springer, 248 pp.
## http://www.springer.com/statistics/computational+statistics/book/978-3-642-28069-6.
##
## An ordinary differential equation model with forcing functions
## The ozone model  - needs file "light.rda"
## implemented by Karline Soetaert
## =============================================================================

library (deSolve)

load(file = "light.rda")   # contains data.frame Light
head(Light, n = 4)

irradiance(seq(from = 0, to = 1, by = 0.25))
k3 <- 1e-11; k2 <- 1e10; k1a <- 1e-30
k1b <- 1; sigma <- 1e11
yini <- c(O = 0, NO = 1.3e8, NO2 = 5e11, O3 = 8e11)

chemistry <- function(t, y, parms) {
with(as.list(y), {

k1  <- k1a + k1b*radiation

dO   <-  k1*NO2 - k2*O
dNO  <-  k1*NO2        - k3*NO*O3 + sigma
dNO2 <- -k1*NO2        + k3*NO*O3
dO3  <-           k2*O - k3*NO*O3
list(c(dO, dNO, dNO2, dO3), radiation = radiation)
})
}

times <- seq(from = 8/24, to = 5, by = 0.01)
out <- ode(func = chemistry, parms = NULL, y = yini,
times = times, method = "bdf")

plot(out, type = "l", lwd = 2, las = 1)``````

Go back