Codigo_Parcial_III.R

# Punto 2c, usamos los valores de alpha y beta que calculamos y el resto de parámetros en Gause (1932)
# El tiempo siempre lo modificamos hasta observar bien los efectos de la competencia

library(deSolve)

## Warning: package 'deSolve' was built under R version 3.3.2

Lvcomp2 <- function (t, n, parms) {
  with(as.list(parms), {
    dn1dt = r1*n[1]*((k1-n[1]-alpha*n[2])/k1)
    dn2dt = r2*n[2]*((k2-n[2]-beta*n[1])/k2)
    list(c(dn1dt, dn2dt))
  })
}

parms <- c(k1=13, k2=5.8, r1 = 0.22, r2 = 0.06, alpha = 2.186, beta = 0.457)
initialN=c(1, 1)
out <- ode(func = Lvcomp2, y = initialN, parms = parms, times = 1:4000)
matplot(out[,1], out[,-1], type = "l")

# Punto 2d, alteramos el valor de alpha según el experimento de Gause (1932)
parms <- c(k1=13, k2=5.8, r1 = 0.22, r2 = 0.06, alpha = 3.15, beta = 0.457)
initialN=c(1, 1)
out <- ode(func = Lvcomp2, y = initialN, parms = parms, times = 1:1000)
matplot(out[,1], out[,-1], type = "l")

# Punto 2e, modificamos alpha y beta para que se cumpla la coexistencia
parms <- c(k1=13, k2=5.8, r1 = 0.22, r2 = 0.06, alpha = 2.186, beta = 0.439)
initialN=c(1, 1)
out <- ode(func = Lvcomp2, y = initialN, parms = parms, times = 1:2000)
matplot(out[,1], out[,-1], type = "l")

# Punto 3c, utilizamos valores de V y P que se desplacen de la intersección de las isoclinas (equilibrios)
LotVmod <- function (Time, State, Pars) {
  with(as.list(c(State, Pars)), {
    dV = (r*V)-(alpha*V*P)
    dP = (beta*V*P)-(q*P)
    return(list(c(dV, dP)))
  })
}

Pars <- c(r = 0.7, alpha = 1, q = 2.088, beta = 0.5)
State <- c(V = 10, P = 0.4)
Time <- seq(0, 100, by = 1)
out <- as.data.frame(ode(func = LotVmod, y = State, parms = Pars, times = Time))
matplot(out[,-1], type = "l", xlab = "time", ylab = "population")
legend("topright", c("conejos", "zorros"), lty = c(1,2), col = c(1,2), box.lwd = 0)

# Punto 4c, modificamos la función de metapoblaciones para incorporar un modelo con colonización interna y efecto rescate
metopob <- function (t, y, parms) {
  f = y[1]
  with(as.list(parms), {
    dp = pi*f*(1-f)-pe*f*(1-f)
    return(list(dp))
  })
}

prms <- c(pi=0.3, pe=0.5)
Initial.f=c(0.4)
out.M <- data.frame(ode(y = Initial.f, parms = prms, times = 1:100, func = metopob))
plot(out.M[,2]~ out.M[,1], type = "l", ylim=c(0,1), ylab = "f", xlab = "tiempo")