####################################################
######## Modelos em tempo continuo lineares ########
####################################################
install.packages("deSolve")

library(deSolve)

rm(list=ls())
########

#### 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

# 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.5, alpha = NULL)
plot(out$N2~out$N1, col=cinzas, ylim= c(0,30), ylab="N2", xlab="N1")

######################
##### Jacobiana: #####
######################


#Valores dos parametros
r1 = 0.05 
r2=0.03
beta1 = 1.2
beta2 = 0.5
K1=80
K2=50

# Valores do equilibrio de interesse

N1bar=(K1-beta1*K2)/(1-beta1*beta2)
N2bar=(K2-beta2*K1)/(1-beta1*beta2)

# Entradas da Jacobiana:
j11=r1/K1 - 2*r1*N1bar/K1 - r1*beta1*N2bar/K1
j12=-r1*N1bar*beta1/K1
j21=r2*N2bar*beta2/K2
j22=r2/K2 - 2*r2*N2bar/K2 - r2*beta2*N1bar/K2

# Montando a matriz jacobiana
jac=matrix(c(j11,j12,
             j21,j22), ncol=2, nrow=2, byrow=T)

# verificando autovalores da matriz jacobiana avaliada no ponto de equilibrio de interesse:
EV=eigen(jac)
EV

##### Caso 2 (uma das spp vai a extincao):
#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")

# 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.5, alpha = NULL)
plot(out$N2~out$N1, col=cinzas, ylim= c(0,30), ylab="N2", xlab="N1")


######################
##### Jacobiana: #####
######################

# Valores do equilibrio de interesse
N1bar=40
N2bar=0

#Valores dos parametros
r1 = 0.05 
r2=0.03
beta1 = 0.2
beta2 = 2
K1=40
K2=50

# Entradas da Jacobiana:
j11=r1/K1 - 2*r1*N1bar/K1 - r1*beta1*N2bar/K1
j12=-r1*N1bar*beta1/K1
j21=r2*N2bar*beta2/K2
j22=r2/K2 - 2*r2*N2bar/K2 - r2*beta2*N1bar/K2

#Montando a matriz jacobiana
jac=matrix(c(j11,j12,
             j21,j22), ncol=2, nrow=2, byrow=T)

# verificando autovalores da matriz jacobiana avaliada no ponto de equilibrio de interesse:
EV=eigen(jac)
EV