Este tutorial pressupõe que você leu o tutorial sobre integração numérica.
Diagramas de bifurcação representam as soluções (de longo prazo) de um modelo em função de algum parâmetro, cujo efeito queremos investigar. A ideia é que, à medida que esse parâmetro muda, as soluções mudam de forma “bem comportada”, e isso nos ajuda a entender melhor o comportamento geral do modelo.
Neste tutorial, vamos estudar um modelo simples de predador-presa (o Rosenzweig-MacArthur), e ver como a quantidade de recursos por presa (\(K\)) muda a dinâmica.
Este modelo é expresso como:
\[ \begin{aligned} \frac{dR}{dt} &= rR \left( 1 - \frac{R}{K} \right) - \frac{a R C}{1+ahR} \\ \frac{dC}{dt} &= \frac{e a R C}{1+ahR} - d C \end{aligned} \]
Usamos o mesmo método já explicado aqui para integrar este modelo numericamente:
##Carregando bibliotecas
library(deSolve)
### Modelo Rosenzweig-MacArthur###
## Sistema de EDOs a ser integrado
RM <- function(t, y, parms){
with(as.list(c(y, parms)),
{
dR = r*R*(1-R/K) - a*R*C/(1+a*h*R)
dC = e*a*R*C/(1 + a*h*R) - d*C
return(list(c(R=dR,C=dC)))
})
}
##parâmetros
parms1 = c(r=1, K=10, a=1, h=0.1, e=0.1, d=0.1)
## condições iniciais
y0 = c(R=1,C=1)
## rodando a integração numérica do tempo 1 a 1000
out = ode(y=y0, times = seq(from=1, to=1000, by=0.5), func = RM, parms = parms1)
##plotando a solução
matplot(x=out[,1], y=out[,2:3], type = "l", lwd=2, lty=1,
col=c("blue", "darkgreen"), xlab="Time", ylab="Population size")
legend("topright", legend = c("Resource","Consumer"),
lty = 1, lwd=2, col=c("blue", "darkgreen"))
Para os parâmetros escolhidos acima, a solução de longo prazo (assintótica) é um ponto fixo. Vamos ver isso no espaço de fase, ou seja, o espaço de Predadores (\(P\)) vs. Presas (\(V\)). Notamos que as setas estão “circulando”, mas sempre apontam para dentro, e assim a trajetória se move em direção ao meio, ao ponto fixo.
## Espaço de fases e pontos fixos
##primeiro precisamos calcular as coordenadas para os vetores
xs = seq(from=0.9, to=1.3, length.out=7)
ys = seq(from=0.94, to=1.04, length.out=7)
coords = expand.grid(R=xs, C=ys)
## E, em seguida, faça um loop sobre todas as coordenadas para calcular as derivadas nesses pontos
dvs = matrix(NA, ncol=2, nrow=nrow(coords))
for(i in 1:nrow(coords))
dvs[i,] = unlist(RM(t=1, y = coords[i,], parms = parms1))
## agora plotamos a trajetória
plot(x=out[,2], y=out[,3], type="l", lwd=2, col="blue",
xlab="Resource",ylab="Consumer", xlim=c(0.9,1.3), ylim=c(0.94,1.04))
## e adicionamos o campo vetorial
arrows(x0=coords[,1], y0=coords[,2], x1=coords[,1]+dvs[,1]*0.5,
y1=coords[,2]+dvs[,2]*0.5, length=0.1, lwd=2)
Aumentando a capacidade de carga \(K\) de \(10\) para \(15\), agora vemos oscilações…
# agora K = 15
parms2 = c(r=1, K=15, a=1, h=0.1, e=0.1, d=0.1)
out2 = ode(y=y0, times = seq(from = 1, to = 1000, by=0.5), func = RM, parms = parms2)
##plotando a solução
matplot(x=out2[,1], y=out2[,2:3], type = "l", lwd=2, lty=1,col=c("blue", "darkgreen"),
xlab="Time", ylab="Population size")
legend("topleft", legend = c("Resource","Consumer"), lty = 1, lwd=2,
col=c("blue", "darkgreen"))
Olhando novamente para o gráfico de espaço de fase, vemos agora que o fluxo (as setas) de dentro apontam para fora, em direção a um ciclo limite, e as setas de fora apontam para dentro. O ciclo limite corresponde à solução periódica que acabamos de ver.
## calculando os vetores novamente
xs = seq(from=0,to=6,length.out=8)
ys = seq(from=0,to=2,length.out=8)
coords = expand.grid(R=xs,C=ys)
dvs = matrix(NA, ncol=2, nrow=nrow(coords))
for(i in 1:nrow(coords))
dvs[i,] = unlist(RM(t=1, y = coords[i,], parms = parms2))
##A trajetoria
plot(x=out2[,2],y=out2[,3], type="l", lwd=2, col="blue",
xlab="Resource", ylab="Consumer", xlim=c(0,6),ylim=c(0,2))
##e os vetores
arrows(x0=coords[,1], y0=coords[,2], x1=coords[,1]+dvs[,1]*0.5,
y1=coords[,2]+dvs[,2]*0.5, length=0.1, lwd=2)
## Warning in arrows(x0 = coords[, 1], y0 = coords[, 2], x1 = coords[, 1] + : zero-
## length arrow is of indeterminate angle and so skipped
Vimos as soluções para dois valores de \(K\), \(10\) e \(15\), então queremos plotá-los como uma função de \(K\). No segundo caso, há oscilações, então ao invés de pegar toda a solução, escolhemos apenas o mínimo e o máximo da solução (depois de muito tempo). Quando a solução é um ponto fixo, o mínimo e o máximo devem coincidir.
##plotando tamanhos de população mínimo e máximo com diferentes valores de K
##objeto com os números de linha que usaremos para o gráfico a seguir
lines = (nrow(out)-500):nrow(out)
##criando um gráfico vazio
plot(0.1, type="n", xlim=c(0,20), ylim=range(out2[lines,2:3]), log="y",
xlab="K", ylab="Min and Max population")
##pontos para o K = 10
points(x = c(10,10), y = range(out[lines,3]), pch=21,bg="darkgreen", cex=3)
points(x = c(10,10), y = range(out[lines,2]), pch=21,bg="blue", cex=3)
##pontos para o K = 15
points(c(15,15), range(out2[lines,3]), pch=21,bg="darkgreen", cex=3)
points(c(15,15), range(out2[lines,2]), pch=21,bg="blue", cex=3)
Este é um diagrama de bifurcação pouco informativo: tem apenas dois pontos em \(K\)! Vamos tentar com muitos valores de \(K\).
O que acontece quando alteramos a capacidade de carga \(K\) de valores muito pequenos para valores muito grandes? Para valores muito pequenos, o recurso não vai sustentar a população consumidora, mas para valores maiores ok \(K\), ambas as espécies devem ser beneficiadas… certo?
## este bloco calcula soluções para muitos K's, deve levar algum tempo
KK = seq(from = 0.5, to=25, by=0.5)
rminmax = matrix(NA, ncol=2, nrow=length(KK))
cminmax = matrix(NA, ncol=2, nrow=length(KK))
## Faça um loop sobre todos os valores de K e obtenha os tamanhos mínimos e máximos da população
for(i in 1:length(KK)){
parmsi = c(r=1, K=KK[i], a=1, h=0.1, e=0.1, d=0.1)
y0 = c(R=1,C=1)
out3 = ode(y=y0, times = seq(from = 1, to = 1000, by=0.5), func = RM, parms = parmsi)
rminmax[i,] = range(out3[(nrow(out3)-500):nrow(out3),2])
cminmax[i,] = range(out3[(nrow(out3)-500):nrow(out3),3])
}
plot(x=KK, y=rminmax[,1], type="l", lwd=2, col="blue",ylim=range(rminmax), log="y",
xlab="K", ylab="Min and Max population")
points(x=KK, y=rminmax[,2], type="l", lwd=2, col="blue")
points(x=KK, y=cminmax[,1], type="l", lwd=2, col="darkgreen",ylim=range(rminmax))
points(x=KK, y=cminmax[,2], type="l", lwd=2, col="darkgreen",ylim=range(rminmax))
Bem, a primeira previsão foi OK (observe que o gráfico acima usa uma escala logarítmica), mas para \(K\) altos, os mínimos da oscilação vão para valores muito baixos, então as populações têm alto risco de extinção. Esse fenômeno é o chamado paradoxo do enriquecimento.
Um tipo especial de diagrama de bifurcação pode ser usado quando temos parâmetros que oscilam com o tempo, e queremos ver como isso interage com o sistema. Vamos considerar novamente as equações de Rosenzweig-MacArthur, mas agora fazemos \(r\), a taxa de crescimento da presa, oscilar senoidalmente no tempo:
\[ \begin{aligned} \frac{dR}{dt} &= r(t) R \left( 1 - \frac{R}{K} \right) - \frac{a R C}{1+ahR} \\ \frac{dC}{dt} &= \frac{e a R C}{1+ahR} - d C \\ r(t) &= r_0 (1+\alpha \sin(2\pi t/T)) \end{aligned} \]
Integramos de forma usual:
### Dinâmica dos recursos do consumidor em um ambiente sazonal ###
## sequência temporal
time <- seq(0, 2000, by = .5)
## paramêtros: vetor
parameters <- c(r0=1, alpha=0.1, T=80, K=10, a=1, h=0.1, e=0.1, d=0.1)
## condições iniciais: vetor
state <- c(R = 1, C = 1)
## Função R para calcular o valor das derivadas em cada valor de tempo
## Use os nomes das variáveis conforme definido nos vetores acima
RM2 <- function(t, state, parameters){
with(as.list(c(state, parameters)), {
r = r0 * (1 + alpha*sin(2*pi*t/T))
dR = R * ( r*(1 - R/K) - a*C / (1 + a*h*R) )
dC = e*a*R*C / (1 + a*h*R) - d*C
return(list(c(dR, dC)))
})
}
## Integrando com ode
out <- ode(y = state, times = time, func = RM2, parms = parameters)
## plotando
matplot(x = out[,1], y = out[,2:3], type="l", lwd=2, lty = 1,
col=c("blue", "darkgreen"), xlab = "Time", ylab = "Population size")
legend("topright", c("Resource", "Consumer"), lty=1, lwd=2, col=c("blue", "darkgreen"))
Observe que, mesmo com \(K\) pequeno, as soluções oscilam devido à oscilação de \(r(t)\).
Agora usamos uma ferramenta que é a favorita de todos os tempos dos físicos: o diagrama de ressonância. Funciona exatamente como um diagrama de bifurcação, mas o parâmetro que é alterado é o período (ou frequência) da oscilação externa.
## Um diagrama de ressonância ##
## Nova sequência de tempo
time <- seq(0, 6000, by = 1)
## Sequência de valores de T
TT <- seq(1, 80, by = 2)
## Uma matriz para armazenar os resultados
results <- matrix(ncol=4, nrow=length(TT),
dimnames=list(NULL, c("R.min","R.max","C.min","C.max")))
## Loop sobre TT
for(i in 1:length(TT)){
parameters <- c(r0=1, alpha=0.1, T=TT[i], K=10, a=1, h=0.1, e=0.1, d=0.1)
tmp1 <- ode(y = state, times = time, func = RM2, parms = parameters)
results[i,1:2] <- range(tmp1[1001:nrow(tmp1), 2])
results[i,3:4] <- range(tmp1[1001:nrow(tmp1), 3])
}
## Plotando
plot(R.min ~ TT , data=results, type="l", lwd=2, lty = 1,
col="blue", xlab = "T", ylab = "Min / Max population size",
log="y", ylim = range(results))
lines(R.max ~ TT, data=results, type="l", lwd=2, lty = 1, col=c("blue"))
lines(C.min ~ TT, data=results, type="l", lwd=2, lty = 1, col=c("darkgreen"))
lines(C.max ~ TT, data=results, type="l", lwd=2, lty = 1, col=c("darkgreen"))
legend("topright", c("Resource", "Consumer"), lty=1, lwd=2, col=c("blue", "darkgreen"))
Vemos um pico forte! (lembre-se que esta é uma escala logarítmica). A frequência em que esse pico ocorre é a frequência de ressonância do sistema, e está relacionada com a frequência natural do sistema (que existe mesmo quando vai para um ponto fixo com parâmetros constantes!). A oscilação externa excita a frequência natural e aciona ciclos de grande amplitude, como quando empurramos uma gangorra (ou gangorra, o balancín).
Primeiro: boa sorte. Em segundo lugar, você provavelmente terá que amostrar o espaço, em vez de passar por tudo. O método recomendado para fazer isso são amostras de hipercubo latino ( [Latin Hypercube samples])(http://en.wikipedia.org/wiki/Latin_hypercube_sampling), que usa uma amostragem aleatória de combinaçãoes de parâmetros, garantindo uma boa varredura do espaço de parâmetros. Observe, no entanto, que esse método é uma maneira de amostrar o espaço de parâmetros e fazer estatísticas úteis com ele, portanto, o resultado só fará sentido se você souber interpretar adequadamente os resultados. Dito isso, existem implementações para R e python: