#########################################################
######## Modelos em tempo continuo nao-lineares ########
########################################################
rm(list=ls())
install.packages("deSolve")

library(deSolve)


#############################
######## Predacao LV ########
#############################

#### Definindo a funcao de predacao:
PredacaoLV <- function (Time, State, Pars) {
  with(as.list(c(State, Pars)), {
    dV = r*V - alpha*V*P
    dP = -q*P + beta*P*V
    return(list(c(dV, dP)))
  })
}


#Parametros do modelo
Pars <- c(r = 0.15, q=15, alpha = 0.01, beta = 1)
#Tamanho inicial das populacoes
State <- c(V = 10, P =10)
#Tempo a ser rodado a solucao numerica
Time <- seq(0, 100, by = 0.1)

out <- as.data.frame(ode(func = PredacaoLV, y = State, parms = Pars, times = Time))
colnames(out)=c("tempo", "V", "P")
out


par(mfrow=c(1,2)) #Plotando 2 resultados na mesma tela

# Plotando os tamanhos de cada especie pelo tempo
plot(out$V~out$tempo , type = "l", xlab = "tempo", ylab = "Numero de individuos", 
     col=1, lty=1, ylim=c(0,150))
lines(out$P~out$tempo , type = "l", col=2, lty=2)
legend("topright", c("V", "P"), lty = c(1,2), col = c(1,2), box.lwd = 0)

# Plotando a densidade de uma sp pela outra (escala de cinza de acordo com o tempo, 
# onde mais claro eh o inicio)
cinzas=  gray.colors(dim(out)[1], start = 1, end = 0, gamma = 0.7, alpha = NULL)
plot(1, type="n", xlab="V", ylab="P", xlim=c(0,30), ylim= c(0,170))
abline(h=Pars[1]/Pars[3],lty=1)
abline(v=Pars[2]/Pars[4],lty=2)
points(out$P~out$V, col=cinzas, ylab="P", xlab="V", pch=16)
#dev.off()

#### Acrescentaremos outro estado inicial no 2o grafico
State2 <- c(V = 13, P =13)

out2 <- as.data.frame(ode(func = PredacaoLV, y = State2, parms = Pars, times = Time))
colnames(out2)=c("tempo", "V", "P")
out2

par(mfrow=c(1,1)) #Plotando 2 resultados na mesma tela

cinzas=  gray.colors(dim(out)[1], start = 1, end = 0, gamma = 0.7, alpha = NULL)
plot(1, type="n", xlab="V", ylab="P", xlim=c(0,30), ylim= c(0,170))
abline(h=Pars[1]/Pars[3],lty=1)
abline(v=Pars[2]/Pars[4],lty=2)
points(out$P~out$V, col=cinzas, ylab="P", xlab="V", pch=16)
points(out2$P~out2$V, col=cinzas, ylab="P", xlab="V", pch=16)

plot(out2$V~out2$tempo , type = "l", xlab = "tempo", ylab = "Numero de individuos", 
     col=1, lty=1, ylim=c(0,150))
lines(out2$P~out$tempo , type = "l", col=2, lty=2)
legend("topright", c("V", "P"), lty = c(1,2), col = c(1,2), box.lwd = 0)


##########################
######## Predacao ########
##########################

#### Definindo a funcao de predacao com resposta funcional tipo I:
Predacao <- function (Time, State, Pars) {
  with(as.list(c(State, Pars)), {
    dV = r*V*((K1 - V)/K1) - alpha*V*P
    dP = -q*P + beta*P*V
    return(list(c(dV, dP)))
  })
}

####################################################################
###### Caso 1: so presas persistem. Predadores vao a extincao ######
####################################################################

#Parametros do modelo
Pars <- c(r = 1, q=150, alpha = 0.01, beta = 1, K1=100)
#Tamanho inicial das populacoes
State <- c(V = 10, P =10)
#Tempo a ser rodado a solucao numerica
Time <- seq(0, 100, by = 1)

out <- as.data.frame(ode(func = Predacao, y = State, parms = Pars, times = Time))
colnames(out)=c("tempo", "V", "P")
out


par(mfrow=c(1,2)) #Plotando 2 resultados na mesma tela

# Plotando os tamanhos de cada especie pelo tempo
plot(out$V~out$tempo , type = "l", xlab = "tempo", ylab = "Numero de individuos", 
     col=1, lty=1, ylim=c(0,110))
lines(out$P~out$tempo , type = "l", col=2, lty=2)
legend("topright", c("V", "P"), lty = c(1,2), col = c(1,2), box.lwd = 0)

# Plotando a densidade de uma sp pela outra (escala de cinza de acordo com o tempo, 
# onde mais claro eh o inicio)
cinzas=  gray.colors(dim(out)[1], start = 1, end = 0, gamma = 0.7, alpha = NULL)
plot(1, type="n", xlab="V", ylab="P", xlim=c(0,150), ylim= c(0,110))
abline(a=Pars[1]/Pars[3],b=-(Pars[1]/Pars[3])/Pars[5],lty=1)
abline(v=Pars[2]/Pars[4],lty=2)
points(out$P~out$V, col=cinzas, ylab="P", xlab="V", pch=16)
dev.off()

###################################################
###### Caso 2: Presas e predadores coexistem ######
###################################################

#Parametros do modelo
Pars <- c(r = 1, q=80, alpha = 0.01, beta = 1, K1=100)
#Tamanho inicial das populacoes
State <- c(V = 10, P =10)
#Tempo a ser rodado a solucao numerica
Time <- seq(0, 50, by = 1)

out <- as.data.frame(ode(func = Predacao, y = State, parms = Pars, times = Time))
colnames(out)=c("tempo", "V", "P")
out

par(mfrow=c(1,2)) #Plotando 2 resultados na mesma tela

# Plotando os tamanhos de cada especie pelo tempo
plot(out$V~out$tempo , type = "l", xlab = "tempo", ylab = "Numero de individuos", 
     col=1, lty=1, ylim=c(0,110), xlim=c(0,20))
lines(out$P~out$tempo , type = "l", col=2, lty=2)
legend("topright", c("V", "P"), lty = c(1,2), col = c(1,2), box.lwd = 0)

# Plotando a densidade de uma sp pela outra (escala de cinza de acordo com o tempo, 
# onde mais claro eh o inicio)
cinzas=  gray.colors(dim(out)[1], start = 0.8, end = 0, gamma = 0.7, alpha = NULL)
plot(1, type="n", xlab="V", ylab="P", xlim=c(0,150), ylim= c(0,110))
abline(a=Pars[1]/Pars[3],b=-(Pars[1]/Pars[3])/Pars[5],lty=1)
abline(v=Pars[2]/Pars[4],lty=2)
legend("topright", c("dV/dt=0", "dP/dt=0"), lty = c(1,2), box.lwd = 0)
points(out$P~out$V, col=cinzas, pch=16)
dev.off()


#################################################################################
###### Caso 3: Predadores persistem e presas vao a extincao. Caso estranho ######
#################################################################################

#Parametros do modelo
Pars <- c(r = 1, q=0.05, alpha = 0.01, beta = 1, K1=0.1)
#Tamanho inicial das populacoes
State <- c(V = 10, P =10)
#Tempo a ser rodado a solucao numerica
Time <- seq(0, 300, by = 1)

out <- as.data.frame(ode(func = Predacao, y = State, parms = Pars, times = Time))
colnames(out)=c("tempo", "V", "P")
out


par(mfrow=c(1,2)) #Plotando 2 resultados na mesma tela

# Plotando os tamanhos de cada especie pelo tempo
plot(out$V~out$tempo , type = "l", xlab = "tempo", ylab = "Numero de individuos", 
     col=1, lty=1, ylim=c(0,110))
lines(out$P~out$tempo , type = "l", col=2, lty=2)
legend("topright", c("V", "P"), lty = c(1,2), col = c(1,2), box.lwd = 0)

# Plotando a densidade de uma sp pela outra (escala de cinza de acordo com o tempo, 
# onde mais claro eh o inicio)
cinzas=  gray.colors(dim(out)[1], start = 1, end = 0, gamma = 0.7, alpha = NULL)
plot(1, type="n", xlab="V", ylab="P", xlim=c(0,5), ylim= c(0,50))
abline(a=Pars[1]/Pars[3],b=-(Pars[1]/Pars[3])/Pars[5],lty=1)
abline(v=Pars[2]/Pars[4],lty=2)
legend("topright", c("dV/dt=0", "dP/dt=0"), lty = c(1,2), box.lwd = 0)
points(out$P~out$V, col=cinzas, pch=16)
dev.off()