Modelo de Levins aprimorado para 2 espécies: mudanças entre as edições

De Física Computacional
Ir para navegação Ir para pesquisar
(Criou página com '{{Ecologia| Modelo de Lotka-Volterra amortecido | Modelo espacialmente explícito }} {{Ecologia| Modelo de Lotka-Volterra amortecido | Modelo espacialmente explícito }}')
 
Sem resumo de edição
 
(5 revisões intermediárias pelo mesmo usuário não estão sendo mostradas)
Linha 1: Linha 1:
{{Ecologia| [[Modelo de Lotka-Volterra amortecido]] | Modelo espacialmente explícito }}
{{Ecologia| [[Modelo de Lotka-Volterra amortecido]] | [[Modelo de Levins aprimorado para 2 espécies II]] }}
{{Ecologia| [[Modelo de Lotka-Volterra amortecido]] | Modelo espacialmente explícito }}
 
=== Comparação com outros modelos ===
Para duas espécies, vamos trabalhar com o seguinte sistema de equações:
 
<math>\begin{align}
\frac{dx}{dt} = & c_{1}x\left(1-D-x\right)-e_1x-\mu xy \\
\frac{dy}{dt} = & yc_{y}\left(x-yx\right)-ye_{y}\\
\end{align}</math>
 
Onde <math>x</math> é a presa e <math>y</math> o predador. Podemos comparar com outros modelos. Como um modelo inspirado no [[modelo de Levins]], podemos ver semelhanças bastante fortes desde já. Lembrando do modelo de Levins rapidamente:
 
<math>\frac{dp}{dt}=c\left(h-p\right)p-ep</math>
 
O termo <math>\left(1-D-x\right)</math> na equação da presa é a quantidade de fragmentos disponíveis para a colonização da presa, uma vez que precisamos descontar os fragmentos destruídos, além dos ocupados. Como todos parâmetros são em termos de proporções, temos probabilidade de selecionarmos um fragmento disponível para a colonização da presa. E o termo <math>\left(x_{1}-yx_{1}\right)</math> na equação do predador tem a mesma função, porém agora precisamos lembrar que os fragmentos que o predador pode ocupar são os que já estão ocupados pela presa, descontando os que já possuem predador. Então considerando propriedades de conjuntos, e olhando em termos de probabilidades esta equação, temos a probabilidade <math>P\left(x\right)</math>  de selecionar um fragmento ocupada pela presa, mas precisamos descontar a probabilidade de selecionar um fragmento que também já esteja ocupada por um predador <math>P\left(x\cap y\right)</math>. A maior diferença se dá porque a presa tem um termo a mais correspondendo ao decréscimo da população devido a predação, evento não considerado no modelo de Levins, mas considerado no modelo de Lotka-Volterra.
 
 
Por isso a próxima comparação que podemos fazer é com o modelo de [[Modelo de Lotka-Volterra amortecido|Lotka-Volterra amortecido]]:
 
 
<math>\begin{align}
\frac{dx}{dt} = & x\left(a-\alpha y\right)-kx^{2} \\
\frac{dy}{dt} = & y\left(-c+\gamma x\right)\\
\end{align}</math>
 
 
Vamos reescrever nosso sistema de equações:
 
<math>\begin{align}
\frac{dx}{dt} = & x\left(c_{1}\left(1-D\right)-\mu y\right)-c_{1}x^{2}-e_{1}x \\
\frac{dy}{dt} = & y\left(-e_{y}+c_{y}x\right)-c_{y}xy^{2}\\
\end{align}</math>
 
 
Olhando agora primeiro para a equação da presa, a maior diferença que encontramos, e é que temos o termo de extinção local que não é devido a predação <math>e_1x</math>,  algo já familiar do modelo de Levins.  E já quando olhamos para a equação do predador, temos o termo <math display="inline">\left(c_{y}x\right)y^{2}</math> , podemos pensar também como um termo logístico, impondo um limite na população de predadores, da mesma forma que ocorre com as presas. Porém agora o termo logístico não é uma constante, mas depende também da população de presas.
 
=== Análises dos pontos de equilíbrio ===
 
==== Sendo todos o parâmetros como variáveis ====
Vamos tentar fazer análises parecidas agora. Primeiro buscar os pontos de equilíbrio pro sistema:
 
<math display="block">\begin{array}{cc}
x\left[\left(c_{1}\left(1-D\right)-\mu y\right)-c_{1}x-e_{1}\right] & =0\\
y\left[\left(-e_{y}+c_{y}x\right)-c_{y}xy\right] & =0
\end{array}</math>
 
<div class="list">
 
<span> </span>Temos então 4 pontos de equilíbrio:
 
* Primeiro: <math display="inline">\left(0,0\right)\rightarrow</math> extinção de ambas as espécies
* Segundo: <math display="inline">\left(\left(1-D\right)-\frac{e_{1}}{c_{1}},0\right)\rightarrow</math> extinção do predador
</div>Os pontos de equilíbrio com ambas as espécies são extremamente mais complicados. Para o terceiro ponto temos:
 
<math display="block">\begin{array}{cc}
x= & -\frac{\sqrt{c_{y}^{2}\mu^{2}+\left(4c_{1}c_{y}e_{y}+2c_{y}^{2}e_{1}+\left(2D-2\right)c_{1}c_{y}^{2}\right)\mu+c_{y}^{2}e_{1}^{2}+\left(2D-2\right)c_{1}c_{y}^{2}e_{1}+\left(D^{2}-2D+1\right)c_{1}^{2}c_{y}^{2}}+c_{y}\mu+c_{y}e_{1}+\left(D-1\right)c_{1}c_{y}}{2c_{1}c_{y}}\\
y= & \frac{\sqrt{c_{y}^{2}\mu^{2}+\left(4c_{1}c_{y}e_{y}+2c_{y}^{2}e_{1}+\left(2D-2\right)c_{1}c_{y}^{2}\right)\mu+c_{y}^{2}e_{1}^{2}+\left(2D-2\right)c_{1}c_{y}^{2}e_{1}+\left(D^{2}-2D+1\right)c_{1}^{2}c_{y}^{2}}+c_{y}\mu-c_{y}e_{1}+\left(1-D\right)c_{1}c_{y}}{2c_{y}\mu}
\end{array}</math>
 
E o quarto e último ponto:
 
<math display="block">\begin{array}{cc}
x= & \frac{\sqrt{c_{y}^{2}\mu^{2}+\left(4c_{1}c_{y}e_{y}+2c_{y}^{2}e_{1}+\left(2D-2\right)c_{1}c_{y}^{2}\right)\mu+c_{y}^{2}e_{1}^{2}+\left(2D-2\right)c_{1}c_{y}^{2}e_{1}+\left(D^{2}-2D+1\right)c_{1}^{2}c_{y}^{2}}-c_{y}\mu-c_{y}e_{1}-\left(D-1\right)c_{1}c_{y}}{2c_{1}c_{y}}\\
y=- & \frac{\sqrt{c_{y}^{2}\mu^{2}+\left(4c_{1}c_{y}e_{y}+2c_{y}^{2}e_{1}+\left(2D-2\right)c_{1}c_{y}^{2}\right)\mu+c_{y}^{2}e_{1}^{2}+\left(2D-2\right)c_{1}c_{y}^{2}e_{1}+\left(D^{2}-2D+1\right)c_{1}^{2}c_{y}^{2}}-c_{y}\mu+c_{y}e_{1}-\left(1-D\right)c_{1}c_{y}}{2c_{y}\mu}
\end{array}</math>Estes últimos pontos (com o equilíbrio entre as espécies) são complicados demais para nos referirmos explicitamente daqui em diante quando necessário utilizamos alguma notação como <math display="inline">\left(x_{j},y_{j}\right)</math> para denotar o <math display="inline">j-\acute{e}simo</math> ponto de equilíbrio, sem escrevê-lo por inteiro. Vamos linearizar e analisar a estabilidade em torno dos pontos de equilíbrio. Para o ponto <math display="inline">\left(x_{1},y_{1}\right)=\left(0,0\right)</math> é muito simples:
 
<math display="block">\frac{dx}{dt}=x\left(c_{1}\left(1-D\right)-\mu y\right)-c_{1}x^{2}-e_{1}x</math><math display="block">\frac{dy}{dt}=y\left(-e_{y}+c_{y}x\right)-c_{y}xy^{2}</math>
 
Separando os termos lineares dos não lineares:
 
<math display="block">\frac{dx}{dt}=\left[x\left(c_{1}\left(1-D\right)-e_{1}\right)\right]-\left[\mu yx+c_{1}x^{2}\right]</math><math display="block">\frac{dy}{dt}=-\left[ye_{y}\right]+\left[c_{y}xy-c_{y}xy^{2}\right]</math>
 
Calculando então os limites:
 
<math display="block">\lim_{\left(x,y\right)\rightarrow\left(0,0\right)}-\frac{\mu yx+c_{1}x^{2}}{x\left(c_{1}\left(1-D\right)-e_{1}\right)}=\lim_{\left(x,y\right)\rightarrow\left(0,0\right)}-\frac{\mu y+c_{1}x}{c_{1}\left(1-D\right)-e_{1}}=0</math><math display="block">\lim_{\left(x,y\right)\rightarrow\left(0,0\right)}-\frac{c_{y}xy-c_{y}xy^{2}}{ye_{y}}=\lim_{\left(x,y\right)\rightarrow\left(0,0\right)}-\frac{c_{y}x-c_{y}xy}{e_{y}}=0</math>
 
Desprezando os termos não lineares:
 
<math display="block">\left(\begin{array}{c}
\dot{x}\\
\dot{y}
\end{array}\right)=\left(\begin{array}{cc}
c_{1}\left(1-D\right)-e_{1} & 0\\
0 & -e_{y}
\end{array}\right)\left(\begin{array}{c}
x\\
y
\end{array}\right)</math>
 
Logo um dos autovalores é real, sempre negativo <math display="inline">\left(-e_{y}\right)</math> e o segundo é <math display="inline">\left(c_{1}\left(1-D\right)-e_{1}\right)</math>. Basicamente se o termo de colonização do guanaco for maior que o da extinção <math display="inline">\left(c_{1}\left(1-D\right)>e_{1}\right)</math>, temos dois autovalores de sinais opostos, ou seja, um ponto de instabilidade. Aumentando a taxa de extinção do guanaco (ou a destruição do habitat) passamos a ter <math display="inline">c_{1}\left(1-D\right)<e_{1}</math>, e a extinção de ambas as espécies. Aqui podemos notar uma diferença dos modelos de Lotka-Volterra: a extinção de ambas as espécies pode ser um ponto de equilíbrio estável. Ou seja nos modelos anteriores o sistema nunca iria para o equilíbrio com extinção total, agora há essa possibilidade dependendo do nível de destruição do sistema.
 
Para o segundo ponto de equilíbrio <math display="inline">\left(x_{2},y_{2}\right)=\left(x_{2},0\right)</math> precisamos fazer a translação: <math display="inline">u=x-x_{2}</math> e <math display="inline">v=y</math>.
 
<div class="list">
 
<span> </span><math display="block">\frac{du}{dt}=\left[\left(-2c_{1}x_{2}+\left(1-D\right)c_{1}-e_{1}\right)u+\left(-\mu x_{2}\right)v+\left(-c_{1}x_{2}^{2}-e_{1}x_{2}+\left(1-D\right)x_{2}c_{1}\right)\right]+\left[-\mu uv-c_{1}u^{2}\right]</math><math display="block">\frac{dv}{dt}=\left[\left({\it c_{y}}x_{2}-{\it e_{y}}\right)v\right]+\left[-{\it c_{y}}v^{2}x_{2}-{\it c_{y}}uv^{2}+{\it c_{y}}uv\right]</math>
</div>Podemos ver facilmente que é semi linear para <math display="inline">\dot{u}</math> uma vez que ao fazermos os:
 
<math display="block">\lim_{\left(x,y\right)\rightarrow\left(0,0\right)}\frac{\text{parte não-linear}}{\text{parte linear}}</math><math display="block">\lim_{\left(u,v\right)\rightarrow\left(0,0\right)}\frac{-\mu uv-c_{1}u^{2}}{\left(-2c_{1}x_{2}+\left(1-D\right)c_{1}-e_{1}\right)u+\left(-\mu x_{2}\right)v}</math>
 
Fazendo a substituição <math display="inline">r^{2}=x^{2}+y^{2}</math> e <math display="inline">\left(x,y\right)=\left(r\cos\theta,r\sin\theta\right)</math>:
 
<math display="block">\lim_{r\rightarrow0}\frac{r^{2}\left(-\mu\cos\theta\sin\theta-c_{1}\cos\theta^{2}\right)}{\left(-2c_{1}x_{2}+\left(1-D\right)c_{1}-e_{1}\right)r\cos\theta+\left(-\mu x_{2}\right)r\sin\theta}</math><math display="block">\lim_{r\rightarrow0}\frac{r\left(-\mu\cos\theta\sin\theta-c_{1}\cos\theta^{2}\right)}{\left(-2c_{1}x_{2}+\left(1-D\right)c_{1}-e_{1}\right)\cos\theta+\left(-\mu x_{2}\right)\sin\theta}=0</math>
 
É interessante notar que a parte constante vai zerar. E para <math display="inline">\dot{v}</math>:
 
<math display="block">\lim_{\left(u,v\right)\rightarrow\left(0,0\right)}\frac{\left[-{\it c_{y}}v^{2}x_{2}-{\it c_{y}}uv^{2}+{\it c_{y}}uv\right]}{\left[\left({\it c_{y}}x_{2}-{\it e_{y}}\right)v\right]}=\frac{\left[-{\it c_{y}}vx_{2}-{\it c_{y}}uv+{\it c_{y}}u\right]}{\left[\left({\it c_{y}}x_{2}-{\it e_{y}}\right)\right]}=0</math>
 
Analisando o comportamento linear então na proximidade do ponto de equilíbrio:
 
<math display="block">\left(\begin{array}{c}
\dot{u}\\
\dot{v}
\end{array}\right)=\left(\begin{array}{cc}
-2c_{1}x_{2}+\left(1-D\right)c_{1}-e_{1} & -\mu x_{2}\\
0 & {\it c_{y}}x_{2}-{\it e_{y}}
\end{array}\right)\left(\begin{array}{c}
u\\
v
\end{array}\right)</math>E os autovalores são então:
 
* <math display="inline">\lambda_{1}=\left(1-D-\frac{e_{1}}{c_{1}}\right)c_{y}-e_{y}</math>
* <math display="inline">\lambda_{2}=e_{1}-\left(1-D\right)c_{1}</math>
 
Podemos ver que o segundo autovalor tem uma relação entre a extinção e a colonização dos guanacos, se tivermos um valor suficientemente alto de colonização então <math display="inline">\left(1-D\right)c_{1}>e_{1}\rightarrow\lambda_{2}<0</math> . E no primeiro autovalor temos uma relação entre a colonização e morte de predadores, se a morte de predadores for suficientemente alta <math display="inline">e_{y}>\left(1-D-\frac{e_{1}}{c_{1}}\right)c_{y}</math> então <math display="inline">\lambda_{1}<0</math> . Com as duas condições satisfeitas,  temos um ponto de equilíbrio onde os predadores são extintos, e só restam presas.
 
E por fim, temos dois pontos que suportam a coexistência entre duas espécies. Para <math display="inline">\left(x_{3},y_{3}\right)</math>:
 
<math display="block">\frac{du}{dt}=\begin{array}{c}
\left[\left(\left(1-D\right)c_{1}-e_{1}-2c_{1}x_{3}-\mu y_{3}\right)u+\left(-\mu x_{3}\right)v+\left(-\mu x_{3}y_{3}-c_{1}x_{3}^{2}-e_{1}x_{3}+\left(1-D\right)c_{1}x_{3}\right)\right]+\\
\left[-\mu uv-c_{1}u^{2}\right]
\end{array}</math><math display="block">\frac{dv}{dt}=\begin{array}{c}
\left[\left(c_{y}y_{3}-y_{3}^{2}c_{y}\right)u+\left(c_{y}x_{3}-e_{y}-2c_{y}x_{3}y_{3}\right)v+\left(c_{y}x_{3}y_{3}-y_{3}^{2}c_{y}x_{3}-e_{y}y_{3}\right)\right]+\\
\left[-2c_{y}uvy_{3}-c_{y}v^{2}x_{3}-c_{y}uv^{2}+c_{y}uv\right]
\end{array}</math>Calculando os limites para verificar o comportamento semi linear:<math display="block">\lim_{\left(x,y\right)\rightarrow\left(0,0\right)}\frac{\text{parte não-linear}}{\text{parte linear}}</math>
 
Para <math display="inline">\dot{u}</math>:
 
<math display="block">\lim_{\left(u,v\right)\rightarrow\left(0,0\right)}\frac{-\mu uv-c_{1}u^{2}}{\left(\left(1-D\right)c_{1}-e_{1}-2c_{1}x_{3}-\mu y_{3}\right)u+\left(-\mu x_{3}\right)v}</math><math display="block">\lim_{r\rightarrow0}\frac{r^{2}\left(-\mu\cos\theta\sin\theta-c_{1}\cos^{2}\theta\right)}{\left(\left(1-D\right)c_{1}-e_{1}-2c_{1}x_{3}-\mu y_{3}\right)r\cos\theta+\left(-\mu x_{3}\right)r\sin\theta}</math><math display="block">\lim_{r\rightarrow0}\frac{r\left(-\mu\cos\theta\sin\theta-c_{1}\cos^{2}\theta\right)}{\left(\left(1-D\right)c_{1}-e_{1}-2c_{1}x_{3}-\mu y_{3}\right)\cos\theta+\left(-\mu x_{3}\right)\sin\theta}=0</math>
 
E para <math display="inline">\dot{v}</math>:
 
<math display="block">\lim_{\left(u,v\right)\rightarrow\left(0,0\right)}\frac{-2c_{y}uvy_{3}-c_{y}v^{2}x_{3}-c_{y}uv^{2}+c_{y}uv}{\left(c_{y}y_{3}-y_{3}^{2}c_{y}\right)u+\left(c_{y}x_{3}-e_{y}-2c_{y}x_{3}y_{3}\right)v}</math><math display="block">\lim_{r\rightarrow0}\frac{r^{2}\left(-2c_{y}\cos\theta\sin\theta y_{3}-c_{y}\sin^{2}\theta x_{3}-c_{y}\cos\theta\sin^{2}\theta+c_{y}\cos\theta\sin\theta\right)}{\left(c_{y}y_{3}-y_{3}^{2}c_{y}\right)r\cos\theta+\left(c_{y}x_{3}-e_{y}-2c_{y}x_{3}y_{3}\right)r\sin\theta}</math><math display="block">\lim_{r\rightarrow0}\frac{r\left(-2c_{y}\cos\theta\sin\theta y_{3}-c_{y}\sin^{2}\theta x_{3}-c_{y}\cos\theta\sin^{2}\theta+c_{y}\cos\theta\sin\theta\right)}{\left(c_{y}y_{3}-y_{3}^{2}c_{y}\right)\cos\theta+\left(c_{y}x_{3}-e_{y}-2c_{y}x_{3}y_{3}\right)\sin\theta}=0</math>
 
Montando a matriz linearizada então:
 
<math display="block">\left(\begin{array}{c}
\dot{u}\\
\dot{v}
\end{array}\right)=\left(\begin{array}{cc}
\left(1-D\right)c_{1}-e_{1}-2c_{1}x_{3}-\mu y_{3} & -\mu x_{3}\\
c_{y}y_{3}-y_{3}^{2}c_{y} & c_{y}x_{3}-e_{y}-2c_{y}x_{3}y_{3}
\end{array}\right)\left(\begin{array}{c}
u\\
v
\end{array}\right)</math>Infelizmente o autovalor é extremamente complexo se vamos considerar todos os parâmetros como variáveis. A partir de agora, vamos trabalhar atribuindo valores numéricos aos parâmetros, exceto o fato de destruição do habitat, que é um parâmetro que estamos interessado em analisar como afeta a dinâmica.
 
==== Caso particular 1: Utilizando os parâmetros do artigo ====
Utilizando o seguinte conjunto de parâmetros
 
* <math display="inline">c_{1}=0.1</math>
* <math display="inline">e_{1}=0.025</math>
* <math display="inline">e_{y}=0.015</math>
* <math display="inline">\mu=0.3</math>
* <math display="inline">c_{y}=0.015</math>
Temos os seguintes pontos de equilíbrio com os respectivos autovalores:
 
* <math display="inline">\left(x_{1},y_{1}\right)=\left(0,0\right)</math>
** <math display="inline">\lambda_{1,k}\in\left\{ -\frac{3}{200},\frac{\left(0.75-D\right)}{10}\right\}</math>
* <math display="inline">\left(x_{2},y_{2}\right)=\left(0.75-D,0\right)</math>
** <math display="inline">\lambda_{2,k}\in\left\{ -\frac{12D+3}{800},\frac{\left(D-0.75\right)}{10}\right\}</math>
* <math display="inline">\left(x_{3},y_{3}\right)=\left(\frac{-4D-9-\sqrt{16D^{2}+72D+273}}{8},\frac{-4D+\left(\sqrt{16D^{2}+72D+273}+15\right)}{24}\right)</math>
** <math display="inline">\lambda_{3,1}=\frac{23\sqrt{16D^{2}+72D+273}+92D+231-\sqrt{2}\sqrt{\left(1156D+1233\right)\sqrt{16D^{2}+72D+273}+4624D^{2}+23016D+33369}}{3200}</math>
** <math display="inline">\lambda_{3,2}=\frac{23\sqrt{16D^{2}+72D+273}+92D+231+\sqrt{2}\sqrt{\left(1156D+1233\right)\sqrt{16D^{2}+72D+273}+4624D^{2}+23016D+33369}}{3200}</math>
* <math display="inline">\left(x_{4},y_{4}\right)=\left(\frac{-4D-9+\sqrt{16D^{2}+72D+273}9}{8},\frac{-4D-\left(\sqrt{16D^{2}+72D+273}+15\right)}{24}\right)</math>
** <math display="inline">\lambda_{4,1}=\frac{-23\sqrt{16D^{2}+72D+273}+92D+231-\sqrt{2}\sqrt{\left(-1156D-1233\right)\sqrt{16D^{2}+72D+273}+4624D^{2}+23016D+33369}}{3200}</math>
** <math display="inline">\lambda_{4,2}=\frac{-23\sqrt{16D^{2}+72D+273}+92D+231+\sqrt{2}\sqrt{\left(-1156D-1233\right)\sqrt{16D^{2}+72D+273}+4624D^{2}+23016D+33369}}{3200}</math>
 
Para o primeiro ponto ele está fixo <math display="inline">\left(0,0\right)</math> e podemos perceber que um autovalor é negativo e o outro depende se <math display="inline">D</math> é maior ou menor que <math display="inline">0.75</math>. Para fatores de destruição maior, ambos autovalores são negativos e temos um ponto de equilíbrio estável com a extinção de ambas  espécies, para fatores de destruição menor, temo um ponto instável.
 
Para o segundo ponto de equilíbrio, temos novamente uma dependência da relação entre <math display="inline">D</math> e <math display="inline">0.75</math>. O próprio pronto vai se aproximando a <math display="inline">\left(0,0\right)</math> conforme <math display="inline">D</math> cresce. Além disso para valores em que <math display="inline">D<0.75</math> temos um ponto de equilíbrio estável e para acima deste valores se torna instável. Lembrando que para <math display="inline">D>0.75</math> o ponto já passa a ter uma população negativa de pumas, o que não faz sentido ecologicamente, enquanto o primeiro ponto de equilíbrio <math display="inline">\left(0,0\right)</math> se torna estável. Então “recuperamos”o resultado do primeiro ponto de equilíbrio.
 
Os próximos como pode-se imaginar são mais complicados de se analisar, por isso vamos utilizar gráficos para trabalhar.
 
<div class="list">[[Ficheiro:3-4P.png|esquerda|280x200px|miniaturadaimagem|Os gráfico superiores fazem referência ao terceiro ponto de equilíbrio, e os inferiores ao quarto ponto. À esquerda a localização do ponto de equilíbrio, e à direita os autovalores.]]<span> </span>
</div>
Os dois pontos de equilíbrio tem uma população negativa, e um auto valor positivo, portanto são pontos que não fazem sentido ecologicamente e matematicamente são pontos de equilíbrio instável. Desta forma, o sistema sempre evolui para a extinção de ambas as espécies, como podemos ver no espaço de fazes. Seja calculando a evolução do sistema para diferentes condições inciais, ou fazendo o rascunho do diagrama de fase.
<div class="list">[[Ficheiro:diagramas_caso1.png|centro|miniaturadaimagem|643x269px|A esquerda temos a evolução do sistema no espaço de estados para diferentes condições iniciais, calculado via Python, e à direta o diagrama de fase plotado via [https://www.wolframalpha.com/widgets/view.jsp?id=9298fea31cf266903b3df7174b95ddd7 WolframAlpha].]]<span> </span>
</div>
 
==== Caso particular 2: Fatores de colonização aumentados ====
Análogo ao problema do modelo de Levins, este modelo ainda carrega a deficiência de não considerar a configuração espacial. Vamos analisar outro cenário que damos multiplicamos por 4 as taxas de colonização, com o objetivo de analisarmos um cenário em que há coexistência entre presas epredadores.
 
Repetindo os cálculos temos que o primeiro ponto de equilíbrio é <math display="inline">\left(0,0\right)</math>, com os autovalores <math display="inline">\left(-\frac{3}{200},\frac{16}{40}\left(0.9375-D\right)\right)</math>, então agora para <math display="inline">D<0.9375</math>, este é um ponto de equilíbrio instável. O segundo ponto de equilíbrio é <math display="inline">\left(D-0.9375,0\right)</math> com os autovalores <math display="inline">\left(\frac{48}{800}\left(0.6875-D\right),\frac{16}{40}\left(D-0.9375\right)\right)</math>, então para <math display="inline">D<0.6875</math> temos temos o par de sinai dos autovalores <math display="inline">\left(+,-\right)</math> , para valores maiores que <math display="inline">D>0.9375</math> <math display="inline">\left(-,+\right)</math> e para valores intermediários <math display="inline">0.6875<D<0.9375</math> temos <math display="inline">\left(-,-\right)</math>. O que isso significa?
 
 
Bom, para valores maiores que <math display="inline">0.9375</math> então temos um ponto de equilíbrio instável,e nessa situação o ponto estável é em <math display="inline">\left(0,0\right)</math> ou seja, esperamos o a extinção de de ambas as espécies. Para valores entre <math display="inline">0.6875<D<0.9375</math> temos um ponto de equilíbrio estável, nessa situação, apenas a presa sobrevive. E para valores menores, temos um ponto de equilíbrio instável. E esperamos que os próximos dois pontos de equilíbrio nos ajude a responder o que acontece. Semelhante à discussão anterior, os próximos dois pontos são melhores analisados através de gráficos.
 
<div class="center">[[Ficheiro:3P.png|miniaturadaimagem|280x200px|Análise do terceiro ponto de equilíbrio para o caso particular 2. À esquerda a psição do ponto e a direita seu auto-valor.]]</div>
Para o terceiro ponto de equilíbrio,  sempre temos população de presas negativas, o que não fez sentido ecologicamente. Esse resultadoé coerente  matematicamente pois  a parte real dos autovalores são sempre positivos, então é um ponto instável. E por fim, o mais interessante é o quarto ponto de equilíbrio. Para este ponto, <math display="inline">D=0.6875</math> é exatamente o valor no qual deixa de ser um ponto que ecologicamente faz sentido e também deixa de ser um ponto de equilíbrio estável, e para valores menores que este, tempos um ponto de equilíbrio estável. Então resumidamente, temos
[[Ficheiro:4.png|direita|miniaturadaimagem|280x200px|Análise do quarto ponto de equilíbrio para o caso particular 2. À esquerda a psição do ponto e a direita seu auto-valor.]]
 
*<math display="inline">D<0.6875</math>: Coexistência de ambas espécies;
*<math display="inline">0.6875<D<0.9375</math>: Sobrevivência da presa;
*<math display="inline">0.9375<D</math>: Extinção de todas as espécies.
 
Plotando então a evolução e diagramas de fases para 3 valores de <math display="inline">D</math> (<math display="inline">0.5,0.7,0.96</math>), temos:
[[Ficheiro:diagramas.png|centro|miniaturadaimagem|367x300px|Na esquerda superior temos a evolução do sistema no espaço de fase plotado via Python para os diferentes valores de <math>D</math>. As outras três imagens são os diagramas de fases pltoas via WolframAlpha também para os três diferentes valores de  <math>D</math>.]]
 
=== Códigos ===
 
Script para o [https://wxmaxima-developers.github.io/wxmaxima/ WxMaxima] utilizado para efetuar os cálculos simbólicos necessários para analisar os pontos de equilíbrio considerando todos parâmetros como variáveis:
Cálculos simbólico do modelo de Levins aprimorado para 2 espécies considerando todos parâmetros como variáveis
Jhordan Silveira de Borba
sbjhordan@gmail.com
  --> ratprint: false$
Declaração do sistema de equações:
  --> dx:(x*(c_1*(1-D)-mu*y)-c_1*x^2-e_1*x);
dy:(y*(-e_y+c_y*x)-c_y*x*y^2);
Soluções e verificação das soluções
  --> sol:(algsys([dx, dy], [x,y]));
N:1;radcan(dx),sol[N][1],sol[N][2]; radcan(dy),sol[N][1],sol[N][2];
Segundo ponto:
- Deslocamento
- Verificação que a solução linear tem o mesmo ponto de equilíbrio e autovalor
    - Aqui se a partir de A2 tivermos x_2=x, então x é a posição do ponto de equilíbrio, e não uma variável.
  --> radcan(dx),x=u+x_2,y=v;radcan(dy),x=u+x_2,y=v;
  --> B:matrix([u],[v]);
A2:matrix([-2*c_1*x_2+(1-D)*c_1-e_1,-mu*x_2],[0,c_y*x_2-e_y]),x_2=x;
C2:matrix([-c_1*x_2^2-e_1*x_2+(1-D)*x_2*c_1],[0]),x_2=x;
radcan(A2.B),u=0,v=0,sol[2][1];
radcan(C2),u=0,v=0,sol[2][1];
radcan(eigenvalues(A2)),sol[2][1];
Terceiro e quarto ponto ponto:
- Deslocamento
- Verificação que a solução linear tem o mesmo ponto de equilíbrio e autovalor
    - Aqui se a partir de A3 tivermos x_3=x e y_3=y, então x e y são a posição do ponto de
equilíbrio, e não uma variável.
  --> radcan(dx),x=u+x_3,y=v+y_3;radcan(dy),x=u+x_3,y=v+y_3;
  --> A3:matrix([(1-D)*c_1-e_1-2*c_1*x_3-mu*y_3,-mu*x_3],[c_y*y_3-y_3^2*c_y,c_y*x_3-e_y-2*c_y*x_3*y_3]),x_3=x,y_3=y;
C3:matrix([-mu*x_3*y_3-c_1*x_3^2-e_1*x_3+(1-D)*c_1*x_3],[c_y*x_3*y_3-y_3^2*c_y*x_3-e_y*y_3]),x_3=x,y_3=y;
radcan(A3.B),u=0,v=0,sol[3][1],sol[3][2];
radcan(C3),u=0,v=0,sol[3][1],sol[3][2];
radcan(C3),u=0,v=0,sol[4][1],sol[4][2];
radcan(eigenvalues(A3)),sol[3][1],sol[3][2];
 
Script para calcularmos os pontos de equilíbrio e autovalores do sistema tendo a fração de destruição do sistema como única variável:
Análise utilizando como variável apenas D
Jhordan Silveira de Borba
sbjhordan@gmail.com
- Repetimos todos os cálculos substituindo as variáveis pelos valores
- Plotamos os gráficos
- Obs.: Para obtermos a parte real usamos realpart(%) e a imaginária
imagpart(%)
--> c_1:0.1;e_1:0.025;e_y:0.015;mu:0.3;c_y:0.015;
dx:(x*(c_1*(1-D)-mu*y)-c_1*x^2-e_1*x);
dy:(y*(-e_y+c_y*x)-c_y*x*y^2);
sol:(algsys([dx, dy], [x,y]));
B:matrix([u],[v]);
A1:matrix([c_1*(1-D)-e_1,0],[0,-e_y]);
A2:matrix([-2*c_1*x_2+(1-D)*c_1-e_1,-mu*x_2],[0,c_y*x_2-e_y]),x_2=x;
A3:matrix([(1-D)*c_1-e_1-2*c_1*x_3-mu*y_3,-mu*x_3],[c_y*y_3-y_3^2*c_y,c_y*x_3-e_y-2*c_y*x_3*y_3]),x_3=x,y_3=y;
l1:radcan(eigenvalues(A1));
l2:radcan(eigenvalues(A2)),sol[2][1];
l3:radcan(eigenvalues(A3)),sol[3][1],sol[3][2];
l4:radcan(eigenvalues(A3)),sol[4][1],sol[4][2];
  --> X:x,sol[3][1];Y:y,sol[3][2];
wxplot2d([X,Y], [D,0,1],[legend,"x","y"]);
wxplot2d([l3[1][1],l3[1][2]], [D,0,1]);
X:x,sol[4][1];Y:y,sol[4][2];
wxplot2d([X,Y], [D,0,1],[legend,"x","y"]);
wxplot2d([l4[1][1],l4[1][2]], [D,0,1]);
Script para calcularmos os vetores em torno de cada ponto de equilíbrio, para que possamos fazer um rascunho do plano de fase.
Calculo de  vetores para o rascunho do plano de fase com D=0.1:
Jhordan Silveira de Borba
  --> D:0.1;c_1:0.1;e_1:0.025;e_y:0.015;mu:0.3;c_y:0.015;
dx:(x*(c_1*(1-D)-mu*y)-c_1*x^2-e_1*x);
dy:(y*(-e_y+c_y*x)-c_y*x*y^2);
sol:(algsys([dx, dy], [x,y])),numer;
B:matrix([u],[v]);
A1:matrix([c_1*(1-D)-e_1,0],[0,-e_y]);
A2:matrix([-2*c_1*x_2+(1-D)*c_1-e_1,-mu*x_2],[0,c_y*x_2-e_y]),x_2=x;
C2:matrix([-c_1*x_2^2-e_1*x_2+(1-D)*x_2*c_1],[0]),x_2=x;
A3:matrix([(1-D)*c_1-e_1-2*c_1*x_3-mu*y_3,-mu*x_3],[c_y*y_3-y_3^2*c_y,c_y*x_3-e_y-2*c_y*x_3*y_3]),x_3=x,y_3=y;
C3:matrix([-mu*x_3*y_3-c_1*x_3^2-e_1*x_3+(1-D)*c_1*x_3],[c_y*x_3*y_3-y_3^2*c_y*x_3-e_y*y_3]),x_3=x,y_3=y;
l1:radcan(eigenvalues(A1)),numer;
l2:radcan(eigenvalues(A2)),sol[2][1],numer;
l3:radcan(eigenvalues(A3)),sol[3][1],sol[3][2],numer;
l4:radcan(eigenvalues(A3)),sol[4][1],sol[4][2],numer;
  --> radcan((A1.B));
radcan((A2.B)),sol[2][1];
radcan((A3.B)),sol[3][1],sol[3][2];
radcan((A3.B)),sol[4][1],sol[4][2];
Script em Python para plotar o rascunho do diagrama de fases, pode ser editado para visualizarmos só a vizinhança de cada ponto.
# -*- coding: utf-8 -*-
# Modelo aprimorado de Levins para 2 espécies: Diagrama de fases
# Jhordan Silveira de Borba
# sbjhordan@gmail.com
import matplotlib.pyplot as plt
import numpy as np
def phase_lot():
    p1=[0.,0.]                                    # Ponto de equilíbrio 1
    p2=[0.65,0.]                                  # Ponto de equilíbrio 2
    p3=[-3.27,1.31]                              # Ponto de equilíbrio 3
    p4=[0.92,-0.09]                              # Ponto de equilíbrio 4
   
    X = np.arange(0,1, 0.1/2)                    # Eixo X
    Y = np.arange(0,1, 0.1/2)                    # Eixo Y
    U,V=np.meshgrid(X,Y)
       
    c=0
    for x in X:
        l=0
        for y in Y:
            #Distâncias
            d=[]
            d.append(np.sqrt((x-p1[0])*(x-p1[0])+(y-p1[1])*(y-p1[1])))
            d.append(np.sqrt((x-p2[0])*(x-p2[0])+(y-p2[1])*(y-p2[1])))
            d.append(np.sqrt((x-p3[0])*(x-p3[0])+(y-p3[1])*(y-p3[1])))
            d.append(np.sqrt((x-p4[0])*(x-p4[0])+(y-p4[1])*(y-p4[1])))
            for i in range(len(d)):
                if (min(d)==d[i]):
                    p=i
           
            #Calculamos o vetor baseado no ponto mais próximo:
            if(p==0):
                u=x-p1[0]
                v=y-p1[1]
                a=(13*u)/200
                b=-(3*v)/200             
            elif (p==1):
                u=x-p2[0]
                v=y-p2[1]
                a=-(70*v+23*u+25)/36
                b=-(21*v)/40
            elif(p==2):
                u=x-p3[0]
                v=y-p3[1]                                     
                a=(285*v+80*u-5)/29
                b=(11*v-10*u)/17
            elif(p==3):
                u=x-p4[0]
                v=y-p4[1]                                     
                a=-(47*v+15*u+12)/17
                b=(20*v-21*u-66)/18
            else:
                print("Algo deu errado")
            m=np.sqrt(a*a+b*b) # Módulo do vetor
            if(m==0):
                m=1
    # Normalizamos e salvamos o vetor: [dx/dt,dy/dt]=[a,b]
            U[l,c]=a/m
            V[l,c]=b/m
            l=l+1
        c=c+1
   
    fig, ax = plt.subplots()
    ax.quiver(X, Y, U, V)
    plt.plot(p0[0],p0[1],'ko') #Plotamos o primeiro ponto de equilíbrio
    plt.plot(p1[0],p1[1],'ko') #Plotamos o segundo ponto de equilíbrio
    plt.plot(p3[0],p2[1],'ko') #Plotamos o terceiro ponto de equilíbrio
    plt.plot(p4[0],p3[1],'ko') #Plotamos o quarto ponto de equilíbrio
    plt.show()
Script em Python para analisar a evolução do sistema no espaço de fases, utilizando o método de Euler:
# -*- coding: utf-8 -*-
# Modelo aprimorado de Levins para 2 espécies: evolução do sistema no espaçode estados
# Jhordan Silveira de Borba
# sbjhordan@gmail.com
def phase2():
    #Parâmetros   
    c1=0.1                  # Fator de colonização do guanaco
    cy=0.015                # Fator de colonização do puma
    e1=0.025                # Fator de extinção local do guanaco
    ey=0.015                # Fator de extinção local do puma
    u1=0.3                  # Fator de predação do guanaco
    N=int(200000)            # Duração da simulação
    d=0.01                  # Tamanho do passo
    D=0.5                    # Destruição do habitat
   
    p1=[0.,0.]      # Ponto de equilíbrio 1
    p2=[0.65,0.]    # Ponto de equilíbrio 2
    p3=[-3.26,1.30] # Ponto de equilíbrio 3
    p4=[0.92,-0.09] # Ponto de equilíbrio 4
   
    #Valores iniciais:
    init=[(0.4,0.4),(0.3,0.3),(0.05,0.05)]
    c=0
    for it in init:
        x=[]
        y=[]
        x.append(it[0])
        y.append(it[1])
        for i in range(N-1):
            x.append(x[i]+d*(c1*x[i]*(1-D-x[i])-e1*x[i]-u1*x[i]*y[i]))
            y.append(y[i]+d*(cy*y[i]*(x[i]-y[i]*(x[i]))-ey*y[i]))
        if(c==0):
            plt.plot(x,y,'g-')
            plt.plot(it[0],it[1],'go')
            D=0.7
        elif (c==1):
            plt.plot(x,y,'b-')
            plt.plot(it[0],it[1],'bo')
            D=0.95
        else:
            plt.plot(x,y,'r-')
            plt.plot(it[0],it[1],'ro')
        c=c+1
       
    plt.plot(p1[0],p1[1],'ko')
    plt.plot(p2[0],p2[1],'ko')
    plt.plot(p3[0],p3[1],'ko')
    plt.plot(p4[0],p4[1],'ko')
    plt.xlim(0,0.5)
    plt.ylim(0,0.5)
    plt.show()
 
=== Principal material utilizado ===
 
# [https://arxiv.org/pdf/1409.0024.pdf Mathematical model of livestock and wildlife: Predationand competition under environmental disturbances] (Fabiana Laguna e outros, Ecological Modelling)
 
 
{{Ecologia| [[Modelo de Lotka-Volterra amortecido]] | [[Modelo de Levins aprimorado para 2 espécies II]] }}

Edição atual tal como às 21h53min de 2 de maio de 2021

Anterior: Modelo de Lotka-Volterra amortecido | Índice: Ecologia | Próximo: Modelo de Levins aprimorado para 2 espécies II

Comparação com outros modelos

Para duas espécies, vamos trabalhar com o seguinte sistema de equações:

Onde é a presa e o predador. Podemos comparar com outros modelos. Como um modelo inspirado no modelo de Levins, podemos ver semelhanças bastante fortes desde já. Lembrando do modelo de Levins rapidamente:

O termo na equação da presa é a quantidade de fragmentos disponíveis para a colonização da presa, uma vez que precisamos descontar os fragmentos destruídos, além dos ocupados. Como todos parâmetros são em termos de proporções, temos probabilidade de selecionarmos um fragmento disponível para a colonização da presa. E o termo na equação do predador tem a mesma função, porém agora precisamos lembrar que os fragmentos que o predador pode ocupar são os que já estão ocupados pela presa, descontando os que já possuem predador. Então considerando propriedades de conjuntos, e olhando em termos de probabilidades esta equação, temos a probabilidade de selecionar um fragmento ocupada pela presa, mas precisamos descontar a probabilidade de selecionar um fragmento que também já esteja ocupada por um predador . A maior diferença se dá porque a presa tem um termo a mais correspondendo ao decréscimo da população devido a predação, evento não considerado no modelo de Levins, mas considerado no modelo de Lotka-Volterra.


Por isso a próxima comparação que podemos fazer é com o modelo de Lotka-Volterra amortecido:



Vamos reescrever nosso sistema de equações:


Olhando agora primeiro para a equação da presa, a maior diferença que encontramos, e é que temos o termo de extinção local que não é devido a predação , algo já familiar do modelo de Levins. E já quando olhamos para a equação do predador, temos o termo , podemos pensar também como um termo logístico, impondo um limite na população de predadores, da mesma forma que ocorre com as presas. Porém agora o termo logístico não é uma constante, mas depende também da população de presas.

Análises dos pontos de equilíbrio

Sendo todos o parâmetros como variáveis

Vamos tentar fazer análises parecidas agora. Primeiro buscar os pontos de equilíbrio pro sistema:

Temos então 4 pontos de equilíbrio:

  • Primeiro: extinção de ambas as espécies
  • Segundo: extinção do predador

Os pontos de equilíbrio com ambas as espécies são extremamente mais complicados. Para o terceiro ponto temos:

E o quarto e último ponto:

Estes últimos pontos (com o equilíbrio entre as espécies) são complicados demais para nos referirmos explicitamente daqui em diante quando necessário utilizamos alguma notação como para denotar o ponto de equilíbrio, sem escrevê-lo por inteiro. Vamos linearizar e analisar a estabilidade em torno dos pontos de equilíbrio. Para o ponto é muito simples:

Separando os termos lineares dos não lineares:

Calculando então os limites:

Desprezando os termos não lineares:

Logo um dos autovalores é real, sempre negativo e o segundo é . Basicamente se o termo de colonização do guanaco for maior que o da extinção , temos dois autovalores de sinais opostos, ou seja, um ponto de instabilidade. Aumentando a taxa de extinção do guanaco (ou a destruição do habitat) passamos a ter , e a extinção de ambas as espécies. Aqui podemos notar uma diferença dos modelos de Lotka-Volterra: a extinção de ambas as espécies pode ser um ponto de equilíbrio estável. Ou seja nos modelos anteriores o sistema nunca iria para o equilíbrio com extinção total, agora há essa possibilidade dependendo do nível de destruição do sistema.

Para o segundo ponto de equilíbrio precisamos fazer a translação: e .

Podemos ver facilmente que é semi linear para uma vez que ao fazermos os:

Fazendo a substituição e :

É interessante notar que a parte constante vai zerar. E para :

Analisando o comportamento linear então na proximidade do ponto de equilíbrio:

E os autovalores são então:

Podemos ver que o segundo autovalor tem uma relação entre a extinção e a colonização dos guanacos, se tivermos um valor suficientemente alto de colonização então . E no primeiro autovalor temos uma relação entre a colonização e morte de predadores, se a morte de predadores for suficientemente alta então . Com as duas condições satisfeitas, temos um ponto de equilíbrio onde os predadores são extintos, e só restam presas.

E por fim, temos dois pontos que suportam a coexistência entre duas espécies. Para :

Calculando os limites para verificar o comportamento semi linear:

Para :

E para :

Montando a matriz linearizada então:

Infelizmente o autovalor é extremamente complexo se vamos considerar todos os parâmetros como variáveis. A partir de agora, vamos trabalhar atribuindo valores numéricos aos parâmetros, exceto o fato de destruição do habitat, que é um parâmetro que estamos interessado em analisar como afeta a dinâmica.

Caso particular 1: Utilizando os parâmetros do artigo

Utilizando o seguinte conjunto de parâmetros

Temos os seguintes pontos de equilíbrio com os respectivos autovalores:

Para o primeiro ponto ele está fixo e podemos perceber que um autovalor é negativo e o outro depende se é maior ou menor que . Para fatores de destruição maior, ambos autovalores são negativos e temos um ponto de equilíbrio estável com a extinção de ambas espécies, para fatores de destruição menor, temo um ponto instável.

Para o segundo ponto de equilíbrio, temos novamente uma dependência da relação entre e . O próprio pronto vai se aproximando a conforme cresce. Além disso para valores em que temos um ponto de equilíbrio estável e para acima deste valores se torna instável. Lembrando que para o ponto já passa a ter uma população negativa de pumas, o que não faz sentido ecologicamente, enquanto o primeiro ponto de equilíbrio se torna estável. Então “recuperamos”o resultado do primeiro ponto de equilíbrio.

Os próximos como pode-se imaginar são mais complicados de se analisar, por isso vamos utilizar gráficos para trabalhar.

Os gráfico superiores fazem referência ao terceiro ponto de equilíbrio, e os inferiores ao quarto ponto. À esquerda a localização do ponto de equilíbrio, e à direita os autovalores.

Os dois pontos de equilíbrio tem uma população negativa, e um auto valor positivo, portanto são pontos que não fazem sentido ecologicamente e matematicamente são pontos de equilíbrio instável. Desta forma, o sistema sempre evolui para a extinção de ambas as espécies, como podemos ver no espaço de fazes. Seja calculando a evolução do sistema para diferentes condições inciais, ou fazendo o rascunho do diagrama de fase.

A esquerda temos a evolução do sistema no espaço de estados para diferentes condições iniciais, calculado via Python, e à direta o diagrama de fase plotado via WolframAlpha.

Caso particular 2: Fatores de colonização aumentados

Análogo ao problema do modelo de Levins, este modelo ainda carrega a deficiência de não considerar a configuração espacial. Vamos analisar outro cenário que damos multiplicamos por 4 as taxas de colonização, com o objetivo de analisarmos um cenário em que há coexistência entre presas epredadores.

Repetindo os cálculos temos que o primeiro ponto de equilíbrio é , com os autovalores , então agora para , este é um ponto de equilíbrio instável. O segundo ponto de equilíbrio é com os autovalores , então para temos temos o par de sinai dos autovalores , para valores maiores que e para valores intermediários temos . O que isso significa?


Bom, para valores maiores que então temos um ponto de equilíbrio instável,e nessa situação o ponto estável é em ou seja, esperamos o a extinção de de ambas as espécies. Para valores entre temos um ponto de equilíbrio estável, nessa situação, apenas a presa sobrevive. E para valores menores, temos um ponto de equilíbrio instável. E esperamos que os próximos dois pontos de equilíbrio nos ajude a responder o que acontece. Semelhante à discussão anterior, os próximos dois pontos são melhores analisados através de gráficos.

Análise do terceiro ponto de equilíbrio para o caso particular 2. À esquerda a psição do ponto e a direita seu auto-valor.

Para o terceiro ponto de equilíbrio, sempre temos população de presas negativas, o que não fez sentido ecologicamente. Esse resultadoé coerente matematicamente pois a parte real dos autovalores são sempre positivos, então é um ponto instável. E por fim, o mais interessante é o quarto ponto de equilíbrio. Para este ponto, é exatamente o valor no qual deixa de ser um ponto que ecologicamente faz sentido e também deixa de ser um ponto de equilíbrio estável, e para valores menores que este, tempos um ponto de equilíbrio estável. Então resumidamente, temos

Análise do quarto ponto de equilíbrio para o caso particular 2. À esquerda a psição do ponto e a direita seu auto-valor.
  • : Coexistência de ambas espécies;
  • : Sobrevivência da presa;
  • : Extinção de todas as espécies.

Plotando então a evolução e diagramas de fases para 3 valores de (), temos:

Na esquerda superior temos a evolução do sistema no espaço de fase plotado via Python para os diferentes valores de . As outras três imagens são os diagramas de fases pltoas via WolframAlpha também para os três diferentes valores de .

Códigos

Script para o WxMaxima utilizado para efetuar os cálculos simbólicos necessários para analisar os pontos de equilíbrio considerando todos parâmetros como variáveis:

Cálculos simbólico do modelo de Levins aprimorado para 2 espécies considerando todos parâmetros como variáveis
Jhordan Silveira de Borba
sbjhordan@gmail.com
 -->	ratprint: false$
Declaração do sistema de equações:
 -->	dx:(x*(c_1*(1-D)-mu*y)-c_1*x^2-e_1*x);
	dy:(y*(-e_y+c_y*x)-c_y*x*y^2);
Soluções e verificação das soluções
 -->	sol:(algsys([dx, dy], [x,y]));
	N:1;radcan(dx),sol[N][1],sol[N][2]; radcan(dy),sol[N][1],sol[N][2];
Segundo ponto:
	- Deslocamento
	- Verificação que a solução linear tem o mesmo ponto de equilíbrio e autovalor
	    - Aqui se a partir de A2 tivermos x_2=x, então x é a posição do ponto de equilíbrio, e não uma variável.
 -->	radcan(dx),x=u+x_2,y=v;radcan(dy),x=u+x_2,y=v;
 -->	B:matrix([u],[v]);
	A2:matrix([-2*c_1*x_2+(1-D)*c_1-e_1,-mu*x_2],[0,c_y*x_2-e_y]),x_2=x;
	C2:matrix([-c_1*x_2^2-e_1*x_2+(1-D)*x_2*c_1],[0]),x_2=x;
	radcan(A2.B),u=0,v=0,sol[2][1];
	radcan(C2),u=0,v=0,sol[2][1];
	radcan(eigenvalues(A2)),sol[2][1];
Terceiro e quarto ponto ponto:
	- Deslocamento
	- Verificação que a solução linear tem o mesmo ponto de equilíbrio e autovalor
	    - Aqui se a partir de A3 tivermos x_3=x e y_3=y, então x e y são a posição do ponto de
	equilíbrio, e não uma variável.
 -->	radcan(dx),x=u+x_3,y=v+y_3;radcan(dy),x=u+x_3,y=v+y_3;
 -->	A3:matrix([(1-D)*c_1-e_1-2*c_1*x_3-mu*y_3,-mu*x_3],[c_y*y_3-y_3^2*c_y,c_y*x_3-e_y-2*c_y*x_3*y_3]),x_3=x,y_3=y;
	C3:matrix([-mu*x_3*y_3-c_1*x_3^2-e_1*x_3+(1-D)*c_1*x_3],[c_y*x_3*y_3-y_3^2*c_y*x_3-e_y*y_3]),x_3=x,y_3=y;
	radcan(A3.B),u=0,v=0,sol[3][1],sol[3][2];
	radcan(C3),u=0,v=0,sol[3][1],sol[3][2];
	radcan(C3),u=0,v=0,sol[4][1],sol[4][2];
	radcan(eigenvalues(A3)),sol[3][1],sol[3][2];

Script para calcularmos os pontos de equilíbrio e autovalores do sistema tendo a fração de destruição do sistema como única variável:

Análise utilizando como variável apenas D
Jhordan Silveira de Borba 
sbjhordan@gmail.com

- Repetimos todos os cálculos substituindo as variáveis pelos valores
- Plotamos os gráficos
- Obs.: Para obtermos a parte real usamos realpart(%) e a imaginária
imagpart(%) 
-->	c_1:0.1;e_1:0.025;e_y:0.015;mu:0.3;c_y:0.015;
	dx:(x*(c_1*(1-D)-mu*y)-c_1*x^2-e_1*x);
	dy:(y*(-e_y+c_y*x)-c_y*x*y^2);
	sol:(algsys([dx, dy], [x,y]));
	B:matrix([u],[v]);
	A1:matrix([c_1*(1-D)-e_1,0],[0,-e_y]);
	A2:matrix([-2*c_1*x_2+(1-D)*c_1-e_1,-mu*x_2],[0,c_y*x_2-e_y]),x_2=x;
	A3:matrix([(1-D)*c_1-e_1-2*c_1*x_3-mu*y_3,-mu*x_3],[c_y*y_3-y_3^2*c_y,c_y*x_3-e_y-2*c_y*x_3*y_3]),x_3=x,y_3=y;
	l1:radcan(eigenvalues(A1));
	l2:radcan(eigenvalues(A2)),sol[2][1];
	l3:radcan(eigenvalues(A3)),sol[3][1],sol[3][2];
	l4:radcan(eigenvalues(A3)),sol[4][1],sol[4][2];
 -->	X:x,sol[3][1];Y:y,sol[3][2];
	wxplot2d([X,Y], [D,0,1],[legend,"x","y"]);
	wxplot2d([l3[1][1],l3[1][2]], [D,0,1]);
	X:x,sol[4][1];Y:y,sol[4][2]; 
	wxplot2d([X,Y], [D,0,1],[legend,"x","y"]); 
	wxplot2d([l4[1][1],l4[1][2]], [D,0,1]);

Script para calcularmos os vetores em torno de cada ponto de equilíbrio, para que possamos fazer um rascunho do plano de fase.

Calculo de  vetores para o rascunho do plano de fase com D=0.1:
Jhordan Silveira de Borba
 -->	D:0.1;c_1:0.1;e_1:0.025;e_y:0.015;mu:0.3;c_y:0.015;
	dx:(x*(c_1*(1-D)-mu*y)-c_1*x^2-e_1*x);
	dy:(y*(-e_y+c_y*x)-c_y*x*y^2);
	sol:(algsys([dx, dy], [x,y])),numer;
	B:matrix([u],[v]);
	A1:matrix([c_1*(1-D)-e_1,0],[0,-e_y]);
	A2:matrix([-2*c_1*x_2+(1-D)*c_1-e_1,-mu*x_2],[0,c_y*x_2-e_y]),x_2=x;
	C2:matrix([-c_1*x_2^2-e_1*x_2+(1-D)*x_2*c_1],[0]),x_2=x;
	A3:matrix([(1-D)*c_1-e_1-2*c_1*x_3-mu*y_3,-mu*x_3],[c_y*y_3-y_3^2*c_y,c_y*x_3-e_y-2*c_y*x_3*y_3]),x_3=x,y_3=y;
	C3:matrix([-mu*x_3*y_3-c_1*x_3^2-e_1*x_3+(1-D)*c_1*x_3],[c_y*x_3*y_3-y_3^2*c_y*x_3-e_y*y_3]),x_3=x,y_3=y;
	l1:radcan(eigenvalues(A1)),numer;
	l2:radcan(eigenvalues(A2)),sol[2][1],numer;
	l3:radcan(eigenvalues(A3)),sol[3][1],sol[3][2],numer;
	l4:radcan(eigenvalues(A3)),sol[4][1],sol[4][2],numer;
 -->	radcan((A1.B));
	radcan((A2.B)),sol[2][1];
	radcan((A3.B)),sol[3][1],sol[3][2];
	radcan((A3.B)),sol[4][1],sol[4][2];

Script em Python para plotar o rascunho do diagrama de fases, pode ser editado para visualizarmos só a vizinhança de cada ponto.

# -*- coding: utf-8 -*-
# Modelo aprimorado de Levins para 2 espécies: Diagrama de fases
# Jhordan Silveira de Borba
# sbjhordan@gmail.com

import matplotlib.pyplot as plt
import numpy as np


def phase_lot():
    p1=[0.,0.]                                    # Ponto de equilíbrio 1
    p2=[0.65,0.]                                  # Ponto de equilíbrio 2 
    p3=[-3.27,1.31]                               # Ponto de equilíbrio 3 
    p4=[0.92,-0.09]                               # Ponto de equilíbrio 4
    
    X = np.arange(0,1, 0.1/2)                     # Eixo X
    Y = np.arange(0,1, 0.1/2)                     # Eixo Y
    U,V=np.meshgrid(X,Y)
        
    c=0
    for x in X:
        l=0
        for y in Y:
            #Distâncias
            d=[]
            d.append(np.sqrt((x-p1[0])*(x-p1[0])+(y-p1[1])*(y-p1[1])))
            d.append(np.sqrt((x-p2[0])*(x-p2[0])+(y-p2[1])*(y-p2[1])))
            d.append(np.sqrt((x-p3[0])*(x-p3[0])+(y-p3[1])*(y-p3[1])))
            d.append(np.sqrt((x-p4[0])*(x-p4[0])+(y-p4[1])*(y-p4[1])))
            for i in range(len(d)):
                if (min(d)==d[i]):
                    p=i
            
            #Calculamos o vetor baseado no ponto mais próximo:
            if(p==0):
                u=x-p1[0]
                v=y-p1[1]
                a=(13*u)/200
                b=-(3*v)/200               
            elif (p==1):
                u=x-p2[0]
                v=y-p2[1]
                a=-(70*v+23*u+25)/36
                b=-(21*v)/40 
            elif(p==2):
                u=x-p3[0]
                v=y-p3[1]                                       
                a=(285*v+80*u-5)/29
                b=(11*v-10*u)/17
            elif(p==3):
                u=x-p4[0]
                v=y-p4[1]                                       
                a=-(47*v+15*u+12)/17
                b=(20*v-21*u-66)/18
            else:
                print("Algo deu errado")

            m=np.sqrt(a*a+b*b) 		 # Módulo do vetor
            if(m==0):
                m=1
	    # Normalizamos e salvamos o vetor: [dx/dt,dy/dt]=[a,b]
            U[l,c]=a/m
            V[l,c]=b/m
            l=l+1
        c=c+1
    
    fig, ax = plt.subplots()
    ax.quiver(X, Y, U, V)
    plt.plot(p0[0],p0[1],'ko') #Plotamos o primeiro ponto de equilíbrio
    plt.plot(p1[0],p1[1],'ko') #Plotamos o segundo ponto de equilíbrio
    plt.plot(p3[0],p2[1],'ko') #Plotamos o terceiro ponto de equilíbrio
    plt.plot(p4[0],p3[1],'ko') #Plotamos o quarto ponto de equilíbrio
    plt.show()

Script em Python para analisar a evolução do sistema no espaço de fases, utilizando o método de Euler:

# -*- coding: utf-8 -*-
# Modelo aprimorado de Levins para 2 espécies: evolução do sistema no espaçode estados
# Jhordan Silveira de Borba
# sbjhordan@gmail.com

def phase2():

    #Parâmetros    
    c1=0.1                   # Fator de colonização do guanaco 
    cy=0.015                 # Fator de colonização do puma
    e1=0.025                 # Fator de extinção local do guanaco
    ey=0.015                 # Fator de extinção local do puma
    u1=0.3                   # Fator de predação do guanaco
    N=int(200000)            # Duração da simulação
    d=0.01                   # Tamanho do passo
    D=0.5                    # Destruição do habitat
    
    p1=[0.,0.]      # Ponto de equilíbrio 1
    p2=[0.65,0.]    # Ponto de equilíbrio 2
    p3=[-3.26,1.30] # Ponto de equilíbrio 3
    p4=[0.92,-0.09] # Ponto de equilíbrio 4
    
    #Valores iniciais:
    init=[(0.4,0.4),(0.3,0.3),(0.05,0.05)]
    c=0
    for it in init:
        x=[]
        y=[]
        x.append(it[0])
        y.append(it[1])
        for i in range(N-1):
            x.append(x[i]+d*(c1*x[i]*(1-D-x[i])-e1*x[i]-u1*x[i]*y[i]))
            y.append(y[i]+d*(cy*y[i]*(x[i]-y[i]*(x[i]))-ey*y[i]))
        if(c==0):
            plt.plot(x,y,'g-')
            plt.plot(it[0],it[1],'go')
            D=0.7
        elif (c==1):
            plt.plot(x,y,'b-')
            plt.plot(it[0],it[1],'bo')
            D=0.95
        else:
            plt.plot(x,y,'r-')
            plt.plot(it[0],it[1],'ro')
        c=c+1
        
    plt.plot(p1[0],p1[1],'ko')
    plt.plot(p2[0],p2[1],'ko')
    plt.plot(p3[0],p3[1],'ko')
    plt.plot(p4[0],p4[1],'ko')
    plt.xlim(0,0.5)
    plt.ylim(0,0.5)
    plt.show()

Principal material utilizado

  1. Mathematical model of livestock and wildlife: Predationand competition under environmental disturbances (Fabiana Laguna e outros, Ecological Modelling)


Anterior: Modelo de Lotka-Volterra amortecido | Índice: Ecologia | Próximo: Modelo de Levins aprimorado para 2 espécies II