Skip to content

[R-sig-dyn-mod] ode with variable dependent on results from previous time step

2 messages · Guadalupe Peralta, Thomas Petzoldt

Original
Guadalupe Peralta
# 
*Dear all,
*

*I'm using the ode function and I wonder how to incorporate a variable
into the model that results from a previous time step.*

*Lets say I have two groups of organisms (P and A) with multiple
species within each group (nP=3, nA=2). Species from group A interact
with species from group P, with an interaction coefficient of 'c'. My
function (funcPA) describing the changes in time is as follows:*


funcPA <- function (t, y, params, nsp) {
  with(as.list(c(y,params)), {
    nP <- nsp[1]
    nA <- nsp[2]
    ntot <- nP + nA
    P <- y[1:nP]
    A <- y[(nP + 1):ntot]
    dP <- array(0, nP)
    dA <- array(0, nA)
    c <- params$c
    for (i in 1:nP) {
      for (j in 1:nA){
        dP[i] <- P[i] + sum(c*A[j]*P[i])
        dP
        dA[j] <- A[j] + sum(c*P[i]*A[j])
        dA }}
    res <-list(c(dP, dA))
  })
}

nsp <- c(nP=3, nA=2)
c <- list()
c[[1]] <- matrix(c(0.1,0.11,0.12,0,0,0.1),nrow=3, ncol=2)
y <- c(P1=2, P2=1, P3=2, A1=1, A2=2)

ode(y=y, times=seq(0,1,by=0.1), func=funcPA, nsp=nsp, parms=c)


*Now, I want to include a variable Ak on the dP[i] equation of the
funcPA function so that it looks like this:*

dP[i] <- P[i] + sum((c*A[j]*P[i]) / (Ak[i]))
*where Ak[i] is the number of individuals from group A (considering
all A species) that interact with each particular P species (P[i]) on
each time step. For the first time step Ak can be calculated from the
initial abundances of species As (i.e. A1=1, A2=2) and a binary
version of the interaction coefficients 'c' (to determine which A
species interact with each P species). However, as the abundances of A
will change from one time step to the next, I wonder how should I
include it the equation? or should I include Ak somehow in the ode
function?*


*Thanks for your help.
Kind regards
Guadalupe
*
Thomas Petzoldt
# 
Hi,

I don't yet completely understand what you want, but maybe the following 
can help you a little bit.

1) your code is not completely correct, as parms$c does not contain the 
interaction matrix. Add a print to your function to see this:

     c <- params$c
     print(c)


2) using for-loops in a model function is quite unusual and, together 
with calculating sums a strong indication to me, that you should 
consider matrix multiplication, see page 124 vs. 123 in our tutorial:

http://desolve.r-forge.r-project.org/user2014/tutorial.pdf#124

That example shows a multi-species Lotka-Volterra model. Could it be 
possible that you want something like that?

3) I wonder also, why you use c*P[i] and not c[i,j] and why all 
interactions are positive.

4) Finally: in ODE models, "time steps" are infinitesimally small (i.e. 
close to zero). This means that interaction happens instantaneously or 
with other words: all variables are a result from previous time. Time is 
continuous, so time "step" and "step before" do not make much sense -- 
except the technical fact that the solver approximates continuous time 
with finite steps; or "time delay" for delay differential equations, but 
this would be another story.

Regards,

Thomas P.


Am 27.02.2017 um 05:02 schrieb Guadalupe Peralta: