# MODELO BÁSICO DE LOTKA VOLTERRA
j <- 500 # iteraciones
alfa <- 0.25
beta <- 0.01
gamma <- 1
sigma <- 0.01
h <- 0.25
H <- P <- seq(1, j, 1)
H[1] <- 80
P[1] <- 60
i <- 2
for (i in 2:j){
HM <- H[i-1]
PM <- P[i-1]
K1 <- h*(alfa*HM - beta*HM*PM)
L1 <- h*(sigma*PM*HM - gamma*PM)
K2 <- h*(alfa*(HM+(K1/2)) - beta*(HM+K1/2)*(PM+L1/2))
L2 <- h*(sigma*(PM+(L1/2))*(HM+K1/2) - gamma*(PM+(L1/2)))
K3 <- h*(alfa*(HM+(K2/2)) - beta*(HM+K2/2)*(PM+L2/2))
L3 <- h*(sigma*(PM+L2/2)*(HM+K2/2) - gamma*(PM+L2/2))
K4 <- h*(alfa*(HM+K3) - beta*(HM+K3)*(PM+L3))
L4 <- h*(( (sigma*(PM+L3)*(HM+K3)) - gamma *(PM+L3)))
H[i] <- HM +(K1+2*K2+2*K3+K4)/6
P[i] <- PM+(L1+2*L2+2*L3+L4)/6
}
plot(H, type= "l", ylim = c(0, max(H)), xla = "", ylab = "N de individuos", col = "blue", lwd = 2)
lines(P, type= "l", col = "red", lwd = 2)
X11()
plot (H,P, type= "l", col = "green", lwd = 2)
lines(c(0, max(H)), c((alfa/beta), (alfa/beta)), type = "l", lty=2, lwd=1)
lines(c((gamma/sigma), (gamma/sigma)), c(0, max(P)), type = "l", lty=2, lwd=1)
points((gamma/sigma), (alfa/beta), type = "p", col = "yellow", lwd = 3)
# MODELO DE PRESA DEPREDADOR CON CRECIMIENTO LOGÍSTICO PARA LA PRESA
j <- 500 # iteraciones
alfa <- 0.25
beta <- 0.01
gamma <- 0.5
sigma <- 0.01
h <- 0.25
Q <- 100
H <- P <- seq(1, j, 1)
H[1] <- 80
P[1] <- 30
i <- 2
for (i in 2:j){
HM <- H[i-1]
PM <- P[i-1]
K1 <- h*(alfa*HM*(1-HM/Q) - beta*HM*PM)
L1 <- h*(sigma*PM*HM - gamma*PM)
K2 <- h*(alfa*(HM+(K1/2)) *(1-(HM+(K1/2))/Q) - beta*(HM+K1/2)*(PM+L1/2))
L2 <- h*(sigma*(PM+(L1/2))*(HM+K1/2) - gamma*(PM+(L1/2)))
K3 <- h*(alfa*(HM+(K2/2)) *(1-(HM+(K2/2))/Q) - beta*(HM+K2/2)*(PM+L2/2))
L3 <- h*(sigma*(PM+L2/2)*(HM+K2/2) - gamma*(PM+L2/2))
K4 <- h*(alfa*(HM+K3) *(1-(HM+K3)/Q) - beta*(HM+K3)*(PM+L3))
L4 <- h*(( (sigma*(PM+L3)*(HM+K3)) - gamma *(PM+L3)))
H[i] <- HM +(K1+2*K2+2*K3+K4)/6
P[i] <- PM+(L1+2*L2+2*L3+L4)/6
}
plot(H, type= "l", ylim = c(0, max(H)), xla = "", ylab = "N de individuos", col = "blue", lwd = 2)
lines(P, type= "l", col = "red", lwd = 2)
X11()
plot (H,P, type= "l", col = "green", lwd = 2)
lines(c(0, Q), c((alfa/beta), 0), type = "l", lty=2, lwd=1)
lines(c((gamma/sigma), (gamma/sigma)), c(0, max(P)), type = "l", lty=2, lwd=1)
# MODELO DE PRESA DEPREDADOR CON CRECIMIENTO LOGÍSTICO PARA AMBOS
j <- 1000 # iteraciones
alfa <- 0.8
alfa.dep <- alfa/4
beta <- 0.07
ce <- 0.7
h <- 0.25
Q <- 100
H <- P <- seq(1, j, 1)
H[1] <- 20
P[1] <- 10
i <- 2
for (i in 2:j){
HM <- H[i-1]
PM <- P[i-1]
K1 <- h*(alfa*HM*(1-HM/Q) - beta*HM*PM)
L1 <- h*(alfa.dep*PM*(1-PM/(ce*HM)))
K2 <- h*(alfa*(HM+(K1/2)) *(1-(HM+(K1/2))/Q) - beta*(HM+K1/2)*(PM+L1/2))
L2 <- h*(alfa.dep*(PM+L1/2)*(1- (PM+L1/2)/(ce*(HM+K1/2))))
K3 <- h*(alfa*(HM+(K2/2)) *(1-(HM+(K2/2))/Q) - beta*(HM+K2/2)*(PM+L2/2))
L3 <- h*(alfa.dep*(PM+L2/2)*(1- (PM+L2/2)/(ce*(HM+K2/2))))
K4 <- h*(alfa*(HM+K3) *(1-(HM+K3)/Q) - beta*(HM+K3)*(PM+L3))
L4 <- h*(alfa.dep*(PM+L3)*(1- (PM+L3)/(ce*(HM+K3))))
H[i] <- HM +(K1+2*K2+2*K3+K4)/6
P[i] <- PM+(L1+2*L2+2*L3+L4)/6
}
plot(H, type= "l", ylim = c(min(P), max(H)), xla = "", ylab = "N de individuos", col = "blue", lwd = 2)
lines(P, type= "l", col = "red", lwd = 2)
X11()
plot (H,P, type= "l", col = "green", lwd = 2)
lines(c(0, Q), c((alfa/beta), 0), type = "l", lty=2, lwd=1)
lines(c(0, H[j]), c(0, P[j]), type = "l", lty=2, lwd=1)
#####################################################################
# PATRONES DE TURING
library(sp)
windows()
# Definiendo la cuadrícula
N <- 256
h <- 0.2 # Tamaño del paso en "x" y "y"
x <- h*(0:N) # coordenadas "x" de la gradilla
y <- h*(0:N) # coordenadas "y"
# Gradilla 2D con coordenadas "x" e "y"
X <- matrix(rep(x,each=length(y)),nrow=length(y));
Y <- matrix(rep(y,length(x)),nrow=length(y))
# Paso del tiempo (muy corto)
dt <- .01*h^2
# Valor de los parámetros según el paper
a <- 22
b <- 84 # Capacidad de carga
c <- 113.33
d <- 27.2 # Relación entre las difusiones
e <- 18
# Datos en t=0
us <- (a*b)/(c*e) # Estado estacionario de u
vs <- (c*e^2)/(b*a) # Estado estacionario de v
u <- matrix(rep(us,length(X)), nrow = length(y), ncol = length(x), byrow = FALSE) # Estado estacionario de u
v <- matrix(rep(vs,length(X)), nrow = length(y), ncol = length(x), byrow = FALSE) # Estado estacionario de v
u <- u + rnorm(length(X), 0.01) # agregamos pequeñas perturbaciones alrededor del estado estacionario
v <- v + rnorm(length(X), 0.01) # agregamos pequeñas perturbaciones alrededor del estado estacionario
# Pasos temporales
t <- 0
tmax <- 100
pasos <- round(tmax/dt) # cantidad de pasos temporales
# Desarrollo de la historia
for(i in 1:pasos){
t <- t+dt
# Difusión hacia los puntos cardinales
u.E <- u[,c((N+1),1:N)]
u.W <- u[, c(2:(N+1),1)]
u.N <- u[c(2:(N+1), 1),]
u.S <- u[c((N+1),1:N),]
v.E <- v[,c((N+1),1:N)]
v.W <- v[, c(2:(N+1),1)]
v.N <- v[c(2:(N+1), 1),]
v.S <- v[c((N+1),1:N),]
# Desarrollo de la cinética no lineal + difusión
u2 <- u^2
v2 <- v^2
u <- u + dt*((a*v*u2)-(e*u2)) + dt*(u.E+u.W+u.N+u.S-(4*u))/h^2
v <- v + dt*((v*b) - (c*u2*v2)) + d*dt*(v.E+v.W+v.N+v.S-(4*v))/h^2
# Mapeamos
vec.u <- as.vector(u)
vec.X <- as.vector(X)
vec.Y <- as.vector(Y)
datos <- cbind(vec.X, vec.Y, vec.u)
colnames(datos) <- c("X", "Y", "u")
datos <- data.frame(datos)
coord<- cbind(datos$X, datos$Y)
datos.sp <-SpatialPoints(coord)
datos.df <- SpatialPointsDataFrame(datos.sp, datos)
datos.grid <- as(datos.df, "SpatialPixelsDataFrame")
if(i%%50 == 0){
graphics.off()
}
print(spplot(datos.grid, "u", xlab = NULL, ylab = NULL, main = paste("Paso =", i), scales = list(draw = TRUE), cuts = 2, at = c(0, 1, 2), col.regions = gray.colors(2, start = 0, end = 1), colorkey= FALSE))
}