Cálculo de pontos fixos e sua estabilidade.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Cálculo de pontos fixos e sua estabilidade"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Este é um guia para calcular os pontos fixos de um modelo, e posteriormente suas condições para estabilidade, utilizando a biblioteca *Sympy*."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Por que usar Sympy?\n",
    "\n",
    "O *Sympy* é um pacote que trabalha com linguagem matemática simbólica. Ou seja, o foco dos cálculos são símbolos e não números. Com o *Sympy* é possível trabalhar com expressões matemáticas sem precisar designar valore às variáveis e parâmetros. Com esta ferramenta, podemos calcular expressões analíticas dos pontos fixos de um modelo, assim como sua matriz jacobiana "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Começamos o algoritmo  importando a biblioteca *Sympy*:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Definindo um modelo\n",
    "\n",
    "O próximo passo é definir o modelo a ser analisado. Nesse guia, escolhemos um modelo de duas espécies que competem pelo mesmo recurso, dado por:\n",
    "\n",
    "$$\\displaystyle \\dfrac{d N_1}{dt} = r_1 N_1 \\left(1- \\dfrac{N_1}{K_1} - b_{12}\\dfrac{N_2}{K_1}\\right)$$\n",
    "$$\\displaystyle \\dfrac{d N_2}{dt} = r_2 N_2 \\left(1- \\dfrac{N_2}{K_2} - b_{21}\\dfrac{N_1}{K_2}\\right)$$\n",
    "\n",
    "Este é conhecido como o modelo Lotka-Volterra para competição de duas espécies. Os parâmetros $r_1$ e $r_2$ correspondem às taxas de crescimento das espécies 1 e 2. Da mesma forma, $K_1$ e $K_2$ são as capacidades suporte de cada espécie. O parâmetro $b_{12}$ descreve o impacto da competição na espécie 1. De maneira análoga, $b_{21}$ é o impacto da competição na espécie 2. Para simplificar as equações, fazemos algumas manipulações para reduzir o número de parâmetros do modelo, conforme realizado nas aulas do curso. Temos como resultado o seguinte modelo:\n",
    "\n",
    "$$\\displaystyle \\dfrac{d u_1}{dt} =  r_1 u_1 \\left(1- u_1 - \\alpha_{12}u_2\\right)$$\n",
    "$$\\displaystyle \\dfrac{d u_2}{dt} = r_2 u_2 \\left(1- u_2 - \\alpha_{21}u_1\\right)$$\n",
    "\n",
    "De modo que $u_1=\\frac{N_1}{K_1}$, $u_2= \\frac{N_2}{K_2}$. Os novos parâmetros são dados por:  $\\alpha_{12}= \\frac{K_2}{K_1}b_{12}$ e $\\alpha_{21}= \\frac{K_1}{K_2}b_{21}$\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Começamos então, com o Sympy, declarando com quais variáveis e parâmetros estamos trabalhando. Fazemos isso a partir da função *symbols*:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Declaramos as variáveis u1 e u2 e as armazenamos em uma lista\"Variaveis\" \n",
    "u_1,u_2=symbols('u_1 u_2')\n",
    "Variaveis=[u_1,u_2]\n",
    "\n",
    "# Declaramos cada parâmetro e os armazenamos em uma lista \"Parametros\" \n",
    "r_1, r_2, alpha_12, alpha_21=symbols('r_1 r_2 alpha_12 alpha_21')\n",
    "Parametros=[rho, alpha_12, alpha_21]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Agora basta escrever as equações do modelo a partir dos símbolos que acabamos de definir. É conveniente guardar o lado direito das equações do modelo em uma lista \"Modelos\":"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Guardamos o lado direito das equações do modelo na lista Modelo\n",
    "# É muito importante preservar a ordem que construímos a lista Variaveis.\n",
    "# Como guardamos primeiro u_1 e depois o u_2, a lista das equações do modelo \n",
    "# deve ter como a primeira entrada a derivada de u_1 em relação ao tempo\n",
    "# e em seguida a derivada de u_2 em relação ao tempo.\n",
    "\n",
    "\n",
    "#                   du_1/dt                                  du_2/dt\n",
    "Modelo=[r_1*u_1*(1- u_1  - alpha_12*u_2), r_2*u_2*(1- u_2  - alpha_21*u_1)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pronto! o modelo está definido e agora podemos manipular as expressões com funções do *Sympy*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pontos Fixos\n",
    "\n",
    "Um ponto fixo é uma configuração de um sistema em que não ha variação das variáveis de interesse. No contexto os modelos ecológicos, seriam pontos em que as derivadas dos tamanhos populacionáis são nulas, ou seja, não há variação das populações estudadas.\n",
    "\n",
    "Para o modelo de competição que usamos como exemplo, um ponto fixo é qualquer ponto $u^* = \\{u_1^*, u_2^*\\}$ que resulte em:\n",
    "\n",
    "$$\\displaystyle \\left.\\dfrac{du_1}{dt}\\right|_{u^*}= \\left.\\dfrac{du_2}{dt}\\right|_{u^*} =0$$\n",
    "\n",
    "Ou seja, as derivadas de $u_1$ e $u_2$ avaliadas em $u^*$ são nulas. Dessa forma, o lado esquerdo das equações do modelo serão iguais a zero quando o sistema se encontra em um ponto fixo:\n",
    "\n",
    "$$\\displaystyle r_1 u_1^* \\left(1- u_1^* - \\alpha_{12}u_2^*\\right)=0$$\n",
    "$$\\displaystyle  r_2 u_2^* \\left(1- u_2^* - \\alpha_{21}u_1^*\\right)=0$$\n",
    "\n",
    "Encontrar os pontos fixos do sistema entã se torna um problema de aritimética. Basta encontrar os valores de $u_1^*$ e $u_2^*$ que satisfazem o sistema acima."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Uma das grandes vantagens de usar o *Sympy* é poder resolver este sistema, que muitas vezes pode ser bastante complicado, com apenas uma função: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left\\{\\left( 0, \\  0\\right), \\left( 0, \\  1\\right), \\left( 1, \\  0\\right), \\left( \\frac{\\alpha_{12} - 1}{\\alpha_{12} \\alpha_{21} - 1}, \\  \\frac{\\alpha_{21} - 1}{\\alpha_{12} \\alpha_{21} - 1}\\right)\\right\\}$"
      ],
      "text/plain": [
       "{(0, 0), (0, 1), (1, 0), ((alpha_12 - 1)/(alpha_12*alpha_21 - 1), (alpha_21 - 1)/(alpha_12*alpha_21 - 1))}"
      ]
     },
     "execution_count": 67,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# A função nonlinsolve iguala o lado direito das equações do modelo a zero  \n",
    "# e nos retorna as expressões de u_1 e u_2 que satisfazem este sistema.\n",
    "# Devemos colocar como input o lado direito das equações e as variáveis de\n",
    "# interesse:\n",
    "Pontos_Fixos=nonlinsolve(Modelo, Variaveis)\n",
    "\n",
    "# As vezes o resultado pode ser sair bastante desorganizado.\n",
    "# A função simplify nos ajuda organizando e simplificando \n",
    "# a expressão matemática:\n",
    "Pontos_Fixos= simplify(Pontos_Fixos)\n",
    "\n",
    "# Por fim, printamos o resultado:\n",
    "Pontos_Fixos"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sucesso! Cada par do output, encapsulado por parênteses, representa um ponto fixo, de modo que o primeiro termo representa o valor de $u_1^*$ e o segundo o valor de $u_2^*$. Dessa forma, a ordem da lista \"Variaveis\" é preservada nos resultados da função. \n",
    "\n",
    "Vamos agora rapidamente analisar os resultados. Temos:\n",
    "\n",
    "- Um ponto fixo em que o tamanho populacional das duas espécies é nulo: $u_1^*=u_2^*=0$\n",
    "- Um ponto fixo onde a espécie 2 se extingue, $u_2^*=0$, e a espécie 1 se encontra em sua capacidade suporte $u_1^*=1$\n",
    "- Um ponto fixo onde a espécie 1 se extingue, $u_1^*=0$, e a espécie 2 se encontra em sua capacidade suporte $u_2^*=1$\n",
    "- Um ponto de coexistência de ambas as espécies"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Matriz Jacobiana"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "De modo a seguir a análise do modelo, devemos escrever a matriz jacobiana modelo. Neste sistemas, de uma maneira informal, podemos pensar que a Jacobiana expressa os efeitos de cada espécie sobre ela e sobre a outra, para uma certa combinação de densidades $u_1$ e $u_2$.  O *Sympy* nos permite realizar esse passo com uma única função!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}r_{1} \\left(- \\alpha_{12} u_{2} - 2 u_{1} + 1\\right) & - \\alpha_{12} r_{1} u_{1}\\\\- \\alpha_{21} r_{2} u_{2} & r_{2} \\left(- \\alpha_{21} u_{1} - 2 u_{2} + 1\\right)\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[r_1*(-alpha_12*u_2 - 2*u_1 + 1),               -alpha_12*r_1*u_1],\n",
       "[              -alpha_21*r_2*u_2, r_2*(-alpha_21*u_1 - 2*u_2 + 1)]])"
      ]
     },
     "execution_count": 68,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Para calcular a matriz jacobiana do modelo, devemos transformar as listas\n",
    "# Variaveis e Modelo em matrizes do Sympy.\n",
    "M=Matrix(Modelo)\n",
    "V=Matrix(Variaveis)\n",
    "\n",
    "#Chamamos a função Jacobian que nos dá a matriz jacobiana do modelo M em função das variáveis V\n",
    "J=M.jacobian(V)\n",
    "\n",
    "#Simplificando temos:\n",
    "simplify(J)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Estabilidade de cada ponto fixo"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "De modo a determinar as condições de estabilidade de um ponto fixo, devemos avaliar a matriz jacobiana do sistema neste ponto. Os autovalores da matriz estabelecem as condições de estabilidade:\n",
    "\n",
    " * Caso a parte real de todos os autovalores sejam negativos, este é um ponto fixo **estável**.\n",
    " * Caso a parte real de algum autovalor for positivo, denominamos este um ponto fixo **instável**.\n",
    " * A existência de algum autovalor com parte real nula torna este ponto fixo **Indeterminado**.\n",
    "\n",
    "As condições acima são para a parte real do autovalor apoenas. Em muitos casos os autovalores terão também uma parte imaginária, o que nos informa sobre oscilações no equilíbrio. Em sistemas computacionais simbólicos como o que estamos usando aqui, a parte imaginária parte estará expressa pelo número $i$, ou por uma raiz quadrada de um número negativo, somada à parte real. Há casos ainda em que você só terá esta parte imaginária. Nos exemplos deste tutorial teremos autovalores sem partes imaginárias.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Vamos obter as condições de estabilidade de cada ponto fixo! Para encontrar os autovalores da Jacobiana com o *sympy* basta utilizar a função \"subs\", que substitui as expressões das pontos fixos obtidas previamente na matriz jacobiana."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Ponto Fixo 1 ($u_1=0$ e $u_2=0$)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left\\{ r_{1} : 1, \\  r_{2} : 1\\right\\}$"
      ],
      "text/plain": [
       "{r_1: 1, r_2: 1}"
      ]
     },
     "execution_count": 69,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Substituindo u_1=0 e u_2=0\n",
    "J0=J.subs([(u_1,0), (u_2,0)])\n",
    "\n",
    "#simplificado\n",
    "simplify(J0.eigenvals())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "O output pode ser lido da seguinte forma:\n",
    "\n",
    "{Autovalor 1 : Multiplicidade do autovalor 1, Autovalor 2 : Multiplicidade do autovalor 2}\n",
    "\n",
    "Não estamos interessados na multiplicidade dos autovalores, apenas em seu sinal. Vemos no output que os autovalores correspondem às taxas de crescimento intrínseco, $r_1$ e $r_2$, sem uma parte imaginária. Então o que precisamos inspecionar é a parte real do autovalor, que neste caso corresponde às taxas de crescimento instríseco das duas espécies. Aplicando os critérios que descrevemos acima, este ponto fixo é estável apenas se os dois valores forem menores que zero:\n",
    "\n",
    " $$r_1 <0$$ e $$r_2 <0$$\n",
    " \n",
    "Ou seja, caso a taxa de crescimento de ambas as espécies sejam negativas, o ponto fixo que representa a extinção de ambas é estável. Isso significa que, sob esta condição, se as duas populações são zero ($u_1=0$, $u_2=0$), retornarão a zero mesmo de uma pequena perturbação. Por outro lado, basta que uma das espécies tenha taxa de crescimento de positivo que este ponto não será estável. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Ponto Fixo 2 ($u_1=0$ e $u_2=1$)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left\\{ - r_{2} : 1, \\  - r_{1} \\left(\\alpha_{12} - 1\\right) : 1\\right\\}$"
      ],
      "text/plain": [
       "{-r_2: 1, -r_1*(alpha_12 - 1): 1}"
      ]
     },
     "execution_count": 70,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Substituindo u_1=0 e u_2=K_2\n",
    "J1=J.subs([(u_1,0), (u_2,1)])\n",
    "\n",
    "#simplificado\n",
    "simplify(J1.eigenvals())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Aqui vemos que os autovalores são o negativo da taxa de crescimento intríseco de uma população, e uma função da taxa de crescimento da outra espécie e do efeito da espécie da espécie 1 sobre a 2. Dessa forma, este ponto fixo é estável apenas se:\n",
    "\n",
    " $$r_1\\left(1-\\alpha_{12}\\right)<0$$ e $$r_2 >0$$\n",
    " \n",
    "Ou seja, a situação em que apenas a espécie 2 persiste só é estável se a taxa de crescimento desta espécie é positiva **e** se o efeito competitivo da espécie 1 sobre a 2 é fraco, menor do que o efeito dela sobre ela mesma (ou seja, $\\alpha_{12}<1$)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Ponto fixo 3 ($u_1 =1$, $u_2 = 0$)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left\\{ - r_{1} : 1, \\  - r_{2} \\left(\\alpha_{21} - 1\\right) : 1\\right\\}$"
      ],
      "text/plain": [
       "{-r_1: 1, -r_2*(alpha_21 - 1): 1}"
      ]
     },
     "execution_count": 71,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Substituindo u_1=K_1 e u_2=0\n",
    "J2=J.subs([(u_1,1), (u_2,0)])\n",
    "\n",
    "#simplificado\n",
    "simplify(J2.eigenvals())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dessa forma, este ponto fixo é estável apenas se:\n",
    "\n",
    " $$r_2\\left(1-\\alpha_{21}\\right)<0$$ e $$r_1 >0$$\n",
    " \n",
    "Ou seja, o caso oposto ao anterior, que descreve a estabilidade da situação em que apenas a espécie 2 persiste."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Ponto fixo de coexistência"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left\\{ - \\frac{\\alpha_{12} r_{1} + \\alpha_{21} r_{2} - r_{1} - r_{2}}{2 \\left(\\alpha_{12} \\alpha_{21} - 1\\right)} - \\frac{\\sqrt{4 \\alpha_{12}^{2} \\alpha_{21}^{2} r_{1} r_{2} - 4 \\alpha_{12}^{2} \\alpha_{21} r_{1} r_{2} + \\alpha_{12}^{2} r_{1}^{2} - 4 \\alpha_{12} \\alpha_{21}^{2} r_{1} r_{2} + 2 \\alpha_{12} \\alpha_{21} r_{1} r_{2} - 2 \\alpha_{12} r_{1}^{2} + 2 \\alpha_{12} r_{1} r_{2} + \\alpha_{21}^{2} r_{2}^{2} + 2 \\alpha_{21} r_{1} r_{2} - 2 \\alpha_{21} r_{2}^{2} + r_{1}^{2} - 2 r_{1} r_{2} + r_{2}^{2}}}{2 \\left(\\alpha_{12} \\alpha_{21} - 1\\right)} : 1, \\  - \\frac{\\alpha_{12} r_{1} + \\alpha_{21} r_{2} - r_{1} - r_{2}}{2 \\left(\\alpha_{12} \\alpha_{21} - 1\\right)} + \\frac{\\sqrt{4 \\alpha_{12}^{2} \\alpha_{21}^{2} r_{1} r_{2} - 4 \\alpha_{12}^{2} \\alpha_{21} r_{1} r_{2} + \\alpha_{12}^{2} r_{1}^{2} - 4 \\alpha_{12} \\alpha_{21}^{2} r_{1} r_{2} + 2 \\alpha_{12} \\alpha_{21} r_{1} r_{2} - 2 \\alpha_{12} r_{1}^{2} + 2 \\alpha_{12} r_{1} r_{2} + \\alpha_{21}^{2} r_{2}^{2} + 2 \\alpha_{21} r_{1} r_{2} - 2 \\alpha_{21} r_{2}^{2} + r_{1}^{2} - 2 r_{1} r_{2} + r_{2}^{2}}}{2 \\left(\\alpha_{12} \\alpha_{21} - 1\\right)} : 1\\right\\}$"
      ],
      "text/plain": [
       "{-(alpha_12*r_1 + alpha_21*r_2 - r_1 - r_2)/(2*(alpha_12*alpha_21 - 1)) - sqrt(4*alpha_12**2*alpha_21**2*r_1*r_2 - 4*alpha_12**2*alpha_21*r_1*r_2 + alpha_12**2*r_1**2 - 4*alpha_12*alpha_21**2*r_1*r_2 + 2*alpha_12*alpha_21*r_1*r_2 - 2*alpha_12*r_1**2 + 2*alpha_12*r_1*r_2 + alpha_21**2*r_2**2 + 2*alpha_21*r_1*r_2 - 2*alpha_21*r_2**2 + r_1**2 - 2*r_1*r_2 + r_2**2)/(2*(alpha_12*alpha_21 - 1)): 1, -(alpha_12*r_1 + alpha_21*r_2 - r_1 - r_2)/(2*(alpha_12*alpha_21 - 1)) + sqrt(4*alpha_12**2*alpha_21**2*r_1*r_2 - 4*alpha_12**2*alpha_21*r_1*r_2 + alpha_12**2*r_1**2 - 4*alpha_12*alpha_21**2*r_1*r_2 + 2*alpha_12*alpha_21*r_1*r_2 - 2*alpha_12*r_1**2 + 2*alpha_12*r_1*r_2 + alpha_21**2*r_2**2 + 2*alpha_21*r_1*r_2 - 2*alpha_21*r_2**2 + r_1**2 - 2*r_1*r_2 + r_2**2)/(2*(alpha_12*alpha_21 - 1)): 1}"
      ]
     },
     "execution_count": 72,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Ponto de coexistência\n",
    "J3=J.subs([(u_1,(alpha_12-1)/(alpha_12*alpha_21 -1)), (u_2,(alpha_21-1)/(alpha_12*alpha_21 -1))])\n",
    "\n",
    "#simplificando\n",
    "simplify(J3.eigenvals())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Opa! Aqui percebemos uma das vantagens do *Sympy*. Uma expressão como essa é bastante complicada e realizar estas contas na mão até chegar às expressões acima seria muito trabalhoso.\n",
    "\n",
    "Claro que não é muito prático trabalhar analiticamente com uma expressão desta, mas sempre podemos apenas introduzir os valores dos parâmetros na expressão e verificar as condições de estabilidade graficamente. Por exemplo, você pode investigar um gráfico do valor de cada expressão acima com alguns valores de parâmetros fixos, e um deles apenas variando.\n",
    "\n",
    "Para facilitar seu trabalho, você por obter as expressões acima na sintaxe de linguagens computacionais, para criar uma função em seu ambiente de programação. para isso, crie uma lista com os autovalores, e use a função \"print\" sobre cada elemento da lista:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-(alpha_12*r_1 + alpha_21*r_2 - r_1 - r_2)/(2*(alpha_12*alpha_21 - 1)) + sqrt(4*alpha_12**2*alpha_21**2*r_1*r_2 - 4*alpha_12**2*alpha_21*r_1*r_2 + alpha_12**2*r_1**2 - 4*alpha_12*alpha_21**2*r_1*r_2 + 2*alpha_12*alpha_21*r_1*r_2 - 2*alpha_12*r_1**2 + 2*alpha_12*r_1*r_2 + alpha_21**2*r_2**2 + 2*alpha_21*r_1*r_2 - 2*alpha_21*r_2**2 + r_1**2 - 2*r_1*r_2 + r_2**2)/(2*(alpha_12*alpha_21 - 1))\n"
     ]
    }
   ],
   "source": [
    "Auto_vals=simplify(J3.eigenvals())\n",
    "Auto_vals=list(Auto_vals)\n",
    "print(Auto_vals[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Também podemos utilizar outros métodos, como integração numérica do modelo para realizar este tipo de análise."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}