Modelo de Lotka-Volterra amortecido

De Física Computacional
Revisão de 17h44min de 2 de maio de 2021 por Jhordan (discussão | contribs)
(dif) ← Edição anterior | Revisão atual (dif) | Versão posterior → (dif)
Ir para navegação Ir para pesquisar


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.

Onde é uma constante positiva. Agora temos três pontos de equilíbrio, novamente temos . Mas temos outro ponto quando apenas :

Então 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:

Isolando então na segunda equação e substituindo na primeira:

Então:

Então nosso outro ponto de equilíbrio é dado por. 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 :

Para e fazendo a mudança de variáveis , e , temos:

E:

Então:

Podemos fazer uma substituição de variável novamente usando coordenadas polares[1] e , , então:

E o limite:

E por fim, vamos estudar os pontos de estabilidade em torno de , tendo agora e , então, manipulando novamente, primeiro trabalhando com :

E com :

Calculando então os limites, relacionado a :

E novamente fazendo a substituição e , , então:

E a :

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

Então:

Novamente temos duas raízes reais de sinais opostos, então temos um ponto de sela, uma instabilidade. Para o segundo ponto:

Então:

Temos duas raízes reais uma é negativa . Porém a classificação do ponto depende dos parâmetros escolhidos. Se ou seja , temos uma instabilidade, uma sela, mas se , então temos um nó hiperbólico, ou seja, estabilidade. E por último:

Então:

Então:

Então o comportamento do sistema vai depender da escolha de parâmetros, porém como todas constantes são positivas vai ser sempre negativo, e a única forma de ser instável é se , para garantir que o número seja real, e ainda para que o autovalor seja positivo. Analisando então essa última desigualdade:

Elevando ao quadrado:

Dividindo por

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

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

Como é necessariamente um termo positivo:

Ou seja, podemos restringir a condição de instabilidade para pois é uma desigualdade válida para qualquer escolha de parâmetros válidos. Mas lembrando que para o outro ponto é estável se , 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 (termo de crescimento dos predadores) ou (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 :

Ou seja, temos apenas a parte imaginária, e retornamos à estabilidade do centro. Dessa forma ou 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: , , , ,. Então os pontos de equilíbrio são:

  • , uma sela, instável;
  • , outro ponto instável, outra sela;
  • , 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 , desprezando então os termos não-lineares:

Próximo a :

De :

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 .


Por curiosidade, se fazemos : Para :

E para :

Plotando o resultado:

A direita o rascunho do diagrama de fase, e a esquerda a evolução para as condições iniciais .
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
  1. Limits and Continuity (Rebecca Noonan-Heale ,Universidade de Utah)


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