[R-sig-dyn-mod] forcing variables (dummy ?) problem
Fred Jean
fjean at univ-brest.fr
Thu Oct 20 23:24:03 CEST 2011
Dear Thomas and Karline
Many thanks for your avices and for giving attention to my newbie's
problems.
I didn't progress a lot those last days (had many courses).
At this point, before those last two emails, I followed the advice of
Karline and here is a piece of code. As you will see, I commented the
code from what I understood. I will now try to take into account advice
from Thomas.
Many thanks again to both of you
Fred
##============================================
require(simecol)
rm(list=ls())
# external parameters
# time for constructing functions and for length of simulations
t1 = seq(from = 0, to = 800, by = 1)
# temperature forcing function parameters
Tm <- 12 # mean temperature
TAmpl <- 6 # temperature variation amplitude
TPer <- 365 # temperature sine period
TFaz <- 365/6 # temperature sine phase
#substrate forcing
Sm <- 0.5 # mean substrate input
SAmpl <- 0.5 # substrate input amplitude
SPer <- 365/2 # substrate input sine period
SFaz <- 365/12 # substrate input sine phase
#### The model
lvfred <- odeModel(
main = function(time, init, parms, inputs) {
# interpolate forcings - output of those is one interpolated value
# at each integration step
s.in <- approxTime1(inputs, time, rule = 2)["s.in"]
t.in <- approxTime1(inputs, time, rule = 2)["t.in"]
# model equations
with(as.list(c(init, parms)),{
ds <- s.in - t.in*b*s*p + g*k
dp <- c*s*p - d*k*p
dk <- e*p*k - f*k
list(c(ds, dp, dk), s.in, t.in)
})
},
# values of parameters
parms = c(b = 0.1, c = 0.1, d = 0.1, e = 0.1, f = 0.1, g = 0.05),
# duration of simulations and integration timestep
times = c(from = 0, to = 800, by = 0.1),
# One may add sub-equations in a file
# "equations.file.R" -- see daphnia example in simecolModel
# equations = equations.file,
inputs = as.matrix(
data.frame(
#time = c(0, 99, 100, 150, 200, 365),
#s.in = c(0.1, 0.2, 0.2, 0.2, 0.2, 0.1)
time = t1,
s.in = Sm - SAmpl * sin((2*pi*t1 / SPer) - SFaz),
t.in = Tm - TAmpl * sin((2*pi*t1 / TPer) - TFaz)
)
),
init = c(s = 1, p = 1, k = 1), # substrate, producer, consumer
solver = "rk4"
)
##============================================
More information about the R-sig-dynamic-models
mailing list