Análise de estabilidade de equações diferenciais lineares atrasadas

Anterior: Establidade em equações diferenciais com atraso | Índice: Ecologia | Próximo: Exemplo: sistema de equações diferenciais com retardo (SIRS)

Analisar a estabilidade local de equações diferenciais atrasadas é mais desafiador que realizar a mesma análise para equações diferenciais ordinárias. Isto ocorre devido a dimensionalidade infinita do sistema. Por exemplo considerando uma equação diferencial linear atrasada simples:

$$\frac{dy}{dt}=\alpha y_{\tau }$$

Onde $y_t\equiv y(t-1)$, só há um único ponto de equilíbrio em $y=0$. Para equações diferenciais ordinárias do tipo:

$$a_n{y}^{\left(n\right)}+a_{n-1}{y}^{(n-1)}+\dots +a_1{y}^{\text{'}}+a_{0}y=0$$

Assumindo que as soluções vão ser da forma $y\left(t\right)=C{e}^{\lambda t}$, pode-se substituir:

$$\begin{array}{rl}C{e}^{\lambda t}(a_n{\lambda }^n+a_{n-1}{\lambda }^{n-1}+\dots +a_1\lambda +a_0)& =0\\ a_n{\lambda }^n+a_{n-1}{\lambda }^{n-1}+\dots +a_1\lambda +a_0\end{array}$$

Então este é o polinômio característico. E sua solução nos dá a solução para a equação diferencial. Por exemplo:

$$y^{\prime }-\alpha y=0\to \lambda =\alpha$$

Logo $y\left(t\right)=C{e}^{\alpha t}$. Esta solução pode ser conferida resolvendo diretamente esta EDO simples:

$$\int \frac{dy}{y}=\int \alpha dt\to y=C{e}^{\alpha t}$$

Agora supondo uma solução análoga $y\left(t\right)=C{e}^{\lambda t}$ para a equação com atraso:

$$\frac{dy}{dt}=\frac{d}{dt}\left(C{e}^{\lambda t}\right)=\lambda C{e}^{\lambda t}$$E$$y(t-1)=C{e}^{\lambda (t-1)}=C{e}^{\lambda t}{e}^{-\lambda }$$

Então a equação característica é:

$$\begin{array}{rl}{y}^{\prime }-\alpha y_{\tau }& =0\\ \lambda C{e}^{\lambda t}-\alpha C{e}^{\lambda t}{e}^{-\lambda }& =0\\ \lambda -\alpha e^{-\lambda }& =0\end{array}$$

Para as equações diferenciais ordinárias havia um polinômio e o teorema fundamental da álgebra permitia dizer quantas raízes esperar. Porém para as equações diferenciais atrasadas, não há teorema algum sobre a quantidade de raízes, este número poderia inclusive ser infinito. Como exemplo, pode-se considerar soluções reais e complexas separadamente conforme será visto na sequência.

Solução real

Supondo que a solução é real, pode-se plotar $z_1=\lambda$ e $z_2=\alpha e^{-\lambda }$ separadamente, e então e procurar por intersecções. Isto é, quando os dois termos possuem o mesmo valor para um mesmo $\lambda$, pois consequentemente então $\lambda -\alpha e^{-\lambda }=0$. Se $\alpha >0$, há uma única intersecção, onde $\lambda >0$, então $y=C{e}^{\lambda t}$ aumenta exponencialmente ao infinito quando $t\to \mathrm{\infty }$. Desta forma, o ponto de equilíbrio é instável.

A esquerda o gráfico gerado para um valor qualquer em que $\alpha >0$, e a direita para três valores em que $\alpha <0$, especificamente ${\alpha }_c<\alpha <0$ em preto, $\alpha ={\alpha }_c$ em verde, e $\alpha <{\alpha }_c$ em azul.

Se $\alpha <0$ pode haver 2 intersecções $({\alpha }_c<\alpha )$, 1 intersecção o $({\alpha }_c=\alpha )$ ou nenhuma intersecção $(\alpha <{\alpha }_c)$. Para identificar qual é este ponto crítico ${\alpha }_c$, basta perceber que neste ponto a reta $z_1$ é tangente à curva $z_2$ no ponto $({\lambda }_c,z_c)$. Logo, a inclinação de ambos deve ser a mesma, $m_1=m_2$. A inclinação da reta $z_1$ é simplesmente $m_1=1$. Então a inclinação da curva também deve ser:$$m_2=\frac{d{z}_2}{d\lambda }|{\lambda }_c=-{\alpha }_c{e}^{-{\lambda }_c}=1$$$${\alpha }_c=-e^{{\lambda }_c}$$

Substituindo a constante ${\alpha }_c$ em $z_2$:

$$z_c=-e^{{\lambda }_c}{e}^{-{\lambda }_c}=-1$$

Pode-se obter agora ${\lambda }_c$ a partir da reta $z_1$, $z_c=-1={\lambda }_c$. Dessa forma o valor crítico é então ${\alpha }_c=-e^{-1}$. Logo, se $\alpha \in \left[{\alpha }_{c}0\right]$, então as raízes da equação característica são reais e negativas e a solução exponencial associada decai para $0$ com o tempo. As raízes $\lambda$ também são chamadas de autovalores, termo que será empregado daqui em diante.

A esquerda a solução para $\alpha =1>0$ e a direita para ${\alpha }_c\le \alpha =-1<0$, considerando como condição inicial $y_0=2$.

Para $\alpha <{\alpha }_c$, não há decaimento ou crescimento exponencial nos componentes da solução, para isto vamos analisar as soluções complexas. As soluções podem ser obtidas numericamente via Mathematica:

sol = NDSolve[{y'[t] == -2.0*y[t - 1], y[t /; t <= 1] == 2}, y, {t, 0, 30];
Plot[y[t] /. sol, {t, 0, 30}, PlotRange -> All]

Solução Complexa

Substituindo então $\lambda ={\lambda }_r+i{\lambda }_i$, na equação característica, obtém-se:

$$\begin{array}{rl}\lambda -\alpha e^{-\lambda }& =0\\ \lambda & =\alpha e^{-\lambda }\\ {\lambda }_r+i{\lambda }_i& =\alpha e^{[-{\lambda }_r-i{\lambda }_i]}\\ {\lambda }_r+i{\lambda }_i& =\alpha e^{-{\lambda }_r}{e}^{-i{\lambda }_i}\\ \left[{\lambda }_r\right]+i\left[{\lambda }_i\right]& =[\alpha e^{-{\lambda }_r}\cos {\lambda }_i]+i[\alpha e^{-{\lambda }_r}\sin (-{\lambda }_i)]\end{array}$$

Logo:

$$\begin{array}{rl}{\lambda }_r& =+\alpha e^{-{\lambda }_r}\cos {\lambda }_i\\ {\lambda }_i& =-\alpha e^{-{\lambda }_r}\sin {\lambda }_i\end{array}$$

Calculando a razão entre os termos:

$$\frac{{\lambda }_r}{{\lambda }_i}=\frac{+\alpha e^{-{\lambda }_r}\cos {\lambda }_i}{-\alpha e^{-{\lambda }_r}\sin {\lambda }_i}=-\cot \left({\lambda }_i\right)$$

Então utilizando ${\lambda }_i$ como parâmetro, pode-se obter equações paramétricas para ${\lambda }_r$ e $\alpha$:

$$\begin{array}{rl}{\lambda }_r& =-\cot \left({\lambda }_i\right){\lambda }_i\\ \alpha & =-\frac{{\lambda }_i}{\exp ({\lambda }_i\cot \left({\lambda }_i\right))\sin {\lambda }_i}\end{array}$$

Soluções das equações paramétricas. Em vermelho o caso especial quando ${\lambda }_i=0$.

Por conta da periodicidade das funções trigonométricas, muitas curva são traçadas quando varía-se ${\lambda }_i$ entre $-\mathrm{\infty }$ e $+\mathrm{\infty }$. Além, disto para $\lambda ={\lambda }_r$ então ${\lambda }_i=0$, desta forma pode-se plotar diretamente $\alpha \left({\lambda }_r\right)$:

$$\alpha ={\lambda }_r{e}^{{\lambda }_r}$$

Utilizando o Geogebra, isto pode ser feito realizando cada entrada manualmente:

-x cotg(x)
d(x) = (-x)/(e^(x cotg(x)) sen(x))
Curva(d(u), c(u), u, -60 , 60)
x = y e^y

A solução geral é $y(t;a)=\sum _c{e}^{{\lambda }_{n}t}$, onde o somatório é sobre todos os valores de ${\lambda }_n$ para um dado parâmetro $\alpha$, então o estado de equilíbrio é estável para valores de $\alpha$ em que todos os autovalores tem valores reais negativos. Ou seja, para valores entre ${\alpha }_{c2}$ e o eixo $\alpha =0$. pois para $\alpha >0$, então ${\lambda }_r>0$ quando ${\lambda }_i=0$. Então o próximo passo é identificar segundo pronto crítico $a_{c2}$.

Ele ocorre quando

${\lambda }_r=0$

e

$a\ne 0$

e

${\lambda }_i\ne 0$

. Então substituindo:

$$\begin{array}{rl}{\lambda }_{r^{\prime }}& ={\alpha }^{\prime }{e}^{-\lambda {\text{'}}_r}\cos {\lambda }_{i^{\prime }}=0\end{array}$$Logo $${\lambda }_{i^{\prime }}=(\pm \frac{\pi }{2}+n\pi ),n\in \mathbb{\mathbb{Z}}$$

E substituindo:$$\begin{array}{rl}{\lambda }_{i^{\prime }}& =-{\alpha }^{\prime }{e}^{-\lambda {\text{'}}_r}\sin {\lambda }_{i^{\prime }}\\ (\pm \frac{\pi }{2}+n\pi )& =-{\alpha }^{\prime }\sin (\pm \frac{\pi }{2}+n\pi )\end{array}$$

Usando a propriedade $\sin (a+b)=\sin a\cos b+\sin b\cos a$ então:

$$\begin{array}{rl}\sin (\pm \frac{\pi }{2}+n\pi )& =\sin (\pm \frac{\pi }{2})\cos \left(n\pi \right)+\sin \left(n\pi \right)\cos (\pm \frac{\pi }{2})\\ & =\pm \sin \left(\frac{\pi }{2}\right)\cos \left(n\pi \right)\\ & \pm {(-1)}^n\end{array}$$Então:$$\begin{array}{rl}(\pm \frac{\pi }{2}+n\pi )& =\mp {(-1)}^n{\alpha }^{\prime }\\ {\alpha }^{\prime }=& {(-1)}^n(\mp n\pi -\frac{\pi }{2})\end{array}$$Então estes são os valores de $\alpha$ em que temos ${\lambda }_r=0$. Novamente pode ser visualizado via Geogebra:

((-1)^n (-n π - π / 2),0)
((-1)^n (+n π - π / 2),0)

Para $n=0$, obtém-se então o ponto crítico desejado:$${\alpha }_{c2}=-\frac{\pi }{2}$$

O primeiro ponto no eixo negativo de $\alpha$ em que tem-se ${\lambda }_r$. Logo, o ponto de equilíbrio é estável se $\alpha \in (-\frac{\pi }{2},0)$, se não, é instável. Pode-se observador que os resultados obtidos para o caso em que os autovalores eram apenas reais está contido neste resultado.

A esquerda a solução para $\alpha =-2<{\alpha }_{c2}$ e a direita para ${\alpha }_{c2}<\alpha =-2<{\alpha }_c<0$, considerando como condição inicial $y_0=2$.

Antes de concluir, pode-se prestar uma atenção especial para o ponto em que as soluções complexas encontram a solução real, este é exatamente o primeiro ponto crítico. Lembrando que o primeiro ponto crítico era $(a_c,\lambda )=(\frac{1}{e},-1)$, este resultado concorda com a equação obtida para ${\lambda }_i=0$:$$\alpha ={\lambda }_r{e}^{{\lambda }_r}=-e^{-1}$$

E também no limite das equações paramétricas:

$$\begin{array}{rl}\underset{{\lambda }_i\to 0}{\lim }{\lambda }_r& =-\underset{{\lambda }_i\to 0}{\lim }\cot \left({\lambda }_i\right){\lambda }_i=-1\\ \underset{{\lambda }_i\to 0}{\lim }\alpha & =\underset{{\lambda }_i\to 0}{\lim }-\frac{{\lambda }_i}{\exp ({\lambda }_i\cot \left({\lambda }_i\right))\sin {\lambda }_i}=-\frac{1}{e}\end{array}$$No Mathematica:

{Limit[-Cot[x]*x, x -> 0], Limit[-x/(Exp[x*Cot[x]]*Sin[x]), x -> 0]}

Principais materiais utilizados

1. Delay-Differential Equations (Richard Bertram, Universidade Estadual da Flórida)
2. Homogeneous Differential Equations (Paul Dawkins, Universidade de Lamar)

Anterior: Establidade em equações diferenciais com atraso | Índice: Ecologia | Próximo: Exemplo: sistema de equações diferenciais com retardo (SIRS)