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:
![{\displaystyle {\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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/cd3347106e5a4f0c4a7b12713f9d73b8c2cc4654)
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:
![{\displaystyle {\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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/65b4804aa8983ee34a9004fa1535e924c76a8b34)
E o quarto e último ponto:
![{\displaystyle {\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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/93cc2caccd95bae7e160f4e52a1b723f4f21e6f2)
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
![{\textstyle \left(x_{j},y_{j}\right)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/520da3fd262227052fa72a64cf10f139ec89e5c0)
para denotar o
![{\textstyle j-{\acute {e}}simo}](https://wikimedia.org/api/rest_v1/media/math/render/svg/4329597bd3a0354666b0ecf544aaf985a775ef2d)
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
![{\textstyle \left(x_{1},y_{1}\right)=\left(0,0\right)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/7373851a055d78205a6feed98aa7fd0fc0e956a0)
é muito simples:
![{\displaystyle {\frac {dx}{dt}}=x\left(c_{1}\left(1-D\right)-\mu y\right)-c_{1}x^{2}-e_{1}x}](https://wikimedia.org/api/rest_v1/media/math/render/svg/51fae24d831f79b3733d888ee5178ab04aad24a5)
![{\displaystyle {\frac {dy}{dt}}=y\left(-e_{y}+c_{y}x\right)-c_{y}xy^{2}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/41816cc591c565fcebee644c820cee33f0b3a385)
Separando os termos lineares dos não lineares:
![{\displaystyle {\frac {dx}{dt}}=\left[x\left(c_{1}\left(1-D\right)-e_{1}\right)\right]-\left[\mu yx+c_{1}x^{2}\right]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/1d5957eaca2437b74ff78180bbd697a60bcfdae2)
![{\displaystyle {\frac {dy}{dt}}=-\left[ye_{y}\right]+\left[c_{y}xy-c_{y}xy^{2}\right]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/cace48f23ca8a797b83e29c197f64978f3d3b95b)
Calculando então os limites:
![{\displaystyle \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}](https://wikimedia.org/api/rest_v1/media/math/render/svg/631f1eb3f38c05f76d35ad9a47c5e0fe0a1f3811)
![{\displaystyle \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}](https://wikimedia.org/api/rest_v1/media/math/render/svg/6007273ea890c25baa7a0d159883afdeeb334fac)
Desprezando os termos não lineares:
![{\displaystyle \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)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/feddd8d999d1569cdc3449bb5341904e9bf11eaa)
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
.
![{\displaystyle {\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]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/9b37392d787aeef0179f0fe1b4d75ffef7d23193)
![{\displaystyle {\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]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/b7280994bb97938fa3a845c0b9dc9bf515d97581)
Podemos ver facilmente que é semi linear para
uma vez que ao fazermos os:
![{\displaystyle \lim _{\left(x,y\right)\rightarrow \left(0,0\right)}{\frac {\text{parte não-linear}}{\text{parte linear}}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/3856b45a529d1c03c79d2935e463a3ad7a864780)
![{\displaystyle \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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/9cbb5961d8156c5116f504e381634d0a527d24a7)
Fazendo a substituição
e
:
![{\displaystyle \lim _{r\rightarrow 0}{\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 }}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/f369fefd8d2429b7ab7f06dc9b0db0301d657e9a)
![{\displaystyle \lim _{r\rightarrow 0}{\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}](https://wikimedia.org/api/rest_v1/media/math/render/svg/db137524126a3ca3d4ea82047681b2bc5d385d9b)
É interessante notar que a parte constante vai zerar. E para
:
![{\displaystyle \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}](https://wikimedia.org/api/rest_v1/media/math/render/svg/e99aa9c08e302639167c227107a26030d2ed9f0c)
Analisando o comportamento linear então na proximidade do ponto de equilíbrio:
![{\displaystyle \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)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/84345616cfe2939e0cc0e448bba2b9f27271b425)
E os autovalores são então:
![{\textstyle \lambda _{1}=\left(1-D-{\frac {e_{1}}{c_{1}}}\right)c_{y}-e_{y}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/5f2e5d3cc47dc4912e334d06bf24c1b76a3a1db4)
![{\textstyle \lambda _{2}=e_{1}-\left(1-D\right)c_{1}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/4ddfbb8ae36dfed6feb5c59890a7fbd2a6fe32a2)
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
:
![{\displaystyle {\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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/e30d0442101efaa32f3701459181c9d7efc820eb)
![{\displaystyle {\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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/3ccd5ba3af39144b732bd38ce14b985844ee9222)
Calculando os limites para verificar o comportamento semi linear:
![{\displaystyle \lim _{\left(x,y\right)\rightarrow \left(0,0\right)}{\frac {\text{parte não-linear}}{\text{parte linear}}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/3856b45a529d1c03c79d2935e463a3ad7a864780)
Para
:
![{\displaystyle \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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/63631ab5c2190deb42eb5f35742be11aec4d9264)
![{\displaystyle \lim _{r\rightarrow 0}{\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 }}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/3900b233724ed3ccaab6cc638439666422032744)
![{\displaystyle \lim _{r\rightarrow 0}{\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}](https://wikimedia.org/api/rest_v1/media/math/render/svg/4b6d02b33e3bc84495bcdb092691df3b683a1393)
E para
:
![{\displaystyle \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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/e1669232d44611545b67738c834d8deefbb940a9)
![{\displaystyle \lim _{r\rightarrow 0}{\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 }}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/efdc6c9921bc5e90cde14839fb05379d68dee7cc)
![{\displaystyle \lim _{r\rightarrow 0}{\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}](https://wikimedia.org/api/rest_v1/media/math/render/svg/998fd2b0dba6f9a0c9909293663588709ad4d1aa)
Montando a matriz linearizada então:
![{\displaystyle \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)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/9c3ef51cee792360da4d3bb312416dc078f05e5c)
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
![{\textstyle c_{1}=0.1}](https://wikimedia.org/api/rest_v1/media/math/render/svg/1d1ca74d44e3f5db153008c67255c774587e1438)
![{\textstyle e_{1}=0.025}](https://wikimedia.org/api/rest_v1/media/math/render/svg/8d2607b4240c07dddcd9537babb69c0e69e859cc)
![{\textstyle e_{y}=0.015}](https://wikimedia.org/api/rest_v1/media/math/render/svg/ec697ce46fd453b94e18019d74e08787e018c718)
![{\textstyle \mu =0.3}](https://wikimedia.org/api/rest_v1/media/math/render/svg/0c87206132353f9a0a728ae3ae7f50f17c65e842)
![{\textstyle c_{y}=0.015}](https://wikimedia.org/api/rest_v1/media/math/render/svg/6ed5b942f3332ca6da84138615a053516e474771)
Temos os seguintes pontos de equilíbrio com os respectivos autovalores:
![{\textstyle \lambda _{1,k}\in \left\{-{\frac {3}{200}},{\frac {\left(0.75-D\right)}{10}}\right\}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/a2ee893a3700e03ee3c947245d1fac7dd863621f)
![{\textstyle \lambda _{2,k}\in \left\{-{\frac {12D+3}{800}},{\frac {\left(D-0.75\right)}{10}}\right\}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/3cdaccd31e560383bb7d6114f93a404c0e33878f)
![{\textstyle \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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/11af170e3a4366765d1b3e344d584e1a1ac5e98a)
![{\textstyle \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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/764e7a6c49b88e6005613f6a0086cb3fdf803084)
![{\textstyle \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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/9dfe4fe962158e107eb0b13c884ce40adb24a97c)
![{\textstyle \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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/383af3b05fd1c8be98a3b2fce7fa74e786949324)
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
![{\displaystyle D}](https://wikimedia.org/api/rest_v1/media/math/render/svg/f34a0c600395e5d4345287e21fb26efd386990e6)
. As outras três imagens são os diagramas de fases pltoas via WolframAlpha também para os três diferentes valores de
![{\displaystyle D}](https://wikimedia.org/api/rest_v1/media/math/render/svg/f34a0c600395e5d4345287e21fb26efd386990e6)
.
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
- Mathematical model of livestock and wildlife: Predationand competition under environmental disturbances (Fabiana Laguna e outros, Ecological Modelling)