3.1.Analise Qualitativa de EDOs-Conceitos Basicos.sagews

ANÁLISE QUALITATIVA DE EDOs: CONCEITOS BÁSICOS

Suponha que precisamos achar a solução de uma EDO de ordem 1, $dy/dx=f(x,y)$ mas não somos capazes de descobrir ou inventar um método para resolvê-la analiticamente. À primeira vista parece grave mas, de fato, não é; isso porque a própria EDO pode nos dizer como é o comportamento da solução mesmo que não a tenhamos explicitamente. Em outras palavras, podemos analisar a EDO qualitativamente a partir do esboço gráfico de suas soluções mesmo sem resolvê-la.

Campos de Direção

Veremos que, quando a função $f(x,y)$ na EDO depende apenas de $y$, podemos responder diversas perguntas qualitativas sobre as propriedades da solução: como a solução se comporta na vizinhança de um dado ponto? como ela se comporta quando $x\rightarrow\infty$?, etc Mas antes, relembremos um conceito do Cálculo Diferencial (e mostrado no capítulo anterior): a derivada $dy/dx$ de uma função diferenciável $y=y(x)$ fornece a inclinação das linhas tangentes aos pontos do gráfico dessa função.

Inclinação (Slope)

Como a solução $y=y(x)$ de uma equação diferencial de ordem 1 (expressa na forma normal) $$\frac{dy}{dx}=f(x,y)$$ é necessariamente diferenciável no intervalo de definição $I$, ela também é contínua nesse intervalo; logo, para cada ponto $(x,y(x))$, existe uma linha tangente ao gráfico. A função $f$ que aparece na forma normal é chamada função inclinação. A inclinação da linha tangente à curva que representa a solução, no ponto $(x,y(x))$ é o valor da primeira derivada $dy/dx$ nesse ponto que, pela EDO, é justamente o valor da função inclinação $f(x,y(x))$
Agora, suponha que $(x,y)$ representa um ponto qualquer em uma região do plano $xy$ na qual a função $f$ está definida. Então, o valor $f(x,y)$ correspondente a esse ponto representa a inclinação de uma linha ou de um segmento de linha. Como exemplo concreto, considere a equação $dy/dx=0.2xy$, na qual $f(x,y)=0.2xy$. Tome o ponto $(2,3)$; a inclinação de um segmento nesse ponto é $f(2,3)=0.2(2)(3)=1.2$ A figura abaixo mostra um segmento de reta com inclinação 1.2 passando pelo ponto (2,3). A figura mostra também que, se uma curva solução também passa pelo mesmo ponto, ela tangencia o segmento de reta, ou seja, o segmento de reta é uma mini linha tangente no ponto em questão. Para o traçado da curva, usamos o fato de que a solução geral da EDO em questão é $y(x)=c e^{x^2/10}$ (em breve aprenderemos como obter essa solução).

# Definindo variáveis do problema
x = var('x')
c = var('c')

# Valor da constante c para o ponto (2,3)
c = 3/exp(2^2/10)

# Função y(x): solução particular da EDO para a constante c
y = c*exp(x^2/10)

# Reta tangente à curva y(x) no ponto (2,3)
rt = 1.2*x+0.6

# Gráficos
P = Graphics()
P += plot(y, x, 0, 4, ymin=0, ymax=5)
P += point((2,3),size=40)
P += plot(rt, x, 1.4, 2.6, ymin=0, ymax=5)
P.show()

Campo de Direção

Se, agora, estendermos o procedimento descrito acima, ou seja, se avaliarmos $f$ para todos os pontos de uma grade retangular no plano-$xy$ e desenharmos um segmento de reta em cada um deles, com a inclinação $f(x,y)$, então a coleção de todos esses segmentos é chamada campo de direções da equação diferencial $dy/dx=f(x,y)$. Visualmente, o campo de direções indica a forma de uma família de curvas solução da EDO e, com isso, pode-se ter uma idéia do comportamento qualitativo das soluções. Se olharmos para uma única curva, percebemos que ela segue o padrão de fluxo do campo de direções. Vejamos alguns exemplos para esclarecer estes fatos. O campo de direções pode ser esboçado manualmente (Euler já fazia isso!) mas é muito mais fácil se usarmos recursos computacionais. O Sage possui algumas funções para isso. Primeiro, vamos exibir o campo de inclinações somente com segmentos de reta (somente as inclinações). Observe que, para tanto, a única informação necessária é a função $f(x,y)$.

Exemplo 1

Considere a EDO $\displaystyle\frac{dy}{dx}=0.2xy$. Usando o Sage, traçamos o campo de direções e sobrepomos algumas soluções. Examine a figura e comprove o que afirmamos até aqui.

# Define variável independente e constante de integração
x=var('x')
c=var('c')

# Define variável dependente
y=function('y')(x)

# Solução geral da EDO
sol = c*exp(0.1*x^2)

# Função f(x,y) da EDO
f(x, y) = 0.2*x*y

# Gráficos
P = Graphics()
P += plot_slope_field(f,[x,-pi,pi],[y,-4,4])                           # Campo de Inclinações
P += plot(sol.substitute(c==3), xmin=-pi, xmax=pi, ymin=-4, ymax=4)    # Soluções Particulares (diferentes valores da constante de integração)
P += plot(sol.substitute(c==1.5), xmin=-pi, xmax=pi, ymin=-4, ymax=4)
P += plot(sol.substitute(c==-3), xmin=-pi, xmax=pi, ymin=-4, ymax=4)
P += plot(sol.substitute(c==-1.5), xmin=-pi, xmax=pi, ymin=-4, ymax=4)

show(P, figsize=(6, 4))

Como curiosidade, existe uma maneira alternativa de exibir o campo de direções, entendendo-o como um campo vetorial. Neste gráfico, os segmentos de reta são apresentados como setas (enfatizando, portanto, o sentido do campo) e o comprimento das setas é proporcional à magnitude do vetor em cada ponto. Observe que são fornecidos, como argumentos, dois componentes do vetor, $x'=1$ e $y'=0.2xy$. Esta forma de apresentação será utilizada futuramente quando estudarmos sistemas de equações diferenciais. Então, aprenderemos mais detalhes.

# Idêntico ao exemplo anterior a menos do gráfico do campo de direções ao invés do campo de inclinações
x=var('x')
c=var('c')

y=function('y')(x)
sol = c*exp(0.1*x^2)

f(x, y) = 0.2*x*y
P = Graphics()
P += plot_vector_field([1, 0.2*x*y],[x,-pi,pi],[y,-4,4])                # Campo de Direções
P += plot(sol.substitute(c==3), xmin=-pi, xmax=pi, ymin=-4, ymax=4)
P += plot(sol.substitute(c==1.5), xmin=-pi, xmax=pi, ymin=-4, ymax=4)
P += plot(sol.substitute(c==-3), xmin=-pi, xmax=pi, ymin=-4, ymax=4)
P += plot(sol.substitute(c==-1.5), xmin=-pi, xmax=pi, ymin=-4, ymax=4)
show(P)

Exemplo 2

Considere a EDO $\displaystyle\frac{dy}{dx}=y-x$. Esboce o campo de direções. Aproveitemos este exemplo para apresentar uma Receita Básica de resolução analítica de EDOs de Ordem 1, usando o SAGE. Basta seguir o código fornecido abaixo. Neste procedimento, determinamos a solução geral da EDO, identificamos a constante arbitrária e atribuímos valores a ela, obtendo algumas soluções particulares que são exibidas graficamente.

# Define as variáveis independente e dependente
x = var('x')
y = function('y')(x)

# Define a equação diferencial na forma dy/dx = f(x,y): o comando diff(y,x) é o comando SAGE para a derivada de ordem 1
edo = diff(y, x) == y - x

# Resolve a EDO usando o comando desolve do SAGE e armazena a expressão na variável sol
sol = desolve(edo, y)

# Exibe a solução geral da EDO
sol.show()

# Identifica as variáveis contidas na solução geral
# Neste caso, são a constante _C e a variável independente x, exibidas pelo SAGE nesta ordem
show(sol.variables())

# Obtem uma referência para a variável _C (que é a primeira, com índice 0)
c=sol.variables()[0]
show(c)

# Define a função f(x,y) da EDO para exibir o campo de inclinações
f = function('f')(x,y)
f(x,y) = y - x

# Obtem algumas soluções particulares para graficar
# Usa o comando substitute do SAGE para substituir os diferentes valores da constante C
sol_particular1 = sol.substitute(c==3)
sol_particular2 = sol.substitute(c==0)
sol_particular3 = sol.substitute(c==-1)

P = Graphics()
P += plot_slope_field(f,[x,-4,4],[y,-4,4])                                   # Exibe o campo de inclinações
P += plot(sol_particular1, xmin=-4, xmax=4, ymin=-4, ymax=4, color='red')    # Exibe a solução particular 1
P += plot(sol_particular2, xmin=-4, xmax=4, ymin=-4, ymax=4, color='green')  # Exibe a solução particular 2
P += plot(sol_particular3, xmin=-4, xmax=4, ymin=-4, ymax=4, color='blue')   # Exibe a solução particular 3
show(P)

$\displaystyle {\left({\left(x + 1\right)} e^{\left(-x\right)} + C\right)} e^{x}$

($\displaystyle C$, $\displaystyle x$)

$\displaystyle C$

Existe outra maneira (equivalente) de resolver uma EDO usando SAGE. Neste caso, ao invés de atribuirmos valores à constante arbitrária de integração, $C$, impomos diferentes condições iniciais que, como sabemos, equivale a determinar o valor de $C$ na expressão da solução geral da EDO. Portanto, estamos resolvendo diferentes Problemas a Valor Inicial.

# Define as variáveis independente e dependente
x = var('x')
y = function('y')(x)

# Define a equação diferencial na forma dy/dx = f(x,y): o comando diff(y,x) é o comando SAGE para a derivada de ordem 1
edo = diff(y, x) == y - x

# Resolve a EDO usando o comando desolve do SAGE passando o parâmetro ics, que é a condição inicial
# A condição inicial é fornecida através de ics=[x0,y0], correspondente a y(x0)=y0
sol_particular1 = desolve(edo, y, ics=[-4, -2.9])
sol_particular2 = desolve(edo, y, ics=[-4, -3])
sol_particular3 = desolve(edo, y, ics=[-4, -3.1])

# Define a função f(x,y) da EDO para exibir o campo de inclinações
f = function('f')(x,y)
f(x,y) = y - x

P = Graphics()
P += plot_slope_field(f,[x,-4,4],[y,-4,4])                                   # Exibe o campo de inclinações
P += plot(sol_particular1, xmin=-4, xmax=4, ymin=-4, ymax=4, color='red')    # Exibe a solução particular 1
P += plot(sol_particular2, xmin=-4, xmax=4, ymin=-4, ymax=4, color='green')  # Exibe a solução particular 2
P += plot(sol_particular3, xmin=-4, xmax=4, ymin=-4, ymax=4, color='blue')   # Exibe a solução particular 3
show(P)

EDO de Ordem 1 Autônoma

Outra maneira de classificar uma EDO é se ela é autônoma ou não autônoma. Uma EDO é dita autônoma se a variável independente não aparece explicitamente; caso contrário, ela é não autônoma. Portanto, a EDO é autônoma se for expressa como: $$f(y,y')=0\quad\quad ou\quad\quad \frac{dy}{dx}=f(y)$$ Assim, a equação $\frac{dy}{dx}=1+y^2$ é autônoma enquanto $\frac{dy}{dx}=0.2xy$ é não autônoma. Muitas equações diferenciais que modelam problemas reais são autônomas; já vimos vários exemplos como no decaimento radiativo, na propagação de doenças contagiosas, na lei de Newton do resfriamento, etc.

Pontos Críticos

Os zeros da função $f(y)$ são especialmente importantes. Dizemos que um número real $c$ é um ponto crítico da EDO autônoma se for um zero de $f$, ou seja, se $f(c)=0$. Um ponto crítico também é chamado de ponto de equilíbrio ou de ponto estacionário. Se substituirmos uma função constante $y(x)=c$ na forma normal da EDO ambos os lados da equação ficam iguais a zero. Isto significa que

Se $c$ é um ponto crítico da EDO, então $y(x)=c$ é uma solução constante da EDO autônoma.

Uma solução constante $y(x)=c$ é chamada solução de equilíbrio; as soluções de equilíbrio são as únicas soluções constantes. Nos problemas físicos reais, usualmente, a variável independente é o tempo $t$; assim, uma solução ser constante significa que ela não varia com o tempo, o que justifica a expressão solução de equilíbrio ou estacionária. Por outro lado, para uma solução não constante $y(x)$ podemos descobrir aspectos importantes do comportamento dessa solução, mais precisamente, se ela é crescente ou decrescente, simplesmente determinando o sinal algébrico da derivada $dy/dx$. Para isso, identificamos os intervalos do domínio para os quais $f(y)$ (que é igual a $dy/dx$) é positiva ou negativa, como se pode ver no exemplo a seguir.

Exemplo

A EDO $\displaystyle\frac{dP}{dt}=P(a-bP)$, onde $a$ e $b$ são constantes positivas é uma EDO autônoma uma vez que a variável independente, $t$, não aparece em $f(P)=P(a-bP)$. Fazendo $f(P)=P(a-bP)=0$, vemos que $0$ e $a/b$ são pontos críticos e as correspondentes soluções de equilíbrio são $P(t)=0$ e $P(t)=a/b$. Exibindo esses pontos críticos sobre uma linha horizontal (ou vertical), observamos que a linha fica subdividida em três intervalos distintos: $(-\infty,0)$, $(0,a/b)$ e $(a/b,\infty)$. As setas do gráfico indicam o sinal algébrico de $f(P)=P(a-bP)$ nesses intervalos: negativo para a esquerda e positivo para a direita. Lembrando que $f(P)=dP/dt$ e do que sabemos sobre o sinal das derivadas, concluímos que as setas mostram os intervalos em que a$f(P)$ é crescente ou decrescente. A tabela abaixo contém todas as informações que explicam o gráfico fornecido em seguida. Observe que, para traçar o gráfico, tomamos $a/b=2$, sem perda de generalidade. O gráfico é chamado Retrato de Fase (neste caso, unidimensional). Quando estudarmos EDOs de ordem superior e sistemas de EDOs encontraremos retratos de fase muito mais interessantes!

No Retrato de Fase acima, mostramos apenas o eixo $P$; se agora, porém, levarmos em consideração o eixo $t$, os três intervalos se transformam em três regiões do plano $tP$ definidas por:

• $R1: -\infty\lt P\lt 0$
• $R2: 0\lt P\lt a/b$
• $R3: a/b \lt P\lt \infty$

onde $-\infty \lt t \lt \infty$.
As três regiões são separadas pelas retas horizontais correspondentes às duas soluções de equilíbrio, $P=0$ e $P=a/b$.

# Gráfico da linha de fase (vertical)
PG11=plot([],figsize=(1,4),axes = False,xmin=-0.1,xmax=0.1)
#PG11=plot([],axes = False,figsize=(2,4),)
PG11 += line( [ (0, -4), (0, 4) ])
PG11 +=point((0,0),size=100)
PG11 +=point((0,2),size=100)
PG11 +=point((0,-1),size=100,marker='v',color="red")
PG11 +=point((0,1),size=100,marker='^',color="red")
PG11 +=point((0,3),size=100,marker='v',color="red")
#PG11 +=text("0",(0.01,0))
#T2=text("a/b",(-0.05,2))
#T3=text("Eixo P",(-0.05,3.8))

# Define as variáveis da EDO
# Tomamos valores particulares para os parâmetros a e b
t=var('t')
a=2
b=1

# Gráfico do campo de inclinações
PG5=plot([],figsize=(6,4))
P=function('P')(t)
f(t,P) = P*(a-b*P)
PG = plot_slope_field(f,[t,0,4],[P,-4,4], plot_points=25)

# Gráfico de soluções não constantes nas regiões R1, R2 e R3
p=function('p')(t)
edo = diff(p, t) == p*(a-b*p)

# A EDO é resolvida numericamente com o comando desolve_rk4 do SAGE
# que usa o método de Runge-Kutta de ordem 4
# Optamos pela solução numérica porque a solução analítica do SAGE fornece P(t) implicitamente
# o que dificulta a construção dos gráficos
PG1=desolve_rk4(edo,p, ics=[0,0.1],step=0.1, end_points=4, output='plot',color='green')
PG2=desolve_rk4(edo,p, ics=[0,4],step=0.1, end_points=4, output='plot')
PG3=desolve_rk4(edo,p, ics=[3,-4],step=0.1, end_points=0, output='plot')

# Gráfico das soluções constantes ou de equilíbrio correspondentes aos pontos críticos
LL1=line( [ (0, 0), (4, 0) ],color='red')
LL2=line( [ (0, a/b), (4, a/b) ],color='red')
PG6=PG5+PG+PG1+PG2+PG3+LL1+LL2
PG11.axes(False)
PAA=graphics_array((PG11,PG6))
show(PAA)

Analisando o Retrato de Fase

Além das duas soluções constantes ou de equilíbrio (em vermelho), o gráfico mostra três soluções particulares de forma que $P(t)$ é decrescente em $R1$, crescente em $R2$ e decrescente em $R3$. Se denotarmos a condição inicial $P(0)=P_0$, então teremos os seguintes comportamentos para cada uma das regiões:

• Para $P_0 \lt 0$ (Região $R1$): a solução $P(t)$ é limitada acima e decresce indefinidamente conforme $t$ aumenta; visualmente, P(t) se afasta cada vez mais de P(t)=0, à medida que $t$ aumenta
• Para $0 \lt P_0 \lt a/b$ (Região $R2$): a solução $P(t)$ é limitada tanto acima quanto abaixo; como ela é decrescente, para $t\rightarrow\infty$, $P(t)\rightarrow a/b$; em outras palavras, a solução se afasta da solução de equilíbrio $P=0$ e se aproxima da solução de equilíbrio $P=a/b$
• Para $P_0 \gt a/b$ (Região $R3$): a solução $P(t)$ é limitada abaixo; para $t\rightarrow\infty$, $P(t)\rightarrow a/b$, ou seja, com o passar do tempo, a solução se aproxim da solução de equilíbrio $P=a/b$

Observe como fomos capazes de prever, com base apenas nos gráficos, o comportamento de longo prazo das soluções da EDO; observe, também, que os dois gráficos exbidos lado a lado contêm informações qualitativas comuns. Finalmente, embora tenhamos mostrado apenas uma solução em cada região, basta examinar o gráfico para se convencer de que, por exemplo, todas as soluções cujas condições iniciais pertençam à região $R2$ terão o mesmo comportamento qualitativo da curva verde (e, analogamente, para as outras regiões).

Pontos Críticos Atratores e Repulsores

O que veremos, agora, é uma classificação dos pontos críticos. Observe o retrato de fase novamente. Primeiro, examinemos o ponto crítico $P=a/b$. O que vemos é que todas as soluções que começam "próximos" a esse ponto crítico, em ambos os lados, tendem para ele (assintoticamente): $\displaystyle\lim_{t\rightarrow\infty} P(t)=a/b$. Um ponto crítico deste tipo é dito assintoticamente estável e também é denominado atrator, uma vez que "atrai" para si todas as soluções em uma certa vizinhança. Por outro lado, para o ponto crítico $P=0$, o comportamento é exatamente oposto: todas as soluções que começam "próximo" a ele, tendem a se afastar. Este tipo de ponto crítico é dito instável ou, por razões óbvias, repulsor. Existem pontos críticos para os quais, em um dos lados as soluções são atraídas e no outro lado, elas são repelidas; neste caso, o ponto crítico é chamado semi-estável.

x=var('x')
a,b=var('a','b')
assume(b>0)
assume(a>0)
show(integrate(((b^2+4*a*x^2)^(1/2))/(x^2),x))

$\displaystyle 2 \, \sqrt{a} \operatorname{arsinh}\left(\frac{2 \, \sqrt{a} x}{b}\right) - \frac{\sqrt{4 \, a x^{2} + b^{2}}}{x}$