Modelo de Lotka-Volterra amortecido

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

Uma versão do modelo de Lotka-Volterra aprimorada inclui um termo de saturação na população de presas, isto é, um termo logístico (que inibe o aumento exponencial) visando representar a finitude dos recursos disponíveis para uma espécie. Este modelo é chamado de modelo de Lotka-Volterra amortecido.

• $\frac{dx}{dt}=x(a-\alpha y)-k{x}^2$
• $\frac{dy}{dt}=y(-c+\gamma x)$

Onde $k$ é uma constante positiva. Agora temos três pontos de equilíbrio, novamente temos $(0,0)$. Mas temos outro ponto quando apenas $y=0$:

$$a-\alpha y-kx=a-kx=0\to x=\frac{a}{k}$$

Então $(\frac{a}{k},0)$ ou seja um ponto de equilíbrio com apenas a sobrevivência da presa. E ainda temos também um ponto de coexistência entre ambas as espécies:

$$a-\alpha y-kx=0$$$$-c+\gamma x=0$$

Isolando então $x$ na segunda equação $x=\frac{c}{\gamma }$ e substituindo na primeira:

$$a-\alpha y-k\frac{c}{\gamma }=0$$Então:$$y=\frac{a}{\alpha }-\frac{kc}{\alpha \gamma }=\frac{\gamma a-kc}{\alpha \gamma }$$

Então nosso outro ponto de equilíbrio é dado por$(\frac{c}{\gamma },\frac{\gamma a-kc}{\alpha \gamma })$. Infelizmente as equações não são separáveis como no caso anterior, portanto vamos analisar se a equação é semi-linear nas proximidades dos pontos de equilíbrio. Para $(0,0)$:

$$\underset{(x,y)\to (0,0)}{\lim }-\frac{\alpha xy-k{x}^2}{xa}=\underset{(x,y)\to (0,0)}{\lim }-\frac{\alpha y-kx}{a}=0$$$$\underset{(x,y)\to (0,0)}{\lim }-\frac{\gamma xy}{cy}=\underset{(x,y)\to (0,0)}{\lim }-\frac{\gamma }{c}x=0$$

Para $(\frac{a}{k},0)$ e fazendo a mudança de variáveis $u=x-\frac{a}{k}$, e $v=y$, temos:

$$\frac{du}{dt}=((u+\frac{a}{k})a-k{(u+\frac{a}{k})}^2)-\left(\alpha (u+\frac{a}{k})v\right)$$$$\frac{du}{dt}=(ua+\frac{a^2}{k}-(u^{2}k+\frac{a^2}{k}+2ua))-(\alpha uv+\frac{\alpha a}{k}v)$$$$\frac{du}{dt}=-(ua+\frac{\alpha a}{k}v)-(\alpha uv+u^{2}k)$$

E:

$$\frac{dv}{dt}=(-vc)+\left(\gamma v(u+\frac{a}{k})\right)=v(\gamma \frac{a}{k}-c)+\left(\gamma vu\right)$$

Então:

$$\underset{(u,v)\to (0,0)}{\lim }\left(\frac{\alpha uv+u^{2}k}{ua+\frac{\alpha a}{k}v}\right)=\underset{(u,v)\to (0,0)}{\lim }\left(\frac{\alpha uv+u^{2}k}{u(a+\frac{\alpha a}{k}\frac{v}{u})}\right)=\underset{(u,v)\to (0,0)}{\lim }\left(\frac{\alpha v+uk}{(a+\frac{\alpha a}{k}\frac{v}{u})}\right)$$

Podemos fazer uma substituição de variável novamente usando coordenadas polares $r^2=u^2+v^2$ e $u=r\cos \theta$, $v=r\sin \theta$, então:

$$\underset{r\to 0}{\lim }\left(\frac{\alpha r\sin \theta +kr\cos \theta }{(a+\frac{\alpha a}{k}\frac{r\sin \theta }{r\cos \theta })}\right)=\underset{r\to 0}{\lim }\left(\frac{\alpha kr\sin \theta +k^{2}r\cos \theta }{(ka+\alpha a\tan \theta )}\right)=0$$

E o limite:

$$\underset{(u,v)\to (0,0)}{\lim }\frac{\gamma vu}{v(\gamma \frac{a}{k}-c)}=\underset{(u,v)\to (0,0)}{\lim }\frac{\gamma u}{(\gamma \frac{a}{k}-c)}=0$$

E por fim, vamos estudar os pontos de estabilidade em torno de $(\frac{c}{\gamma },\frac{\gamma a-kc}{\alpha \gamma })$, tendo agora $u=x-\frac{c}{\gamma }$ e $v=\frac{\gamma a-kc}{\alpha \gamma }-y$, então, manipulando novamente, primeiro trabalhando com $\dot{x}$:

$$\frac{dx}{dt}=(xa-k{x}^2)-\left(\alpha xy\right)$$$$\frac{du}{dt}=((\frac{c}{\gamma }+u)a-k{(\frac{c}{\gamma }+u)}^2)-\left(\alpha (\frac{c}{\gamma }+u)(v+\frac{\gamma a-kc}{\alpha \gamma })\right)$$$$\frac{du}{dt}=(\frac{ca}{\gamma }+au-\frac{k{c}^2}{{\gamma }^2}-k{u}^2-\frac{2ck}{\gamma }u)-(\frac{c\alpha }{\gamma }v+\frac{\gamma ac-k{c}^2}{{\gamma }^2}+\alpha uv+\left[\frac{\gamma a-kc}{\gamma }\right]u)$$$$\frac{du}{dt}=((\frac{ca}{\gamma }-\frac{k{c}^2}{{\gamma }^2}-\frac{\gamma ac-k{c}^2}{{\gamma }^2})+(a-\frac{2ck}{\gamma }-\left[\frac{\gamma a-kc}{\gamma }\right]-ku)u-\frac{c\alpha }{\gamma }v)-\left(\alpha uv\right)$$$$\frac{du}{dt}=((-\frac{kc}{\gamma }-ku)u-\frac{c\alpha }{\gamma }v)-\left(\alpha uv\right)$$$$\frac{du}{dt}=-(\frac{kc}{\gamma }u+\frac{c\alpha }{\gamma }v)-(\alpha uv+k{u}^2)$$

E com $\dot{y}$:

$$\frac{dy}{dt}=(-yc)+\left(\gamma yx\right)$$$$\frac{dv}{dt}=(-(v+\frac{\gamma a-kc}{\alpha \gamma })c)+\left(\gamma (v+\frac{\gamma a-kc}{\alpha \gamma })(\frac{c}{\gamma }+u)\right)$$$$\frac{dv}{dt}=(-vc-\frac{ac}{\alpha }+\frac{k{c}^2}{\alpha \gamma })+(vc+\frac{\gamma ac-k{c}^2}{\alpha \gamma }+\gamma uv+\left[\frac{\gamma a-kc}{\alpha }\right]u)$$$$\frac{dv}{dt}=\gamma uv+\left[\frac{\gamma a-kc}{\alpha }\right]u$$

Calculando então os limites, relacionado a $u$:

$$\underset{(u,v)\to (0,0)}{\lim }\left(\frac{\alpha uv+k{u}^2}{\frac{kc}{\gamma }u+\frac{c\alpha }{\gamma }v}\right)=\underset{(u,v)\to (0,0)}{\lim }\left(\frac{\alpha uv+k{u}^2}{u(\frac{kc}{\gamma }+\frac{c\alpha }{\gamma }\frac{v}{u})}\right)=\underset{(u,v)\to (0,0)}{\lim }\left(\frac{\alpha v+ku}{\frac{kc}{\gamma }+\frac{c\alpha }{\gamma }\frac{v}{u}}\right)$$

E novamente fazendo a substituição $r^2=u^2+v^2$ e $u=r\cos \theta$, $v=r\sin \theta$, então:

$$\underset{r\to 0}{\lim }\left(\frac{\alpha r\sin \theta +kr\cos \theta }{\frac{kc}{\gamma }+\frac{c\alpha }{\gamma }\frac{r\sin \theta }{r\cos \theta }}\right)=\underset{r\to 0}{\lim }\left(\frac{r\alpha \gamma \sin \theta +rk\gamma \cos \theta }{kc+c\alpha \tan \theta }\right)=0$$

E a $v$:

$$\underset{(u,v)\to (0,0)}{\lim }\frac{\gamma uv}{\left[\frac{\gamma a-kc}{\alpha }\right]u}=\underset{(u,v)\to (0,0)}{\lim }\frac{\gamma v}{\left[\frac{\gamma a-kc}{\alpha }\right]}=0$$

Então os três pontos são semi-lineares. A partir disto, podemos analisar os tipos de estabilidades. Para $(0,0)$:

$$\left(\begin{array}{c}\dot{x}\\ \dot{y}\end{array}\right)=\left(\begin{array}{cc}a& 0\\ 0& -c\end{array}\right)\left(\begin{array}{c}x\\ y\end{array}\right)$$

Então:

$$(a-\lambda )(c+\lambda )=0$$ Novamente temos duas raízes reais de sinais opostos, então temos um ponto de sela, uma instabilidade. Para o segundo ponto:

$$\left(\begin{array}{c}\dot{u}\\ \dot{v}\end{array}\right)=\left(\begin{array}{cc}-a& -\frac{\alpha a}{k}\\ 0& \gamma \frac{a}{k}-c\end{array}\right)\left(\begin{array}{c}u\\ v\end{array}\right)$$

Então:

$$(-a-\lambda )((\frac{\gamma a}{k}-c)-\lambda )=0$$$$(a+\lambda )((\frac{\gamma a}{k}-c)-\lambda )=0$$

Temos duas raízes reais uma é negativa ${\lambda }_1=-a$. Porém a classificação do ponto depende dos parâmetros escolhidos. Se $(\frac{\gamma a}{k}-c)>0$ ou seja $\gamma a>ck$, temos uma instabilidade, uma sela, mas se $\gamma a<ck$, então temos um nó hiperbólico, ou seja, estabilidade. E por último:

$$\left(\begin{array}{c}\dot{u}\\ \dot{v}\end{array}\right)=\left(\begin{array}{cc}-\frac{kc}{\gamma }& -\frac{c\alpha }{\gamma }\\ \frac{\gamma a-kc}{\alpha }& 0\end{array}\right)\left(\begin{array}{c}u\\ v\end{array}\right)$$

Então:

$$(\frac{kc}{\gamma }+\lambda )\lambda +\frac{c\alpha }{\gamma }\left(\frac{\gamma a-kc}{\alpha }\right)=0$$$$\gamma {\lambda }^2+kc\lambda +c(\gamma a-kc)=0$$

Então:

$$\lambda =-\left(\frac{kc}{2\gamma }\right)\pm \frac{\sqrt{{\left(kc\right)}^2-4\gamma c(\gamma a-kc)}}{2\gamma }$$$$\lambda =-\left(\frac{kc}{2\gamma }\right)\pm \frac{\sqrt{[{\left(kc\right)}^2+4\gamma c^{2}k]-4c{\gamma }^{2}a}}{2\gamma }$$

Então o comportamento do sistema vai depender da escolha de parâmetros, porém como todas constantes são positivas $-\left(\frac{kc}{2\gamma }\right)$ vai ser sempre negativo, e a única forma de ser instável é se ${\left(kc\right)}^2+4\gamma c^{2}k>4c{\gamma }^{2}a$, para garantir que o número seja real, e ainda $\frac{\sqrt{{\left(kc\right)}^2+4\gamma c^{2}k-4c{\gamma }^{2}a}}{2\gamma }>\frac{kc}{2\gamma }$ para que o autovalor seja positivo. Analisando então essa última desigualdade:

$$\frac{\sqrt{{\left(kc\right)}^2+4\gamma c^{2}kc-4c{\gamma }^{2}a}}{2\gamma }>\frac{kc}{2\gamma }$$$$\sqrt{{\left(kc\right)}^2+4\gamma c^{2}k-4c{\gamma }^{2}a}>kc$$

Elevando ao quadrado:

$${\left(kc\right)}^2+4\gamma c^{2}k-4c{\gamma }^{2}a>{\left(kc\right)}^2$$$$4\gamma c^{2}k-4c{\gamma }^{2}a>0$$

Dividindo por $4c\gamma$

$$ck-\gamma a>0$$$$ck>\gamma a$$

Então essa é a condição para que o autovalor seja positivo. E olhando pra primeira desigualdade:

$${\left(kc\right)}^2+4\gamma c^{2}k>4c{\gamma }^{2}a$$$$\frac{k^{2}c}{4\gamma }+ck>\gamma a$$

Para garantir que o número seja real. Então temos duas desigualdades para satisfazer para que seja instável:

• $ck>\gamma a$
• $\frac{k^{2}c}{4\gamma }+ck>\gamma a$

Como $\left(\frac{k^{2}c}{4\gamma }\right)$ é necessariamente um termo positivo:

$$\gamma a<ck<\frac{k^{2}c}{4\gamma }+ck$$

Ou seja, podemos restringir a condição de instabilidade para $ck>\gamma a$ pois $ck<\frac{k^{2}c}{4\gamma }+ck$ é uma desigualdade válida para qualquer escolha de parâmetros válidos. Mas lembrando que para o outro ponto é estável se $ck>\gamma a$ , então sempre temos um ponto de equilíbrio estável e outro instável (além do a instabilidade na origem). Ou seja, por exemplo para valores mais baixos de de $c$ (termo de crescimento dos predadores) ou $k$ (termo de amortecimento da presa) temos coexistência de ambas as espécies, se aumentar estes parâmetros, somente a presa sobrevive. Além disso, quando temos $k=0$:

$$\lambda =\pm \frac{\sqrt{-4c{\gamma }^{2}a}}{2\gamma }=\pm i\sqrt{ca}$$

Ou seja, temos apenas a parte imaginária, e retornamos à estabilidade do centro. Dessa forma $k=0$ ou $k\ne 0$ determina se o sistema atinge equilíbrio ou permanece oscilatório, e depois o conjunto de parâmetros, caso atinja estabilidade, determina se ela é atingida com ou sem predadores.

Exemplo

Para estudarmos melhor uma situação, vamos escolher os parâmetros: $a=1$, $\alpha =0.5$, $k=0.5$, $c=0.75$,$\gamma =0.5$. Então os pontos de equilíbrio são:

• $(0,0)\to \lambda =\{1,-0.75\}$, uma sela, instável;
• $(2,0)\to \lambda =\{-1,0.25\}$, outro ponto instável, outra sela;
• $(1.5,0.5)\to \lambda =\{-0.375+0.22i,-0.375-0.22i\}$, agora temos um nó assintoticamente estável, pois as parte reais são negativas e as imaginárias são os conjugados.

E fazendo um rascunho do plano de fases, próximo de $(0,0)$, desprezando então os termos não-lineares:

$$\left(\begin{array}{c}\dot{x}\\ \dot{y}\end{array}\right)=\left(\begin{array}{cc}a& 0\\ 0& -c\end{array}\right)\left(\begin{array}{c}x\\ y\end{array}\right)=\left(\begin{array}{cc}1& 0\\ 0& -0.75\end{array}\right)\left(\begin{array}{c}x\\ y\end{array}\right)$$$$\left(\begin{array}{c}\dot{x}\\ \dot{y}\end{array}\right)=\left(\begin{array}{c}x\\ -0.75y\end{array}\right)$$

Próximo a $(\frac{a}{k},0)=(2,0)$:

$$\left(\begin{array}{c}\dot{u}\\ \dot{v}\end{array}\right)=\left(\begin{array}{cc}-a& -\frac{\alpha a}{k}\\ & \gamma \frac{a}{k}-c\end{array}\right)\left(\begin{array}{c}u\\ v\end{array}\right)$$$$\left(\begin{array}{c}\dot{u}\\ \dot{v}\end{array}\right)=\left(\begin{array}{cc}-1& -1\\ & 0.25\end{array}\right)\left(\begin{array}{c}u\\ v\end{array}\right)=\left(\begin{array}{cc}-1& -1\\ & 0.25\end{array}\right)\left(\begin{array}{c}x-2\\ y\end{array}\right)=\left(\begin{array}{c}\dot{x}\\ \dot{y}\end{array}\right)$$$$\left(\begin{array}{c}\dot{x}\\ \dot{y}\end{array}\right)=\left(\begin{array}{c}(2-x)-y\\ 0.25y\end{array}\right)$$

De $(\frac{c}{\gamma },\frac{\gamma a-kc}{\alpha \gamma })=(1.5,0.5)$:

$$\left(\begin{array}{c}\dot{u}\\ \dot{v}\end{array}\right)=\left(\begin{array}{cc}-\frac{kc}{\gamma }& -\frac{c\alpha }{\gamma }\\ \frac{\gamma a-kc}{\alpha }& 0\end{array}\right)\left(\begin{array}{c}u\\ v\end{array}\right)$$$$\left(\begin{array}{c}\dot{u}\\ \dot{v}\end{array}\right)=\left(\begin{array}{cc}-0.75& -0.75\\ 0.25& 0\end{array}\right)\left(\begin{array}{c}u\\ v\end{array}\right)=\left(\begin{array}{cc}-0.75& -0.75\\ 0.25& 0\end{array}\right)\left(\begin{array}{c}x-1.5\\ y-0.5\end{array}\right)=\left(\begin{array}{c}\dot{x}\\ \dot{y}\end{array}\right)$$$$\left(\begin{array}{c}\dot{x}\\ \dot{y}\end{array}\right)=\left(\begin{array}{c}0.75[(1.5-x)+(0.5-y)]\\ 0.25(x-1.5)\end{array}\right)$$

Plotando temos então:

A esquerda o rascunho do diagrama de fases, e a direita a evolução do sistema para as condições iniciais $(x_0,y_0)=(2,0.05)$.

Por curiosidade, se fazemos $k=0$: Para $(0,0)$:

$$\left(\begin{array}{c}\dot{x}\\ \dot{y}\end{array}\right)=\left(\begin{array}{cc}a& 0\\ 0& -c\end{array}\right)\left(\begin{array}{c}x\\ y\end{array}\right)=\left(\begin{array}{cc}1& 0\\ 0& -0.75\end{array}\right)\left(\begin{array}{c}x\\ y\end{array}\right)$$

E para $(\frac{c}{\gamma },\frac{a}{\alpha })=(1.5,2.0)$:

$$\left(\begin{array}{c}\dot{u}\\ \dot{v}\end{array}\right)=\left(\begin{array}{cc}0& -\frac{\alpha c}{\gamma }\\ \frac{\gamma a}{\alpha }& 0\end{array}\right)\left(\begin{array}{c}u\\ v\end{array}\right)$$$$\left(\begin{array}{c}\dot{x}\\ \dot{y}\end{array}\right)=\left(\begin{array}{cc}0& -0.75\\ 1& 0\end{array}\right)\left(\begin{array}{c}x-1.5\\ y-2\end{array}\right)=\left(\begin{array}{c}-0.75(y-2)\\ x-1.5\end{array}\right)$$Plotando o resultado:

A direita o rascunho do diagrama de fase, e a esquerda a evolução para as condições iniciais $(x_0,y_0)=(1,0.5)$.

Códigos

Os seguintes códigos escritos em Python foram utilizados para obter a solução numérica e plotar o rascunho do diagrama de fases próximo aos pontos de equilíbrio.

import matplotlib.pyplot as plt
import numpy as np

# LOTKA-VOLTERRA: Solução numérica
# Jhordan Silveira de Borba
# sbjhordan@gmail.com

def sol_lot():
    x=[]
    y=[]
    x.append(1)      # População inicial de presas
    y.append(0.5)    # População inicial de predadores
    N=400000         # Quantidade de passos
    d=0.0001         # Tamanhodos passos
    a=0              # 1: Amortecido, 0: Sem amortecimento
    
    for i in range(N-1):
        x.append(x[i]+d*(x[i]*(1-0.5*y[i])-0.5*x[i]*x[i]*a))
        y.append(y[i]+d*(y[i]*(-0.75+0.5*x[i])))
    #Plotamos a evolução temporal das frações de população
    X=np.arange(len(x)) #Eixo x
    plt.plot(X,x,'b-')
    plt.plot(X,y,'k-')
    plt.xlabel('Passo')
    plt.ylim(0,6)
    plt.show()

A função abaixo tem como finalidade plotar um rascunho do plano de fase próximo ao pontos de equilíbrio do sistema:

import matplotlib.pyplot as plt
import numpy as np

# LOTKA-VOLTERRA: Plano de fase
# Jhordan Silveira de Borba
# sbjhordan@gmail.com

def phase_lot():
    X = np.arange(0, 4, 0.2)  # Eixo x
    Y = np.arange(0, 4, 0.2)  # Eixo Y
    U,V=np.meshgrid(X,Y)
    
    p1=[0.,0.]           # Ponto de equilíbrio 1
    p2=[1.5,2.0]         # Ponto de equilíbrio 2
    p3=[1.5*100,0.5*100] # Ponto de equilíbrio 3
    
    c=0
    for x in X:
        l=0
        for y in Y:
            #Distâncias
            d1=np.sqrt((x-p1[0])*(x-p1[0])+(y-p1[1])*(y-p1[1]))
            d2=np.sqrt((x-p2[0])*(x-p2[0])+(y-p2[1])*(y-p2[1]))
            d3=np.sqrt((x-p3[0])*(x-p3[0])+(y-p3[1])*(y-p3[1]))
            
            #encontrar o ponto de equilíbrio mais próximo
            p=1
            if(d2<d1):
                p=2
                if(d3<d2):
                    p=3
            elif(d3<d1):
                p=3
                
            # Calculamos o vetor de variação do estado baseado no ponto mais próximo:
            # [dx/dt,dy/dt]=[a,b]
            if(p==1):
                a=x  
                b=-0.75*y
            elif(p==2):
                a=-0.75*(y-2)
                b=x-1.5
            elif(p==3):
                a=0.75*((1.5-x)+(0.5-y))
                b=0.25*(x-1.5)
            else:
                print("Algo deu errado")
            m=np.sqrt(a*a+b*b) # Módulo do vetor para normalizar
            if(m==0):
                m=1
            U[l,c]=a/m
            V[l,c]=b/m
            l=l+1
        c=c+1
    
    # Plotamos o resultado
    fig, ax = plt.subplots()
    ax.quiver(X, Y, U, V)      # Os vetores
    plt.plot(p1[0],p1[1],'ro') # O ponto de equilíbrio 1
    plt.plot(p2[0],p2[1],'go') # O ponto de equilíbrio 2
    plt.plot(p3[0],p3[1],'go') # O ponto de equilíbrio 3
    plt.show()

Principais materiais utilizados

1. Numerical Methods (Instituto de Tecnologia de Massachusetts)
2. Numerical Solution of Ordinary Differential Equations (R. Sureshkumar,Instituto de Tecnologia de Massachusetts)
3. A estes materiais, somam-se os vistos em Modelo de Lotka-Volterra

Citações

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