##########################################
#### Plantas anuais com duas geracoes ####
##########################################

rm(list=ls())

S=20 #numero medio de smts por ind.
alpha=0.5 #germinacao
sigma=0.2 # sobrevivencia (1-m)

A=alpha*sigma*S
A
B=(sigma^2)*alpha*(1-alpha)*S
B

A^2 #condicao que A^2 > 4(-B) (raizes reais)
4*(-B)

R1=A/2+(sqrt(A^2-4*(-B)))/2 #calculo das raizes da equacao de segundo grau 
R2=A/2-(sqrt(A^2-4*(-B)))/2

R1 #veja como sao os valores de r1 e r2 e entenda como essa populacao cresce no tempo (grafico abaixo)
R2

p0=3 #condicoes iniciais
p1=5

tempo.Max=20
p=rep(0,tempo.Max)
p[1]=p0
p[2]=p1

for(n in 2:(tempo.Max-1)){ #calculando iterativamente cada tamanho populacional a cada instante de tempo
  p[n+1]=A*p[n]+B*p[n-1]
}

tempo=seq(0,(tempo.Max-1))

cbind(tempo,p)

plot(p~tempo, pch=19, xlab="n", ylab="p(n)")

# Calculando os parametros C1 e C2 (que dependem das condicoes iniciais) e que acompanham as raizes na solucao da equacao
C1=(p1-p0*R2)/(R1-R2)

C2=(-p1+p0*R1)/(R1-R2)


# Colocando a solucao geral na forma de funcao para poder calcular qq tempo sem precisar dos anteriores
Eq.Solucao <- function(c1,c2,r1,r2,n){ 
  pn=C1*(R1^n) + C2*(R2^n) 
  return(pn)
}

Eq.Solucao(C1,C2,R1,R2,8)


# Colocando a solucao geral com todas as contas dentro da funcao

Eq.Solucao2 <- function(a,b,P0,P1,n){ 
  r1=b/2+(sqrt(a^2-4*(-b)))/2
  r2=b/2-(sqrt(a^2-4*(-b)))/2
  c1=(P1-P0*r2)/(r1-r2)
  c2=(-P1+P0*r1)/(r1-r2)
  
  pn=C1*(R1^n) + C2*(R2^n)
  return(pn)
}

A
B
p0
p1

Eq.Solucao2(2,0.2,3,5,8)
##########################################################################
#### Plantas anuais com duas geracoes - mudando valores de parametros ####
##########################################################################

rm(list=ls())

S=20 #numero medio de smts por ind.
alpha=0.25 #germinacao
sigma=0.1 #sobrevivencia (1-m)

A=alpha*sigma*S
A
B=(sigma^2)*alpha*(1-alpha)*S
B

A^2 #condicao que A^2 > 4(-B) (raizes reais)
4*(-B)

R1=A/2+(sqrt(A^2-4*(-B)))/2 #calculo das raizes da equacao de segundo grau 
R2=A/2-(sqrt(A^2-4*(-B)))/2

R1 #veja como sao os valores de r1 e r2 e entenda como essa populacao cresce no tempo (grafico abaixo)
R2

p0=3 #condicoes iniciais
p1=5

tempo.Max=20
p=rep(0,tempo.Max)
p[1]=p0
p[2]=p1

for(n in 2:(tempo.Max-1)){ #calculando iterativamente cada tamanho populacional a cada instante de tempo
  p[n+1]=A*p[n]+B*p[n-1]
}

tempo=seq(0,(tempo.Max-1))

cbind(tempo,p)


plot(p~tempo, pch=19, xlab="n", ylab="p(n)")


C1=(p1-p0*R2)/(R1-R2)

C2=(-p1+p0*R1)/(R1-R2)


Eq.Solucao(C1,C2,R1,R2,8)

#################################
#### Exemplo 1 de oscilacoes ####
#################################

A=-0.5
B=0.1

#Notem que esses valores nao fazem sentido biologico, ja que A e B devem ser maiores que 0

A^2 #condicao que A^2 > 4(-B) (raizes reais)
4*(-B)

R1=A/2+(sqrt(A^2-4*(-B)))/2 #calculo das raizes da equacao de segundo grau 
R2=A/2-(sqrt(A^2-4*(-B)))/2

R1 #veja como sao os valores de r1 e r2 e entenda como essa populacao cresce no tempo (grafico abaixo)
R2

p0=3
p1=5

tempo.Max=20
p=rep(0,tempo.Max)
p[1]=p0
p[2]=p1

for(n in 2:(tempo.Max-1)){ #calculando iterativamente cada tamanho populacional a cada instante de tempo
  p[n+1]=A*p[n]+B*p[n-1]
}

tempo=seq(0,(tempo.Max-1))

cbind(tempo,p)


plot(p~tempo, pch=19, xlab="n", ylab="p(n)")
lines(p~tempo)

#################################
#### Exemplo 2 de oscilacoes ####
#################################


A=-2
B=0.5

A^2 #condicao que beta^2 > 4gama (raizes reais)
4*(-B)

R1=A/2+(sqrt(A^2-4*(-B)))/2 #calculo das raizes da equacao de segundo grau 
R2=A/2-(sqrt(A^2-4*(-B)))/2

R1 #veja como sao os valores de r1 e r2 e entenda como essa populacao cresce no tempo (grafico abaixo)
R2

p0=3
p1=5

tempo.Max=10
p=rep(0,tempo.Max)
p[1]=p0
p[2]=p1

for(n in 2:(tempo.Max-1)){ #calculando iterativamente cada tamanho populacional a cada instante de tempo
  p[n+1]=A*p[n]+B*p[n-1]
}

tempo=seq(0,(tempo.Max-1))

cbind(tempo,p)


plot(p~tempo, pch=19, xlab="n", ylab="p(n)")
lines(p~tempo)


plot(p~tempo, pch=19, xlab="n", ylab="p(n)", xlim=c(0,7), ylim=c(-215,480))
lines(p~tempo)