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

library(deSolve)


########

#### Definindo a funcao de competicao:
Competicao <- function (Time, State, Pars) {
  with(as.list(c(State, Pars)), {
    dN1 = r1*N1*((K1 - N1 - beta1*N2)/K1)
    dN2 = r2*N2*((K2 - N2 - beta2*N1)/K2)
    return(list(c(dN1, dN2)))
  })
}

#####
##### Caso 1 (spp coexistem):
#####

# Parametros do modelo
Pars <- c(r1 = 0.05, r2=0.03, beta1 = 1.2, beta2 = 0.5, K1=80, K2=50)
# Tamanho inicial das populacoes
State <- c(N1 = 10, N2 =10)
# Tempo a ser rodado a solucao numerica
Time <- seq(0, 300, by = 1)

# Saidas gravadas no arquivo out
out <- as.data.frame(ode(func = Competicao, y = State, parms = Pars, times = Time))
colnames(out)=c("tempo", "N1", "N2")
out

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

# Plotando os tamanhos de cada especie pelo tempo
plot(out$N1~out$tempo , type = "l", xlab = "tempo", 
     ylab = "Numero de individuos", col=1, lty=1,ylim=c(0,50))
lines(out$N2~out$tempo , type = "l", col=2, lty=2)
legend("topright", c("N1", "N2"), 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(out$N2~out$N1, col=cinzas, xlim=c(0,100),ylim= c(0,100), 
     ylab="N2", xlab="N1", pch=16)
abline(a=Pars[5]/Pars[3],b=-(Pars[5]/Pars[3])/Pars[5],lty=1)
abline(a=Pars[6],b=-Pars[6]/(Pars[6]/Pars[4]),lty=2)
legend("topright", c("dN1/dt=0", "dN2/dt=0"), lty = c(1,2), box.lwd = 0)

dev.off()

#####
##### Caso 2 (coexistencia instavel):
#####

#Parametros do modelo
Pars <- c(r1 = 0.05, r2=0.03, beta1 = 2, beta2 = 2, K1=40, K2=50)
#Tamanho inicial das populacoes
State <- c(N1 = 10, N2 =10)
#Tempo a ser rodado a solucao numerica
Time <- seq(0, 300, by = 1)
par=as.data.frame(Pars, byrow=TRUE)


#Saidas gravadas no arquivo out
out <- as.data.frame(ode(func = Competicao, y = State, parms = Pars, times = Time))
colnames(out)=c("tempo", "N1", "N2")
out

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

# Plotando os tamanhos de cada especie pelo tempo
plot(out$N1~out$tempo , type = "l", xlab = "tempo", ylab = "Numero de individuos", 
     col=1, lty=1, ylim=c(0,50))
lines(out$N2~out$tempo , type = "l", col=2, lty=2)
legend("topright", c("N1", "N2"), 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(out$N2~out$N1, col=cinzas, xlim=c(0,60),ylim= c(0,60), 
     ylab="N2", xlab="N1",  pch=16)
abline(a=Pars[5]/Pars[3],b=-(Pars[5]/Pars[3])/Pars[5],lty=1)
abline(a=Pars[6],b=-Pars[6]/(Pars[6]/Pars[4]),lty=2)
legend("topright", c("dN1/dt=0", "dN2/dt=0"), lty = c(1,2), box.lwd = 0)

dev.off()

#############
############# Outras condicoes iniciais para mesmos parametros do modelo: 
#############
Pars <- c(r1 = 0.05, r2=0.03, beta1 = 2, beta2 = 2, K1=40, K2=50)
#Tamanho inicial das populacoes
State <- c(N1 = 25, N2 =10)
#Tempo a ser rodado a solucao numerica
Time <- seq(0, 300, by = 1)
par=as.data.frame(Pars, byrow=TRUE)


#Saidas gravadas no arquivo out
out <- as.data.frame(ode(func = Competicao, y = State, parms = Pars, times = Time))
colnames(out)=c("tempo", "N1", "N2")
out

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

# Plotando os tamanhos de cada especie pelo tempo
plot(out$N1~out$tempo , type = "l", xlab = "tempo", ylab = "Numero de individuos", 
     col=1, lty=1, ylim=c(0,50))
lines(out$N2~out$tempo , type = "l", col=2, lty=2)
legend("topright", c("N1", "N2"), 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(out$N2~out$N1, col=cinzas, xlim=c(0,60),ylim= c(0,60), 
     ylab="N2", xlab="N1",  pch=16)
abline(a=Pars[5]/Pars[3],b=-(Pars[5]/Pars[3])/Pars[5],lty=1)
abline(a=Pars[6],b=-Pars[6]/(Pars[6]/Pars[4]),lty=2)
legend("topright", c("dN1/dt=0", "dN2/dt=0"), lty = c(1,2), box.lwd = 0)

dev.off()

##### Caso 3 (sp 1 sempre exclui a sp2):
#Parametros do modelo
Pars <- c(r1 = 0.05, r2=0.03, beta1 = 0.2, beta2 = 2, K1=40, K2=50)
#Tamanho inicial das populacoes
State <- c(N1 = 10, N2 =10)
#Tempo a ser rodado a solucao numerica
Time <- seq(0, 300, by = 1)

#Saidas gravadas no arquivo out
out <- as.data.frame(ode(func = Competicao, y = State, parms = Pars, times = Time))
colnames(out)=c("tempo", "N1", "N2")
out

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

# Plotando os tamanhos de cada especie pelo tempo
plot(out$N1~out$tempo , type = "l", xlab = "tempo", ylab = "Numero de individuos", 
     col=1, lty=1, ylim=c(0,50))
lines(out$N2~out$tempo , type = "l", col=2, lty=2)
legend("topright", c("N1", "N2"), 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(out$N2~out$N1, col=cinzas,xlim=c(0,100), ylim= c(0,210), 
     ylab="N2", xlab="N1", pch=16)
abline(a=Pars[5]/Pars[3],b=-(Pars[5]/Pars[3])/Pars[5],lty=1)
abline(a=Pars[6],b=-Pars[6]/(Pars[6]/Pars[4]),lty=2)
legend("topright", c("dN1/dt=0", "dN2/dt=0"), lty = c(1,2), box.lwd = 0)

dev.off()


#############
############# Outras condicoes iniciais para mesmos parametros do modelo: 
#############
Pars <- c(r1 = 0.05, r2=0.03, beta1 = 0.2, beta2 = 2, K1=40, K2=50)
#Tamanho inicial das populacoes
State <- c(N1 = 1, N2 =50)
#Tempo a ser rodado a solucao numerica
Time <- seq(0, 300, by = 1)


#Saidas gravadas no arquivo out
out <- as.data.frame(ode(func = Competicao, y = State, parms = Pars, times = Time))
colnames(out)=c("tempo", "N1", "N2")
out

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

# Plotando os tamanhos de cada especie pelo tempo
plot(out$N1~out$tempo , type = "l", xlab = "tempo", ylab = "Numero de individuos", 
     col=1, lty=1,ylim=c(0,50))
lines(out$N2~out$tempo , type = "l", col=2, lty=2)
legend("topright", c("N1", "N2"), 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(out$N2~out$N1, col=cinzas,xlim=c(0,100), ylim= c(0,210), 
     ylab="N2", xlab="N1",  pch=16)
abline(a=Pars[5]/Pars[3],b=-(Pars[5]/Pars[3])/Pars[5],lty=1)
abline(a=Pars[6],b=-Pars[6]/(Pars[6]/Pars[4]),lty=2)
legend("topright", c("dN1/dt=0", "dN2/dt=0"), lty = c(1,2), box.lwd = 0)

dev.off()

##### Caso 4 (sp 2 sempre exclui a sp 1):

#Parametros do modelo
Pars <- c(r1 = 0.05, r2=0.03, beta1 = 2, beta2 = 0.25, K1=40, K2=50)
#Tamanho inicial das populacoes
State <- c(N1 = 10, N2 =10)
#Tempo a ser rodado a solucao numerica
Time <- seq(0, 300, by = 1)
par=as.data.frame(Pars, byrow=TRUE)


#Saidas gravadas no arquivo out
out <- as.data.frame(ode(func = Competicao, y = State, parms = Pars, times = Time))
colnames(out)=c("tempo", "N1", "N2")
out

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

# Plotando os tamanhos de cada especie pelo tempo
plot(out$N1~out$tempo , type = "l", xlab = "tempo", ylab = "Numero de individuos", 
     col=1, lty=1,ylim=c(0,50))
lines(out$N2~out$tempo , type = "l", col=2, lty=2)
legend("topright", c("N1", "N2"), 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(out$N2~out$N1, col=cinzas, xlim=c(0,210),ylim= c(0,100), 
     ylab="N2", xlab="N1",  pch=16)
abline(a=Pars[5]/Pars[3],b=-(Pars[5]/Pars[3])/Pars[5],lty=1)
abline(a=Pars[6],b=-Pars[6]/(Pars[6]/Pars[4]),lty=2)
legend("topright", c("dN1/dt=0", "dN2/dt=0"), lty = c(1,2), box.lwd = 0)

dev.off()