dinâmica

Exercícios Resolvidos - Sistemas não lineares


Problema 1

Uma partícula com massa m , desloca-se ao longo do eixo dos x sob a ação de uma força resultante F x que depende da posição x e da componente da velocidade v x . 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:
(a) F x = m x (1 + v x )
(b) F x = m x ( x 2 + v x 1)

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

˙ x = v x ˙ v x = F x m

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

J ( x , v x ) = v x x v x v x ( x (1 + v x )) x ( x (1 + v x )) v x = 01 1 v x x

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

A = J (0 , 0) = 01 10

Que tem traço nulo e determinante (positivo) igual a 1. Como tal, o ponto de equilíbrio poderá ser centro ou foco (na aproximação linear não há dúvida que é centro, mas devido aos termos não lineares um centro pode tornar-se foco; deverá ser conferido no retrato de fase). 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

E comprova-se que o ponto de equilíbrio é um centro. Qualquer curva de evolução com v x maior que -1 é um ciclo.

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

J ( x , v x ) = v x x v x v x x ( x 2 + v x 1) x x ( x 2 + v x 1) v x = 01 3 x 2 v x + 1 x

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

A 1 = J (0 , 0) = 01 1 0

Com determinante negativo e, como tal, o ponto (0, 0) é ponto de sela.

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

A 2 = J (1 , 0) = 01 2 1

Com traço t = 1 e determinante d = 2 . Como d é maior que t 2 /4 , o ponto (1, 0) é foco atrativo.

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

A 3 = J ( 1 , 0) = 01 21

Com traço t = 1 e determinante d = 2 . Como d é maior que t 2 /4 , o ponto (-1, 0) é 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]);

Retrato de fase

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 l = 50  cm e massa m = 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

¨ θ = g l sin θ C l m | ˙ θ | ˙ θ
Se a massa m estiver concentrada numa esfera de raio R = 2  cm, a expressão para a constante C é dada pela equação 4.14: C = π ρ R 2 /4 , onde ρ = 1 . 2  kg/m3 é a massa volúmica do ar. Trace os gráficos de θ ( t ) , ω ( t ) 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:
(a) O pêndulo parte do repouso com um ângulo inicial θ = 120 .
(b) 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 t = 0 até t = 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 ]

Conclui-se 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 é um resultado convergente com 4 algarismos significativos. 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 roda três voltas completas, 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 ¨ x + ˙ x 2 + 4 x 2 = 4 ,
(a) Escreva as equações de evolução do sistema dinâmico associado à equação.
(b) Encontre os pontos de equilíbrio do sistema.
(c) Determine a matriz jacobiana.
(d) Caracterize cada um dos pontos de equilíbrio.
(e) Se em t = 0 os valores da variável x e da sua derivada são x 0 = 1 e ˙ x 0 = 1 , determine (numericamente) os valores da variável e da sua derivada em t = 2 .

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

v = ˙ x

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

˙ v + v 2 + 4 x 2 = 4

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

˙ x = v ˙ v = 4 v 2 4 x 2

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

4 x 2 = 4 x =± 1

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

J ( x , v ) = v x v v (4 v 2 4 x 2 ) x (4 v 2 4 x 2 ) v = 01 8 x 2 v

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

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

J (1 , 0) = 01 80

Como o traço dessa matriz é nulo e o determinante é 8, os valores próprios são números imaginários e o ponto x = 1 , v = 0 é um centro ou foco (o retrato de fase mostra que é centro).

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

J ( 1 , 0) = 01 8 0

Como o traço dessa matriz é nulo e o determinante é -8, os valores próprios são reais, com sinais opostos. O ponto x = 1 , v = 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 x e ˙ x , em t = 2 , são: 0.5869 e 0.8277.

Problema 8

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

˙ x = 2 x y 3 x 4 ˙ y = y 4 2 x 3 y
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:
(a) Determine o versor na direção da velocidade de fase em qualquer ponto do eixo dos x e em qualquer ponto do eixo dos y .
(b) Determine o versor na direção da velocidade de fase em qualquer ponto das duas retas y = x e y = x .
(c) 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?
(d) Diga se existem ciclos, órbitas homoclínicas ou heteroclínicas e, caso a resposta seja afirmativa, quantas.

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

u = x 4 ˆ ı = ˆ u = ˆ ı

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

u = y 4 ˆ = ˆ u u = ˆ

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

u = x 4 ˆ ı x 4 ˆ

com módulo igual a 2 x 4 e versor:

e u = x 4 ˆ ı x 4 ˆ 2 x 4 = 1 2( ˆ ı ˆ )

Na reta y = x ,

u = 3 x 4 ˆ ı + 3 x 4 ˆ = e u = 1 2( ˆ ı + ˆ )

(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 é ponto de sela.

Folha de Descartes

(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, asteroide, sonda espacial, etc) de massa m 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 é,

U = G M m x 2 + y 2
onde G é a constante de gravitação universal, M é a massa do Sol, e as coordenadas x e y 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 G M será igual a 4 π 2 .
(a) Encontre as equações de movimento do corpo celeste, em unidades de anos para o tempo e UA para as distâncias.
(b) O cometa Halley aproxima-se 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 0 . 587ˆ ı , com velocidade inicial 11 . 50ˆ , com intervalos de tempo t = 0 . 05  anos. Trace a órbita desde t = 0 até t = 100  anos. Que pode concluir acerca do erro numérico?
(c) Repita o procedimento da alínea anterior com t = 0 . 02  anos e trace a órbita desde t = 0 até t = 150  anos. Que pode concluir acerca do erro numérico?
(d) 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: x , y , ˙ x e ˙ y . 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 ˙ x e ˙ y . Para aplicar as equações de Lagrange é necessário definir xp e yp como derivadas x e y 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)   4 π 2 m x y 2 + x 2 3/2 + m x p p = 0
(%i23) eq1: solve(%,xpp)[1];
(%o23)   x p p = 4 π 2 x y 2 + x 2 3/2
(%i24) diff(diff(Ec,yp),t) - diff(Ec,y) + diff(U,y) = 0;
(%o24)   m y p p + 4 π 2 m y y 2 + x 2 3/2 = 0
(%i25) eq2: solve(%,ypp)[1];
(%o25)   y p p = 4 π 2 y y 2 + x 2 3/2

As equações de movimento são:

¨ x = 4 π 2 x ( x 2 + y 2 ) 3/2 ¨ y = 4 π 2 y ( x 2 + y 2 ) 3/2

(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 ( t , x , y , ˙ x , ˙ y ) em diferentes instantes entre 0 e 100. Como tal, o gráfico de trajetória do cometa ( y vs x ) 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 t = 0 . 05 não é suficientemente pequeno e o resultado obtido têm erro numérico muito elevado.

(c) Reduzindo o valor dos incrementos de 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 t 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 t aproximadamente 36 anos. Como foram usados incrementos de t iguais a 0.02 = 1/50, 36 anos aparecerá na posição 1801 da lista. Observando a lista de valores de x nessa parte da lista:

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

Conclui-se que o valor mínimo de x (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.