Considere o seguinte sistema de equações diferenciais abaixo, que representa o modelo das interações entre predador e presa em um determinado ecossistema:
\begin{aligned} & \frac{d y_{1}}{d x}=\alpha\left(1-y_{2}\right) y_{1} \\\ & \frac{d y_{2}}{d x}=-\beta\left(1-y_{1}\right) y_{2} \end{aligned}
onde y_{1} e y_{2} representam, respectivamente, a população de predadores e a população de presas, com y_{1}(0)=1,5 e y_{2}(0)=0,75. Os parâmetros podem ser tomados como \alpha=\beta=1. Para este sistema, pede-se:
a) Determine as soluções estacionárias deste modelo e a sua estabilidade.
b) Encontre o tamanho máximo do passo para o método de Euler explícito, de modo que a solução seja estável e não-oscilatória.
c) Determine a equação de iteração do método de Euler implícito.

Passo 1

O enunciado fornece o sistema
\begin{aligned} & \frac{d y_{1}}{d x}=\alpha\left(1-y_{2}\right) y_{1},\\[4pt] & \frac{d y_{2}}{d x}=-\beta\left(1-y_{1}\right) y_{2}, \end{aligned}
com condições iniciais \;y_{1}(0)=1{,}5,\;y_{2}(0)=0{,}75\; e parâmetros \alpha=\beta=1. As perguntas são:
a) Determinar as soluções estacionárias (pontos de equilíbrio) e a estabilidade de cada uma.
b) Encontrar o tamanho máximo do passo para o método de Euler explícito, de modo que a solução seja estável e não-oscillatória (interprete que a análise se faz pela linearização e critério de estabilidade do método sobre autovalores reais).
c) Determinar a equação de iteração do método de Euler implícito (Backward Euler) aplicada ao sistema.
As informações relevantes dadas: as expressões do lado direito do sistema (com \alpha=\beta=1) e as condições iniciais (que servem só como exemplo numérico, não são necessárias para a análise dos equilíbrios e das fórmulas do método numérico).]

Passo 2

Encontrar os pontos de equilíbrio: são os pontos (y_1,y_2) tais que os termos à direita zero. Com \alpha=\beta=1 o sistema fica
\begin{aligned} & \frac{d y_{1}}{d x}=(1-y_{2})\,y_{1},\\[4pt] & \frac{d y_{2}}{d x}=-(1-y_{1})\,y_{2}. \end{aligned}
Para equilíbrio precisamos de (1-y_{2})y_{1}=0 e -(1-y_{1})y_{2}=0. As possibilidades são y_1=0 ou y_2=1 da primeira equação, e y_2=0 ou y_1=1 da segunda. Combinando temos quatro pontos candidatos:
(0,0),\quad (0,1),\quad (1,0),\quad (1,1).
Esses são todos os pontos estacionários.]

Passo 3

Para estudar estabilidade linearizamos o sistema em torno de um ponto de equilíbrio usando a matriz Jacobiana. Escrevendo f_1(y_1,y_2)=(1-y_2)y_1 e f_2(y_1,y_2)=-(1-y_1)y_2, os parciais são:
\frac{\partial f_1}{\partial y_1}=1-y_2,\qquad \frac{\partial f_1}{\partial y_2}=-y_1,
\frac{\partial f_2}{\partial y_1}=y_2,\qquad \frac{\partial f_2}{\partial y_2}=y_1-1.
Logo a Jacobiana é
J(y_1,y_2)= \begin{pmatrix} 1-y_2 & -y_1\\[4pt] y_2 & y_1-1 \end{pmatrix}.
A estabilidade linear de um equilíbrio é dada pelos autovalores de J nesse ponto: autovalores com parte real negativa → assintoticamente estável (linearmente); parte real positiva → instável; autovalores puramente imaginários ou zeros → caso não hiperbólico (necessita análise não linear).]

Passo 4

Avaliar J em cada equilíbrio e calcular autovalores.
No ponto (0,0):
J(0,0)=\begin{pmatrix}1 & 0\\[4pt]0 & -1\end{pmatrix}.
Os autovalores são \lambda_1=1 e \lambda_2=-1. Há um autovalor positivo e um negativo, portanto (0,0) é uma sela (saddle) e é instável (trajetórias são atraídas numa direção e repelidas em outra).]

Passo 5

No ponto (1,0):
J(1,0)=\begin{pmatrix}1 & -1\\[4pt]0 & 0\end{pmatrix}.
Característica: determinante =0. Os autovalores obtêm-se resolvendo \det(J-\lambda I)=0, isto dá \lambda(\lambda-1)=0, então \lambda=1 e \lambda=0. Como existe um autovalor positivo \lambda=1, o ponto (1,0) é instável (não é um nó estável); o autovalor zero indica não-hiperbolicidade numa direção.]

Passo 6

No ponto (0,1):
J(0,1)=\begin{pmatrix}0 & 0\\[4pt]1 & -1\end{pmatrix}.
Os autovalores satisfazem \lambda(\lambda+1)=0, logo \lambda=0 e \lambda=-1. Há um autovalor negativo e um nulo; isso significa que linearmente o ponto tem uma direção estável (autovalor -1) e uma direção neutra (autovalor 0). Portanto (0,1) é não-hiperbólico: não é assintoticamente estável no sentido linearmente uniforme (a presença de \lambda=0 impede concluir estabilidade assintótica pelo linearizado), mas tem estabilidade em uma direção e neutralidade em outra.]

Passo 7

No ponto (1,1):
J(1,1)=\begin{pmatrix}0 & -1\\[4pt]1 & 0\end{pmatrix}.
A matriz é a de uma rotação; os autovalores são as raízes de \lambda^2+1=0, ou seja \lambda=\pm i (puro imaginário). Logo a linearização indica um centro (autovalores puramente imaginários). Isso é não-dissipativo linearmente: não há amortecimento nem crescimento exponencial linear; precisamos verificar a natureza não linear para concluir estabilidade mais precisa.]

Passo 8

Explicação não-linear para (1,1): procurar uma quantidade conservada (integral primeira) que prove que as órbitas são fechadas ao redor de (1,1). Vamos construir H(y_1,y_2) tal que \dfrac{dH}{dx}=0 ao longo das soluções. Tentar a forma
H(y_1,y_2)=\ln(y_1y_2)-y_1-y_2.
Calculamos a derivada ao longo das soluções:
\frac{dH}{dx}=\frac{1}{y_1}\frac{dy_1}{dx}+\frac{1}{y_2}\frac{dy_2}{dx}-\frac{dy_1}{dx}-\frac{dy_2}{dx}.
Substituindo \dfrac{dy_1}{dx}=(1-y_2)y_1 e \dfrac{dy_2}{dx}=-(1-y_1)y_2 e simplificando (é um cálculo direto, feito termo a termo), obtém-se \dfrac{dH}{dx}=0. Logo H é conservada e as curvas de nível de H são trajetórias do sistema. Em particular sistemas com uma integral primeira desse tipo têm órbitas fechadas em volta do ponto crítico central. Portanto (1,1) é um centro não linear: estável no sentido de Lyapunov (soluções próximas permanecem próximas) mas não assintoticamente estável — as soluções próximas são periódicas (órbitas fechadas) e não convergem para (1,1).

Passo 9

Resumo das classificações (conclusão da parte a), com linguagem direta):
  • Ponto de equilíbrio (0,0): autovalores \lambda=1,-1. Tipo: sela (saddle). Instável.
  • Ponto de equilíbrio (1,0): autovalores \lambda=1,0. Instável (existe uma direção de crescimento).
  • Ponto de equilíbrio (0,1): autovalores \lambda=0,-1. Não-hiperbólico; possui uma direção estável e uma neutra (não assintoticamente estável pelo linearizado).
  • Ponto de equilíbrio (1,1): autovalores \lambda=\pm i. Centro não linear com integral primeira \;H=\ln(y_1y_2)-y_1-y_2\;: órbitas fechadas em torno de (1,1). Estável no sentido de Lyapunov, não assintoticamente estável.

Passo 10

Passando para a parte b) (critérios de estabilidade do método de Euler explícito). Vamos lembrar um resultado básico e muito usado em provas: ao aplicar o método de Euler explícito a uma equação linear escalar teste \dfrac{dz}{dx}=\lambda z, o fator de amplificação por passo h é
\mu(h)=1+h\lambda,
ou seja z_{n+1}=\mu(h) z_n. Para que a discretização seja estável e não-oscillatória numa componente associada a um autovalor real \lambda devemos ter \mu(h) real e 0<\mu(h)\le 1. (Explicação: \mu com módulo menor que 1 dá decaimento; \mu>0 evita alternância de sinal — ou seja oscilações de sinal a cada passo.)

Passo 11

Aplicar esse critério aos autovalores reais do Jacobiano do sistema. Os autovalores reais que aparecem nos equilíbrios são \lambda=1,\;0,\;-1 (ver passos anteriores). Para cada \lambda a condição 0<1+h\lambda\le 1 traduz-se em:
  • se \lambda=-1: 0<1-h\le 1 \Rightarrow 0<h\le 1.
  • se \lambda=0: 1+h\cdot 0=1 (marginalmente neutro; Euler preserva exatamente esse modo sem decaimento).
  • se \lambda=1: 0<1+h\le 1\Rightarrow 0<-h impossível para h>0. Ou seja, uma componente com \lambda=1 não pode ser tornada estável pelo Euler explícito (o fator \mu=1+h>1 cresce qualquer h>0).

Passo 12

Interpretação final para a pergunta b): como o sistema possui um autovalor positivo \lambda=1 em alguns pontos (por exemplo em (0,0) e (1,0)), não existe passo h>0 que torne globalmente estável o método explícito em torno desses equilíbrios instáveis (não se pode numericamente estabilizar um modo que é linearmente instável). Contudo, se a intenção é garantir que as componentes correspondentes a autovalores negativos não oscillatem (por exemplo a direção com \lambda=-1 em (0,1)), a condição fornece:
h_{\max}=1.
Ou seja, para a componente com \lambda=-1 devemos escolher 0<h\le 1 para ter fator de amplificação no intervalo [0,1] (não-oscilatório e não-expansivo). Se se exige estrita estabilidade (decrescimento estrito) então 0<h<1. Em resumo: o maior passo que ainda dá fator não-negativo e não-expansivo para o autovalor -1 é h=1; mas ressalto que a presença do modo \lambda=1 impede estabilidade global por Euler explícito.

Passo 13

Agora a parte c): escrever a equação de iteração do método de Euler implícito (Backward Euler). A regra geral do método implícito aplicada a um sistema \dfrac{dy}{dx}=f(y) é
y_{n+1}=y_n+h\,f(y_{n+1}).
Aplicando ao nosso sistema (com \alpha=\beta=1) obtemos, escrevendo índices n e n+1:
\begin{aligned} & y_{1}^{\,n+1}=y_{1}^{\,n}+h\,(1-y_{2}^{\,n+1})\,y_{1}^{\,n+1},\\[6pt] & y_{2}^{\,n+1}=y_{2}^{\,n}-h\,(1-y_{1}^{\,n+1})\,y_{2}^{\,n+1}. \end{aligned}]

Passo 14

Manipular essas equações para expor as formas de solução implícita (isolando y_{1}^{\,n+1} e y_{2}^{\,n+1} em termos um do outro). Desenvolvendo os termos que envolvem y^{n+1}:
\begin{aligned} & y_{1}^{\,n+1}=y_{1}^{\,n}+h\,y_{1}^{\,n+1}-h\,y_{1}^{\,n+1}y_{2}^{\,n+1},\\[6pt] & y_{2}^{\,n+1}=y_{2}^{\,n}-h\,y_{2}^{\,n+1}+h\,y_{1}^{\,n+1}y_{2}^{\,n+1}. \end{aligned}
Agrupando y^{n+1} em cada equação:
\begin{aligned} & y_{1}^{\,n+1}\big(1-h+h\,y_{2}^{\,n+1}\big)=y_{1}^{\,n},\\[6pt] & y_{2}^{\,n+1}\big(1+h-h\,y_{1}^{\,n+1}\big)=y_{2}^{\,n}. \end{aligned}
Portanto as fórmulas implícitas equivalentes são
\boxed{\,y_{1}^{\,n+1}=\dfrac{y_{1}^{\,n}}{1-h+h\,y_{2}^{\,n+1}}\,,\qquad y_{2}^{\,n+1}=\dfrac{y_{2}^{\,n}}{1+h-h\,y_{1}^{\,n+1}}\,}
que constituem um sistema acoplado de equações não-lineares para (y_{1}^{\,n+1},y_{2}^{\,n+1}). Na prática resolve-se esse sistema por iteração (por exemplo método de Newton ou ponto fixo alternado) a cada passo.]

Resposta

(0,0) — sela (instável).
(1,0) — autovalores \{1,0\}, instável (não-hiperbólico).
(0,1) — autovalores \{0,-1\}, não-hiperbólico (uma direção estável e uma neutra).
(1,1) — autovalores \{\pm i\}, centro não-linear (órbitas periódicas; estável no sentido de Lyapunov, não assintoticamente estável).
h_{\max}=1 (para a componente associada a \lambda=-1 garantir fator 0\le 1+h\lambda\le 1; observar que a presença de \lambda=1 impede estabilidade global pelo Euler explícito).
Equações de iteração do Euler implícito:
y_{1}^{\,n+1}=\dfrac{y_{1}^{\,n}}{1-h+h\,y_{2}^{\,n+1}},\qquad y_{2}^{\,n+1}=\dfrac{y_{2}^{\,n}}{1+h-h\,y_{1}^{\,n+1}}.

Exercícios relacionados