Профильная смена "Экология"

# устанавливаю нужную библиотеку
install.packages("deSolve")

# включаю нужную библиотеку
library(deSolve)

МОДЕЛЬ "ХИЩНИК - ЖЕРТВА"

# создаю функцию, решающую систему уравнений для модели Вольтерры и Лотки

fun_vl <- function(t, y, parms){
  X <- y[1]
  Y <- y[2]
  with(as.list(p),{
    dX.dt <- a1 * X - b21 * X * Y - c1 * X^2
    dY.dt <- b12 * X * Y - d2 * Y - c2 * Y^2
    return(list(c(dX.dt, dY.dt)))
  })
}

# задаю параметры модели - МЕНЯТЬ ЗДЕСЬ!
a1 <- 0.5 # коэффициент естественного прироста жертвы
b21 <- 0.01 # коэффициент, характеризующий скорость поедания жертвы хищником
b12 <- 0.01 # коэффициент изменения плотности популяции хищника, обусловленный поеданием жертвы
d2 <- 0.2 # коэффициент смертности хищника
c1 <- 0 # коэффициент внутривидовой конкуренции жертвы
c2 <- 0 # коэффициент внутривидовой конкуренции хищника

# ввожу начальные параметры системы и собираю параметры в вектор
y0 <- c(X = 25, Y = 5) # задаю начальные плотности популяций жертвы (X) и хищника (Y)
p <- c(a1 = a1, b21 = b21, b12 = b12, d2 = d2, c1 = c1, c2 = c2) # собираю параметры
t0 <- seq(0, 2000, 0.1) # задаю время наблюдения - 2000 условных лет (диапазон от 0 до 200 с шагом 0,1)

# вызываю пользовательскую функцию для расчёта модели
out <- ode(y = y0, t = t0, fun_vl, parms = p)

# рисую графики распределения функции
plot(t0, out[, 2], type = "l", xlab = "t", ylab = "x, y", ylim = c(0, 300), col = "red", lwd = 4) # строю первую линию,
# указав, что она красного цвета, а по оси Y будет диапазон от 0 до 300
# для оси X данные берутся из переменной t0, для оси Y - из второго столбика таблицы out
# подпись для оси X - t, для оси Y - X, Y
lines(t0, out[, 3], col = "blue", lwd = 3) # добавляю вторую линию,
# указав, что она синего цвета, данные из переменной t0 для оси X и из третьего столбца таблицы out
legend("topright", legend = c("вид-жертва", "вид-хищник"), col = c("red", "blue"),
       lty = c(1, 1), lwd = 4, xjust = 1, yjust = 1) # добавляю легенду, указав её место, типы линий и подписи
title("Графическое выражение модели системы 'хищник - жертва'") # добавляю заголовок

ОБОБЩЁННАЯ МОДЕЛЬ ВЗАИМОДЕЙСТВИЙ ПОПУЛЯЦИЙ

# создаю функцию, решающую систему уравнений для модели Вольтерры и Лотки
fun_vl <- function(t, y, parms){
  x1 <- y[1]
  x2 <- y[2]
  with(as.list(p),{
    dX.dt <- a1 * x1 + b12 * x1 * x2 - c1 * x1^2
    dY.dt <- a2 * x2 + b21 * x1 * x2 - c2 * x2^2
    return(list(c(dX.dt, dY.dt)))
  })
}

# задаю параметры модели
a1 <- 0.5 # коэффициент размножения жертвы
a2 <- 0.2 # коэффициент размножения хищника
b12 <- -0.01 # коэффициент влияния хищника на жертву
b21 <- 0.01 # коэффициент влияния жертвы на хищника
c1 <- 0.005 # коэффициент внутривидовой конкуренции жертвы
c2 <- 0.001 # коэффициент внутривидовой конкуренции хищника

# ввожу начальные параметры системы и собираю параметры в вектор
y0 <- c(x1 = 25, x2 = 5)
p <- c(a1 = a1, a2 = a2, b12 = b12, b21 = b21, c1 = c1, c2 = c2)
t0 <- seq(0, 200, 0.1)

# вызываю пользовательскую функцию для расчёта модели, перевожу переменную out в формат дата фрейма и добавляю в неё столбик с нумерацией строк
out <- ode(y = y0, t = t0, fun_vl, parms = p)
out <- as.data.frame(out)
out$id <- 1:nrow(out)

# рисую графики распределения функции
plot(out$id, out$x1, type = "l", xlab = "steps", ylab = "x, y", ylim = c(0, 200), col = "red", lwd = 4)
lines(out$id, out$x2, col = "blue", lwd = 3)
legend("topright", legend = c("вид-жертва", "вид-хищник"), col = c("red", "blue"), lty = c(1, 1), lwd = 4, xjust = 1, yjust = 1)

title("Графическое выражение модели системы 'хищник - жертва'")