Exibir código-fonte para Análise de estabilidade de equações diferenciais lineares atrasadas

{{Ecologia| [[Probabilidade básica]] |[[Linearização de sistemas de equações não lineares]]}}

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:

<math display="block">\frac{dy}{dt}=\alpha y_{\tau}</math>

Onde <math display="inline">y_{t}\equiv y\left(t-1\right)</math>, só há um único ponto de equilíbrio em <math display="inline">y=0</math>. Para equações diferenciais ordinárias do tipo:

<math display="block">a_{n}y^{\left(n\right)}+a_{n-1}y^{\left(n-1\right)}+\dots+a_{1}y^{'}+a_{0}y=0</math>

Assumindo que as soluções vão ser da forma <math display="inline">y\left(t\right)=Ce^{\lambda t}</math>, pode-se substituir:

<math display="block">\begin{align}
Ce^{\lambda t}\left(a_{n}\lambda^{n}+a_{n-1}\lambda^{n-1}+\dots+a_{1}\lambda+a_{0}\right) & =0\\
a_{n}\lambda^{n}+a_{n-1}\lambda^{n-1}+\dots+a_{1}\lambda+a_{0}\end{align}</math>

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

<math display="block">y'-\alpha y=0\rightarrow\lambda=\alpha</math>

Logo <math display="inline">y\left(t\right)=Ce^{\alpha t}</math>. Esta solução pode ser conferida resolvendo diretamente esta EDO simples:

<math display="block">\int\frac{dy}{y}=\int\alpha dt\rightarrow y=Ce^{\alpha t}</math>

Agora supondo uma solução análoga <math display="inline">y\left(t\right)=Ce^{\lambda t}</math> para a equação com atraso:

<math display="block">\frac{dy}{dt}=\frac{d}{dt}\left(Ce^{\lambda t}\right)=\lambda Ce^{\lambda t}</math>E<math display="block">y\left(t-1\right)=Ce^{\lambda\left(t-1\right)}=Ce^{\lambda t}e^{-\lambda}</math>

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

<math display="block">\begin{align}
y'-\alpha y_{\tau} & =0\\
\lambda Ce^{\lambda t}-\alpha Ce^{\lambda t}e^{-\lambda} & =0\\
\lambda-\alpha e^{-\lambda} & =0\end{align}</math>

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 <math display="inline">z_{1}=\lambda</math> e <math display="inline">z_{2}=\alpha e^{-\lambda}</math> separadamente, e então e procurar por intersecções. Isto é, quando os dois termos possuem o mesmo valor para um mesmo <math display="inline">\lambda</math>, pois consequentemente então <math display="inline">\lambda-\alpha e^{-\lambda}=0</math>. Se <math display="inline">\alpha>0</math>, há uma única intersecção, onde <math display="inline">\lambda>0</math>, então <math display="inline">y=Ce^{\lambda t}</math> aumenta exponencialmente ao infinito quando <math display="inline">t\rightarrow\infty</math>. Desta forma, o ponto de equilíbrio é instável.
[[Ficheiro:delay1.png|esquerda|miniaturadaimagem|A esquerda o gráfico gerado para um valor qualquer em que <math>\alpha>0</math>, e a direita para três valores em que <math>\alpha <0</math>, especificamente <math>\alpha_c<\alpha<0</math> em preto, <math>\alpha = \alpha_c</math> em verde, e <math>\alpha < \alpha_c</math> em azul.]]


Se <math display="inline">\alpha<0</math> pode haver 2 intersecções <math display="inline">\left(\alpha_{c}<\alpha\right)</math>, 1 intersecção o <math display="inline">\left(\alpha_{c}=\alpha\right)</math> ou nenhuma intersecção <math display="inline">\left(\alpha<\alpha_{c}\right)</math>. Para identificar qual é este ponto crítico <math>\alpha_c</math>, basta perceber que neste ponto a reta <math display="inline">z_{1}</math> é tangente à curva <math display="inline">z_{2}</math> no ponto <math display="inline">\left(\lambda_{c},z_{c}\right)</math>. Logo, a inclinação de ambos deve ser a mesma, <math display="inline">m_{1}=m_{2}</math>. A inclinação da reta <math display="inline">z_{1}</math> é simplesmente <math display="inline">m_{1}=1</math>. Então a inclinação da curva também deve ser:<math display="block">m_{2}=\frac{dz_{2}}{d\lambda}|\lambda_{c}=-\alpha_{c}e^{-\lambda_{c}}=1</math><math display="block">\alpha_{c}=-e^{\lambda_{c}}</math>


Substituindo a constante <math>\alpha_c</math> em <math display="inline">z_{2}</math>:

<math display="block">z_{c}=-e^{\lambda_{c}}e^{-\lambda_{c}}=-1</math>

Pode-se obter agora <math display="inline">\lambda_{c}</math> a partir da reta <math display="inline">z_{1}</math>, <math display="inline">z_{c}=-1=\lambda_{c}</math>. Dessa forma o valor crítico é então <math display="inline">\alpha_{c}=-e^{-1}</math>. Logo, se <math display="inline">\alpha\in\left[\alpha_{c}0\right]</math>, então as raízes da equação característica são reais e negativas e a solução exponencial associada decai para <math display="inline">0</math> com o tempo. As raízes <math display="inline">\lambda</math> também são chamadas de autovalores, termo que será empregado daqui em diante.
[[Ficheiro:delay2.png|esquerda|miniaturadaimagem|A esquerda a solução para <math display="inline">\alpha=1>0</math> e a direita para <math display="inline">\alpha_c \leq \alpha=-1<0</math>, considerando como condição inicial <math display="inline">y_{0}=2</math>.]]


Para <math display="inline">\alpha<\alpha_{c}</math>, 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 [https://www.wolfram.com/mathematica/index.html.pt-br?footer=lang 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 <math display="inline">\lambda=\lambda_{r}+i\lambda_{i}</math>, na equação característica, obtém-se:

<math display="block">\begin{align}
\lambda-\alpha e^{-\lambda} & =0\\
\lambda & =\alpha e^{-\lambda}\\
\lambda_{r}+i\lambda_{i} & =\alpha e^{\left[-\lambda_{r}-i\lambda_{i}\right]}\\
\lambda_{r}+i\lambda_{i} & =\alpha e^{-\lambda_{r}}e^{-i\lambda_{i}}\\
\left[\lambda_{r}\right]+i\left[\lambda_{i}\right] & =\left[\alpha e^{-\lambda_{r}}\cos\lambda_{i}\right]+i\left[\alpha e^{-\lambda_{r}}\sin\left(-\lambda_{i}\right)\right]\end{align}</math>

Logo:

<math display="block">\begin{align}
\lambda_{r} & =+\alpha e^{-\lambda_{r}}\cos\lambda_{i}\\
\lambda_{i} & =-\alpha e^{-\lambda_{r}}\sin\lambda_{i}\end{align}</math>

Calculando a razão entre os termos:

<math display="block">\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)</math>

Então utilizando <math display="inline">\lambda_{i}</math> como parâmetro, pode-se obter equações paramétricas para <math display="inline">\lambda_{r}</math> e <math display="inline">\alpha</math>:

<math display="block">\begin{align}
\lambda_{r} & =-\cot\left(\lambda_{i}\right)\lambda_{i}\\
\alpha & =-\frac{\lambda_{i}}{\exp\left(\lambda_{i}\cot\left(\lambda_{i}\right)\right)\sin\lambda_{i}}\end{align}</math>
[[Ficheiro:delay3.png|miniaturadaimagem|Soluções das equações paramétricas. Em vermelho o caso especial quando <math>\lambda_i=0</math>.]]
Por conta da periodicidade das funções trigonométricas, muitas curva são traçadas quando varía-se <math display="inline">\lambda_{i}</math> entre <math display="inline">-\infty</math> e <math display="inline">+\infty</math>. Além, disto para <math display="inline">\lambda=\lambda_{r}</math> então <math display="inline">\lambda_{i}=0</math>, desta forma pode-se plotar diretamente <math display="inline">\alpha\left(\lambda_{r}\right)</math>:

<math display="block">\alpha=\lambda_{r}e^{\lambda_{r}}</math>

<div class="list">

<span> Utilizando o [https://www.geogebra.org/ Geogebra], isto pode ser feito realizando cada entrada manualmente:</span>
 -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 é <math display="inline">y\left(t;a\right)=\sum_{c}e^{\lambda_{n}t}</math>, onde o somatório é sobre todos os valores de <math display="inline">\lambda_{n}</math> para um dado parâmetro <math display="inline">\alpha</math>, então o estado de equilíbrio é estável para valores de <math display="inline">\alpha</math> em que todos os autovalores tem valores reais negativos. Ou seja, para valores entre  <math display="inline">\alpha_{c2}</math> e o eixo <math display="inline">\alpha=0</math>. pois para <math display="inline">\alpha>0</math>, então <math display="inline">\lambda_{r}>0</math> quando <math display="inline">\lambda_{i}=0</math>. Então o próximo passo é identificar  segundo pronto crítico <math display="inline">a_{c2}</math>.
</div>Ele ocorre quando <math display="inline">\lambda_{r}=0</math> e <math display="inline">a\neq0</math> e <math display="inline">\lambda_{i}\neq0</math>. Então substituindo:

<math display="block">\begin{align}
\lambda_{r}' & =\alpha'e^{-\lambda'_{r}}\cos\lambda_{i}'=0\end{align}</math>Logo
<math display="block">\lambda_{i}'=\left(\pm\frac{\pi}{2}+n\pi\right),n\in\mathbb{Z}</math>


E substituindo:<math display="block">\begin{align}
\lambda_{i}' & =-\alpha'e^{-\lambda'_{r}}\sin\lambda_{i}'\\
\left(\pm\frac{\pi}{2}+n\pi\right) & =-\alpha'\sin\left(\pm\frac{\pi}{2}+n\pi\right)\end{align}</math>


Usando a propriedade <math display="inline">\sin\left(a+b\right)=\sin a\cos b+\sin b\cos a</math> então:

<math display="block">\begin{align}
\sin\left(\pm\frac{\pi}{2}+n\pi\right) & =\sin\left(\pm\frac{\pi}{2}\right)\cos\left(n\pi\right)+\sin\left(n\pi\right)\cos\left(\pm\frac{\pi}{2}\right)\\
 & =\pm\sin\left(\frac{\pi}{2}\right)\cos\left(n\pi\right)\\
 & \pm\left(-1\right)^{n}\end{align}</math>Então:<math display="block">\begin{align}
\left(\pm\frac{\pi}{2}+n\pi\right) & =\mp\left(-1\right)^{n}\alpha'\\
\alpha'= & \left(-1\right)^{n}\left(\mp n\pi-\frac{\pi}{2}\right)\end{align}</math>Então estes são os valores de <math display="inline">\alpha</math> em que temos <math display="inline">\lambda_{r}=0</math>. Novamente pode ser visualizado via Geogebra:
 ((-1)^n (-n π - π / 2),0)
 ((-1)^n (+n π - π / 2),0)
Para <math display="inline">n=0</math>, obtém-se então o ponto crítico desejado:<math display="block">\alpha_{c2}=-\frac{\pi}{2}</math>

O primeiro ponto no eixo negativo de <math display="inline">\alpha</math> em que tem-se <math display="inline">\lambda_{r}</math>. Logo, o ponto de equilíbrio é estável se <math display="inline">\alpha\in\left(-\frac{\pi}{2},0\right)</math>, 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.
[[Ficheiro:delay4.png|esquerda|miniaturadaimagem|A esquerda a solução para <math display="inline">\alpha=-2< \alpha_{c2}</math> e a direita para <math display="inline">\alpha_{c2}< \alpha=-2< \alpha_c<0</math>, considerando como condição inicial <math display="inline">y_{0}=2</math>.]]
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 <math display="inline">\left(a_{c},\lambda\right)=\left(\frac{1}{e},-1\right)</math>, este resultado concorda com a equação obtida para <math>\lambda_i=0</math>:<math display="block">\alpha=\lambda_{r}e^{\lambda_{r}}=-e^{-1}</math>

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

<math display="block">\begin{align}
\lim_{\lambda_{i}\rightarrow0}\lambda_{r} & =-\lim_{\lambda_{i}\rightarrow0}\cot\left(\lambda_{i}\right)\lambda_{i}=-1\\
\lim_{\lambda_{i}\rightarrow0}\alpha & =\lim_{\lambda_{i}\rightarrow0}-\frac{\lambda_{i}}{\exp\left(\lambda_{i}\cot\left(\lambda_{i}\right)\right)\sin\lambda_{i}}=-\frac{1}{e}\end{align}</math>No Mathematica:
 {Limit[-Cot[x]*x, x -> 0], Limit[-x/(Exp[x*Cot[x]]*Sin[x]), x -> 0]}

====== Principais materiais utilizados ======
#[https://www.math.fsu.edu/~bertram/lectures/delay.pdf Delay-Differential Equations] (Richard Bertram, Universidade Estadual da Flórida)
#[https://tutorial.math.lamar.edu/classes/de/hohomogeneousde.aspx Homogeneous Differential Equations]  (Paul Dawkins, Universidade de Lamar)


{{Ecologia| [[Probabilidade básica]] |[[Linearização de sistemas de equações não lineares]]}}