[R] coupled ODE population model

[spencerg]

Hi, Justin: 


      In addition to the example below, Thomas Petzoldt also has two 
vignettes included with his "simecol" package.  These might help you 
more easily use his package. 


      Have you looked for what you want using the "RSiteSearch.function" 
in the "RSiteSearch" package?  Consider, for example, the following: 


library(RSiteSearch)
popec <- RSiteSearch.function('population ecology')
hits(popec)  #  100 help pages
summary(popec) #   29 packages
HTML(popec) # displays a table in a web browser


      This finds 100 different help pages in 29 different packages 
including this search term.  The last line opens a web browser with a 
table of the results sorted to put the package with the most hits first 
and providing links to HTML versions of those 100 help pages.  The 
RSiteSearch package also includes the capability to combine different 
searches using union and intersection. 


      Hope this helps. 
      Spencer

Thomas Petzoldt wrote:
> Justin Frank wrote:
>>
>> I'm fairly new to R, and I'm trying to write out a population model that
>> satisfies the following;
>>
>> the system consists of s species, i= 1, 2,...,s
>> network of interactions between species is specified by a (s x s) 
>> real matrix,
>> C[i,j]
>>
>> x[i] being the relative population of the "ith" species (0 =< x[i] =< 1,
>> sum(x[i]=1)
>>
>> the evolution rule being considered is as follows;
>>
>> xprime[i] = f[i] if x[i] > 0 or f[i] >= 0
>>
>> xprime[i] = 0 if x[i] = 0 and f[i] < 0
>>
>> where f[i] = sum(C[i,j]*x[j]) - x[i]*sum(C[k,j]*x[j])
>>
>>
>> I have a bit of attempted code written out, but are there any tricks 
>> or tips
>> that would condense or make this mess look nicer?
>>
>> -Justin
>
> Hi Justin,
>
> you may consider using package simecol (http://www.simecol.de). You 
> find several examples on the help pages and also some PDF files (so 
> called vignettes). If you have many species, you should consider using 
> matrix formulations. The example below demonstrates a multi-species 
> Lotka-Volterra model with an interaction matrix A.
>
> While the example is not yet your model it may serve as a starting 
> point. If you have any further questions, please let me know.
>
>
> Thomas Petzoldt
>
>
>
>
>
> library(simecol)
> ##########################################################
> ## Basic Multi Species Predator Prey Model
> ##  1) parameters for pairwise single species interaction
> ##########################################################
>
> LVPP <- new("odeModel",
>   main = function(t, n, parms) {
>     with(parms, {
>       dn <- r * n  + n * (A %*% n)
>       return(list(c(dn)))
>     })
>   },
>   parms = list(
>     r = c(r1 = 0.1, r2 = 0.1, r3 = -0.1, r4 = -0.1),
>     ## only pairwise interactions:
>     A = matrix(c(0.0, 0.0, -0.2, 0.0,      # prey 1
>                  0.0, 0.0, 0.0, -0.1,      # prey 2
>                  0.2, 0.0, 0.0, 0.0,       # predator 1; eats prey 1
>                  0.0, 0.1, 0.0, 0.0),      # predator 2; eats prey 2
>                  nrow = 4, ncol = 4, byrow=TRUE)
>   ),
>   times = seq(from=0, to=500, by = 0.1),
>   init  = c(prey1=1, prey2=1, pred1=2, pred2=2),
>   solver = "lsoda"
> )
>
> plot(sim(LVPP))
>
> ##########################################################
> ## 2) multi species interactions
> ##########################################################
>
> ## make two clones of LVPP
> LVPPweak <- LVPPstrong <- LVPP
>
> ## a helper function
>
> ## this copies the negative of the upper triangular to the lower
> makeSym <- function(A, f = -1) {
>   ind <- lower.tri(A)
>   A[ind] <- f * t(A)[ind]
>   A
> }
>
> ## weak coupling
> A <- matrix(c(0.0, 0.0, -0.1, -0.001,       # prey 1
>               NA,  0.0, -0.002, -0.2,       # prey 2
>               NA,  NA,   0.0,   0.0,        # predator 1
>               NA,  NA,   NA,    0.0),       # predator 2
>               nrow = 4, ncol = 4, byrow=TRUE)
>
> parms(LVPPweak)$A <- makeSym(A)
>
> ## stronger coupling
> A <- matrix(c(0.0, 0.0, -0.1, -0.05,        # prey 1
>               NA,  0.0, -0.02, -0.2,        # prey 2
>               NA,  NA,   0.0,   0.0,        # predator 1
>               NA,  NA,   NA,    0.0),       # predator 2
>               nrow = 4, ncol = 4, byrow=TRUE)
>
> parms(LVPPstrong)$A <- makeSym(A)
>
>
> LVPPweak <- sim(LVPPweak)
> LVPPstrong <- sim(LVPPstrong)
>
>
> plot(LVPPweak)
> plot(LVPPstrong)
>
>
> o <- out(LVPPweak)
> par(mfrow=c(2,2))
> plot(o$prey1, o$pred1)
> plot(o$prey2, o$pred2)
>
> plot(o$prey1, o$pred2)
> plot(o$prey2, o$pred1)
>
>
> o <- out(LVPPstrong)
> par(mfrow=c(2,2))
> plot(o$prey1, o$pred1)
> plot(o$prey2, o$pred2)
>
> plot(o$prey1, o$pred2)
> plot(o$prey2, o$pred1)
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>