## https://strimas.com/post/lotka-volterra/ library(tidyverse) library(deSolve) library(FME) params1 <- c(alpha=1.2,beta=0.2, delta =0.5, gamma=0.2) params2 <- c(alpha=1.0,beta=0.2, delta =0.5, gamma=0.2) params3 <- c(alpha=0.8,beta=0.2, delta =0.5, gamma=0.2) init <- c(x=1, y=1) times <- seq(0,100,by=1) deriv <- function(t, state, pars) { with(as.list(c(state, pars)), { d_x <- alpha * x - beta * x * y d_y <- delta * beta * x * y - gamma * y return(list(c(x = d_x, y = d_y))) }) } lv_results1 <- ode(init, times, deriv, params1) lv_results2 <- ode(init, times, deriv, params2) lv_results3 <- ode(init, times, deriv, params3) lv_model <- function(pars, times = seq(0, 50, by = 1)) { # initial state state <- c(x = 1, y = 2) # derivative deriv <- function(t, state, pars) { with(as.list(c(state, pars)), { d_x <- alpha * x - beta * x * y d_y <- delta * beta * x * y - gamma * y return(list(c(x = d_x, y = d_y))) }) } # solve ode(y = state, times = times, func = deriv, parms = pars) } lv_results1 <- lv_model(pars = params1, times = seq(0, 50, by = 0.25)) lv_results2 <- lv_model(pars = params2, times = seq(0, 50, by = 0.25)) lv_results3 <- lv_model(pars = params3, times = seq(0, 50, by = 0.25)) lv_results1 %>% data.frame() %>% gather(var, pop, -time) %>% mutate(var = if_else(var == "x", "Prey", "Predator")) %>% ggplot(aes(x = time, y = pop)) + geom_line(aes(color = var)) + scale_color_brewer(NULL, palette = "Set1") + labs(title = "Lotka-Volterra predator prey model", subtitle = paste(names(params1), params1, sep = " = ", collapse = "; "), x = "Time", y = "Population density") plot(lv_results3[,2],lv_results3[,3],xlab="Prey",ylab="Predator",ylim=c(0,20),xlim=c(0,20), type="l",col="blue",cex.axis=1.5,cex.lab=1.5) lines(lv_results2[,2],lv_results2[,3],col="red") lines(lv_results1[,2],lv_results1[,3],col="green")