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

[Soetaert, Karline]

Fred,

I you type in R:
library(deSolve)
?forcings

A help page will be opened showing how to use forcing functions in the deSolve package. deSolve is the differential equation solver package used in simecol.


Karline

-----Original Message-----
From: r-sig-dynamic-models-bounces at r-project.org [mailto:r-sig-dynamic-models-bounces at r-project.org] On Behalf Of Fred Jean
Sent: 17 October 2011 22:23
To: R-sig-dynamic-models at r-project.org
Subject: [R-sig-dyn-mod] forcing variables (dummy ?) problem

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.)

_______________________________________________
R-sig-dynamic-models mailing list
R-sig-dynamic-models at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models