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