[R-sig-dyn-mod] forcing variables (dummy ?) problem

Fred Jean fjean at univ-brest.fr
Mon Oct 17 22:22:37 CEST 2011


Dear list
I'm testing the use of the simecol package to get involved in 
mechanistic modelling in biology with R. (I'm an octave user for 
differential equations and R user for statistics, but I'd like to shift 
to get everything in R)

Since one week (I'm a stubborn dummy) I'm trying to understand how 
things work with forcing (external) data, but I'm getting nut because I 
don't understand how to get it work when having forcing data /at each 
timestep/.
As a sort of minimal example, here is an attempt to get such forcing 
inspired from the lv3 model example :

##============================================
require(simecol)

lvfred <- odeModel(
   main = function(time, init, parms, inputs) {

     # calcul par interpolation de la valeur de s.in à chaque valeur de
     # time
     s.in <- approxTime1(inputs, time, rule = 2)["s.in"]
     #cat(length(s.in)) # = 1 à chaque time

     with(as.list(c(init, parms)),{
         ds <- s.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)
       })
   },

   parms = c(b = 0.1, c = 0.1, d = 0.1, e = 0.1, f = 0.1, g = 0.05),
   times  = seq(from = 0, to = 365, by = 1),
   #
   # On peut rajouter des sous equations au modèle dans un fichier avec
   # un nom du type "equations.file.R"
   #
   # 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 = times,
       s.in = c(rep(0.1,99), rep(0.2,101), rep(0.1, 166))
       )
     ),
   init = c(s = 1, p = 1, k = 1), # substrate, producer, consumer
   solver = "lsoda"
   )

this code works perfectly but I don't understand how and where in the 
code I should  include e.g. a forcing temperature function depending on 
the vector /times/ as in the short code her after :

#
Tm <- 12 	# Température moyenne
Ampl <- 6 	# Amplitude de variation
Per <- 365	# Période de la sinusoïde
Faz <- 365/12	# Phase de la sinusoïde

tempC <-  Tm - Ampl * sin((2*pi*times / Per) - Faz)
#

And once I get this /tempC/ signal, how to include it ? Should I use the 
/equations slot/ ?

Usin tempC, I tried to get a /t.in/ forcing variable build with 
approxTime1() as for /s.in/ but didn't succeed so far... I tried to find 
examples with such forcings in simecolModels package, but none of them 
seemed to have such tricks.

Could you please tell me where to find solutions to those problems ?

Best whishes

Fred J.


-- 
Dr Fred JEAN - UMR IRD CNRS 6539 / LEMAR / Univ. Brest
Pho:+33 (0)2 98 49 86 38 // Fax:+33 (0)2 98 49 86 45
http://www.univ-brest.fr/IUEM/UMR6539/
== C'est curieux chez les marins, ce besoin de faire des phrases (M.A.)



More information about the R-sig-dynamic-models mailing list