Problemas Resolvidos

10. Sistemas não lineares

Problema 1

Uma partícula com massa , desloca-se ao longo do eixo dos sob a ação de uma força resultante que depende da posição e da componente da velocidade . Para cada um dos casos seguintes encontre os pontos de equilíbrio, diga que tipo de ponto equilíbrio é cada um (estável ou instável; centro, foco, nó ou ponto de sela) e desenhe o retrato de fase mostrando as órbitas mais importantes:

As equações de evolução são:

(a) Nos pontos de equilíbrio, = 0 e = 0, ou seja, existe um único ponto de equilíbrio em ( , ) = (0, 0). A matriz jacobiana é:

E a matriz da aproximação linear na vizinhança do ponto de equilíbrio é:

Que tem traço negativo e determinante (positivo) igual a 1. Como tal, o ponto de equilíbrio é um centro. O retrato de fase traça-se com o comando:

(%i1) plotdf ([vx,-x*(1+vx)], [x,vx], [x,-5,5], [vx,-5,5]);

Retrato de fase do problema 10.1a

(b) Nos pontos de equilíbrio, = 0 e = 0, ou seja, existem três pontos de equilíbrio em ( , ) = (0, 0), (-1, 0) e (1, 0). A matriz jacobiana é:

A matriz da aproximação linear na vizinhança do ponto de equilíbrio (0, 0) é:

Com determinante negativo e, como tal, o ponto de equilíbrio é um ponto de sela.

No ponto de equilíbrio (1, 0) a matriz da aproximação linear é:

Com traço = -1 e determinante = 2. Como é maior que , o ponto (1, 0) é um foco atrativo.

No ponto de equilíbrio (−1, 0) a matriz da aproximação linear é:

Com traço = 1 e determinante = 2. Como é maior que , o ponto (−1, 0) é um foco repulsivo.

O retrato de fase traça-se com o comando:

(%i2) plotdf ([vx,-x*(x^21+vx)], [x,vx], [x,-2.5,2.5], [vx,-2,5]);

(b) Um ponto de sela em ( , ) = (0, 0), um foco instável em ( , ) = (−1, 0) e um foco estável em ( , ) = (1, 0).

Retrato de fase do problema 10.1b

Problema 4

A amplitude de oscilação de um pêndulo decresce, devido à força de resistência do ar e ao atrito no eixo. Admita um pêndulo de comprimento = 50 cm e massa = 0.150 kg, em que o atrito no eixo é desprezável mas a resistência do ar não. A equação de movimento é a equação 8.8

Se a massa estiver concentrada numa esfera de raio = 2 cm, a expressão para a constante é dada pela equação 4.14: , onde = 1.2 kg/m é a massa volúmica do ar. Trace os gráficos de , e da curva de evolução no espaço de fase e explique o significado físico da solução, para os dois casos seguintes:

  1. O pêndulo parte do repouso com um ângulo inicial  = 120°.
  2. O pêndulo é lançado desde = 60°, com velocidade angular inicial = −7.8 s−1.

(a) Usando o programa rk, com intervalos de tempo de 0.1, desde = 0 até = 50,

(%i3) [g, l, m]: [9.8, 0.5, 0.15]$
(%i4) C: %pi*1.2*0.02^2/4$
(%i5) s: rk ([w,-g*sin(q)/l-C*l*abs(w)*w/m], [q,w], [2*%pi/3,0], [t,0,50,0.1])$
(%i6) last (s);
(%o6)  [ 50.0, 0.408596821360162, 6.193790347574476 ]

Executando novamente o programa rk com intervalos de tempo dez vezes menores,

(%i7) s: rk ([w,-g*sin(q)/l-C*l*abs(w)*w/m], [q,w], [2*%pi/3,0], [t,0,50,0.01])$
(%i8) last (s);
(%o8)  [ 50.0, - 0.8184365726225941, 5.503739362621793 ]

Mostrando que é necessário reduzir ainda mais o valor dos intervalos de tempo, para obter uma solução convergente:

(%i9) s: rk ([w,-g*sin(q)/l-C*l*abs(w)*w/m], [q,w], [2*%pi/3,0], [t,0,50,0.005])$
(%i10) last (s);
(%o10)  [ 50.0, - 0.8184437883132009, 5.503721542035767 ]

Que mostra que o resultado já é convergente. O gráfico do ângulo e da velocidade angular, em função do tempo, obtém-se com o comando:

(%i11) plot2d ([[discrete, makelist([p[1],p[2]],p,s)], [discrete, makelist([p[1],p[3]],p,s)]], [legend, "angulo", "vel. angular"], [xlabel, "t"]);

Ângulo e velocidade angular do pêndulo amortecido

E a curva de evolução no espaço de fase é o gráfico da velocidade angular em função do ângulo, obtido com o seguinte comando:

(%i12) plot2d ([discrete, makelist([p[2],p[3]],p,s)], [xlabel, "angulo"], [ylabel, "vel. angular"]);

Curva de evolução do pêndulo amortecido

Os dois gráficos mostram que pêndulo oscila com amplitude que decresce lentamente.

(b) Usando o programa rk, com os mesmos intervalos de tempo usados para obter os gráficos na alínea anterior,

(%i13) s: rk ([w,-g*sin(q)/l-C*l*abs(w)*w/m], [q,w], [%pi/3,-7.8], [t,0,50,0.005])$

Os gráficos do ângulo e da velocidade angular, em função do tempo, e da curva de evolução no espaço de fase, obtêm-se repetindo os mesmos comandos da alínea anterior:

Ângulo e velocidade angular do pêndulo amortecido
Curva de evolução do pêndulo amortecido

O pêndulo faz três voltas completas, rodando no sentido horário, e quando passa a quarta vez pela posição de equilíbrio estável, começa a oscilar com amplitude que decresce lentamente.

Problema 7

Para analisar a equação diferencial não linear ,

  1. Escreva as equações de evolução do sistema dinâmico associado à equação.
  2. Encontre os pontos de equilíbrio do sistema.
  3. Determine a matriz jacobiana.
  4. Caracterize cada um dos pontos de equilíbrio.
  5. Se em = 0 os valores da variável e da sua derivada são = 1 e = 1, determine (numericamente) os valores da variável e da sua derivada em = 2.

(a) Define-se uma segunda variável de estado:

e substitui-se na equação do sistema:

Como tal, as duas equações de evolução —expressões das derivadas das duas variáveis de estado— são:

(b) Para resolver esta alínea não é necessário ter resolvido a alínea anterior. Basta observar que nos pontos de equilíbrio permanece constante e, assim sendo, = = 0. Substituindo na equação do sistema,

(c) Usando as equações obtidas na alínea (a),

(Também pode usar-se a função jacobian do Maxima, para determinar a matriz).

(d) Substituindo = 1 e = 0 na matriz jacobiana obtém-se:

Como o traço dessa matriz é nulo e o determinante é 8, os valores próprios são números imaginários e o ponto = 1, = 0 é um centro.

Substituindo = −1 e = 0 na matriz jacobiana obtém-se:

Como o traço dessa matriz é nulo e o determinante é -8, os valores próprios são reais, com sinais opostos. O ponto = −1, = 0 é então ponto de sela.

(e) Usa-se a função rk do Maxima várias vezes, com valores decrescentes dos intervalos de tempo, até se obterem valores convergentes do resultado:

(%i14) last( rk( [v,4-v^2-4*x^2], [x,v], [1,1], [t,0,2,0.1]));
(%o14)    [ 2.0, 0.58688, 0.82753 ]
(%i15) last( rk( [v,4-v^2-4*x^2], [x,v], [1,1], [t,0,2,0.05]));
(%o15)    [ 2.0, 0.58687, 0.82768 ]

Ou seja, os valores aproximados de e , em = 2, são: 0.5869 e 0.8277

Problema 8

O sistema dinâmico com equações de evolução:

tem um único ponto de equilíbrio na origem. A matriz jacobiana nesse ponto é igual a zero e, portanto, os valores próprios (nulos) não podem ser usados para caraterizar o ponto de equilíbrio. Use o seguinte método para analisar o retrato de fase do sistema:

  1. Determine o versor na direção da velocidade de fase em qualquer ponto do eixo dos e em qualquer ponto do eixo dos .
  2. Determine o versor na direção da velocidade de fase em qualquer ponto das duas retas e .
  3. Faça a mão um gráfico mostrando os versores que encontrou nas alíneas a e b, em vários pontos nos 4 quadrantes do espaço de fase, e trace algumas curvas de evolução seguindo as direções da velocidade de fase. Com base nesse gráfico, que tipo de ponto de equilíbrio julga que é a origem?
  4. Diga se existem ciclos, órbitas homoclínicas ou heteroclínicas e no caso afirmativo quantas.

(a) No eixo dos , é igual a zero e a velocidade de fase é,

No eixo dos , é igual a zero e a velocidade de fase é,

(b) Na reta = , a velocidade de fase é,

com módulo igual a e versor:

Na reta = − ,

(c) A figura seguinte mostra os versores encontrados nas duas alíneas anteriores e algumas curvas de evolução. Como há curvas que se aproximam da origem e curvas que se afastam dele, a origem é um ponto de sela.

Retrato de fase do problema 10.8

(d) Não existem ciclos nem órbitas heteroclínicas. Existem um número infinito de órbitas homoclínicas: todas as curvas de evolução no primeiro e terceiro quadrantes são órbitas homoclínicas.

Problema 10

Qualquer corpo celeste (planeta, cometa, asteróide, sonda espacial, etc) de massa no sistema solar tem uma energia potencial gravítica produzida pelo Sol, que é responsável pelas órbitas elípticas desses corpos. A expressão para a energia potencial é,

onde é a constante de gravitação universal, é a massa do Sol, e as coordenadas e são medidas no plano da órbita do corpo celeste, com origem no Sol. Se as distâncias forem medidas em unidades astronómicas, UA, e os tempos em anos, o produto será igual a .

  1. Encontre as equações de movimento do corpo celeste, em unidades de anos para o tempo e UA para as distâncias.
  2. O cometa Halley chega até uma distância mínima do Sol igual a 0.587 UA. Nesse ponto, a sua velocidade é máxima, igual a 11.50 UA/ano, e perpendicular à sua distância até o Sol. Determine numericamente a órbita do cometa Halley, a partir da posição inicial , com velocidade inicial , com intervalos de tempo Δ  = 0.05 anos. Trace a órbita desde = 0 até = 100 anos. Que pode concluir acerca do erro numérico?
  3. Repita o procedimento da alínea anterior com Δ  = 0.02 anos e trace a órbita desde = 0 até = 150 anos. Que pode concluir acerca do erro numérico?
  4. Diga qual é, aproximadamente, a distância máxima que o cometa Halley se afasta do Sol, e compare a órbita do cometa com as órbitas do planeta mais distante, Neptuno (órbita entre 29.77 UA e 30.44 UA) e do planeta mais próximo do Sol, Mercúrio (órbita entre 0.31 UA e 0.39 UA) (Plutão já não é considerado um planeta).

(a) Há quatro variáveis de estado: , , e . As expressões das energias cinética e potencial são:

(%i16) Ec: m*(xp^2 + yp^2)/2$
(%i17) U: -4*%pi^2*m/sqrt(x^2 + y^2)$

Onde xp e yp representam as velocidades generalizadas e . Para aplicar as equações de Lagrange é necessário definir xp e yp como derivadas e em ordem ao tempo, e definir também xpp e ypp como derivadas de xp e yp:

(%i18) gradef (x,t,xp)$
(%i19) gradef (y,t,yp)$
(%i20) gradef (xp,t,xpp)$
(%i21) gradef (yp,t,ypp)$

As duas equações de Lagrange conduzem às duas equações de movimento:

(%i22) diff(diff(Ec,xp),t) - diff(Ec,x) + diff(U,x) = 0;
(%o22)  
(%i23) eq1: solve(%,xpp)[1];
(%o23)  
(%i24) diff(diff(Ec,yp),t) - diff(Ec,y) + diff(U,y) = 0;
(%o24)  
(%i25) eq1: solve(%,ypp)[1];
(%o25)  

As equações de movimento são:

(b) Usando as condições iniciais dadas e o intervalo de tempo desde 0 até 100, com incrementos de 0.05, a solução numérica do problema obtém-se com o programa rk:

(%i26) o: rk ([xp,yp,rhs(eq1),rhs(eq2)], [x,y,xp,yp], [0.587,0,0,11.5], [t,0,100,0.05])$

onde o é uma lista com várias listas de cinco elementos, com os valores de ( , , , , ) em diferentes instantes entre 0 e 100. Como tal, o gráfico de trajetória do cometa ( vs ) pode ser obtido com o seguinte comando:

(%i27) plot2d ([discrete, makelist([p[2],p[3]], p, o)], [xlabel,"x"], [ylabel,"y"], [y,-10,10], same_xy);

Usou-se a opção same_xy para que a escala nos dois eixos seja igual, mostrando a forma real da trajetória. O resultado é o gráfico seguinte:

Órbita do cometa Halley, com erro numérico

O facto de que o satélite não repete a mesma trajetória, mas aproxima-se cada vez mais do Sol, indica que a sua energia mecânica diminui, em vez de permanecer constante, como era suposto acontecer. Conclui-se então que o intervalo Δ  = 0.05 não é suficientemente pequeno e os dados obtidos têm um erro numérico muito elevado.

(c) Reduzindo o valor dos incrementos do tempo:

(%i28) o: rk ([xp,yp,rhs(eq1),rhs(eq2)], [x,y,xp,yp], [0.587,0,0,11.5], [t,0,100,0.02])$
(%i29) plot2d ([discrete, makelist([p[2],p[3]], p, o)], [xlabel,"x"], [ylabel,"y"], [y,-10,10], same_xy);

Órbita do cometa Halley

O erro numérico é muito menor, mas o cometa continua a perder energia; seria necessário reduzir ainda mais o valor de Δ  para diminuir o erro.

(d) O comando

(%i30) plot2d ([discrete,makelist([p[1],p[2]],p,o)]);

Mostra que o cometa está mais afastado do Sol em aproximadamente 36 anos. Como foram usados incrementos de iguais a 0.02 = 1/50, 36 anos aparecerá na posição 1801 da lista. Observando a lista de valores de nessa parte da lista:

(%i31) makelist (o[i][2], i, 1780, 1820);

Conclui-se que o valor mínimo de (distância máxima ao Sol) é aproximadamente 34.14 UA. Essa distância máxima é maior do que a órbita de Neptuno e a distância mínima, 0.587 UA, está entre as órbitas de Mercúrio e Vénus.