######################################
######## Modelos nao lineares ########
######################################

###################################################
#### Modelo de Saturacao x(n+1)=Kx(n)/(b+x(n)) ####
###################################################

# Removendo todos os objetos previamente guardados
rm(list=ls())

K=10
b=0.5

# Ponto 1 de equilibrio calculado em aula
xbar.1=0
 
tempo.Max=200 # definindo tempo maximo a ser rodado
x=rep(0,tempo.Max) # vetor onde guardaremos os valores
x[1]=xbar.1 # comeco exatamente no pto de eq 1

for(n in 1:(tempo.Max-1)){ #calculando iterativamente
  x[n+1]=K*x[n]/(b+x[n])
}

tempo=seq(0,(tempo.Max-1))
dados=cbind(tempo,x)

dados=as.data.frame(dados)
colnames(dados)=c("tempo", "x")
tail(dados) # fim da tabela com os dados juntos (tempo, x(n))

# Plot da populaçao com saturacao comecando no ponto de eq 1 sem perturbacao
plot(dados$x~dados$tempo, type="o",pch=16, xlab="tempo", ylab="x(n)", 
     ylim=c(0,10), xlim=c(0,50), main="Ponto eq 1 sem perturbacao")


###### Fazendo a mesma coisa pro segundo ponto de equilibrio

# Ponto 2 de equilibrio calculado em aula
xbar.2=K-b

tempo.Max=200 # definindo tempo maximo a ser rodado
x=rep(0,tempo.Max) # vetor onde guardaremos os valores
x[1]=xbar.2 # comecamndo exatamente no pto de eq 2

for(n in 1:(tempo.Max-1)){ #calculando iterativamente
  x[n+1]=K*x[n]/(b+x[n])
}

tempo=seq(0,(tempo.Max-1))
dados=cbind(tempo,x)

dados=as.data.frame(dados)
colnames(dados)=c("tempo", "x")
tail(dados) # fim da tabela de dados

# Plot da populaçao com saturacao comecando no ponto de eq 2 sem perturbacao
plot(dados$x~dados$tempo, type="o",pch=16, xlab="tempo", ylab="x(n)", 
     ylim=c(0,10), xlim=c(0,50), main="Ponto eq 2 sem perturbacao")

#########################################
######### Efetuando perturbacoes ########
########################################

K=10
b=0.5

# Ponto 1 de equilibrio calculado em aula
xbar.1=0

eps=0.5 # perturbacao

tempo.Max=200 # definindo tempo maximo a ser rodado
x=rep(0,tempo.Max) # vetor onde guardaremos os valores
x[1]=xbar.1 + eps #

for(n in 1:(tempo.Max-1)){ #calculando iterativamente
  x[n+1]=K*x[n]/(b+x[n])
}

tempo=seq(0,(tempo.Max-1))
dados=cbind(tempo,x)

dados=as.data.frame(dados)
colnames(dados)=c("tempo", "x")
head(dados)

# Plot da populaçao com saturacao comecando no ponto de eq 1 com perturbacao
plot(dados$x~dados$tempo, type="o",pch=16, xlab="tempo", ylab="x(n)", 
     ylim=c(0,10), xlim=c(0,50), main="Ponto eq 1 com perturbacao")


###### Fazendoperturbacao no 2o ponto de equilibrio

# Ponto 2 de equilibrio calculado em aula
xbar.2=K-b

eps=2 # perturbacao um pouco maior que antes

tempo.Max=200 # definindo tempo maximo a ser rodado
x=rep(0,tempo.Max) # vetor onde guardaremos os valores
x[1]=xbar.2 + eps # aqui poderiamos ter feito a perturbacao para baixo tbm

for(n in 1:(tempo.Max-1)){ #calculando iterativamente
  x[n+1]=K*x[n]/(b+x[n])
}

tempo=seq(0,(tempo.Max-1))
dados=cbind(tempo,x)

dados=as.data.frame(dados)
colnames(dados)=c("tempo", "x")
head(dados)

# Plot da populaçao com saturacao comecando no ponto de eq 2 com perturbacao
plot(dados$x~dados$tempo, type="o",pch=16, xlab="tempo", ylab="x(n)", 
     ylim=c(5,15), xlim=c(0,50), main="Ponto eq 2 com perturbacao")



################ Como fica a taxa de crescimento
# Vamos ver como ficou a taxa de crescimento para o ultimo caso calculado

r=rep(0,tempo.Max)

tempo.Max=200

for(n in 1:(tempo.Max)){
  r[n]=K/(b+dados$x[n])
}


tempo=seq(0,(tempo.Max-1))
dados.r=cbind(tempo,r)

dados.r=as.data.frame(dados.r)
colnames(dados.r)=c("tempo", "r")
head(dados.r)

plot(dados.r$r~dados.r$tempo, type="o",pch=16, xlab="tempo", ylab="r(n)", 
     ylim=c(0.8,1.2), xlim=c(0,50), main="taxa de crescimento")







##########################################################
#### Modelo da Eq Logistica x(n+1)=r0 x(n)*(1-x(n)/K) ####
##########################################################

### Faremos ja com perturbacoes no ponto de eq ja que comecar 
### diretamente no ponto de eq nao altera nada ao longo do tempo

K=100
r0=1.1 # altere para 0.8 para ver o que acontece

# Ponto 1 de equilibrio calculado em aula
xbar.1=0

eps=1 # perturbacao

tempo.Max=200 # definindo tempo maximo a ser rodado
x=rep(0,tempo.Max) # vetor onde guardaremos os valores
x[1]=xbar.1 + eps # aqui poderiamos ter feito a perturbacao para baixo tbm (mas pop ficaria < 0)

for(n in 1:(tempo.Max-1)){ #calculando iterativamente
  x[n+1]=r0*x[n]*(1-x[n]/K)
}

tempo=seq(0,(tempo.Max-1))
dados.1=cbind(tempo,x)

dados.1=as.data.frame(dados.1)
colnames(dados.1)=c("tempo", "x")
head(dados.1)

# Plot da populaçao cresc. logistico comecando no ponto de eq 1 com perturbacao
plot(dados.1$x~dados.1$tempo, type="o",pch=16, xlab="tempo", ylab="x(n)", 
     ylim=c(0,11), xlim=c(0,70), main="Ponto eq 1 com perturbacao")



###### Fazendo perturbacao no 2o ponto de equilibrio

K=100
r0=1.1 # altere para 0.8 para ver o que acontece

# Ponto 2 de equilibrio calculado em aula
xbar.2=K*(1-1/r0)       # note que para r0<1 este ponto tem valor negativo

eps=2 # perturbacao

tempo.Max=200 # definindo tempo maximo a ser rodado
x=rep(0,tempo.Max) # vetor onde guardaremos os valores
x[1]=xbar.2 + eps # aqui poderiamos ter feito a perturbacao para baixo tbm

for(n in 1:(tempo.Max-1)){ #calculando iterativamente
  x[n+1]=r0*x[n]*(1-x[n]/K)
}

tempo=seq(0,(tempo.Max-1))
dados.2=cbind(tempo,x)

dados.2=as.data.frame(dados.2)
colnames(dados.2)=c("tempo", "x")
head(dados.2)

# Plot da populaçao cresc. logistico comecando no ponto de eq 2 com perturbacao
plot(dados.2$x~dados.2$tempo, type="o",pch=16, xlab="tempo", ylab="x(n)", 
     ylim=c(5,15), xlim=c(0,50), main="Ponto eq 2 com perturbacao")


################ Como ficam as taxas de crescimento
# Vamos ver como ficou a taxa de crescimento para o os dois casos calculados

#### 1o caso:
r=rep(0,tempo.Max)

tempo.Max=200

for(n in 1:(tempo.Max)){
  r[n]=r0*(1-dados.1$x[n]/K)
}


tempo=seq(0,(tempo.Max-1))
dados1.r=cbind(tempo,r)

dados1.r=as.data.frame(dados1.r)
colnames(dados1.r)=c("tempo", "r")
head(dados1.r)

plot(dados1.r$r~dados1.r$tempo, type="o",pch=16, xlab="tempo", ylab="r(n)", 
     ylim=c(0.8,1.2), xlim=c(0,100), main="taxa de crescimento")



#### 2o caso:

r=rep(0,tempo.Max)

tempo.Max=200

for(n in 1:(tempo.Max)){
  r[n]=r0*(1-dados.2$x[n]/K)
}


tempo=seq(0,(tempo.Max-1))
dados2.r=cbind(tempo,r)

dados2.r=as.data.frame(dados2.r)
colnames(dados2.r)=c("tempo", "r")
head(dados2.r)

plot(dados2.r$r~dados2.r$tempo, type="o",pch=16, xlab="tempo", ylab="r(n)", 
     ylim=c(0.9,1.1), xlim=c(0,100), main="taxa de crescimento")