dinâmica

10. Sistemas não lineares

Introdução

Hoverboard

Um hoverboard tem apenas um eixo e duas rodas. A pessoa no hoverboard pode rodar com ele à volta do eixo, tal como um pêndulo. O pêndulo tem duas posições de equilíbrio, quando o seu centro de gravidade se encontra na mesma linha vertical que passa pelo eixo. O ponto de equilíbrio por baixo do eixo é ponto de equilíbrio estável e se o centro de gravidade estiver próximo desse ponto, o pêndulo oscila. No caso do hoverboard, o centro de gravidade encontra-se próximo do ponto de equilíbrio por cima do eixo, que é um ponto de equilíbrio instável; como tal, a pessoa rodará, caindo para o chão, ou para a frente ou para trás. O sistema de controlo automático do motor desloca o hoverboard e a pessoa na direção necessária para evitar essa queda. Um pêndulo com o seu centro de gravidade próximo do ponto de equilíbrio instável (por cima do eixo), chama-se pêndulo invertido. Um hoverboard, um segway e um monociclo são exemplos de pêndulos invertidos.

10.1. Aproximação linear

Nos sistemas dinâmicos com duas variáveis de estado:

(10.1)
˙ x 1 = f 1 ( x 1 , x 2 )˙ x 2 = f 2 ( x 1 , x 2 )

cada uma das funções f 1 e f 2 podem ser escritas na forma de uma série de Taylor, na vizinhança de um ponto qualquer ( a , b ) do espaço de fase:

(10.2)
f i ( x 1 , x 2 ) = f i ( a , b ) + ( x 1 a ) f i x 1 ( a , b ) + ( x 2 b ) f i x 2 ( a , b ) + . . .

onde o índice i pode ser 1 ou 2. Se o ponto ( a , b ) é um ponto de equilíbrio, então f 1 ( a , b ) = 0 = f 2 ( a , b ) e, portanto, o primeiro termo das duas séries é nulo. Mudando a origem de coordenadas para o ponto de equilíbrio ( a , b ), isto é, num novo sistema de coordenadas: x = x 1 a , y = x 2 b , as funções são, aproximadamente,

(10.3)
f i ( x , y ) = x f i x 1 ( a , b ) + y f i x 2 ( a , b )    ( i = 1 , 2)

Ou seja, uma combinação linear das novas variáveis x e y , onde as constantes são os valores das derivadas parciais no ponto de equilíbrio ( a , b ). Substituindo essas aproximações no sistema 10.1, obtém-se um sistema linear ( ˙ x = ˙ x 1 e ˙ y = ˙ x 2 , porque a e b são constantes).

(10.4)
˙ x ˙ y = f 1 x 1 f 1 x 2 f 2 x 1 f 2 x 2 ( a , b ) x y

esta aproximação linear é válida apenas numa vizinhança da origem ( x = 0 , y = 0 ), ou seja, quando x 1 e x 2 estejam próximas de a e b .

A matriz quadrada na equação 10.4 chama-se matriz jacobiana e representa-se por J ( x 1 , x 2 ) . Substituindo as coordenadas ( a , b ) do ponto de equilíbrio na matriz jacobiana, obtém-se uma matriz constante. Por cada ponto de equilíbrio existe uma matriz de coeficientes constantes, que define o sistema linear que aproxima bem o sistema não linear na vizinhança do ponto de equilíbrio. Os valores e vetores próprios de cada uma dessas matrizes permitem analisar a estabilidade do sistema, na vizinhança de cada ponto de equilíbrio, da mesma forma que é feito para os sistemas lineares.

Exemplo 10.1

Classifique os pontos de equilíbrio e represente o retrato de fase do sistema:  ˙ x 1 = 4 x 21 4 x 22 ˙ x 2 = x 22 x 21 + 1

Resolução. Já foi demonstrado no exemplo 7.2 do capítulo 7, que este sistema tem quatro pontos de equilíbrio. As funções f 1 e f 2 e os pontos de equilíbrio são armazenados em duas listas assim:

(%i1) f: [4-x1^2-4*x2^2, x2^2-x1^2+1]$
(%i2) equilibrio: solve(f)$

Convém definir também outra lista com os nomes das variáveis de estado:

(%i3) v: [x1, x2]$

A matriz jacobiana, com duas linhas e duas colunas, obtem-se com o comando jacobian do Maxima, que precisa de duas listas: uma lista com as funções e outra lista com os nomes das variáveis

(%i4) J: jacobian (f,v);
(%o4)   2 x 1 8 x 2 2 x 12 x 2

Substituindo as coordenadas dos pontos de equilíbrio, obtêm-se as matrizes dos sistemas lineares que aproximam o sistema na vizinhança de cada um desses pontos. Por exemplo, no primeiro ponto de equilíbrio,

(%i5) subst (equilibrio[1], J);
(%o5)   2 5/2 58 3 52 5/2 5 2 3 5

Para estudar a estabilidade do sistema na vizinhança desse ponto de equilíbrio, calculam-se os valores próprios dessa matriz.

(%i6) eigenvectors (%)$
(%i7) float (%);
(%o7) [[ [ 3 . 963 , 4 . 944] , [1 . 0 , 1 . 0] ] , [[ [ 1 . 0 , 1 . 048] ] , [[ 1 . 0 , 0 . 3896] ] ] ]

O resultado mostra 4 listas; a primeira lista são os valores próprios, a segunda lista são as multiplicidades de cada valor próprio, e as últimas duas listas são os vetores próprios.

Nesse ponto de equilíbrio os valores próprios são reais, com sinais opostos; conclui-se que é um ponto de sela. O quarto ponto de equilíbrio também é ponto de sela:

(%i8) subst (equilibrio[4], J);
(%o8)   2 5/2 5 8 3 5 2 5/2 52 3 5
(%i9) eigenvectors (%)$
(%i10) float (%);
(%o10) [[ [ 4 . 944 , 3 . 963] , [1 . 0 , 1 . 0] ] , [[ [ 1 . 0 , 0 . 3896] ] , [[ 1 . 0 , 1 . 048] ] ] ]

No segundo ponto de equilíbrio:

(%i11) subst (equilibrio[2], J);
(%o11)      2 5/2 58 3 5 2 5/2 5 2 3 5
(%i12) eigenvectors (%)$
(%i13) float (map (rectform, %));
(%o13) [[ [ 3 . 929i 2 . 04 , 3 . 929i 2 . 04] , [1 . 0 , 1 . 0] ] , [[ [ 1 . 0 , 0 . 07912 0 . 634i ] ] , [[ 1 . 0 , 0 . 634i + 0 . 07912] ] ] ]

Como os valores próprios são complexos, com parte real negativa, o ponto de equilíbrio é um foco atrativo (estável). Cálculos semelhantes para o terceiro ponto de equilíbrio mostram que também é um foco, mas repulsivo (instável), porque os valores próprios são complexos, com parte real positiva. O retrato de fase constrói-se usando o comando:

(%i14) plotdf (f, v, [x1,-3,3], [x2,-3,3])$

Na figura 10.1 mostra-se o resultado. Existe um único ponto de equilíbrio estável, um foco atrativo, em ( x 1 , x 2 ) = (1.265, -0.7746). Os outros 3 pontos de equilíbrio, dois pontos de sela e um foco repulsivo, são instáveis. As duas curvas de evolução que foram traçadas a sair do foco repulsivo em ( x 1 , x 2 ) = (-1.265, 0.7746) e a continuação dessas curvas passando pelos pontos de sela, delimitam a região de estabilidade, em que se o estado inicial do sistema estiver nessa região, o estado final aproximar-se-á do ponto de equilíbrio estável.

Retrato de fase dum sistema não linear
Figura 10.1: Retrato de fase do sistema ˙ x 1 = 4 x 21 4 x 22 , ˙ x 2 = x 22 x 21 + 1 .

10.2. O pêndulo

Pendulo com barra e disco

Figura 10.2: Pêndulo.

O tipo de pêndulo estudado nesta secção é formado por um objeto ligado a uma barra rígida atravessada por um eixo horizontal fixo (figura 10.2). Esse tipo de pêndulo pode rodar num plano vertical dando voltas completas. O sistema tem um único grau de liberdade, θ , que é o ângulo que a barra faz com a vertical. Seja θ = 0 quando o pêndulo está na posição mais baixa e θ = π na posição mais alta. A velocidade angular é ˙ θ e a velocidade do centro de massa é r ˙ θ onde r é a distância desde o centro de massa até o eixo.

A energia cinética é:

(10.5)
E c = 1 2 m r 2 ˙ θ 2 + 1 2 I cm ˙ θ 2

Onde m é a massa total e I cm o momento de inércia em relação ao centro de massa. De acordo com o teorema dos eixos paralelos 5.28, o momento de inércia em relação ao eixo do pêndulo é I e = m r 2 + I cm , que pode ser escrito I e = m r 2g , onde r g é o raio de giração em relação ao eixo. Como tal, a energia cinética é

(10.6)
E c = 1 2 m r 2g ˙ θ 2

A energia potencial gravítica é (arbitrando energia nula em θ = π /2 )

(10.7)
U = m g r cos θ

Ignorando a resistência do ar, a equação de Lagrange conduz à equação de movimento:

(10.8)
¨ θ = g l sin θ

onde l = r 2g / r define o comprimento eficaz do pêndulo. No caso particular dum pêndulo simples, em que a massa da barra é desprezável e o objeto é pequeno, l é a distância desde o objeto até o eixo (ver exemplo 8.5 do capítulo 8).

As equações de evolução obtêm-se substituindo a derivada de θ pela velocidade angular ω :

(10.9)
˙ θ = ω ˙ ω = g l sin θ

Estas equações não lineares não podem ser resolvidas analiticamente, mas podem ser resolvidas por aproximação numérica. O comando rk do Maxima usa-se para obter a solução numérica pelo método de Runge-Kutta de quarta ordem; é necessário dar 4 argumentos ao comando: uma lista de expressões para as componentes da velocidade de fase, uma lista com os nomes das variáveis de estado, uma lista com valores iniciais para essas variáveis e um intervalo de valores para a variável independente, incluindo o nome dessa variável, valor inicial, valor final e valor dos incrementos nesse intervalo. O comando rk produz uma lista de pontos que aproximam a solução; cada ponto terá as coordenadas da variável independente, seguida pelas variáveis de estado.

Por exemplo, para um pêndulo com l igual a 50 cm, largado do repouso com ângulo inicial de 30°, a solução aproximada é obtida com ( q e w representam θ e ω ):

(%i15) s: rk([w,-(9.8/0.5)*sin(q)],[q,w], [%pi/6,0], [t,0,5,0.01])$

Os gráficos de θ e ω em função do tempo e a curva de evolução no espaço de fase θω obtêm-se com os seguintes comandos:

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

Os dois gráficos são apresentados na figura 10.3.

Solução do pêndulo com amplitude de 30 graus, em função de t Solução do pêndulo com amplitude de 30 graus, no plano de fase
Figura 10.3: Oscilações de um pêndulo de 50 cm com amplitude de 30°.

A lista de dados numéricos obtida permite concluir que o período de oscilação está entre 1.44 s e 1.45 s. Os gráficos na figura 10.3 são muito parecidos com os gráficos de um oscilador harmónico simples. Se o ângulo inicial for maior, essa semelhança começa a desaparecer. Por exemplo, a figura 10.4 mostra os resultados obtidos com ângulo inicial de 120°.

Solução do pêndulo com amplitude de 120 graus, em função de t Solução do pêndulo com amplitude de 120 graus, no plano de fase
Figura 10.4: Oscilações de um pêndulo de 50 cm com amplitude de 120°.

Nesse caso conclui-se a partir dos dados numéricos que o período de oscilação aumenta, em relação à amplitude de 30° e está entre 1.94 s e 1.95 s.

Nos dois casos apresentados nas figuras 10.3 e 10.4, a curva de evolução é um ciclo, indicando que existe um ponto de equilíbrio estável na região interna do ciclo.

Os pontos de equilíbrio do pêndulo, onde os lados direitos das equações 10.9 são nulos, encontram-se em θ = 0 , ± π , ± 2 π … e ω = 0 .

Os pontos em θ = 0 , ± 2 π , ± 4 π … são realmente o mesmo ponto físico, na posição mais baixa do pêndulo, correspondentes à passagem do pêndulo por essa posição, após um número qualquer de voltas. Os pontos em θ =± π , ± 3 π … são também o mesmo ponto físico, na posição mais alta do pêndulo.

10.3. Aproximação linear do pêndulo

A matriz jacobiana correspondente às equações 10.9 do pêndulo é

(10.10)
01 g l cos θ 0

No ponto de equilíbrio em θ = 0 (em geral, 0, ± 2 π , ± 4 π ,…), a matriz é:

(10.11)
01 g l 0

que é a matriz de um oscilador harmónico simples, analisada no exemplo  9.4 do capítulo 9. Os dois valores próprios são ± i g / l , o ponto de equilíbrio θ = ω = 0 é um centro e se o estado inicial do sistema está próximo desse ponto, o pêndulo oscila com frequência angular = g / l . No caso do pêndulo de 50 cm analisado na secção anterior, essa expressão conduz ao período 1.42 s. Lembre-se que esse valor é apenas uma aproximação, que é melhor quanto menor for a amplitude; os valores do período calculados numericamente na secção anterior são mais realistas.

Na vizinhança do ponto de equilíbrio θ = π (em geral, ± π , ± 3 π ,…), a matriz jacobiana é

(10.12)
01 g l 0

que é a matriz de um oscilador invertido, analisada no exemplo 9.3 do capítulo 9. Os dois valores próprios são ± g / l e o ponto de equilíbrio é ponto de sela (equilíbrio instável).

O retrato de fase no intervalo 10 < θ < 10 , mostrará 3 centros ( 2 π , 0 e 2 π ) e 4 pontos de sela ( 3 π , π , π e 3 π ). No caso l = 50  cm considerado na secção anterior, usa-se o comando:

(%i18) plotdf([w,-(9.8/0.5)*sin(q)],[q,w],[q,-10,10], [w,-20,20]);

O resultado é a figura 10.5. No eixo das abcissas está representado o ângulo θ e no eixo das ordenadas a velocidade angular ω . As duas curvas identificadas com as letras A e B formam parte de uma órbita heteroclínica.

Retrato de fase dum pêndulo
Figura 10.5: Retrato de fase de um pêndulo de 50 cm.

A órbita heteroclínica do pêndulo corresponde ao caso em que a energia mecânica do pêndulo é exatamente igual à energia potencial gravítica no ponto de altura máxima. Usando como referência U = 0 , na posição em que a barra do pêndulo está na horizontal ( θ = π /2 ), a energia potencial no ponto mais alto é U = m g l . Cada uma das curvas A e B corresponde ao movimento em que inicialmente o pêndulo está parado na posição mais alta, desce completando uma oscilação completa e para novamente na posição mais alta, sem voltar a oscilar mais. A diferença entre a órbita heteroclínica e os ciclos, é que nos ciclos as oscilações repetem-se indefinidamente, enquanto que na órbita heteroclínica há apenas uma oscilação.

Dentro da órbita heteroclínica, os ciclos na sua vizinhança correspondem a oscilações em que o pêndulo chega quase até o ponto mais alto, parece ficar parado nesse ponto por alguns instantes e logo desce novamente até o ponto mais baixo, repetindo o movimento no outro lado da vertical.

As órbitas heteroclínicas também são separatrizes no retrato de fase, porque delimitam a região onde existe movimento oscilatório: região sombreada na figura 10.6. Se o estado inicial está dentro dessa região, o pêndulo oscila; caso contrário, o pêndulo descreve movimento circular não uniforme.

Separatrizes no retrato de fase do pêndulo
Figura 10.6: As órbitas heteroclínicas delimitam a região de movimento oscilatório.

As figuras 10.3 e 10.4 mostram que com amplitude 30° a aproximação linear é bastante boa, pois a curva de evolução é muito parecida à do oscilador harmónico simples e o período é próximo do período obtido com a aproximação linear, mas com amplitude de 120°, a aproximação linear já não é muito boa.

10.4. Espaços de fase com várias dimensões

Nos sistemas mecânicos autónomos, por cada grau de liberdade há uma equação de movimento, que implica duas variáveis de estado. Assim sendo, a dimensão do espaço de fase é o dobro do número de graus de liberdade. Se um sistema não é autónomo é necessário acrescentar mais uma dimensão ao espaço de fase, como se mostra na seguinte secção. Existem então sistemas mecânicos com espaços de fase de dimensão 2, 3, 4, 5, …

Nos casos em que o espaço de fase tem mais do que duas dimensões o programa plotdf não pode ser utilizado para esboçar o retrato de fase. É necessário resolver as equações de evolução para alguns valores iniciais específicos e construir gráficos mostrando apenas algumas das variáveis de estado.

10.4.1. Sistemas de equações não autónomas

A forma geral de um sistema com n equações diferenciais não autónomas é:

˙ x 1 = f 1 ( x 1 , x 2 , . . . , x n , t )˙ x 2 = f 1 ( x 1 , x 2 , . . . , x n , t ). . . ˙ x n = f 1 ( x 1 , x 2 , . . . , x n , t )

Para caraterizar cada possível estado do sistema é necessário saber os valores das n variáveis x i e o valor do tempo; ou seja, cada estado é um ponto com n + 1 coordenadas ( x 1 , x 2 ,…, x n , t ) e o espaço de fase tem n + 1 dimensões e o tempo é também variável de estado. As componentes da velocidade de fase são as derivadas das variáveis de estado: ( ˙ x 1 , ˙ x 2 ,…, ˙ x n , ˙ t ). As expressões para as primeiras n componentes são dadas pelo sistema de n equações diferenciais acima e a última componente ˙ t é sempre igual a 1 (derivada de t em ordem a t ). Como tal, o sistema de n equações não autónomas considera-se é equivalente a um sistema autónomo com n + 1 variáveis de estado.

Esse tipo de sistemas de equações podem ser resolvidos também com o comando rk, sem ser necessário indicar t como variável de estado, nem a última componente da velocidade de fase, ˙ t = 1 ; o valor inicial de t dá-se no intervalo de integração e não na lista de valores iniciais das variáveis de estado. No entanto, há que ter em conta que se a velocidade de fase depende da variável independente t , essa variável é também variável de estado.

Exemplo 10.2

A equação diferencial:

t 2 ¨ x + t ˙ x + t 2 1 9 x = 0

é uma equação de Bessel. Escreva a equação como sistema dinâmico e identifique o espaço de fase.

Resolução. Define-se uma variável auxiliar y igual a ˙ x :

(10.13)
˙ x = y

assim sendo, a segunda derivada ¨ x é igual à primeira derivada de y e a equação de Bessel é:

t 2 ˙ y + t y + t 2 1 9 x = 0

resolvendo para ˙ y , obtém-se:

(10.14)
˙ y = 1 9 t 2 1 x y t

Como esta equação não é autónoma, considera-se a variável independente t como mais uma variável de estado, com a equação de evolução trivial:

(10.15)
d t d t = 1

O espaço de fase tem três dimensões e cada estado tem coordenadas ( t , x , y ). O sistema dinâmico é definido pelas 3 equações 10.13, 10.14 e 10.15.

10.4.2. Lançamento de projéteis

Forças num projétil no ar

Figura 10.7: Projétil no ar.

No caso do lançamento de um projétil com velocidade oblíqua, sobre o corpo atuam três forças externas: o peso, m p g , a resistência do ar, F r e a impulsão m a g , onde m p é a massa do projétil e m a a massa do ar que ocupava o mesmo volume do projétil. O problema é semelhante ao problema da queda livre, estudado na secção 4.3.3 do capítulo 4, mas a força de resistência do ar deixa de ser vertical (ver figura 10.7). O peso e a impulsão são verticais, em sentidos opostos, podendo ser combinados numa única força vertical (peso eficaz) de módulo ( m p m a ) g .

Admite-se que a massa volúmica do projétil é muito maior que a massa volúmica do ar e, portanto, o peso eficaz aponta para baixo e m p m a é quase igual a m p . De qualquer modo, a massa do projétil costuma medir-se medindo o seu peso eficaz no ar, assim que o valor medido ( m ) da massa do projétil é realmente m p m a e o peso eficaz é m g .

A força de resistência do ar muda constantemente de sentido, porque é sempre tangente à trajetória e no sentido oposto à velocidade. Como foi explicado no capítulo 4, no caso do ar o número de Reynolds costuma ser elevado e admite-se que a resistência do ar é proporcional ao quadrado da velocidade. Se o projétil é uma esfera de raio R , a expressão do módulo de F r é dada pela equação 4.14 e a força é:

(10.16)
F r = π 4 ρ R 2 v 2 e t

onde ρ é a massa volúmica do ar e e t é o versor tangencial, na direção e sentido do vetor velocidade:

(10.17)
e t = v v

Escolhendo um sistema de eixos em que a gravidade aponta no sentido negativo do eixo dos y e a velocidade inicial v 0 com que é lançado o projétil está no plano x y , o peso e a força de resistência do ar estão sempre no plano x y e o movimento do projétil dá-se nesse plano. Assim sendo, o vetor velocidade é ( v x ˆ ı + v y ˆ ) e a força de resistência do ar é:

(10.18)
F r = π 4 ρ R 2 v 2 x + v 2 y ( v x ˆ ı + v y ˆ )

O vetor do peso é m g ˆ . Aplicando a segunda lei de Newton, obtêm-se as componentes da aceleração:

(10.19)
a x = πρ R 2 4 m v x v 2 x + v 2 y a y = g πρ R 2 4 m v y v 2 x + v 2 y

Estas equações devem ser resolvidas em simultâneo porque as duas componentes v x e v y aparecem nas duas equações. É impossível encontrar a solução exata do problema, mas pode obter-se uma aproximação numérica.

A seguir vão-se comparar as trajetórias de duas esferas diferentes, lançadas com a mesma velocidade inicial para compará-las com a trajetória parabólica que teriam se pudessem ser lançadas no vácuo, sem resistência do ar. Considere-se o caso em que a velocidade inicial é 12 m/s, fazendo um ângulo de 45° com o plano horizontal; as componentes da velocidade inicial são,

(%i19) [vx0, vy0]: float (12*[cos(%pi/4),sin(%pi/4)])$

Começando pelo caso mais fácil, o lançamento dos projéteis no vácuo, as componentes da aceleração são a x = 0 e a y = 9 . 8 . O estado do projétil é ( x , y , v x , v y ) e a velocidade de fase ( v x , v y , a x , a y ). Os valores iniciais da velocidade já foram calculados em (%i19) e arbitre-se que o projétil parte da origem com valores iniciais nulos para x e y . Para integrar as equações de movimento desde t = 0 até t = 2  s, com incrementos de 0.01 s, usa-se o comando:

(%i20) tr1: rk ([vx,vy,0,-9.8], [x,y,vx,vy], [0,0,vx0,vy0], [t,0,2,0.01])$

e o último ponto calculado na lista tr1 é,

(%i21) last (tr1);
(%o21)   [2 . 0 , 16 . 97 , 2 . 629 , 8 . 485 , 11 . 11]

As 5 componentes do ponto são o tempo, as coordenadas da posição e as componentes da velocidade. Este resultado mostra que em t = 2 a bola já está a cair, porque v y é negativa e que já desceu debaixo da altura inicial, porque y é negativa.

Como pretendemos obter a trajetória até a bola regressar à altura y = 0 , é necessário extrair unicamente os pontos da lista tr1 com terceira componente ( y ) positiva. Percorre-se a lista toda, comparando o terceiro elemento de cada ponto com 0, até encontrar o primeiro ponto em que o terceiro elemento é negativo. Isso consegue-se usando o comando sublist_indices do Maxima:

(%i22) first (sublist_indices (tr1, lambda([p],p[3] < 0)));
(%o22)  175

Usou-se lambda que define um operador; neste caso esse operador compara o terceiro elemento da entrada que lhe for dada com zero. O comando sublist_indices percorre a lista tr1 passando cada elemento como entrada para esse operador e, nos casos em que o operador produz o resultado "true", o índice do respetivo elemento da lista é acrescentado a uma sub lista. O comando first seleciona apenas o primeiro elemento nessa sub lista, neste caso, o índice do primeiro ponto em que y é negativo. Como tal, só interessam os primeiros 174 pontos na lista; se o objetivo é construir o gráfico da trajetória, cria-se outra lista com as coordenadas x e y dos primeiros 174 pontos:

(%i23) r1: makelist ([tr1[i][2], tr1[i][3]], i, 1, 174)$

A seguir vai repetir-se o mesmo procedimento para uma bola de ténis e uma bola de ténis de mesa, tendo em conta a resistência do ar. A massa volúmica do ar é aproximadamente 1.2 kg/m3. É conveniente definir uma função que calcula a constante que aparece nas equações de movimento 10.19, em função do raio e a massa de cada uma das bolas; também é conveniente definir a expressão do módulo da velocidade para não ter que escrevê-la várias vezes:

(%i24) c(R,m) := -%pi*1.2*R^2/4/m$
(%i25) v: sqrt(vx^2+vy^2)$

Uma bola de ténis típica tem raio de aproximadamente 3.25 cm e massa 62 gramas. No comando (%i20) é necessário substituir a aceleração da gravidade pelas duas componentes da aceleração (equações 10.19)

(%i26) tr2: rk ([vx, vy, c(0.0325,0.062)*vx*v, -9.8+c(0.0325,0.062)*vy*v], [x,y,vx,vy], [0,0,vx0,vy0], [t,0,2,0.01])$

O primeiro ponto com altura negativa é

(%i27) first (sublist_indices (tr2, lambda([p],p[3] < 0)));
(%o27)  167

e a trajetória da bola de ténis armazena-se noutra variável:

(%i28) r2: makelist ([tr2[i][2],tr2[i][3]],i,1,166)$

Repetem-se os mesmos cálculos para uma bola de ténis de mesa típica, com raio 1.9 cm e massa 2.4 g

(%i29) tr3: rk ([vx, vy, c(0.019,0.0024)*vx*v, -9.8+c(0.019,0.0024)*vy*v], [x,y,vx,vy], [0,0,vx0,vy0], [t,0,2,0.01])$
(%i30) first (sublist_indices (tr3, lambda([p],p[3] < 0)));
(%o30)  133
(%i31) r3: makelist ([tr3[i][2],tr3[i][3]],i,1,132)$

O gráfico das 3 trajetórias constrói-se com o seguinte comando:

(%i32) plot2d ([[discrete, r1], [discrete, r2], [discrete, r3]], [xlabel, "x (m)"], [ylabel, "y (m)"], [y, 0, 12], [legend, "vacuo", "tenis", "tenis de mesa"])$

O resultado é apresentado na figura 10.8.

Trajetórias de 3 projéteis
Figura 10.8: Trajetórias de uma bola no vácuo e bolas de ténis e ténis de mesa no ar.

As trajetórias das bolas no ar não são parábolas, mas no fim curvam-se mais e terminam com uma queda mais vertical. O efeito da resistência do ar é mais visível na bola de ténis de mesa; apesar de ser mais pequena que a bola de ténis, a força de resistência do ar produz nela maior aceleração tangencial negativa, devido à sua menor massa volúmica. Lançadas com a mesma velocidade, o alcance horizontal da bola de ténis de mesa é 6.2 m e o da bola de ténis 12.4 m. O alcance horizontal hipotético das duas bolas, se a resistência do ar pudesse ser ignorada, seria 14.7 m.

10.4.3. Pêndulo de Wilberforce

Pêndulo de Wilberforce

Figura 10.9: Pêndulo
de Wilberforce.

O pêndulo de Wilberforce (figura 10.9) é constituído por um cilindro pendurado de uma mola vertical muito comprida. Quando uma mola é esticada ou comprimida, cada espira muda ligeiramente de tamanho; no pêndulo de Wilberforce, o número elevado de espiras na mola faz com que seja mais visível essa mudança, de forma que enquanto a mola oscila, também se enrola ou desenrola, fazendo rodar o cilindro em relação ao eixo vertical.

O sistema tem dois graus de liberdade, a altura z do centro de massa do cilindro e o ângulo de rotação do cilindro à volta do eixo vertical, θ . Se z = 0 e θ = 0 são escolhidos na posição de equilíbrio, é possível ignorar a energia potencial gravítica que poderá ser eliminada das equações com uma mudança de variáveis (problema 4 do capítulo 9). A energia potencial elástica tem 3 termos, que dependem da elongação da mola z e do seu ângulo de rotação θ ; as energias cinética e potencial são,

(10.20)
E c = 1 2 m ˙ z 2 + 1 2 I cm ˙ θ 2 U = 1 2 k z 2 + 1 2 a θ 2 + b z θ

em que k , a e b são constantes elásticas da mola. As equações de Lagrange, ignorando a resistência do ar e outras forças dissipativas, conduzem às seguintes equações de movimento:

(10.21)
¨ z = k m z b m θ ¨ θ = a I cm θ b I cm z

Para resolver as equações de evolução numericamente, é necessário dar alguns valores típicos para a massa, o momento de inércia e as constantes elásticas,

(%i33) [m, I, k, a, b]: [0.5, 1e-4, 5, 1e-3, 0.5e-2]$

A solução no intervalo de tempo desde 0 até 40, com condição inicial z = 10  cm e as outras variáveis iguais a 0, obtém-se com o seguinte comando:

(%i34) sol: rk(['v,w,-(k*z+b*ang)/m,-(a*ang+b*z)/I], [z,ang,'v,w],[0.1,0,0,0],[t,0,40,0.01])$

A figura 10.10 mostra o gráfico obtido para o ângulo θ e a elongação z , multiplicada por um fator de 100 para que seja visível na mesma escala do ângulo.

Elongação e ângulo  vs tempo no pêndulo de Wilberforce
Figura 10.10: Elongação e ângulo de rotação no pêndulo de Wilberforce.

O gráfico mostra uma caraterística interessante do pêndulo de Wilberforce: se o pêndulo é posto a oscilar, sem rodar, a amplitude das oscilações lineares decresce gradualmente, enquanto que o cilindro começa a rodar com oscilações de torção que atingem uma amplitude máxima quando o cilindro deixa de se deslocar na vertical. A amplitude das oscilações de torção começa logo a diminuir à medida que a oscilação linear cresce novamente. Essa intermitência entre deslocamento vertical e rotação repete-se indefinidamente.

A projeção do retrato de fase nas variáveis z e θ é apresentada na figura 10.11.

Ãngulo vs elongação no pêndulo de Wilberforce
Figura 10.11: Retrato de fase no plano formado pela elongação e o ângulo.

Neste sistema existem duas frequências angulares. A frequência angular longitudinal e a frequência angular de torção,

(10.22)
2 z = k m 2 θ = a I cm

O cilindro num pêndulo de Wilberforce costuma ter quatro porcas que podem ser deslocadas, aumentando ou diminuindo o momento de inércia, para conseguir que as duas frequências sejam muito parecidas e o efeito de alternância entre oscilações lineares e rotacionais seja mais visível. Os valores dos parâmetros usados no exemplo acima, foram escolhidos de forma a garantir duas frequências iguais.

Perguntas

(Para conferir a sua resposta, clique nela.)

  1. O valor aproximado do período de um pêndulo com comprimento l é 2 π l / g , onde g é a aceleração da gravidade. Essa expressão é uma boa aproximação unicamente em algumas situações. Se o ângulo θ é zero no ponto de equilíbrio estável, qual das condições seguintes garante que essa expressão seja uma boa aproximação do seu valor real?

    1. valor máximo da velocidade angular pequeno.
    2. aceleração da gravidade pequena.
    3. comprimento l pequeno.
    4. valores máximos do ângulo e da velocidade angular pequenos.
    5. valor máximo do ângulo pequeno.
  2. A força tangencial numa partícula com velocidade v e posição na trajetória s é: F t = 4 s ( s v 2 ) . Quantos pontos de equilíbrio tem esse sistema?
    1. 1
    2. 2
    3. 3
    4. 4
    5. 0
  3. Qual é a matriz jacobiana do sistema ˙ x = y 2 , ˙ y = x y ?
    1. y 2 11 x y
    2. 02 y 11
    3. 02 y y x
    4. y x 02 y
    5. 11 0 2 y
  4. As equações de evolução de um sistema dinâmico no espaço de fase ( x , y ), são ˙ x = x y , ˙ y = y + 1 . Qual dos seguintes vetores aponta na direção e sentido da velocidade de fase em (1, 2)?
    1. 4ˆ ı + 2ˆ
    2. 2ˆ ı + 4ˆ
    3. 6ˆ ı + 4ˆ
    4. 4ˆ ı + 6ˆ
    5. 2ˆ ı 3ˆ
  5. No retrato de fase na figura, que tipo de ponto de equilíbrio é o ponto (1,0)? Retrato de fase de um sistema com nó e ponto de sela
    1. nó atrativo
    2. foco repulsivo
    3. ponto de sela
    4. foco atrativo
    5. nó repulsivo

Problemas

  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)
  2. Em cada um dos casos seguintes encontre os pontos de equilíbrio e os valores próprios da matriz jacobiana nesses pontos e identifique os tipos de pontos de equilíbrio:
    (a) ˙ x = y 2 + 3 y 10    ˙ y = x y + x + 12
    (b) ˙ x = 3 x y 2 2 y    ˙ y = x y 2
    (c) ˙ x = y 2 + 2 x y + 2    ˙ y = x 2 y 2 2
    (d) ˙ x = x + 4 y y 3    ˙ y = y + 4 x x 3
  3. O diagrama mostra o retrato de fase de um sistema com unicamente 3 pontos de equilíbrio, no caso idealizado em que não existe atrito. Faça (a mão) um esboço da energia potencial e de como seria o retrato de fase do sistema real, considerando as forças de atrito. Retrato de fase de um potencial quártico
  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.
  5. A base do pêndulo da figura 10.2 roda no plano horizontal, com velocidade angular constante ω b , enquanto o pêndulo oscila.
    (a) Demonstre que a equação de movimento é:
    ¨ θ = 1 l sin θ r ω 2 b cos θ g
    onde r é a distância do centro de massa até o eixo e o comprimento eficaz l é o raio de giração ao quadrado, sobre r .
    (b) Trace o gráfico de sin θ r ω 2 b cos θ g em função de θ , entre π e π , para um pêndulo com r = 0 . 3  m e ω b = 2  s−1. Repita o gráfico para ω b = 8  s−1. Com base nos dois gráficos, identifique em cada caso os pontos de equilíbrio estável e instável.
    (c) Demonstre que quando ω b < g / r , existe um único ponto de equilíbrio estável em θ = 0 e um único ponto de equilíbrio instável em θ =± π .
    (d) Se ω b > g / r , demostre que os pontos de equilíbrio em θ = 0 e θ =± π são ambos instáveis e aparecem dois pontos de equilíbrio estável em ± θ 0 , onde θ 0 é um ângulo entre zero e π /2 .
  6. Na trajetória da bola de ténis de mesa calculada na secção 10.4.2, o alcance horizontal da bola é aproximadamente o valor da coordenada x do último ponto da lista de pontos r1. Repita os cálculos, com diferentes valores do ângulo de lançamento, para determinar os valores do alcance com ângulos de 35°, 36°, 37°, 38°, 39° e 40°. Registe numa tabela os valores obtidos para o alcance horizontal, em função do ângulo, com precisão até os milímetros. Com base na tabela, qual é o ângulo de lançamento que produz o maior alcance horizontal? Usando o resultado do problema 12 do capítulo 6, mostre que no vácuo o ângulo que produz o alcance máximo é 45°.
  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 .
  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.
  9. Uma partícula de massa m desloca-se no plano x y sob a ação de uma força conservativa com energia potencial,
    U = k x 2 x 2 + k y 2 y 2
    onde k x e k y são duas constantes positivas. As trajetórias da partícula obtidas com diferentes valores dessas constantes chamam-se figuras de Lissajous.
    (a) Encontre as duas equações de movimento para ¨ x e ¨ y
    (b) Resolva numericamente as equações de movimento, no caso m = 0 . 3 , k x = 2 e k y = 8 (unidades SI), entre t = 0 e t = 2 . 43 , se a partícula partir do ponto (1, 0) com velocidade inicial v = 0 . 6ˆ . Trace o gráfico da trajetória da partícula no plano x y .
    (c) Repita a alínea anterior, mas admitindo que a partícula parte do ponto (1, 0) com velocidade inicial v = 0 . 3ˆ ı + 0 . 6ˆ .
    (d) Observe que o sistema pode ser considerado como um conjunto de dois osciladores harmónicos independentes, nas direções x e y . Calcule o período de oscilação para cada um dos dois osciladores e diga qual é a relação entre os dois períodos.
    (e) Repita os cálculos da alínea c, mudando o valor de k y para 18. Que relação observa entre os gráficos da trajetória e k y / k x ?
  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).

Respostas

Perguntas: 1. D. 2. A. 3. C. 4. D. 5. E.

Problemas

  1. (a) Unicamente um centro em ( x , v x ) = (0, 0).
    Retrato de fase (b) Um ponto de sela em ( x , v x )=(0, 0), um foco instável em ( x , v x )=(-1, 0) e um foco estável em ( x , v x )=(1, 0).
    Retrato de fase
  2. (a) Dois pontos de equilíbrio: (3, -5), com valores próprios 7 e -4, é ponto de sela; (-4, 2), com valores próprios 3 e -7 é ponto de sela.
    (b) Dois pontos de equilíbrio: (0, 0), com valores próprios ± i 2 é centro; (0.763, 0.874), com valores próprios -2.193 e 2.736 é ponto de sela.
    (c) Dois pontos de equilíbrio: ( 2 6/3 , 6/3 ) e ( 2 6/3 , 6/3 ), ambos pontos de sela com valores próprios ± 2 2 .
    (d) Nove pontos de equilíbrio. Um ponto de sela em (0, 0), com valores próprios 3 e -5, outros dois pontos de sela em ( 5 , 5 ) e ( 5 , 5 ), com valores próprios 10 e -12, outros dois pontos de sela em ( 3 , 3 ) e ( 3 , 3 ) com valores próprios 4 e -6 e quatro focos atrativos em ( b a , a ), ( b a , a ), ( a b , b ) e ( a b , b ), com valores próprios 1 +± i 23 , onde a = 2 + 3 e b = 2 3 .
  3. Os pontos de sela continuam sendo pontos de sela e o centro passa a ser foco estável. Potencial quártico e retrato de fase
  4. (a) O pêndulo oscila com amplitude que decresce lentamente:
    Ângulo e velocidade angular do pêndulo amortecido Curva de evolução do pêndulo amortecido
    (b) 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:
    Ângulo e velocidade angular do pêndulo amortecido Curva de evolução do pêndulo amortecido
  5. (b)
    Forca do pêndulo com velocidade angular igual a 2 Forca do pêndulo com velocidade angular igual a 8
    Com ω b = 2  s−1, há um ponto de equilíbrio estável em θ = 0 e um ponto de equilíbrio instável em θ =± π . Com ω b = 8  s−1, há dois pontos de equilíbrio instável em θ = 0 e θ =± π e dois pontos de equilíbrio estável em θ 1 e θ 1 .
  6. ÂnguloAlcance (m)
    35° 6.293
    36° 6.299
    37° 6.301
    38° 6.299
    39° 6.325
    40° 6.314
    O ângulo de 37° produz o alcance máximo. No problema 12 do capítulo 6, o valor máximo do seno é 1, quando 2 θ = 90 ; e, portanto, θ = 45 .
  7. (a) ˙ x = v , ˙ v = 4 v 2 4 x 2  (b) ( x , ˙ x ) = (1, 0) e ( x , ˙ x ) = (-1, 0)   (c) J = 01 8 x 2˙ x   (d) (1, 0) é um centro e (-1, 0) é ponto de sela.  (e) x = 0 . 5869 , ˙ x = 0 . 8277 .
  8. (a) No eixo dos x , ˆ ı . No eixo dos y , ˆ . (b) Na reta y = x , (ˆ ı ˆ )/ 2 . Na reta y = x , ( ˆ ı + ˆ )/ 2 . (c) Ver figura; a origem é ponto de sela. (d) Nenhum ciclo nem órbita heteroclínica; número infinito de órbitas homoclínicas (todas as curvas de evolução no primeiro e terceiro quadrantes). Folha de Descartes
  9. (a) ¨ x = k x m x ¨ y = k y m y
    (b) e (c)
    Figura de Lissajous Figura de Lissajous
    (d) Na direção x , 2.433 s. Na direção y , 1.217 s. O período na direção x é o dobro do período na direção de y . (e) Se k y / k x for um número inteiro, o estado da partícula regressa ao estado inicial depois de descrever uma figura de Lissajous com k y / k x loops segundo o eixo dos x . Figura de Lissajous
  10. (a) ¨ x = 4 π 2 x ( x 2 + y 2 ) 3/2 ¨ y = 4 π 2 y ( x 2 + y 2 ) 3/2
    (b) e (c)
    Órbita do cometa Halley, com erro numérico Órbita do cometa Halley
    Na alínea b o erro numérico é muito elevado; a energia do cometa não permanece constante mais diminui. Na alínea c 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) 34.4 UA. A órbita sai por fora da órbita de Neptuno, e entra até um ponto entre órbitas de Mercúrio e Vénus.
Pergunta 1, resposta A: Errada

Velocidade máxima pequena não implica que o pêndulo esteja próximo do ponto de equilíbrio estável, porque num pêndulo com l elevado,

(clique para continuar)

Pergunta 1, resposta B: Errada

(clique para continuar)

Pergunta 1, resposta C: Errada

(clique para continuar)

Pergunta 1, resposta D: Certa

(clique para continuar)

Pergunta 1, resposta E: Errada

(clique para continuar)

Pergunta 2, resposta A: Certa

(clique para continuar)

Pergunta 2, resposta B: Errada

(clique para continuar)

Pergunta 2, resposta C: Errada

(clique para continuar)

Pergunta 2, resposta D: Errada

(clique para continuar)

Pergunta 2, resposta E: Errada

(clique para continuar)

Pergunta 3, resposta A: Errada

(clique para continuar)

Pergunta 3, resposta B: Errada

(clique para continuar)

Pergunta 3, resposta C: Certa

(clique para continuar)

Pergunta 3, resposta D: Errada

(clique para continuar)

Pergunta 3, resposta E: Errada

(clique para continuar)

Pergunta 4, resposta A: Errada

(clique para continuar)

Pergunta 4, resposta B: Errada

(clique para continuar)

Pergunta 4, resposta C: Errada

(clique para continuar)

Pergunta 4, resposta D: Certa

(clique para continuar)

Pergunta 4, resposta E: Errada

(clique para continuar)

Pergunta 5, resposta A: Errada

(clique para continuar)

Pergunta 5, resposta B: Errada

(clique para continuar)

Pergunta 5, resposta C: Errada

(clique para continuar)

Pergunta 5, resposta D: Errada

(clique para continuar)

Pergunta 5, resposta E: Certa

(clique para continuar)