Trouble interpretting R code for Lotka-Volterra

user_17576976 0
Viewed 65

I was given this code to create a Lotka Volterra predator-prey model. This code runs fine and produces the type of graph I am hoping to get, however, I don't understand what some of the code is doing and therefore don't know what I need to modify when using it on other data at other scales.

The two portions I don't understand are:

    for (n1 in seq((0.1*k[1]),(1.2*k[1]),100)){
  for (n2 in seq(.1*k[2],1.2*k[2],.1*k[2])){

and the if statement in this line:

    (n1==100&n2==100){plot(y[,2],y[,3],type="l",xlim=c(0,1000)

Any help interpreting this so I have a better understanding of how the code works would be greatly appreciated!

    library(deSolve)
    r<<-matrix(c(1,1),ncol=1)
    k<<-matrix(c(1000,1000),ncol=1)
    alpha<<-matrix(c(1,1.5,.5,1),ncol=2)
    t_stop<-100
    color="blue"
    for (n1 in seq((0.1*k[1]),(1.2*k[1]),100)){
      for (n2 in seq(.1*k[2],1.2*k[2],.1*k[2])){
    lvcomp2<-function(t,n,params){ w=n[1]; z=n[2]
    return(list(c(r[1]*n[1]*(k[1]-n[1]-alpha[1,2]*n[2])/k[1],r[2]*n[2]*(k[2]-n[2]- 
    alpha[2,1]*n[1])/k[2])))}
    y<-lsoda(c(n1,n2),c(seq(0,t_stop,1)),lvcomp2,NA)
    if (n1==100&n2==100){plot(y[,2],y[,3],type="l",xlim=c(0,1000), 
    ylim=c(0,1000),col=color,xlab="Abundance of species 1",ylab="Abundance of species 2")}
    if (n2==1200) {lines(y[,2],y[,3], col=color)}
    if (n1==1200){lines(y[,2],y[,3], col=color)}
    if (n2==100) {lines(y[,2],y[,3], col=color)}
    if (n1==100) {lines(y[,2],y[,3], col=color)}}}
modeling r
1 Answers

That code does a good job of obfuscating the process. Some questionable practices in it, such as superassignment <<-, unnecessary actions in loop, et cetera. Here is a simplified version of the code, which we can get by going step-by-step and taking things outside of the loop, filling in constant values where they don't change, and applying a style guide.

library(deSolve)

# in the sample code, both loops iterate over the same thing:
s <- seq(from = 100, to = 1200, by = 100)
lsoda_times <- seq(0, 100, 1)

# all combinations of s to solve an Ordinary Differential Equation for
grid <- expand.grid(s, s)
grid <- subset(grid, Var1 %in% c(1200, 100) | Var2 %in% c(1200, 100))

lvcomp2 <- function(t, n, parms, ...) { # use lsoda convention ..., see ?lsoda
  list(c( # list output to satisfy lsoda requirements
    n[1] * (1000-n[1] - 0.5*n[2]) / 1000, # derivatives
    n[2] * (1000-n[2] - 1.5*n[1]) / 1000
  ))
}

y <- lsoda(c(100, 100), seq(0, 100, 1), lvcomp2, NA)
# make base plot
plot(y[, 2], y[, 3], type = "l", xlim = c(0, 1000), ylim = c(0, 1000),
     col = "blue", xlab = "Abundance of species 1", ylab = "Abundance of species 2"
)

# add lines to the plot
for (i in 1:nrow(grid)) {
  y <- lsoda(c(grid[i, 1], grid[i, 2]), lsoda_times, lvcomp2, NA)
  lines(y[, 2], y[, 3], col = "blue")
}

per @tpetzoldt comment - updated to a more efficient approach first removing 100 unnecessary calls to lsoda

user_148539070
Related