###########################################
##### 2 especies -  Modelo Nao Linear #####
###########################################

####### 
####### Caso que a extincao eh estavel:
####### 

rm(list=ls())

T.max=200

P=matrix(0,T.max,1)
H=matrix(0,T.max,1)

# Primeiro valor dos vetores das spp
P0=3
H0=5
P[1]<-P0
H[1]<-H0

a=1
f=0.5

#Calculando para cada  instante de tempo
for (n in 1:T.max){
  H[n+1]= f*H[n]/(1+a*P[n])
  P[n+1]= H[n]*(1-1/(1+a*P[n]))
}

# Juntando os dados
tempo<-seq(0,T.max)
dados=data.frame(tempo,H,P)
head(dados)

# Fazendo o grafico de Numero de individuos pelo tempo:
plot(dados$tempo, dados$H, type="b", pch=1,xlim=c(0,10), ylim=c(0,10),col="red", 
     xlab="tempo", ylab="Numero de Individuos") 
lines(dados$tempo, dados$P, pch=18, col="blue", type="b", lty=2)
legend(5, 10, legend=c("H", "P"),
       col=c("red", "blue"),pch=c(1,18), cex=0.8)


# Fazendo o grafico de Hospedeiros x Parasitoides (retrato de fases):
plot(dados$H, dados$P, type="b", pch=16,xlim=c(0,7), ylim=c(0,7), 
     xlab="Hospedeiro", ylab="Parasitoide") 



#### Pela jacobiana:

# Ponto de equilibrio 1:
P.eq=0
H.eq=0

j11=f/(1+a*P.eq)
j12=-a*f*H.eq/(1+a*P.eq^2)
j21=1-1/(1+a*P.eq)
j22=-a*H.eq/(1+a*P.eq^2)

Jac=matrix(c(j11,j12,j21,j22), ncol=2,nrow=2, byrow=T)

eigen(Jac)

# Ponto de equilibrio 2:
P.eq=(f-1)/a
H.eq=f/a

P.eq
H.eq

####### 
####### Caso que a extincao nao eh estavel
####### 

rm(list=ls())

T.max=200

P=matrix(0,T.max,1)
H=matrix(0,T.max,1)

# Primeiro valor dos vetores das spp
P0=3
H0=5
P[1]<-P0
H[1]<-H0

a=1
f=2

#Calculando para cada  instante de tempo
for (n in 1:T.max){
  H[n+1]= f*H[n]/(1+a*P[n])
  P[n+1]= H[n]*(1-1/(1+a*P[n]))
}

# Juntando os dados
tempo<-seq(0,T.max)
dados=data.frame(tempo,H,P)
head(dados)

# Fazendo o grafico de Numero de individuos pelo tempo:
plot(dados$tempo, dados$H, type="b", pch=1,xlim=c(0,50), ylim=c(0,10),col="red", 
     xlab="tempo", ylab="Numero de Individuos") 
lines(dados$tempo, dados$P, pch=18, col="blue", type="b", lty=2)
legend(5, 10, legend=c("H", "P"),
       col=c("red", "blue"), pch=c(1,18), cex=0.8)


# Fazendo o grafico de Hospedeiros x Parasitoides (retrato de fases):
plot(dados$H, dados$P, type="b", pch=16,xlim=c(0,7), ylim=c(0,7),
     xlab="Hospedeiro", ylab="Parasitoide") 


#### Pela jacobiana:

# Ponto de equilibrio 1: 
P.eq=0
H.eq=0

j11=f/(1+a*P.eq)
j12=-a*f*H.eq/(1+a*P.eq)^2
j21=1-1/(1+a*P.eq)
j22=a*H.eq/(1+a*P.eq)^2

Jac=matrix(c(j11,j12,j21,j22), ncol=2,nrow=2, byrow=T)

eigen(Jac)

# Ponto de equilibrio 2: 
P.eq=(f-1)/a
H.eq=f/a

P.eq
H.eq

j11=f/(1+a*P.eq)
j12=-a*f*H.eq/(1+a*P.eq)^2
j21=1-1/(1+a*P.eq)
j22=a*H.eq/(1+a*P.eq)^2

Jac2=matrix(c(j11,j12,j21,j22), ncol=2,nrow=2, byrow=T)

eigen(Jac2)



####### Caso que a extincao nao eh estavel e aumentando o valor de f

rm(list=ls())

T.max=200

P=matrix(0,T.max,1)
H=matrix(0,T.max,1)

# Primeiro valor dos vetores das spp
P0=3
H0=5

a=1

#Calculando para cada  instante de tempo
Par.Hosp<-function(par.a,par.f,H.0,P.0,t.max){
  P=matrix(0,T.max,1)
  H=matrix(0,T.max,1)
  P[1]=P.0
  H[1]=H.0
  a=par.a
  f=par.f
    for (n in 1:t.max){
    H[n+1]= f*H[n]/(1+a*P[n])
    P[n+1]= H[n]*(1-1/(1+a*P[n]))
    }
  tempo=seq(0,t.max)
  dados=data.frame(tempo,H,P)
  return(dados)
}

#Juntando os dados
dados.f2=Par.Hosp(a,2,P0,H0,T.max)

dados.f5=Par.Hosp(a,5,P0,H0,T.max)

dados.f10=Par.Hosp(a,10,P0,H0,T.max)



plot(dados.f2$tempo, dados.f2$H, type="b", pch=16,xlim=c(0,20), ylim=c(0,40),col="#DCDCDC", 
     xlab="tempo", ylab="Numero de Hospedeiros") 
lines(dados.f5$tempo, dados.f5$H, pch=16, type="b", lty=2, col="#A9A9A9")
lines(dados.f10$tempo, dados.f10$H, pch=16, type="b", lty=2, col="#778899")
legend(18, 40, legend=c("f=2", "f=5", "f=10"),
       col=c("#DCDCDC", "#A9A9A9", "#778899"), pch=16, cex=0.8)




plot(dados.f2$tempo, dados.f2$P, type="b", pch=16,xlim=c(0,20), ylim=c(0,40),col="#DCDCDC", 
     xlab="tempo", ylab="Numero de Parasitoides") 
lines(dados.f5$tempo, dados.f5$P, pch=16, type="b", lty=2, col="#A9A9A9")
lines(dados.f10$tempo, dados.f10$P, pch=16, type="b", lty=2, col="#778899")
legend(18, 40, legend=c("f=2", "f=5", "f=10"),
       col=c("#DCDCDC", "#A9A9A9", "#778899"), pch=16, cex=0.8)