unesp UNIVERSIDADE ESTADUAL PAULISTA Instituto de Biociências, Letras e Ciências Exatas DEPARTAMENTO DE CIÊNCIAS DE COMPUTAÇÃO E ESTATÍSTICA ANÁLISE DA DINÂMICA DE UM SISTEMA VIBRANTE NÃO IDEAL DE DOIS GRAUS DE LIBERDADE Luiz Oreste Cauz Dissertação de Mestrado Pós-Graduação em Matemática Aplicada MAP - 097 Rua Cristóvão Colombo, 2265 15054-000 - São José do Rio Preto - SP - Brasil Telefone: (017) 221-2444 Fax: (017) 221-2445 Análise da Dinâmica de um Sistema Vibrante não Ideal de Dois Graus de Liberdade Luiz Oreste Cauz 1 Dissertação apresentada ao Instituto de Biociências, Letras e Ciências Exatas da Universidade Estadual Paulista “Júlio de Mesquita Filho”, Campus de São José do Rio Preto, São Paulo, para a obtenção do t́ıtulo de Mestre em Matemática Aplicada. Orientador: Prof. Dr. Masayoshi Tsuchida São José do Rio Preto 2005 1contato: lcauz@yahoo.com.br i À Deus. À minha famı́lia. Aos meus amigos. Dedico ii Agradecimentos À Deus, por tudo. À minha famı́lia que sempre me apoiou. Em especial ao Prof. Dr. Masayoshi Tsuchida, pela amizade, orientação , incentivo e paciência na elaboração deste trabalho. Aos meus professores de graduação e pós-graduação. Aos meus amigos de república e da minha cidade pelos momentos de descontração. Aos meus amigos pós-graduandos do DCCE, do Departamento de Matemática do IBILCE. À todas as pessoas e funcionários do IBILCE que, direta ou indiretamente, contribúıram para a elaboração deste trabalho. Ao CNPq pelo aux́ılio finaceiro. iii Resumo Neste trabalho apresentamos um estudo da dinâmica de um sistema vibrante não ideal, composto por um motor e uma mola, conhecido como vibrador centŕıfugo. O objetivo deste estudo é mostrar a diferença de comportamento do sistema, quando consideramos molas duras (coeficiente de elasticidade cúbica positivo) ou molas suaves (coeficiente de elasticidade cúbica negativo). Para mola dura foi analisada a estabilidade dos pontos de equiĺıbrio, e mostrada por meio da teoria de variedade central e do teorema de Bezout a existência da bifurcação de Hopf. Para mola suave, é mostrada a existência de uma órbita heterocĺınica conectando dois pontos de sela. Usando o método clássico de Melnikov, é discutida a existência ou não do comportamento caótico para um determinado ńıvel de energia e para certos valores do coeficiente de amortecimento. Toda a análise é acompanhada de simulações numéricas para a confirmação dos resultados. Palavras-chave: vibrador centŕıfugo, bifurcação de Hopf, função de Melnikov, efeito Sommerfeld. iv Abstract In this work we present a study of the dynamics of a non-ideal vibrating system, composed by a motor and a spring, which is known as centrifugal vibrator. The purpose of this study is to show the difference of behavior of the system when we consider hard springs (positive coefficient of cubical elasticity) or soft springs (negative coefficient of cubical elasticity). For hard spring the stability of the fixed point was analyzed, and by means of the Central Manifolds Theory and the Bezout theorem the existence of the Hopf Bifurcation is shown. For soft spring, it is shown the existence of a heteroclinic orbit connecting two saddle points. Using the classical Melnikov method it is discussed the existence, or not, of the chaotic behavior for some energy level and certain values of the damping coefficient. All the analysis is followed by numerical simulations to confirm the results. Keywords: centrifugal vibrator, Hopf bifurcation, Melnikov function, Sommerfeld effect. v Sumário Lista de Tabelas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii Lista de Figuras . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xi 1 Introdução 1 1.1 Objetivos do Trabalho . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.2 Apresentação dos Caṕıtulos . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2 Conceitos preliminares: Tópicos de Sistemas Dinâmicos 6 2.1 Sistemas lineares . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.2 Classificação do equiĺıbrio quanto à topologia e à estabilidade . . . . . . . . . 8 2.3 Critério de Routh-Hurwitz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.4 Sistemas não-lineares . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.5 Equivalência topológica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.6 Teoremas locais para sistemas não-lineares . . . . . . . . . . . . . . . . . . . . 12 2.6.1 Teorema de Hartman-Grobman . . . . . . . . . . . . . . . . . . . . . . 12 2.6.2 Teorema das variedades hiperbólicas . . . . . . . . . . . . . . . . . . . 12 2.7 Teoria da variedade central . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.7.1 Teoremas da variedade central . . . . . . . . . . . . . . . . . . . . . . . 14 2.7.2 Teoremas de Carr . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.8 Bifurcações em sistemas dinâmicos de tempo cont́ınuo . . . . . . . . . . . . . . . 16 2.8.1 Bifurcação de Hopf . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.9 Mapas e seção de Poincaré . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.10 Método de Melnikov . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.10.1 Caso autônomo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.10.2 Caso não-autônomo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.11 O mapa de ferradura de Smale . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 vi Sumário vii 2.12 Teorema de Bezout . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3 O Sistema Mecânico e Algumas Simulações Numéricas. 26 3.1 O Vibrador Centŕıfugo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 3.2 Ponto de Equiĺıbrio do Sistema . . . . . . . . . . . . . . . . . . . . . . . . . . 29 3.3 Equações de Movimento Usando a Formulação de Hamilton . . . . . . . . . . . 32 3.4 Resultados Numéricos: Pontos de equiĺıbrio . . . . . . . . . . . . . . . . . . . 34 3.5 Resultados Numéricos: Estudo da ressonância 1:1 . . . . . . . . . . . . . . . . 43 4 Bifurcação de Hopf no Vibrador Centŕıfugo Não Ideal 51 4.1 Estudo do Caso Cŕıtico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 4.2 Estudo de um Caso Particular . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 4.3 Resultados Numéricos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 4.3.1 Bifurcação de Hopf para um caso particular . . . . . . . . . . . . . . . 60 5 Análise da dinâmica do sistema através do método de Melnikov 64 5.1 Formulação do Problema num Campo Hamiltoniano Perturbado . . . . . . . . 64 5.2 Resultados Numéricos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 6 Conclusões 79 A O Método de Redução de Holmes e Marsden 81 A.1 O Método de Redução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 B Teoremas de Wiggins e Shaw 86 Referências Bibliográficas 90 Lista de Tabelas 3.1 Pontos de equiĺıbrio . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 3.2 Parâmetros f́ısicos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 4.1 Autovalores para diferentes valores de B. . . . . . . . . . . . . . . . . . . . . . 61 viii Lista de Figuras 2.1 Variedades de um ponto de sela de um sistema não-linear . . . . . . . . . . . . . . 13 2.2 Significado geométrico da aplicação de Poincaré P : x0 é uma órbita periódica, x∗ 0 é o ponto onde ela intercepta Σ; x1 é uma órbita não periódica que corta Σ em P (x′ 1) 6= x′ 1. 19 2.3 Seções de Poincaré : (a) unidimensional; (b) bidimensional . . . . . . . . . . . . . 20 2.4 Caso ε = 0: linha tracejada; caso ε 6= 0: linha cheia. . . . . . . . . . . . . . . . . . 21 2.5 Emaranhado homocĺınico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.6 Esquema do mapa de ferradura . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.7 Iterações de f−1 sobre o quadrado D . . . . . . . . . . . . . . . . . . . . . . . . 25 3.1 Vibrador centŕıfugo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 3.2 Modelo de torque linear com A= 0.01 e condições inicial (y1, y2, y3, y4) = (0, 0, 0.001, 0). (a) Oscilação da mola; (b) velocidade de vibração da mola; (c) ângulo da massa m; (d) velocidade angular da massa m. . . . . . . . . . . . . . . . . . . . . . . . . . 36 3.3 Modelo de torque exponencial, com A= 0.01 e condição inicial (y1, y2, y3, y4) = (0, 0, 3.001, 0). (a) Oscilação da mola; (b) velocidade de vibração da mola; (c) ângulo da massa m; (d) velocidade angular da massa m. . . . . . . . . . . . . . . . . . . 37 3.4 Modelo de torque exponencial, com A= 0.01 e condição inicial (y1, y2, y3, y4) = (0, 0, 0.001, 0). (a) Deslocamento da mola; (b) velocidade de vibração da mola; (c) ângulo da massa m; (d) velocidade angular da massa m. . . . . . . . . . . . . . . 38 3.5 Comportamento do movimento da massa m para A= 0.1. (a) Modelo de torque linear e condição inicial (y1, y2, y3, y4) = (0, 0, 0.001, 0). (b) Modelo de torque exponencial e condição inicial (y1, y2, y3, y4) = (0, 0, 3.001, 0). . . . . . . . . . . . . . . . . . . . 39 3.6 Comportamento do movimento da massa m para A= 0.2. (a) Modelo de torque linear e condição inicial (y1, y2, y3, y4) = (0, 0, 0.001, 0). (b) Modelo de torque exponencial e condição inicial (y1, y2, y3, y4) = (0, 0, 3.2, 0). . . . . . . . . . . . . . . . . . . . . 39 ix Lista de Figuras x 3.7 Comportamento do movimento da massa m para A= 0.3. (a) Modelo de torque linear e condição inicial (y1, y2, y3, y4) = (0, 0, 3, 0). (b) Modelo de torque exponencial e condição inicial (y1, y2, y3, y4) = (0, 0, 3.4, 0). . . . . . . . . . . . . . . . . . . . . . 40 3.8 Modelo de torque linear. A= 0.980665. Condições iniciais:(y1, y2, y3, y4) = (0, 0, 4, 0). (a) Deslocamento da mola. (b) velocidade de vibração da mola. (c) Ângulo da massa. (d) Velocidade angular . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 3.9 Modelo de torque linear. A= 1.2. Condições iniciais:(y1, y2, y3, y4) = (0, 0, 4, 0). (a) Deslocamento da mola. (b) velocidade de vibração da mola. (c) Ângulo da massa. (d) Velocidade angular . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 3.10 Caracteŕısticas de movimento para parâmetro de controle a = 0.61. Caso ϕ̇ < ω . . 44 3.11 Caracteŕısticas de movimento para parâmetro de controle a = 0.67. Caso ϕ̇ ≈ ω . . 45 3.12 Caracteŕısticas de movimento para parâmetro de controle a = 0.72. Caso ϕ̇ > ω . . 46 3.13 Caracteŕısticas de movimento para parâmetro de controle a = 0.679. Aparecimento de modulações logo após a passagem pela ressonância. . . . . . . . . . . . . . . . . 47 3.14 Caracteŕısticas de movimento para parâmetro de controle a = 0.685. As modulações persistem mais à medida que se afasta da passagem pela ressonância. . . . . . . . . 48 3.15 Caracteŕısticas de movimento para parâmetro de controle a = 0.69. As modulações ocorrem praticamente em todo o intervalo de tempo de simulação. . . . . . . . . . 49 3.16 Amplitude máxima de oscilação da mola em função da freqüência do rotor. Efeito Sommerfeld. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 4.1 Espaços de fases (y1, y2) para diferentes valores de B . . . . . . . . . . . . . . . . 61 4.2 Espaços de fases (y3, y4) para diferentes valores de B . . . . . . . . . . . . . . . . 62 4.3 Diagrama de bifurcação da variável y1 . . . . . . . . . . . . . . . . . . . . . . . 62 4.4 Ampliação da figura 4.3 em torno da origem. . . . . . . . . . . . . . . . . . . . . 63 4.5 Função torque para o caso particular M(y4) = −By4. As retas 1 (B = 0.00005) e 2 (B = 0.00002) correspondem a M ′(y4) < 0, a reta 3 é o caso cŕıtico M ′(y4) = 0 e as retas 4 (B = −0.00002) e 5 (B = −0.00005) referem-se a M ′(y4) > 0. . . . . . . . . 63 5.1 Espaço de fases e a variedade invariante M. . . . . . . . . . . . . . . . . . . . . 70 5.2 Retrato de fases do sistema não perturbado (5.7). Os pontos de sela são (±3.162, 0). 75 5.3 Variedade instável do ponto de sela (−3.162, 0). (a) Caso ε = 0. (b) Caso ε = 0.00001. 75 5.4 Função ∆(ϕ0) para β = 0.1 e s0 = 0.5. Os Zeros são simples. . . . . . . . . . . . . 76 Lista de Figuras xi 5.5 Função ∆(ϕ0) para β = 0.1 e s0 = 0.2325. Os Zeros são duplos. . . . . . . . . . . 76 5.6 Função ∆(ϕ0) para β = 0.1 e s0 = 0.9636. Os Zeros são duplos. . . . . . . . . . . 77 5.7 Função ∆(ϕ0) para β = 0.1 e s0 = 1.0. Não há zeros. . . . . . . . . . . . . . . . . 77 5.8 Função ∆(ϕ0) para β = 0.2006780233 e s0 = 0.5203690288. Os Zeros são duplos. . 78 5.9 Função ∆(ϕ0) para β = 0.21 e s0 = 0.5. Não há zeros. . . . . . . . . . . . . . . . . 78 Caṕıtulo 1 Introdução Nesta dissertação, estudamos a dinâmica de um sistema mecânico conhecido como vi- brador centŕıfugo, que foi inicialmente estudado por Kononenko [10] e recentemente por Dantas e Balthazar [3, 4]. Esse sistema é constitúıdo por um motor que gira uma massa desbalanceada, e é sustentado por uma mola e um amortecedor fig(3.1). As equações de movimento são dadas em [10] pelas expressões    m0ẍ + βẋ + cx = mrϕ̇2 cos ϕ + mrϕ̈senϕ, Iϕ̈ = M(ϕ̇) + mrẍsenϕ + mgrsenϕ, (1.1) onde m0 = m1 +m com m1 a massa do motor e m a massa excêntrica. Ainda, β é o coeficiente de amortecimento, c é a constante linear da mola, I é o momento de inércia da parte giratória do motor e M(ϕ̇) é a diferença entre o torque aplicado ao rotor e o torque de resistência à rotação do rotor. Nos seus estudos, Kononenko considera também que a mola possua uma elasticidade não linear ([10], pág. 66). A variável x representa a oscilação do motor, e ϕ é o ângulo de rotação da massa m. Das equações (1.1) notamos que as variáveis x e ϕ estão acopladas, e isso indica que a rotação do motor e a vibração da mola se interagem mutuamente. O efeito da vibração da mola sobre a rotação do motor se deve ao fato da potência do motor ser limitada, e nesse caso dizemos que o sistema dinâmico é não ideal. No caso ideal, o agente perturbador possui potência ilimitada e não é afetado pelo comportamento do sistema dinâmico. Um dos objetivos de Kononenko é mostrar a interação do sistema oscilante com a fonte de perturbação quando estes estão em condição de ressonância, ou seja, quando a freqüência de rotação do motor está muito próxima da freqüência natural do sistema. Para isso, Kononenko impôs algumas restrições ao sistema, considerando que a força de amortecimento βẋ, bem 1 Caṕıtulo 1 2 como o torque resultante M(ϕ̇) e os termos mrϕ̇2 cos ϕ, mrϕ̈senϕ, mrẍsenϕ e mgrsenϕ sejam pequenos quando comparados com as outras forças que agem sobre o sistema. Com estas restrições o sistema (1.1) foi reescrito como    ẍ + ω2x = ε(q2ϕ̇ 2 cos ϕ + q2ϕ̈senϕ − hẋ), ϕ̈ = ε[M1(ϕ̇) + q3ẍsenϕ], (1.2) onde ω2 = c m0 , εq2 = mr m0 , εq3 = mr I , εM1(ϕ̇) = M(ϕ̇) I e εh = β m0 . Note que se ε = 0 em (1.2), então ϕ = Ωt + ϕ(0) e assim temos a freqüência do motor dϕ dt = Ω = constante. O termo mgrsenϕ é desconsiderado por Kononenko, pois ele propõe um modelo para estudar os efeitos de ressonância no vibrador centŕıfugo, e a constante g acarreta uma posição de equiĺıbrio na variável angular ϕ. Pelas restrições impostas temos que m m0 � 1, m I � 1, M(ϕ̇) I � 1, β m0 � 1 e conseqüen- temente podemos fazer um estudo do sistema quando ω − ϕ̇ = εα0. Introduzindo a mudança de variável x = A cos(ϕ + Ξ), dx dt = −Aωsen(ϕ + Ξ), dϕ dt = Θ o sistema (1.2) é reescrito na forma dΘ dt = ε[M1(Θ) − Aq3ωΘ cos(ϕ + Ξ)senϕ], dA dt = − ε ω [q2Θ 2 cos ϕ + Aωhsen(ϕ + Ξ)]sen(ϕ + Ξ) + ε2 · · · , dΞ dt = ε{α0 − 1 Aω [q2Θ 2 cos ϕ + Aωhsen(ϕ + Ξ)] cos(ϕ + Ξ)} + ε2 · · · , dϕ dt = Θ. (1.3) Para resolver o sistema (1.3) foi usada uma técnica conhecida como método da média de Bogoliubov [10], considerando que as soluções em primeira aproximação das três primeiras equações do sistema (1.3) são dadas na forma Θ = Ω + εU1(t, Ω, a, ξ), A = a + εU2(t, Ω, a, ξ), Ξ = ξ + εU3(t, Ω, a, ξ), onde Ui, i = 1, 2, 3 são funções periódicas no tempo. Caṕıtulo 1. Objetivos do Trabalho 3 A integração da quarta equação de (1.3) em primeira aproximação tem a forma ϕ = Ωt + εΦ(t). Assim, as três primeiras equações são substitúıdas por dΘ dt = ε[M1(Ω) + 1 2 q3aωΩsenξ] da dt = − ε 2 ( ha + q2 Ω2 ω senξ ) dξ dt = ε ( α − q2Ω 2 2aω cos ξ ) (1.4) onde εα = ω − Ω. Deste modo, sob a condição de regime estacionário do sistema anterior (1.4), pode-se observar a amplitude a, fase de oscilação ξ e a freqüência de rotação Θ do motor. Em resumo, Kononenko estabeleceu condições para se estudar a interação entre o sistema oscilante e a fonte de excitação (motor DC) quando estes estão em condição de ressonância (ϕ̇ ≈ ω). É importante ressaltar que Kononenko não leva em consideração a aceleração da gravi- dade g atuando sobre o sistema. Em [3], esta constante toma um papel fundamental na obtenção de um caso particular (bifurcação de Hopf) por meio da eq. (4.26), onde a6 = mgr I . Resultados de análise da dinâmica do sistema vibrante (vibrador centŕıfugo) realizados em [3] e [4] serão apresentados nos caṕıtulos (3), (4) e (5). Alguns sistemas similares são estudados em [5] e [16]. Em [5] apresenta-se uma análise de um sistema dinâmico composto por um motor que gira uma massa desbalanceada, e está montado sobre uma mesa que pode oscilar horizontal- mente. Esta mesa é acoplada lateralmente por amortecedor e mola. É considerado que a fonte de energia é não ideal. Foi analisada a interação entre a dinâmica da estrutura oscilante e a dinâmica do motor. Este trabalho baseou-se nas caracteŕısticas do motor já que na maioria dos trabalhos encontrados na literatura leva-se em conta somente as caracteŕısticas das curvas do torque aplicado ao rotor e o torque de resistência à rotação. Em [16] é apresentado um estudo de um sistema composto por um bloco com um motor e um pêndulo acoplados, e que está suspenso por um amortecedor e uma mola. Neste trabalho é assumido que a caracteŕıstica da curva de torque é linear (M = u1 + u2ϕ̇). Um estudo numérico foi realizado tomando u2 fixo, e u1 como um parâmetro de controle. 1.1. Objetivos do Trabalho 4 1.1 Objetivos do Trabalho O objetivo deste trabalho é o estudo numérico da dinâmica do sistema vibrante composto por um motor DC (motor elétrico de corrente cont́ınua) e uma mola com não linearidade cúbica fig.(3.1). A mola sustenta o motor que gira uma pequena massa a uma distância fixa do eixo de rotação. É assumido que esse motor opera com potência limitada, de modo que temos uma influência do comportamento da mola sobre a rotação do motor, caracterizando assim um sistema mecânico não ideal. Este sistema é conhecido como vibrador centŕıfugo. Os principais resultados que serão apresentados é devido a [3] e [4], onde é mostrado um estudo anaĺıtico sobre a dinâmica do vibrador centŕıfugo com molas duras e molas com menor rigidez. Algumas simulações numéricas foram realizadas com o objetivo de constatar estes resul- tados. Outras simulações também foram realizadas. 1.2 Apresentação dos Caṕıtulos Esta dissertação foi dividida em 6 caṕıtulos. No caṕıtulo 2 apresentamos alguns tópicos em sistemas dinâmicos, os quais julgamos importantes para a compreensão dos caṕıtulos seguintes. O texto neste caṕıtulo é apresentado de forma introdutória. Não nos preocupamos em enunciar teoremas e definições, e sim discorrer sobre o assunto de maneira informal. Enun- ciamos apenas alguns teoremas que foram utilizados nos caṕıtulos posteriores. No caṕıtulo 3 é apresentado o sistema vibrante estudado neste trabalho, bem como suas equações de movimento e algumas simulações numéricas. As equações de movimento são apresentadas de duas formas, uma usando o formalismo de Lagrange [3] e outra usando o formalismo Hamiltoniano [4]. Nas simulações numéricas é mostrada a estabilidade dos pontos de equiĺıbrio, e alguns resultados sobre o comportamento do sistema quando a mola e o motor estão em condição de ressonância, isto é, quando a freqüência de vibração da mola e a freqüência do motor estão próximas. No caṕıtulo 4 apresentamos mais detalhadamente os estudos anaĺıticos feitos por Dantas e Balthazar [3], e contribúımos com o estudo numérico sobre a bifurcação de Hopf. Procedimento análogo é feito no caṕıtulo 5, onde detalhamos outro trabalho de Dantas e Balthazar. Neste artigo é assumido que o coeficiente de não linearidade da mola seja negativo, e desta forma surge um par de pontos de equiĺıbrio do tipo sela e uma órbita heterocĺınica Caṕıtulo 1. Apresentação dos Caṕıtulos 5 conectando esses pontos. Usando um sistema Hamiltoniano perturbado associado à dinâmica do sistema e utilizando-se um método de redução de graus de liberdade [9] e o método de Melnikov [19], os autores mostram a existência de vibrações caóticas no sistema dinâmico. O caṕıtulo 6 contém as conclusões e as considerações finais sobre o trabalho desenvol- vido. Em seguida apresentamos os apêndices que complementam o trabalho e finalmente as referências utilizadas nesta dissertação. Caṕıtulo 2 Conceitos preliminares: Tópicos de Sistemas Dinâmicos Apresentamos neste caṕıtulo, de forma introdutória, alguns tópicos de sistemas dinâmicos que são utilizados nesta dissertação. Como não objetivamos uma apresentação detalhada dos conceitos, sempre que posśıvel, não nos preocupamos em enunciar teoremas e definições, mas sim discorrer sobre o assunto de maneira informal. Enunciamos somente alguns dos teoremas que serão usados posteriormente. Estes tópicos e teoremas que são apresentados, são resul- tados clássicos da teoria de sistemas dinâmicos, e por isso podem ser encontrados facilmente na literatura. Algumas das referências que utilizamos neste caṕıtulo são [1], [2], [6], [8], [13], [14], [17] e [20]. Alguns textos que não foram utilizados neste caṕıtulo mas, que dão suporte para alguns estudos espećıficos que não foram apresentados, serão referenciados no decorrer do texto. 2.1 Sistemas lineares Seja o sistema linear dx dt def = ẋ = Ax x ∈ R n, (2.1) onde A é uma matriz n× n com coeficientes constantes. Uma solução do sistema (2.1) é uma função vetorial x(x0, t) que depende do tempo t e da condição inicial x(0) = x0 (2.2) tais que x(x0, t) é solução de (2.1) e (2.2). 6 Caṕıtulo 2. Sistemas lineares 7 A solução do sistema linear (2.1) é simplesmente x(x0, t) = eAtx0, (2.3) onde eAt é uma matriz exponencial n×n. Isso pode ser facilmente deduzido do desenvolvimento em série da matriz eAt, ou seja eAt = [I + tA + t2 2! A2 + · · · + tn (n)! An + · · ·] (2.4) onde I é a matriz unitária. Derivando (2.4) termo a termo em relação a t, temos d dt eAt = [A + tA2 + · · · + tn−1 (n − 1)! An + · · ·] = AeAt. (2.5) Portanto, pode-se usar (2.4) e (2.5) para mostrar que (2.3) é solução de (2.1) e (2.2). A solução geral do sistema linear é a superposição de n soluções linearmente indepen- dentes {x1(t), · · · , xn(t)} x(t) = n ∑ j=1 cjx j(t), onde as n constantes cj são determinadas pela condição inicial. Se os autovalores {λ1, · · · , λn} de A forem reais e distintos , então A tem n autovetores linearmente independentes {v1, · · · , vn}. Além disso, podemos diagonalizar A por meio de uma transformação linear. A partir de (2.4), é fácil de ver a matriz eAt também será diagonal com autovalores {eλ1t, · · · , eλnt}. Deste modo, a solução de cada xj(t) será xj(t) = eλjtvj. No caso em que A tem autovalores complexos conjugados αj±iβj e autovetores vRe±ivIm, temos as soluções xj Re = eαjt(vRe cos βjt − vImsenβjt), xj Im = eαjt(vResenβjt + vIm cos βjt), onde xj Re e xj Im são respectivamente a parte real e a parte imaginária de xj(t). De maneira simplificada escrevemos xj(t) = eRe(λj)teiIm(λj)tvj. Como eiIm(λj)t é uma função limitada, a estabilidade de xj(t) vai depender essencialmente de Re(λj). Se Re(λj) > 0, eRe(λj)t cresce continuamente com o tempo e xj(t) → ∞ quando Caṕıtulo 2. Classificação do equiĺıbrio quanto à topologia e à estabilidade 8 t → ∞. Isso significa que as trajetórias xj(t) deixam a vizinhança de um ponto de equiĺıbrio P ∗. Inversamente, se Re(λj) < 0, xj(t) → P ∗ quando t → ∞ e nesse caso o ponto de equiĺıbrio é estável. As diferentes possibilidades de combinação dos autovalores, que podem ser reais, ima- ginários puros, todos com parte real positiva, etc..., vão definir não só a estabilidade do ponto de equiĺıbrio mas também a forma das soluções em sua vizinhança. Isso nos leva a classificar os pontos de equiĺıbrio de acordo com a sua natureza. 2.2 Classificação do equiĺıbrio quanto à topologia e à estabilidade O ponto de equiĺıbrio de um sistema linear pode ser classificado de acordo com a topologia do seu retrato de fases e de acordo com sua estabilidade. Essa classificação é realizada em função dos autovalores. Considere o sistema de n equações diferenciais (2.1). O polinômio caracteŕıstico é obtido através de det(A − λI) = 0. Quando todos os autovalores da matriz A tiverem a parte real diferente de zero, o ponto de equiĺıbrio correspondente P ∗ é chamado de hiperbólico, independente do valor da parte imaginária. Quando pelo menos um autovalor tem a parte real nula, o ponto de equiĺıbrio é denominado de não-hiperbólico. Os pontos de equiĺıbrio hiperbólicos podem ser classificados de três formas quanto à estabilidade: atratores, repulsores, e selas. • Se todos os autovalores de A tem a parte real negativa, o ponto de equiĺıbrio é chamado de atrator, sendo que neste caso o equiĺıbrio é assintoticamente estável. Se todos os autovalores de A são complexos, então o atrator é chamado de foco estável, se todos os autovalores de A são reais, o atrator é chamado de nó estável. • Se todos os autovalores da matriz A tem a parte real positiva o ponto de equiĺıbrio é chamado de repulsor ou fonte. Se os autovalores são complexos, então a fonte é chamada de foco instável e, se todos os autovalores de A são reais, a fonte é chamada de nó instável. • Quando alguns autovalores (mas não todos) têm a parte real positiva e o restante tem a parte real negativa, então o ponto de equiĺıbrio é chamado de sela. Caṕıtulo 2. Critério de Routh-Hurwitz 9 Quanto à estabilidade de pontos de equiĺıbrio não-hiperbólicos, pode-se dizer que: • Um ponto de equiĺıbrio não-hiperbólico é instável se um ou mais autovalores de A tem a parte real positiva. • Se alguns autovalores da matriz A tem a parte real negativa, enquanto que os outros autovalores tem a parte real nula, o ponto de equiĺıbrio é chamado de marginalmente estável. • Se todos os autovalores da matriz A são imaginários puros e não-nulos, o ponto de equiĺıbrio é chamado de centro. 2.3 Critério de Routh-Hurwitz Como foi visto, a estabilidade de um ponto de equiĺıbrio é estabelecida pelo sinal da parte real de seus autovalores λj. Portanto, na determinação da estabilidade não é necessário calcular explicitamente os valores de λj, bastando conhecer o sinal das suas partes reais. Num sistema linear, os autovalores da matriz A são ráızes do polinômio caracteŕıstico que é obtido a partir de det(A − λI) = 0. O problema de descobrir se todas as ráızes de um polinômio têm parte real negativa, sem calcular explicitamente essas ráızes, foi solucionado em 1874 por E.J. Routh (1831 − 1907). Em 1895 A. Hurwitz (1859 − 1919) encontrou, independentemente, uma solução equivalente. Considere um sistema dinâmico composto por n equações diferenciais de primeira ordem. O polinômio caracteŕıstico de grau n deste sistema tem a forma geral λn + a1λ n−1 + a2λ n−2 + · · · + an−1λ + an = 0 (2.6) O critério de Hurwitz garante que, se todas as ráızes desse polinômio têm parte real negativa, então os coeficientes aj são positivos. Essa é uma condição necessária mas não suficiente, isto é, o fato dos coeficientes aj serem positivos não garante que as ráızes terão a parte real negativa. As condições necessárias e suficientes para que todas as ráızes do polinômio (2.6) tenham parte real negativa são dadas em termos da matriz de Hurwitz, Hn×n, onde H é constrúıda Caṕıtulo 2. Sistemas não-lineares 10 de seguinte forma                  a1 a3 a5 a7 · · · 0 0 1 a2 a4 a6 · · · 0 0 0 a1 a3 a5 · · · 0 0 0 1 a2 a4 · · · 0 0 ... ... ... ... . . . ... ... 0 0 0 0 · · · an−1 0 0 0 0 0 · · · an−2 an                  . Observe que a diagonal principal dessa matriz contem, sem repetição, os coeficientes a1, a2, · · · , an. O critério de Hurwitz estabelece que Re(λj) < 0 para todo j = 1, 2 · · · , n se são todos positivos os coeficientes aj e se são positivos os determinantes ∆j onde ∆1 = |a1|, ∆2 = ∣ ∣ ∣ ∣ ∣ ∣ a1 a3 1 a2 ∣ ∣ ∣ ∣ ∣ ∣ , ∆3 = ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣ a1 a3 a5 1 a2 a4 0 a1 a3 ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣ ∣ , · · · , ∆n = |H|. 2.4 Sistemas não-lineares Considere o sistema de equações não-lineares ẋ = f(x), x ∈ R n e f(x) = (f1(x), f2(x), · · · , fn(x))T (2.7) para o qual existe um ponto de equiĺıbrio x̄, isto é f(x̄) = 0. Em geral, o estudo global de sistemas não lineares é muito dif́ıcil, então procuramos determinar o comportamento das soluções x(t) na vizinhança de x̄. Assim, podemos linearizar (2.7) em torno de x̄, e deste modo estudar o sistema linear ξ̇ = Aξ, ξ ∈ R n onde A = [∂fi/∂xj]x=x̄ é a matriz Jacobiana da função f(x) = (f1(x), f2(x), · · · , fn(x))T e x = x̄ + ξ, |ξ| � 1. Se certas condições são satisfeitas, as novas variáveis ξ descrevem o comportamento das soluções de (2.7) na vizinhança de x̄. Assim, o estudo da estabilidade de um ponto de equiĺıbrio de um sistema dinâmico não-linear reduz-se ao estudo de um sistema linear Caṕıtulo 2. Equivalência topológica 11 correspondente. Alguns teoremas garantem as condições necessárias e suficientes para que os resultados da aproximação linear sejam válidos, ao menos localmente. Antes de apresentarmos esses teoremas, é necessário introduzir o conceito de equivalência topológica. 2.5 Equivalência topológica Assuma que a função h(x) = y estabeleça uma relação entre o conjunto x e o conjunto y:, isto é h(x) = y ⇐⇒                h1(x1, x2, · · · , xn) = y1 h2(x1, x2, · · · , xn) = y2 ... ... hn(x1, x2, · · · , xn) = yn Suponha que h seja uma função injetora e sobrejetora. Uma função com essas pro- priedades é invert́ıvel, ou seja, existe a função inversa h−1(y) = x. Se h é cont́ınua, invert́ıvel e sua inversa h−1 é cont́ınua, então h é chamada de homeomorfismo, e o domı́nio x e imagem y são considerados homeomorfos. Um homeomorfismo é chamado difeomorfismo se h(x) e h−1(y) são diferenciáveis em todos os pontos. Por exemplo, a função y = h(x) = x2 não é injetora para x ∈] −∞, +∞[, portanto, não representa um homeomorfismo. Entretanto, y = h(x) = x2 é injetora para x ∈]0, +∞[ e, para esse domı́nio, constitui um homeomorfismo, já que é também sobrejetora e tem inversa h−1(y) = √ y cont́ınua. A função h(x) = x3 é um homeomorfismo para x ∈ ] −∞, +∞[, mas não é um difeomorfismo, pois sua inversa não é diferenciável na origem. A função h(x) = ex é injetora, mas não é sobrejetora para x ∈] −∞, +∞[. Quando os retratos de fases dos sistemas dinâmicos: dx dt = f(x) e dy dt = g(x) podem ser relacionados por um homeomorfismo h(x) = y que preserva o sentido do movimento (a orientação) no espaço de fases, então esses sistemas são topologicamente orbitalmente equivalentes. Isso significa que as trajetórias de um sistema podem ser continuamente defor- madas até se tornarem iguais às trajetórias do outro sistema. Deformações cont́ınuas envolvem esticamentos e alongamentos, mas não cortes ou emendas. Dois retratos de fases que apre- sentam a mesma estrutura orbital são qualitativamente equivalentes. Conseqüentemente, eles representam comportamentos dinâmicos similares. Caṕıtulo 2. Teoremas locais para sistemas não-lineares 12 2.6 Teoremas locais para sistemas não-lineares 2.6.1 Teorema de Hartman-Grobman D.M. Grobman, em 1959, e P. Hartman, em 1963, provaram independentemente que, na vizinhança de um ponto de equiĺıbrio hiperbólico, um sistema não-linear de dimensão n ap- resenta um comportamento qualitativamente equivalente ao do sistema linear correspondente [14]. Portanto, o teorema de Hartman-Grobmam garante que a estabilidade de um ponto de equiĺıbrio hiperbólico é preservada quando se lineariza o sistema em torno desse ponto, de modo que o retrato de fases, na sua vizinhança, é topologicamente orbitalmente equivalente ao retrato de fases do sistema linear associado. Dois retratos de fases são topologicamente or- bitalmente equivalentes quando um é uma versão distorcida do outro. Se o ponto de equiĺıbrio é não-hiperbólico, ou seja, se há algum autovalor com parte real nula, então a linearização não permite predizer sua estabilidade. 2.6.2 Teorema das variedades hiperbólicas Seja um sistema de equações diferenciais não-lineares dx dt = f(x), com o campo vetorial f sendo de classe r ( f r vezes diferenciável). Seja P ∗ um ponto de equiĺıbrio de f e considere a matriz jacobiana calculada nesse ponto, a partir da versão linear. Os autovalores corres- pondentes a essa matriz podem ser separados em 3 grupos: σe, σi e σc, dependendo do sinal da parte real desses autovalores. As letras e, i, c são, respectivamente, as iniciais de estável, instável e central. Assim λ ∈ σe se Re(λ) < 0 λ ∈ σi se Re(λ) > 0 λ ∈ σc se Re(λ) = 0 sendo que λ ∈ C. O sub-espaço gerado pelos autovetores cujos autovalores pertencem σe é chamado de auto-espaço estável Ee; aquele gerado pelos autovetores cujos autovalores per- tencem a σi é chamado de auto-espaço instável Ei, e aquele correspondente a σc, de auto- espaço central Ec. Note que esses auto-espaços referem-se à versão linear do sistema de equações diferenciais. Caṕıtulo 2. Teoria da variedade central 13 Se há ne autovalores com parte real negativa, ni autovalores com parte real positiva e nc com parte real nula, então ne + ni + nc = n, sendo n a dimensão do sistema. Soluções pertencentes a Ee apresentam decaimento exponencial, soluções pertencentes a Ei exibem crescimento exponencial, e soluções pertencentes a Ec têm estabilidade neutra. O teorema de Hartman-Grobman é válido somente quando nc = 0. Seja um sistema dinâmico descrito por n equações diferenciais autônomas. Um conjunto S de pontos do espaço de fases n-dimensional é uma variedade invariante local se, para x0 pertencente a esse conjunto S, a solução x(t) com condição inicial x0 está em S para t < T , com T > 0. Se isso é válido para T → ∞, então S é uma variedade invariante. O teorema das variedades hiperbólicas, provado de forma completa por A. Kelley em 1967, afirma que existe uma variedade estável W e tangente em P ∗ ao auto-espaço Ee. Essa variedade é única e é de classe r. Existe também uma variedade instável W i, única e de classe r, tangente em P ∗ ao auto-espaço Ei. A figura (2.1) ilustra esses resultados [14]. Figura 2.1: Variedades de um ponto de sela de um sistema não-linear 2.7 Teoria da variedade central O teorema de Hartman-Grobman e o teorema das variedades hiperbólicas estabelecem que, se um ponto de equiĺıbrio P ∗ é hiperbólico, então W e e Ee se tangenciam em P ∗, assim como W i e Ei se tangenciam em P ∗ e que, na vizinhança desse ponto, o sistema dinâmico não-linear e a versão linearizada são topologicamente orbitalmente equivalentes. O teorema da variedade central estende esses resultados para o caso em que o ponto de equiĺıbrio apresenta Caṕıtulo 2. Teoria da variedade central 14 algum autovalor que é nulo ou tem parte real nula. Mostra-se que o comportamento do sistema dinâmico ao redor de um ponto de equiĺıbrio não-hiperbólico pode ser reduzido ao estudo do comportamento ao longo de uma variedade central W c, que é tangente ao sub-espaço central Ec no ponto de equiĺıbrio em questão. 2.7.1 Teoremas da variedade central Seja um sistema de equações diferenciais não-lineares dx dt = f(x), com o campo vetorial f sendo de classe r. Seja P ∗ um ponto de equiĺıbrio de f e considere a matriz jacobiana calculada nesse ponto, a partir da versão linear. O teorema da variedade central estabelece que se há autovalor da matriz jacobiana λ ∈ σc, então existe uma variedade central W c tangente ao sub-espaço Ec em P ∗. Essa variedade, entretanto, é de classe r − 1 e não é necessariamente única. Isso dificulta a caracterização do comportamento assintótico em torno de P ∗ [6, 14]. 2.7.2 Teoremas de Carr Considere que um sistema dinâmico de dimensão n apresente um ponto de equiĺıbrio. Por conveniência, esse ponto pode ser deslocado para a origem de um novo sistema de coordenadas. A partir dáı, encontram-se seus autovalores e descobre-se que o ponto é não-hiperbólico. Então, determinam-se seus autovetores e reescreve-se o sistema original de modo que a parte linear esteja na forma canônica de Jordan. Dessa maneira, o sistema original passa a ser escrito como dxj dt = c ∑ i=1 Ajixi + fj(x, y) j = 1, 2, · · · , c dyk dt = e ∑ i=1 Bkiyi + gk(x, y) k = 1, 2, · · · , e (2.8) sendo x = (x1, x2, · · · , xc), y = (y1, y2, · · · , ye) e c + e = n. Na notação vetorial, escreve-se dx dt = Ax + f(x, y) dy dt = By + g(x, y) (2.9) onde Aji são os elementos da matriz A e Bki os elementos da matriz B, e f e g são funções de classe C2 com f(0, 0) = g(0, 0) = 0, f ′(0, 0) = g′(0, 0) = 0. Além disso, suponha que todos os autovalores de B tenha a parte real negativa, e que todos os autovalores de A tenha a parte real nula. Caṕıtulo 2. Teoria da variedade central 15 Uma variedade invariante é definida como uma variedade central, se ela pode ser repre- sentada localmente pelas curvas cont́ınuas e diferenciáveis yk = hk(x) para k = 1, 2, · · · , e. Além disso, as funções hk(x) devem ser tais que hk(0) = 0 e (∂hk/∂xi)|x=0 = 0 para todo i. Agora, substituindo yk = hk(x) na equação (2.8) para dxj/dt, obtém-se dxj dt = c ∑ i=1 Ajixi + fj(x, h(x)) j = 1, 2, · · · , c. (2.10) J. Carr, em 1981, provou, de maneira completa, três teoremas que possibilitam calcular uma variedade central, ao menos aproximadamente. Teorema 2.1. Existe uma variedade central y = h(x), |x| < δ, para o sistema (2.9) sendo que h é C2 Esse teorema estabelece a existência de uma variedade central para o sistema original (2.8). Teorema 2.2. [8] Se a origem x = 0 de (2.10) é localmente assintoticamente estável(instável), então a origem de (2.9) é localmente assintoticamente estável(instável). O segundo teorema prova que o sistema original (2.9) e o sistema reduzido (2.10) são equivalentes quanto ao comportamento assintótico das soluções em torno da origem. Se, em um deles, a origem é assintoticamente estável (instável), no outro, a origem também é assintoticamente estável (instável). Isso possibilita uma redução na dimensão do sistema original, já que se descartam as e equações referentes à variedade estável. A dificuldade de se usar esse teorema é a deter- minação da variedade central. Agora mostraremos como h(x) pode ser calculada ao menos aproximadamente. Como y(t) = h(x(t)) temos que Dh(x)[Ax + f(x, h(x))] = Bh(x) + g(x, h(x)) ou N (h(x)) = Dh(x)[Ax + f(x, h(x))] − Bh(x) − g(x, h(x)) = 0 com h(0) = Dh(0) = 0 Caṕıtulo 2. Bifurcações em sistemas dinâmicos de tempo cont́ınuo 16 Teorema 2.3. [8] Se uma função φ(x), com φ(0) = Dφ(0) = 0, puder ser estabelecida tal que N (φ(x)) = O(|x|p) para algum p > 1 quando |x| → 0 então segue que h(x) = φ(x) + O(|x|p) quando |x| → 0. 2.8 Bifurcações em sistemas dinâmicos de tempo cont́ınuo O termo bifurcação, introduzido por Poincaré em 1885, refere-se a mudança qualitativa do espaço de fases de um sistema dinâmico, conforme algum parâmetro do sistema passa por um valor cŕıtico. A idéia de bifurcação está intimamente ligada ao conceito de estabilidade estrutural. Estabilidade estrutural relaciona-se com a preservação (ou não) da topologia do espaço de fases, quando as equações que originam esse espaço são perturbadas através da variação dos valores dos parâmetros das equações. Assim, um sistema dinâmico é estruturalmente estável, se ele é orbitalmente topologicamente equivalente a uma versão perturbada. Se, no entanto, ao se variar o valor do parâmetro em torno de um valor cŕıtico ocorre uma mudança qualitativa no seu espaço de fases, então o sistema dinâmico é estruturalmente instável, para aquele valor cŕıtico do parâmetro. Essa mudança na topologia do espaço de fases é chamada de bifurcação. Bifurcações em sistemas dinâmicos podem ser classificadas como bifurcações locais e bifurcações globais. Bifurcações locais são aquelas que podem ser previstas estudando-se o campo vetorial na vizinhança de um ponto de equiĺıbrio ou de uma órbita fechada. Normalmente, esse estudo é feito pelo cálculo de autovalores. Bifurcações globais são aquelas que não podem ser deduzidas a partir de uma análise local. Apresentaremos nesta seção, apenas o teorema da bifurcação de Hopf. Para maiores detalhes sobre teoria de bifurcações, veja as referências citadas no ińıcio do caṕıtulo. No final deste caṕıtulo, apresentaremos um estudo sobre quebra de órbitas homocĺınicas e heterocĺıni- cas em sistemas bidimensionais, em particular o método clássico de Melnikov. Este assunto será apresentado nesta ordem em virtude de que precisamos primeiramente de conceitos de mapas de Poincaré. Caṕıtulo 2. Mapas e seção de Poincaré 17 2.8.1 Bifurcação de Hopf Teorema 2.4. Suponha que o sistema ẋ = fµ(x), x ∈ R n, µ ∈ R tenha um ponto de equiĺıbrio (x0, µ0), no qual as seguintes propriedades são satisfeitas: (a) Dxfµ0(x0) tem um par de autovalores puramente imaginários e outros com a parte real diferente de zero. Então (a) implica que existe uma curva de equiĺıbrio suave (x(µ), µ) com x(µ0) = µ0. Os autovalores λ(µ), λ̄(µ) de Dxfµ0(x(µ)) que são imaginários em µ = µ0 variam suavemente com µ. (b) Se, além disso d dµ (Re(λ(µ))) |µ=µ0= d 6= 0 então existe uma única variedade central tri-dimensional passando por (x0, µ0) em R n×R, e um sistema de coordenadas suave (preser- vando o plano µ=const.) para o qual a expansão de Taylor até terceira ordem na variedade central é dada por    u̇ = (dµ + k(u2 + v2))u − (ω + cµ + b(u2 + v2))v + O(4) v̇ = (ω + cµ + b(u2 + v2))u + (dµ + k(u2 + v2))v + O(4) Se k 6= 0, então existe uma superf́ıcie de soluções periódicas na variedade central a qual tem tangência quadrática com o auto-espaço de λ(µ0), λ̄(µ0) concordando em segunda ordem com o parabolóide µ = −(a/b)(u2 + v2). Se k < 0, então estas soluções periódicas são ciclos limite estáveis (bifurcação de Hopf supercŕıtica), enquanto que se k > 0, as soluções periódicas são instáveis (bifurcação de Hopf subcŕıtica). Para ver mais detalhes sobre teorema da bifurcação de Hopf consulte [8] e [12]. 2.9 Mapas e seção de Poincaré Denominamos mapa um sistema dinâmico que evolui no tempo de uma forma discreta. Seja um sistema dinâmico cont́ınuo não-linear e o fluxo φt a ele associado. Esse fluxo pode dar origem a um mapa xi+1 = Fµ(xi), (2.11) onde x é um vetor n-dimensional, Fµ é uma função não-linear (µ é o parâmetro de controle), e i representa os passos temporais fixos e discretos, ou passagens sucessivas por uma superf́ıcie de seção do fluxo. Se o fluxo φt é liso (r-vezes continuamente diferenciável) então F é um mapa liso com uma inversa lisa, i.e., um difeomorfismo. A órbita do mapa será então a seqüência de Caṕıtulo 2. Mapas e seção de Poincaré 18 pontos (xi) −∞ +∞ definida pela equação (2.11), que é genericamente denominada uma equação de diferenças. Equações de diferenças podem também ser lineares, e nesse caso xi+1 = Bxi e qualquer ponto inicial gera uma única órbita, desde que B não tenha autovalores nulos (caso degenerado). Uma das maneiras pela qual um fluxo cont́ınuo dá origem a um mapa discreto é pela utilização de seções de Poincaré. A seção de Poincaré é uma maneira de reduzir o estudo de um fluxo num espaço de fases com n dimensões a uma aplicação (difeomorfismo), chamada mapa de Poincaré ou mapa de retorno (“return map”), num espaço de fases com (n − 1) dimensões. Considere por exemplo a equação ẍ + g(x, ẋ) = f(t), (2.12) onde f(t) é uma função periódica de peŕıodo T . O diagrama de fases de (2.12) é tridimensional, e cada estado é representado por (x, ẋ, t). Nesse caso o mapa de Poincaré é obtido simples- mente considerando-se a intersecção da trajetória com o plano (x, ẋ) (seção de Poincaré) toda vez que t for igual a um múltiplo de T . Introduzimos agora o conceito de mapa de Poincaré. Seja o sistema dinâmico autônomo n-dimensional, com soluções periódicas ẋ = f(x), x = (x1, x2, · · · , xn), (2.13) onde f(x) é um campo vetorial não linear. Seja x0 uma órbita periódica (peŕıodo T ) associada ao fluxo φ(t) gerado por (2.13). Tomemos uma hipersuperf́ıcie Σ ⊂ R n de dimensão (n − 1) de tal maneira que o fluxo seja transversal a ela. Seja x∗ 0 o ponto onde a órbita x0 intercepta Σ (fig. 2.2) e U ⊆ Σ uma vizinhança de x∗ 0 . Então o mapa de Poincaré P : U → Σ é definido para um ponto x′ 1 ∈ U por P (x′ 1) = φ(x′ 1, τ), onde τ = τ(x′ 1) é o tempo necessário para que a órbita φ(x1, t) que parte de x′ 1 retorne pela primeira vez a Σ. Em geral τ depende de x′ 1, mas τ → T quando x′ 1 → x∗ 0 . A hipersuperf́ıcie Σ é chamada seção de Poincaré. Uma órbita periódica corresponde a um ponto fixo x∗ 0 da aplicação P , isto é, P (x∗ 0) = x∗ 0. Caṕıtulo 2. Método de Melnikov 19 Figura 2.2: Significado geométrico da aplicação de Poincaré P : x0 é uma órbita periódica, x∗ 0 é o ponto onde ela intercepta Σ; x1 é uma órbita não periódica que corta Σ em P (x′ 1) 6= x′ 1. Em particular, a estabilidade de x∗ 0 relativamente à aplicação P reflete a estabilidade da órbita periódica x0 com relação ao fluxo φ(t). Na fig.(2.2) é representada a órbita periódica x0 (que passa por x∗ 0) e aquela não periódica x1. Dessa figura pode-se depreender o significado geométrico de P . P tem valores em Σ, que é (n − 1) dimensional. Mas x só pode variar em Σ; então P é um mapa (n − 1) dimensional. Na Fig.(2.3) mostra-se como fluxos em duas e três dimensões dão origem a mapas uni- dimensionais e bidimensionais. Definido sobre a superf́ıcie transversal ao fluxo, o mapa de Poincaré relaciona um ponto do fluxo ao seu primeiro ponto de cruzamento com essa mesma superf́ıcie. O resultado das iterações do mapa de Poincaré é portanto a sequência de pontos nos quais o fluxo intercepta a seção de Poincaré. 2.10 Método de Melnikov V. K. Melnikov, em 1963, desenvolveu um método pelo qual se pode provar, analitica- mente, a existência de bifurcações homocĺınicas ou heterocĺınicas em sistemas Hamiltonianos perturbados. 2.10.1 Caso autônomo Considere um sistema dinâmico autônomo, do tipo dx1 dt = f1(x) + εg1(x), dx2 dt = f2(x) + εg2(x), (2.14) Caṕıtulo 2. Método de Melnikov 20 Figura 2.3: Seções de Poincaré : (a) unidimensional; (b) bidimensional sendo 0 < ε � 1 um parâmetro de valor pequeno. Assuma que o campo vetorial f ∈ R 2 é Hamiltoniano e que esse campo é perturbado pelo campo vetorial g ∈ R 2. Assuma que para qualquer valor de ε, o sistema tem um ponto de equiĺıbrio do tipo sela, e quando ε = 0, o sistema (2.14) possui uma órbita homocĺınica x0(t) ligando o ponto de sela a ele mesmo ( o caso onde sistema possui uma órbita heterocĺınica pode ser tratado de forma similar). Melnikov provou que, para ε 6= 0 a órbita homocĺınica somente existe se a função de Melnikov M̄ = ∫ +∞ −∞ [f1(x)g2(x) − f2(x)g1(x)]|x(t)=x0(t)dt se anula para alguma combinação de parâmetros do sistema (Na literatura normalmente usa- se M para denotar a função de Melnikov. Escolhemos usar M̄ aqui, pois M será usado posteriormente para denotar outra função). A função de Melnikov relaciona-se com a distância entre as variedades estável e instável de um ponto de sela do sistema perturbado, sendo que essa distância é medida em relação a um ponto P sobre a órbita homocĺınica do sistema não-perturbado. Caṕıtulo 2. Método de Melnikov 21 A distância d vale d = ε M̄ |f(P )| + O(ε2). Para ε = 0, essas variedades são tangentes. Para ε 6= 0, as variedades somente se tangenciam se a função de Melnikov se anula para alguma combinação de valores do sistema. Além disso, se M̄ 6= 0, então a posição das variedades estável e instável em relação à órbita x0(t) do sistema não-perturbado depende do sinal de M̄ . Se M̄ > 0 a variável estável está “fora” e a variedade instável está “dentro”, e se M̄ < 0 ocorre situação inversa (fig. 2.4) [14]. Figura 2.4: Caso ε = 0: linha tracejada; caso ε 6= 0: linha cheia. 2.10.2 Caso não-autônomo A quebra de uma órbita homocĺınica ou heterocĺınica pode resultar em comportamento caótico. Considere um sistema dinâmico não-autônomo, do tipo dx1 dt = f1(x) + εg1(x, t), dx2 dt = f2(x) + εg2(x, t), (2.15) no qual f = (f1, f2) é um campo Hamiltoniano e g = (g1, g2) é um campo perturbativo periódico de peŕıodo fixo T em t, e 0 < ε � 1. Assuma que f tem um ponto de sela com uma órbita homocĺınica associada. Equivalentemente, se ε 6= 0 podemos escrever o sistema (2.15) na forma estendida dx1 dt = f1(x) + εg1(x, φ), dx2 dt = f2(x) + εg2(x, φ), dφ dt = 1. Caṕıtulo 2. Método de Melnikov 22 Assim, no espaço de fase estendido, pode-se analisar o sistema numa secção de Poincaré definida por φ = φ0 com 0 ≤ φ0 < T . Constrói-se o mapa de Poincaré tomando x1(0), x2(0), φ0 como condições iniciais e integrando o sistema de t = φ0 e t = φ0 + T . Para ε = 0, as variedades estável e instável do ponto de sela são tangentes. Para ε 6= 0, no caso autônomo essas variedades ou são tangentes ou não se tocam. Para ε 6= 0, no caso não- autônomo, essas variedades podem se interseccionar transversalmente. Ao estudar o problema da estabilidade do sistema solar, Poincaré percebeu que se elas se interseccionam de maneira não-tangente, então o número dessas intersecções é infinito. A figura (2.5) mostra como essas intersecções ocorrem, e a figura resultante é chamada de emaranhado homocĺınico. A existência de emaranhado implica caos. Figura 2.5: Emaranhado homocĺınico Em 1892, 1893 e 1899, foram publicados os 3 volumes do livro Les Methodes Nouvelles de la Mecanique Celeste (“Os Novos Métodos da Mecânica Celeste”) de Poincaré. No ter- ceiro livro ele escreveu que não se atreveria a desenhar figura de tal complexidade. Em [1] encontra-se uma citação de Poincaré da obra de sua autoria citada acima no que se refere aos emaranhados homocĺınicos. Segue na ı́ntegra: Que se procure representar tal figura formada por estas duas curvas e suas infinitas inter- secções, onde cada uma corresponde a uma solução duplamente assintótica, suas intersecções formando um tipo de trelissa, de tecido, uma rede de malhas infinitamente apertadas; cada uma das curvas não deve jamais se auto-recortar, mas elas têm de dobrar sobre si mesmas de uma maneira muito complexa para voltar a cortar um número infinito de vezes a malha da rede. A complexidade dessa figura é tão atordoante que eu nem mesmo procuro traça-la. Nada é mais apropriado a nos dar uma idéia da complicação do problema de três corpos e em geral Caṕıtulo 2. O mapa de ferradura de Smale 23 de todos os problemas da Dinâmica em que não há uma integral uniforme... O primeiro esboço do emaranhado homocĺınico é atribúıdo a Birkhoff, realizado por volta de 1930. De acordo com Melnikov, em primeira aproximação em potências de ε, a distância entre as variedades estável e instável, associadas ao ponto de sela para o caso ε 6= 0, é proporcional à função de Melnikov definida como M̄(φ0) = ∫ +∞ −∞ [f1(x)g2(x, t + φ0) − f2(x)g1(x, t + φ0)]|x(t)=x0(t)dt onde nota-se agora a dependência da seção de Poincaré. Em sistemas autônomos, todas as seções de Poincaré são identicas, já que tais sistemas não dependem exclusivamente do tempo. Seja uma seção de Poincaré particular φ0 = φ01. Se para M̄(φ0) = 0, tem-se dM̄(φ0) d(φ0) |φ0=φ01 6= 0, então as variedades estável e instável do sistema perturbado se interseccio- nam transversalmente na seção de Poincaré. Suponha que para alguma combinação de valores dos parâmetros, as variedades estável e instável se interceptam transversalmente. Como são variedades invariantes, um ponto que pertence a essas curvas deve permanecer nelas indefinidamente. Se o ponto de intersecção não é um ponto de equiĺıbrio, então a ocorrência de uma intersecção entre as variedades implica a ocorrência de infinitas intersecções. Portanto, se a função M̄ se anula numa seção de Poincaré, então essa função tem que se anular em infinitos pontos. Cada um desses pontos é chamado de ponto homocĺınico transversal [14]. 2.11 O mapa de ferradura de Smale O mapa de ferradura foi concebido por Smale por volta de 1960. Esse mapa1 consiste em considerar uma aplicação f sobre um quadrado unitário D ∈ R 2, tal que f : D → R 2 onde D = {(x, y) ∈ R 2| 0 ≤ x ≤ 1 e 0 ≤ y ≤ 1}. 1Nesta seção consideraremos apenas mapa de ferradura bi-dimensional. Uma versão tri-dimensional pode ser encontrada em [18, 19] Caṕıtulo 2. O mapa de ferradura de Smale 24 Assumimos que f−1 contrai-se na direção horizontal x por um fator menor que 1/2 e estica-se na direção vertical por um fator maior que 2 e depois dobra. O mapa é definido apenas no quadrado unitário: pontos que abandonam o quadrado são ignorados. Em particular, considere as faixas horizontais H0 = {(x, y) ∈ R 2| 0 ≤ x ≤ 1 e 0 ≤ y ≤ 1 µ }, H1 = {(x, y) ∈ R 2| 0 ≤ x ≤ 1 e 1 − 1 µ ≤ y ≤ 1}, onde µ é uma constante maior que 2. Figura 2.6: Esquema do mapa de ferradura O mapa de ferradura é definido por f−1(H0) =   x y  →   λ 0 0 µ     x y   , f−1(H1) =   x y  →   −λ 0 0 −µ     x y  +   1 µ   , onde 0 < λ < 1 2 . Deste modo, f−1 leva as faixas horizontais H0 e H1 nas faixas verticais V0 e V1 por f−1(H0) = V0 = {(x, y)| 0 ≤ x ≤ λ, 0 ≤ y ≤ 1} e f−1(H1) = V1 = {(x, y)| 1 − λ ≤ x ≤ 1, 0 ≤ y ≤ 1}. O mapa f−1 e o inverso f estão esquematizados na figura (2.6) Caṕıtulo 2. Teorema de Bezout 25 O conjunto invariante Λ do mapa é a coleção de pontos que permanecem em D após todas as iterações de f : Λ = · · · f−2(D) ∩ f−1(D) ∩ D ∩ f(D) ∩ f 2(D) ∩ · · · = ∞ ⋂ n=−∞ fn(D). Esse conjunto consiste de uma intersecção infinita de retângulos horizontais e verticais. A primeira interação de f−1 produz dois retângulos V0 e V1, V0 à esquerda e V1 à direita. Os retângulos V0 e V1 geram quatro retângulos V00, V01, V11 e V10, conforme a figura (2.7). Em geral a n-ésima iteração produz 2n retângulos, de largura λn e altura 1. Os iterados por f produzem 2n faixas horizontais após a n-ésima iteração. Smale provou que o conjunto Figura 2.7: Iterações de f−1 sobre o quadrado D dos pontos Λ que permanecem no quadrado D após a ação de fn e f−n para n → ∞, isto é, a intersecção das 2n faixas verticais e 2n faixas horizontais sobre o quadrado D para n → ∞, é um conjunto do tipo Cantor [11]. Uma questão interessante é qual a importância do mapa de ferradura de Smale? A resposta vem do teorema homocĺınico de Smale-Birkhoff ,o qual estabelece que a dinâmica na vizinhança de um ponto homocĺınico transversal de um difeomorfismo é similar à dinâmica do mapa de ferradura, baseada em esticamentos e dobras. 2.12 Teorema de Bezout Teorema 2.5. Duas curvas planas de ordem m e n sem componentes em comum tem exata- mente m.n intersecções. Caṕıtulo 3 O Sistema Mecânico e Algumas Simulações Numéricas. O objetivo deste caṕıtulo é apresentar o sistema vibrante, bem como suas equações de movimento e por fim apresentar algumas simulações computacionais. Além disso, apresentare- mos as equações de movimento de duas formas. Uma usando as equações de Lagrange, e outra usando o formalismo de Hamilton. 3.1 O Vibrador Centŕıfugo Considere um sistema vibrante com dois graus de liberdade composto por um motor DC de massa m0 com potência limitada, e uma mola com não linearidade cúbica que sustenta o motor. Além disso, esse sistema é excitado por uma pequena massa m desbalanceada, acoplada à parte giratória do motor a uma distância r do eixo de rotação. Devido às caracteŕısticas desse sistema, temos influência rećıproca no comportamento da mola e do motor, e assim temos um sistema dinâmico não ideal ou simplesmente um problema não ideal. Este sistema mecânico é conhecido como vibrador centŕıfugo e está mostrado na figura (3.1). Denotemos por x o deslocamento vertical do motor e por ϕ o deslocamento angular da pequena massa m. Neste trabalho assumimos que o motor não tenha movimentos laterais, então se intro- duzirmos as coordenadas x1 e z1 como sendo as posições vertical e horizontal da massa m temos as relações x1 = x + r cos ϕ e z1 = rsenϕ. (3.1) 26 Caṕıtulo 3. O Vibrador Centŕıfugo 27 Figura 3.1: Vibrador centŕıfugo Usando-se as relações (3.1) obtemos as velocidades generalizadas dadas na forma ẋ1 = ẋ − rϕ̇senϕ e ż1 = rϕ̇ cos ϕ. (3.2) Em termos das coordenadas generalizadas, podemos escrever a energia cinética desse sistema na forma T = 1 2 [ m0ẋ 2 + m(ẋ2 1 + ż2 1) + Jϕ̇2 ] , (3.3) onde J é o momento de inércia da parte giratória do motor. Assim, substituindo as velocidades dadas pelas eqs.(3.2) na eq.(3.3), reescrevemos a energia cinética T em termos do raio r e das variáveis x e ϕ, ou seja T = 1 2 { m0ẋ 2 + m[(ẋ − rϕ̇senϕ)2 + (rϕ̇ cos ϕ)2] + Jϕ̇2 } = 1 2 { m1ẋ 2 + m(−2ẋrϕ̇senϕ + r2ϕ̇2) + Jϕ̇2 } onde m1 = m0 + m é a massa do sistema. A energia potencial do sistema é dada por V = U − Wc, onde U = 1 2 cx2 + 1 4 dx4 é a energia de deformação da mola, e Wc = −mgr cos ϕ denota o trabalho das forças conservativas (peso). O termo 1 4 dx4 que aparece em U determina uma caracteŕıstica de não linearidade da mola do tipo Duffing, e a constante g em V denota a aceleração da gravidade. Assim a energia potencial V é escrita como V = 1 2 cx2 + 1 4 dx4 + mgr cos ϕ. Portanto, depois de obtidas as energias cinética T e potencial V , e usando o formalismo de Lagrange podemos escrever a Lagrangeana L como sendo L = 1 2 { m1ẋ 2 + m(−2ẋrϕ̇senϕ + r2ϕ̇2) + Jϕ̇2 } − 1 2 cx2 − 1 4 dx4 − mgr cos ϕ. (3.4) Caṕıtulo 3. Ponto de Equiĺıbrio do Sistema 28 Assim, as equações de movimento de Lagrange (ver [7] ou [13] ) ficam d dt ( ∂L ∂ẋ ) − ∂L ∂x = −βẋ, (3.5) d dt ( ∂L ∂ϕ̇ ) − ∂L ∂ϕ = M(ϕ̇), (3.6) onde βẋ é uma força linear de resistência ao movimento oscilatório, e M(ϕ̇) é uma função que define as caracteŕısticas do torque resultante gerado pelo motor. Substituindo a Lagrangeana dada pela eq.(3.4) nas derivadas parciais ∂L ∂ẋ , ∂L ∂x , ∂L ∂ϕ̇ e ∂L ∂ϕ das eqs.(3.5) e (3.6) obtemos ∂L ∂ẋ = m1ẋ − mrϕ̇senϕ, ∂L ∂x = −cx − dx3, ∂L ∂ϕ̇ = −mẋrsenϕ + (J + mr2)ϕ̇, ∂L ∂ϕ = 1 2 (−2mẋrϕ̇ cos ϕ) + mgr cos ϕ. Agora, levando estes resultados nas eqs.(3.5) e (3.6), podemos escrever as equações do sistema vibrante na forma m1ẍ + βẋ + cx + dx3 = mrϕ̇2 cos ϕ + mrϕ̈senϕ, (3.7) Iϕ̈ = M(ϕ̇) + mrẍsenϕ + mgrsenϕ, (3.8) onde mrϕ̈ cos ϕ, mrϕ̈senϕ e mrẍsenϕ representam os termos de acoplamento da parte vi- bratória do sistema com a fonte de energia (motor DC). As constantes c e d que aparecem na eq.(3.7) são respectivamente os coeficientes de elasticidade linear e não linear da mola. Além disso, denotamos por J o momento de inércia do rotor e mr2 é o momento de inércia da massa m em rotação. Conseqüentemente o momento de inércia total do sistema é dado por I = J + mr2. Assumimos neste caṕıtulo e no caṕıtulo 4 que todas as constantes que aparecem nas eqs. (3.7) e (3.8) são tomadas estritamente positivas. 3.2. Ponto de Equiĺıbrio do Sistema 29 3.2 Ponto de Equiĺıbrio do Sistema Nesta seção obteremos os pontos de equiĺıbrio do sistema de eqs.(3.7) e (3.8), reescre- vendo-o como um conjunto de quatro equações de primeira ordem. Para isso, introduzindo a mudança de variável y = (y1, y2, y3, y4) T com y1 = x, y2 = ẋ, y3 = ϕ, y4 = ϕ̇, as eqs.(3.7) e (3.8) se escrevem ẏ1 = y2, ẏ2 − a4(seny3)ẏ4 = a1y1 + a2y2 + a3y 3 1 + a4y 2 4 cos y3, ẏ3 = y4, ẏ4 − a5(seny3)ẏ2 = M(y4) I + a6(seny3), (3.9) onde a1 = −c m1 , a2 = −β m1 , a3 = −d m1 , a4 = mr m1 , a5 = mr I e a6 = mgr I . Podemos então reescrever a eq.(3.9) na forma matricial A0ẏ = f(y) onde A0 =         1 0 0 0 0 1 0 −a4(seny3) 0 0 1 0 0 −a5(seny3) 0 1         , e f(y) =         y2 a1y1 + a2y2 + a3y 3 1 + a4y 2 4 cos y3 y4 M(y4) I + a6(seny3)         . Agora, obtendo a matriz inversa de A0, escrevemos as equações do sistema vibrante, eqs.(3.7 e 3.8), como um conjunto de equações de primeira ordem na forma matricial ẏ = A−1 0 f(y), (3.10) sendo A−1 0 =         1 0 0 0 0 α1 0 α1a4(seny3) 0 0 1 0 0 α1a5(seny3) 0 α1         com α1 = 1 1 − a4a5sen2y3 . Caṕıtulo 3. Ponto de Equiĺıbrio do Sistema 30 Note que a4a5 = m2r2 m1(J+mr2) < m2r2 m1mr2 = m m1 < 1, então temos que 1 − a4a5sen2y3 6= 0 para qualquer y3 ∈ R, e desta forma A−1 0 é uma matriz com coeficientes bem definidos para todo y3. Finalmente reescrevemos a eq.(3.10) na forma         ẏ1 ẏ2 ẏ3 ẏ4         =         y2 f1(y1, y2, y3, y4) y4 f2(y1, y2, y3, y4)         (3.11) onde f1 = α1K1 + K2K3, f2 = K4K1 + α1K3, K1 = a1y1 + a2y2 + a3y 3 1 + a4y 2 4 cos y3 K2 = α1a4seny3, K3 = M(y4) I + a6seny3 e K4 = α1a5seny3 Procuremos agora, os pontos de equiĺıbrio da equação acima, isto é, os pontos y∗ = (y∗ 1, y ∗ 2, y ∗ 3, y ∗ 4) que representam as soluções estacionárias dy dt = 0 . Dessa forma temos 0 = y2, 0 = α1(a1y1 + a2y2 + a3y 3 1 + a4y 2 4 cos y3) + α1a4seny3 ( M(y4) I + a6seny3 ) , 0 = y4, 0 = α1a5seny3(a1y1 + a2y2 + a3y 3 1 + a4y 2 4 cos y3) + α1 ( M(y4) I + a6seny3 ) e o ponto de equiĺıbrio é dado por                  y1 = 0, y2 = 0, seny3 = senϕ0 = −M(0) mgr , y4 = 0. (3.12) Note que esse ponto de equiĺıbrio somente existe se o torque inicial M(0) satisfaz a desigualdade |M(0)| ≤ mgr. Em análise de estabilidade de pontos de equiĺıbrio faz-se cons- tantemente o uso de transformação de coordenadas a fim de transladar os pontos de equiĺıbrio para a origem. Com esse objetivo fazemos a transformação de variável y3 → y3 + ϕ0, e a eq. Caṕıtulo 3. Ponto de Equiĺıbrio do Sistema 31 (3.11) é reescrita como         ẏ1 ẏ2 ẏ3 ẏ4         =         y2 g1(y1, y2, y3, y4) y4 g2(y1, y2, y3, y4)         (3.13) onde g1 = α2V1 + V2V3, g2 = V4V1 + α2V3, V1 = a1y1 + a2y2 + a3y 3 1 + a4y 2 4 cos(y3 + ϕ0), V2 = α2a4sen(y3 + ϕ0), V3 = M(y4) I + a6sen(y3 + ϕ0) V4 = α2a5sen(y3 + ϕ0) e α2 = 1 1 − a4a5sen2(y3 + ϕ0) . Agora, este sistema possui um único ponto de equiĺıbrio (y1, y2, y3, y4) = (0, 0, 0, 0). O estudo anaĺıtico de sistemas dinâmicos não lineares é bastante complexo, e portanto limitamos a análise dos mesmos na vizinhança dos pontos de equiĺıbrio utilizando os sistema linear associado. Assim, expandindo a eq.(3.13) em série de Taylor em torno de (0, 0, 0, 0) obtemos         ẏ1 ẏ2 ẏ3 ẏ4         =         0 1 0 0 b1 b2 b3 b4 0 0 0 1 b5 b6 b7 b8                 y1 y2 y3 y4         + O(2), (3.14) onde b1 = κa1, b2 = κa2, b3 = κa4a6senϕ0 cos ϕ0, b4 = κ a4senϕ0M ′(0) I , b5 = κa1a5senϕ0, b6 = κa2a5senϕ0, b7 = κa6 cos ϕ0, b8 = κ M ′(0) I , κ = 1 1 − a4a5sen2ϕ0 e O(2) representam os termos de grau maior ou igual a 2. Temos que o polinômio caracteŕıstico associado à parte linear de (3.14) tem a forma p(X) = 1 1 − a4a5sen2ϕ0 ((1 − a4a5sen2ϕ0)X 4 + ( −M ′(0) I − a2 ) X3 + ( a2 M ′(0) I − a6 cos ϕ0 − a1 ) X2 + ( a2a6 cos ϕ0 + a1 M ′(0) I ) X (3.15) + a1a6 cos ϕ0). Sabe-se que a estabilidade de um ponto de equiĺıbrio é estabelecida pelo sinal da parte real dos autovalores da matriz jacobiana associada ao sistema. No entanto , para se determinar a estabilidade de um ponto de equiĺıbrio não é necessário calcular explicitamente os valores da parte real dos autovalores, mas apenas conhecer o sinal de suas partes reais. O critério Caṕıtulo 3. Equações de Movimento Usando a Formulação de Hamilton 32 estabelecido por Routh-Hurwitz (ver caṕıtulo 2 ), estabelece as condições necessárias para que todos os zeros do polinômio caracteŕıstico associado à matriz jacobiana do sistema linea- rizado tenham a parte real negativa. Da eq.(3.15) temos que p(0) = κa1a6 cos ϕ0 é positivo se cos ϕ0 < 0, então não consideraremos o caso cos ϕ0 > 0, pois isso torna o sistema dinâmico instável. 3.3 Equações de Movimento Usando a Formulação de Hamilton Nesta seção escreveremos as equações de movimento usando o formalismo Hamiltoni- ano. Não faremos nenhuma análise sobre estas equações neste caṕıtulo, mas o estudo destas equações será realizado no caṕıtulo 5. Considere então os momentos Px e Pϕ os quais são definidos em termos das velocidades ẋ e ϕ̇ pelas relações Px = ∂L ∂ẋ e Pϕ = ∂L ∂ϕ̇ (3.16) onde L = L(x, ẋ, ϕ, ϕ̇) é Lagrangeana do sistema dada por (3.4). Assumindo que r = 1 em (3.7) e (3.8) a função Hamiltoniana (ver [7] ou [13]) do sistema vibrante se escreve como sendo H(x, ϕ, Px, Pϕ) = ẋPx + ϕ̇Pϕ − 1 2 (m1ẋ 2 − 2mẋϕ̇senϕ + I1ϕ̇ 2) + 1 2 cx2 + 1 4 dx4 + mg cos ϕ, = ẋPx + ϕ̇Pϕ − 1 2 ẋ(m1ẋ − mϕ̇senϕ) − 1 2 ϕ̇(I1ϕ̇ − mẋsenϕ) + 1 2 cx2 + 1 4 dx4 + mg cos ϕ, (3.17) onde I1 é o momento de inércia do sistema motor-massa quando r = 1 (raio de giração da massa excêntrica). Usando as definições de momento e a Lagrangeana dadas respectivamente pelas eqs. (3.16) e (3.4) obtemos Px = m1ẋ − mϕ̇senϕ e Pϕ = I1ϕ̇ − mẋsenϕ. (3.18) Substituindo (3.18) na eq. (3.17) obtemos H(x, ϕ, Px, Pϕ) = ẋPx + ϕ̇Pϕ 2 + 1 2 cx2 + 1 4 dx4 + mg cos ϕ. (3.19) Caṕıtulo 3. Equações de Movimento Usando a Formulação de Hamilton 33 Agora, isolando ẋ e ϕ̇ em (3.18) obtemos ẋ = (m1I1 − m2sen2ϕ)Px + msenϕ(m1Pϕ + mPxsenϕ) m1(m1I1 − m2sen2ϕ) , (3.20) ϕ̇ = m1Pϕ + mPxsenϕ m1I1 − m2sen2ϕ . (3.21) Assim, substituindo as eqs.(3.20) e (3.21) na Hamiltoniana (3.19) obtemos H(x, ϕ, Px, Pϕ) = I1P 2 x + 2m(senϕ)PxPϕ + m1P 2 ϕ 2(m1I1 − m2sen2ϕ) + 1 2 cx2 + 1 4 dx4 + mg cos ϕ. (3.22) A fim de obter as equações de Hamilton para este problema, consideramos a função Hamiltoniana H(x, ϕ, Px, Pϕ) = Pxẋ + Pϕϕ̇ − L(x, ẋ, ϕ, ϕ̇) e assim, diferenciando-a obtemos dH = ∂H ∂Px dPx + ∂H ∂Pϕ dPϕ + ∂H ∂x dx + ∂H ∂ϕ dϕ, (3.23) dH = (dPx)ẋ + Px(dẋ) + (dPϕ)ϕ̇ + Pϕ(dϕ̇) − dL, (3.24) onde dL = ∂L ∂x dx + ∂L ∂ϕ dϕ + ∂L ∂ẋ dẋ + ∂L ∂ϕ̇ dϕ̇. (3.25) Além disso, usando as equações na forma de Lagrange dadas pelas eqs. (3.5) e (3.6) temos ∂L ∂x = Ṗx + βẋ e ∂L ∂ϕ = Ṗϕ − M(ϕ̇) (3.26) onde Ṗx = d dt ( ∂L ∂ẋ ) e Ṗϕ = d dt ( ∂L ∂ϕ̇ ) . Deste modo, comparando as equações (3.23) e (3.24) e usando conjuntamente as equações (3.25) e (3.26) obtemos ∂H ∂Px dPx + ∂H ∂Pϕ dPϕ + ∂H ∂x dx + ∂H ∂ϕ dϕ = (dPx)ẋ + Px(dẋ) + (dPϕ)ϕ̇ + Pϕ(dϕ̇) − (Ṗx + βẋ)dx − (Ṗϕ − M(ϕ̇))dϕ − Px(dẋ) − Pϕ(dϕ̇). Rearranjando os termos da equação anterior adequadamente obtemos ( ∂H ∂Px − ẋ ) dPx + ( ∂H ∂Pϕ − ϕ̇ ) dPϕ + ( ∂H ∂x + Ṗx + βẋ ) dx + ( ∂H ∂ϕ + Ṗϕ − M(ϕ̇) ) dϕ = 0. (3.27) Caṕıtulo 3. Resultados Numéricos: Pontos de equiĺıbrio 34 Portanto, da equação (3.27) obtemos                      ẋ = ∂H ∂Px , Ṗx = −∂H ∂x − β ∂H ∂Px , ϕ̇ = ∂H ∂Pϕ , Ṗϕ = −∂H ∂ϕ + M ( ∂H ∂Pϕ ) , (3.28) que são as equações Hamiltonianas de movimento do sistema vibrante (vibrador centŕıfugo), sendo a Hamiltoniana H dada por (3.22). 3.4 Resultados Numéricos: Pontos de equiĺıbrio Nesta seção simularemos as equações de movimento (3.11) adotando os seguintes valores para os parâmetros f́ısicos que aparecem no sistema - massa do motor: m0 = 1, 8kg; - coeficiente de elasticidade linear da mola: c = 0, 3N/m; - coeficiente de elasticidade não-linear da mola: d = 0, 06N/m3; - coeficiente de amortecimento: β = 0, 2Ns/m; - momento de inércia do rotor: J = 0, 001Kgm2; - excentricidade da massa: r = 0, 5m; - massa excêntrica: m = 0, 2kg; - aceleração da gravidade: g = 9, 80665m/s2. Para simular estas equações, adotamos primeiramente curvas caracteŕısticas do tipo linear (M(y4) = A−By4), e posteriormente do tipo exponencial (M(y4) = Ae−By4), onde A e B são constantes positivas. Da análise do ponto de equiĺıbrio e das caracteŕısticas das curvas temos que ϕ0 ∈ [π, 2π]. Além disso, no estudo sobre a estabilidade do ponto de equiĺıbrio, conclui-se que o ponto é estável para cos ϕ0 < 0. Como o ponto de equiĺıbrio (3.12) depende de M(0) = A, fixando os valores para as constantes que definem as caracteŕısticas das curvas, obtemos um conjunto de valores de equiĺıbrio estáveis (tab.3.1). Para M(0) = 0.01 o sistema apresenta um ponto de equiĺıbrio estável mostrado na tabela 3.1 e os comportamentos de y1 (oscilação da mola), y2 (velocidade de oscilação da mola), y3 (deslocamento angular da massa m) e y4 (velocidade angular da massa m) estão mostrados nas figuras (3.2) a (3.4). Na figura (3.2), o modelo de torque M(y4) é o linear, e apesar de Caṕıtulo 3. Resultados Numéricos: Pontos de equiĺıbrio 35 Tabela 3.1: Pontos de equiĺıbrio A B y1 y2 y3 y4 0.01 0.05 0 0 3.151789992 0 0.1 0.05 0 0 3.243741827 0 0.2 0.05 0 0 3.346976795 0 0.3 0.05 0 0 3.452491867 0 adotarmos a condição inicial (y1, y2, y3, y4) = (0, 0, 0.001, 0) a trajetória é estável. Quando o modelo de torque é o não linear, a trajetória correspondente a essa condição inicial é instável fig.(3.4). Uma trajetória estável só é posśıvel adotando-se uma condição inicial mais próxima do ponto de equiĺıbrio, como mostra a figura (3.3), onde a condição inicial é (0, 0, 3.001, 0). Uma conclusão imediata dessa simulação é que a bacia de atração do torque não linear é menor que a bacia de atração do torque linear. Nas figuras (3.5) a (3.7) são omitidos os gráficos de y1, y2 e y4, e mostramos apenas o comportamento do deslocamento angular da massa m para A = 0.1, 0.2 e 0.3. Para estes valores de A, as trajetória resultantes são estáveis para o torque linear e não linear. Caṕıtulo 3. Resultados Numéricos: Pontos de equiĺıbrio 36 Figura 3.2: Modelo de torque linear com A= 0.01 e condições inicial (y1, y2, y3, y4) = (0, 0, 0.001, 0). (a) Oscilação da mola; (b) velocidade de vibração da mola; (c) ângulo da massa m; (d) velocidade angular da massa m. Caṕıtulo 3. Resultados Numéricos: Pontos de equiĺıbrio 37 Figura 3.3: Modelo de torque exponencial, com A= 0.01 e condição inicial (y1, y2, y3, y4) = (0, 0, 3.001, 0). (a) Oscilação da mola; (b) velocidade de vibração da mola; (c) ângulo da massa m; (d) velocidade angular da massa m. Caṕıtulo 3. Resultados Numéricos: Pontos de equiĺıbrio 38 Figura 3.4: Modelo de torque exponencial, com A= 0.01 e condição inicial (y1, y2, y3, y4) = (0, 0, 0.001, 0). (a) Deslocamento da mola; (b) velocidade de vibração da mola; (c) ângulo da massa m; (d) velocidade angular da massa m. Caṕıtulo 3. Resultados Numéricos: Pontos de equiĺıbrio 39 Figura 3.5: Comportamento do movimento da massa m para A= 0.1. (a) Modelo de torque linear e condição inicial (y1, y2, y3, y4) = (0, 0, 0.001, 0). (b) Modelo de torque exponencial e condição inicial (y1, y2, y3, y4) = (0, 0, 3.001, 0). Figura 3.6: Comportamento do movimento da massa m para A= 0.2. (a) Modelo de torque linear e condição inicial (y1, y2, y3, y4) = (0, 0, 0.001, 0). (b) Modelo de torque exponencial e condição inicial (y1, y2, y3, y4) = (0, 0, 3.2, 0). Caṕıtulo 3. Resultados Numéricos: Pontos de equiĺıbrio 40 Figura 3.7: Comportamento do movimento da massa m para A= 0.3. (a) Modelo de torque linear e condição inicial (y1, y2, y3, y4) = (0, 0, 3, 0). (b) Modelo de torque exponencial e condição inicial (y1, y2, y3, y4) = (0, 0, 3.4, 0). O ponto de equiĺıbrio do sistema dado pela relação (3.12) é estável somente se |M(0)| ≤ mgr, e quando |M(0)| = mgr o ponto de equiĺıbrio (0, 0, ϕ0, 0) torna-se instável. Para cons- tatar isto numericamente fixamos para B o valor dado pela tabela 1, e assumimos para A os valores 0.980665 e 1.2. No caso do valor limite A = 0.980665, o ponto de equiĺıbrio (0, 0, ϕ0, 0), com ϕ0 = 3π 2 já torna-se instável conforme podemos ver na fig. (3.8). Nesta figura mostramos os comportamentos de yi, i = 1, 2, 3 e 4 no caso em que o torque é linear. Após a fase de transição (t ≈ 5seg), o motor executa movimento circulatório tornando instável o ponto de equiĺıbrio. Na fig. (3.9) mostramos o comportamento do sistema dinâmico para A = 1.2. Caṕıtulo 3. Resultados Numéricos: Pontos de equiĺıbrio 41 Figura 3.8: Modelo de torque linear. A= 0.980665. Condições iniciais:(y1, y2, y3, y4) = (0, 0, 4, 0). (a) Deslocamento da mola. (b) velocidade de vibração da mola. (c) Ângulo da massa. (d) Velocidade angular Caṕıtulo 3. Resultados Numéricos: Pontos de equiĺıbrio 42 Figura 3.9: Modelo de torque linear. A= 1.2. Condições iniciais:(y1, y2, y3, y4) = (0, 0, 4, 0). (a) Deslocamento da mola. (b) velocidade de vibração da mola. (c) Ângulo da massa. (d) Velocidade angular Caṕıtulo 3. Resultados Numéricos: Estudo da ressonância 1:1 43 3.5 Resultados Numéricos: Estudo da ressonância 1:1 Nesta seção apresentaremos alguns resultados numéricos adotando valores para os parâmetros tais que a mola e o motor estejam em condição de ressonância 1 : 1, ou seja, quando a freqüência natural da mola ω e a velocidade de rotação do motor ϕ̇ satisfaçam a relação ϕ̇ ≈ ω. Os valores escolhidos para os parâmetros f́ısicos do sistema estão listados na tabela (3.2). massa do motor: m0 = 1, 2kg; coeficiente de elasticidade linear da mola: c = 1, 0N/m; coeficiente de elasticidade não linear da mola: d = 0, 001N/m3; coeficiente de amortecimento: β = 0, 005Ns/m; momento de inércia do rotor: J = 0, 000017Kgm2; excentricidade da massa: r = 0, 5m; massa excêntrica: m = 0, 1kg; aceleração da gravidade: g = 9, 80665m/s2. Tabela 3.2: Parâmetros f́ısicos Note que com os valores dos parâmetros acima, temos um sistema com mola com ca- racteŕıstica quase linear. Desta forma, a freqüência da mola é aproximadamente ω = √ c m1 que vale ω ≈ 0.88. Consideramos nestas simulações, curvas caracteŕısticas para o torque do tipo linear (M(y4) = a − by4). Nos resultados que apresentaremos a seguir, fixamos b = 0.5 e ado- tamos a como um parâmetro de controle. Os resultados apresentados nas figuras (3.10), (3.11) e (3.12) correspondem respectivamente aos valores dos parâmetros a = 0.61, a = 0.67 e a = 0.72. Em cada figura é apresentado o histórico no tempo de y1, o espaço de fase associado e o histórico no tempo de y4. Estas figuras mostram o comportamento do sistema antes da ressonância (ϕ̇ < ω) fig.(3.10), na ressonância (ϕ̇ ≈ ω ≈ 0.88) fig (3.11) e depois da passagem pela ressonância (ϕ̇ > ω) fig. (3.12). Nas figuras (3.13) e (3.14) estão mostrados o comportamento do sistema quando está próximo à região de ressonância porém, após a passagem por esta região. Na figura (3.13) foi adotado a = 0.679 e na figura (3.14) a = 0.685. Observe que na figura (3.13) o histórico no tempo da mola apresenta inicialmente um espécie de “pulso” e depois um aumento de amplitude e, percebe-se neste caso, interferência rećıproca entre o comportamento da mola e Caṕıtulo 3. Resultados Numéricos: Estudo da ressonância 1:1 44 o motor. Na figura (3.14) ocorre o mesmo comportamento, mas a modulação inicial se estende por um peŕıodo de tempo maior. Este comportamento de modulação constitúı o mecanismo que diminui a amplitude de vibração da mola, à medida que afastamos da região de ressonância depois da passagem. A figura (3.15) reforça essa conclusão, onde foi utilizado a = 69. Na figura (3.16) está mostrada a amplitude máxima de oscilação da mola em função da freqüência do rotor, quando assumimos valores para o parâmetro a entre 0.56 e 0.85 com variação ∆a = 0.01. Esse é o efeito Sommerfeld que motivou os estudos sobre sistemas dinâmicos não ideais. Figura 3.10: Caracteŕısticas de movimento para parâmetro de controle a = 0.61. Caso ϕ̇ < ω Caṕıtulo 3. Resultados Numéricos: Estudo da ressonância 1:1 45 Figura 3.11: Caracteŕısticas de movimento para parâmetro de controle a = 0.67. Caso ϕ̇ ≈ ω Caṕıtulo 3. Resultados Numéricos: Estudo da ressonância 1:1 46 Figura 3.12: Caracteŕısticas de movimento para parâmetro de controle a = 0.72. Caso ϕ̇ > ω Caṕıtulo 3. Resultados Numéricos: Estudo da ressonância 1:1 47 Figura 3.13: Caracteŕısticas de movimento para parâmetro de controle a = 0.679. Aparecimento de modulações logo após a passagem pela ressonância. Caṕıtulo 3. Resultados Numéricos: Estudo da ressonância 1:1 48 Figura 3.14: Caracteŕısticas de movimento para parâmetro de controle a = 0.685. As modulações persistem mais à medida que se afasta da passagem pela ressonância. Caṕıtulo 3. Resultados Numéricos: Estudo da ressonância 1:1 49 Figura 3.15: Caracteŕısticas de movimento para parâmetro de controle a = 0.69. As modulações ocorrem praticamente em todo o intervalo de tempo de simulação. Caṕıtulo 3. Resultados Numéricos: Estudo da ressonância 1:1 50 Figura 3.16: Amplitude máxima de oscilação da mola em função da freqüência do rotor. Efeito Sommerfeld. Caṕıtulo 4 Bifurcação de Hopf no Vibrador Centŕıfugo Não Ideal O propósito deste caṕıtulo é obter as condições fundamentais para o surgimento da bifurcação de Hopf no vibrador centŕıfugo não ideal, ou seja, estamos interessados em um caso cŕıtico no qual o polinômio (3.15) tenha um par de autovalores puramente imaginários e os outros autovalores com a parte real negativa. Além disso, devemos obter todas as hipóteses do teorema da bifurcação de Hopf descritas no caṕıtulo 2. Em seguida realizamos alguns testes numéricos com o objetivo de constatar este resultado. 4.1 Estudo do Caso Cŕıtico Concentremo-nos agora, em obter as condições fundamentais para a ocorrência da bifurcação de Hopf na dinâmica do vibrador centŕıfugo. A condição necessária e suficiente para ocorrência da hipótese (a) do teorema de Hopf (teorema 2.4) para este pro- blema está na existência de números reais t, r e s tais que t 6= 0, r > 0 e s > 0 de forma que o polinômio caracteŕıstico p(X) seja reescrito como p(X) = (X2 + t2)(X2 + rX + s) ou p(X) = X4 + rX3 + (s + t2)X2 + t2rX + st2. (4.1) Note que o polinômio p(X) definido acima tem um par zeros puramente imaginários e os outros com a parte real negativa. 51 Caṕıtulo 4. Estudo do Caso Cŕıtico 52 Fazendo                      A1 = 1 1 − a4a5sen2ϕ0 ( −M ′(0) I − a2 ) , A2 = 1 1 − a4a5sen2ϕ0 ( a2 M ′(0) I − a6 cos ϕ0 − a1 ) , A3 = 1 1 − a4a5sen2ϕ0 ( a2a6 cos ϕ0 + a1 M ′(0) I ) , A4 = 1 1 − a4a5sen2ϕ0 a1a6 cos ϕ0, podemos reescrever (3.15) como p(X) = X4 + A1X 3 + A2X 2 + A3X + A4. (4.2) Agora, se compararmos os coeficientes de (4.1) com os coeficientes de (4.2) obtemos    A1 > 0, A2 > 0, A3 > 0, A4 > 0, A1A2A3 − A2 3 = A2 1A4. (4.3) Usando a relação (4.3) temos −M ′(0) I − a2 > 0, a2a6 cos ϕ0 + a1 M ′(0) I > 0 e a1a6 cos ϕ0 > 0. Portanto            M ′(0) < −a2I, M ′(0) < −a2a6(cos ϕ0)I a1 , cos ϕ0 < 0. Além disso, temos também que t = ± √ A3 A1 , r = A1, s = A2 − A3 A1 = A1A4 A3 . (4.4) Fazendo γ = M ′(0) I e δ = cos ϕ0 em (4.3) temos − a4a5a 2 2a 2 6δ 4 − 2a1a4a5a2a6γδ3 − a2 1a4a5γ 2δ2 − a2a1γ 3 − a2 2a6γ 2δ + a2a 2 6γδ2 + (−a2 2a1 + a4a5a 2 1)γ 2 + (−2a1a2a6 − a3 2a6 + 2a4a5a1a2a6)γδ + a4a5a 2 2a 2 6δ2 + a2a 2 1γ = 0 (4.5) Seja P (γ, δ) o polinômio do lado esquerdo de (4.5). Este polinômio é irredut́ıvel [4]. Agora, com intuito de obter a condição (b) do teorema da bifurcação de Hopf (teorema 2.4), tomemos um par γ0 e δ0 satisfazendo a condição (4.3). Fixando δ = δ0 em (3.15) e tomando γ como um parâmetro de controle podemos reescrever p(X) como f(γ,X) = X4 + A1(γ)X3 + A2(γ)X2 + A3(γ)X + A4(γ). (4.6) Caṕıtulo 4. Estudo do Caso Cŕıtico 53 Assim, usando o teorema da função impĺıcita em (4.6) e levando em consideração as relações (4.1), (4.3) e (4.4), conclúımos que existe uma curva de zeros de p(X), g(γ) = a(γ) + ib(γ) tal que a(γ0) = 0, b(γ0) = √ A3 A1 e f(γ0, ib(γ0)) = 0. Temos também que ∂f ∂X = 4X3 + 3A1(γ)X2 + 2A2(γ) + A3(γ), e então ∂f ∂X (γ0, ib(γ0)) = 2 −A3A 2 1 + i(−2A3 √ A3A1 + A1A2 √ A1A3) A2 1 . (4.7) Note que ib(γ0) não é zeros duplos de p(X), logo ∂f ∂X (γ0, ib(γ0)) 6= 0 e assim temos ∂g ∂γ (γ0) = − ∂f ∂γ (γ, g(γ0)) ∂f ∂X (γ, g(γ0)) , (4.8) onde ∂g ∂γ (γ0) = a′(γ0) + ib′(γ0). (4.9) Além disso, ∂f ∂γ = A′ 1(γ)X3 + A′ 2(γ)X2 + A′ 3(γ)X + A′ 4(γ) = −κX3 + a2κX2 + a1κX com κ = 1 1 − a4a5sen2ϕ0 . Portanto ∂f ∂γ (γ0, ib(γ0)) = κ −a2A3A1 + i(A3 √ A3A1 + a1 √ A3A1A1) A2 1 . (4.10) Agora, substituindo (4.7), (4.9) e (4.10) em (4.8) obtemos a′(γ0) + ib′(γ0) = −κ 2 ( −a2A3A1 + i(A3 √ A3A1 + a1 √ A3A1A1) −A3A2 1 + i(−2A3 √ A3A1 + A2 √ A3A1A1) ) , ou melhor a′(γ0) = −1 2 κ 1 | − A3A2 1 + i(−2( √ A3)3 √ A1 + A2 √ A3( √ A3)3)|2 × A1A3(A 2 1A3a2 + a1A 2 1A2 + A1A2A3 − 2a1A1A3 − 2A2 3), (4.11) onde Ai = Ai(γ0, δ0), i = 1, 2, 3, 4, e γ0, δ0 é um par de valores tais que p(X) tem um par de autovalores puramente imaginários e os outros com a parte real negativa. Caṕıtulo 4. Estudo de um Caso Particular 54 Note que a′(γ0) 6= 0 se e somente se A2 1A3a2 + a1A 2 1A2 + A1A2A3 − 2a1A1A3 − 2A2 3 6= 0. (4.12) A desigualdade (4.12) pode ser expressa como Q(γ, δ) 6= 0 com Q(γ, δ) = a1γ 3 + 3a2a1γ 2 + (−2a1a4a5a6δ 3 + (2a2 1a4a5 + a2 6)δ 2 + (a2 2a6 + 2a1a4a5a6 − 2a1a6)δ + a2 1 − 2a2 1a4a5 + 2a1a 2 2)γ − 2a2a4a5a 2 6δ 4 + 2a1a2a4a5a6δ 3 + (2a2a4a5a 2 6 − a2a 2 6)δ 2 + (−2a1a2a4a5a6 + a3 2a6 + 2a1a2a6)δ − a2 1a2. Observe que P (0, 0) = 0 e Q(0, 0) = −a2 1a2. Disto e da irredutibilidade de P (γ, δ), podemos concluir que o mdc{P (γ, δ), Q(γ, δ)} = 1, pois se assim não fosse, poderia existir um polinômio R(γ, δ) de grau menor que o de Q(γ, δ) tal que Q(γ, δ) = R(γ, δ)P (γ, δ), e desta maneira teŕıamos Q(0, 0) = 0 (contradição). Então, podemos usar o teorema de Bezout sobre intersecções de curvas algébricas (ver [17]) para concluir que as curvas P (γ, δ) = 0 e Q(γ, δ) = 0 se interceptam em um conjunto finito de pontos, o qual denotaremos por S. Desta forma, tomando um par (γ0, δ0) tal que P (γ0, δ0) = 0 e Q(γ0, δ0) 6= 0, e usando o teorema de Hopf temos o seguinte resultado [3]. Se M ′(0) e cos ϕ0 satisfaz em (4.3) e (γ0, δ0) /∈ S, então em uma adequada variedade central o sistema (3.13) pode ser reescrito como    u̇ = (dµ + k(u2 + v2))u − (ω + cµ + b(u2 + v2))v + O(4) v̇ = (ω + cµ + b(u2 + v2))u + (dµ + k(u2 + v2))v + O(4) (4.13) onde d = a′(γ0) e µ = γ − γ0. Observação 4.1. Todos os coeficientes de (4.13) são funções cont́ınuas de γ0, δ0, µ, onde µ varia no intervalo (−ε, ε) cujo tamanho depende de γ0 e δ0. 4.2 Estudo de um Caso Particular Assumindo que M(0) = 0 e cos ϕ0 < 0, temos pelo ponto de equiĺıbrio (3.12) e pelo polinômio P (γ, δ) que senϕ0 = 0, cos ϕ0 = −1 e γ0 = M ′(0) I = 0. Estes valores dos parâmetros satisfazem a condição (4.3) e além disso, substituindo-os na equação (4.11) obtemos a′(γ0) = 1 2 . Caṕıtulo 4. Estudo de um Caso Particular 55 De agora em diante, o objetivo é mostrar que k que aparece em (4.13) é diferente de zero e deste modo obter a bifurcação de Hopf. Fixando estes parâmetros no polinômio caracteŕıstico p(X) obtemos p(X) = X4 − a2X 3 + (a6 − a1)X 2 − a2a6X − a1a6, que possui os zeros X1,2 = ±i √ a6 e X3,4 = a2 ± √ a2 2 + 4a1 2 . Observe que para estes valores, o polinômio caracteŕıstico tem um par de autovalores puramente imaginários e dois autovalores com a parte real negativa. Portanto, sob estas condições temos as hipóteses fundamentais do teorema de bifurcação de Hopf. Além do mais, temos um ponto de equiĺıbrio não-hiperbólico, o qual nos permite estudar o sistema por meio de teoria de variedade central (ver caṕıtulo 2). Assim, expandindo o sistema (3.13) em série de Taylor em torno do ponto de equiĺıbrio (0, 0, 0, 0), e fixando os valores M(0) = 0, cos ϕ0 = −1, γ0 = 0 e senϕ0 = 0, colocamos diretamente a equação (3.13) na forma (2.9), ou seja   ẏ1 ẏ2   =   0 1 a1 a2     y1 y2  +   0 F1(y1, y2, y3, y4)  +   0 O(4)   , (4.14)   ẏ3 ẏ4   =   0 1 −a6 0     y3 y4  +   0 F2(y1, y2, y3, y4)  +   0 O(4)   , (4.15) sendo F1(y1, y2, y3, y4) = a4a6y 2 3 − a4y 2 4 + a3y 3 1 − a4 M (2)(0) 2I y2 4y3 + a1a4a5y1y 2 3 + a2a4a5y2y 2 3 e F2(y1, y2, y3, y4) = −a5(a1y1 + a2y2)y3 + M (2)(0) 2I y2 4 + a6 − 6a4a5a6 6 y3 3 + a5a4y3y 2 4 + M (3)(0) 6I y3 4. As equações de movimento do sistema dado por (4.14) e (4.15) está numa forma adequada para se usar a teoria de variedade central. Então, sejam h : R 2 → R 2, p(., .) e q(., .) funções de classe C2 tais que h(y3, y4) =   p(y3, y4) q(y3, y4)   =   y1 y2   Caṕıtulo 4. Estudo de um Caso Particular 56 e h(0, 0) = h′(0, 0) =   0 0   . Agora, substituindo (y1, y2) T = (h(y3, y4)) T no sistema (4.15) obtemos   ẏ3 ẏ4   =   0 1 −a6 0     y3 y4  +       0    −a5(a1p(y3, y4) + a2q(y3, y4))y3 + M (2)(0) 2I y2 4          +    0 a6 − 6a4a5a6 6 y3 3 + a5a4y3y 2 4 + M (3)(0) 6I y3 4    +   0 O(4)   . (4.16) Pelo segundo teorema de Carr (teorema 2.2), temos que o sistema dado por (4.14) e (4.15) e o sistema (4.16) são equivalentes quanto ao comportamento assintótico em torno da origem, ou seja o sistema bi-dimensional (4.16) contém todas as informações necessárias para se determinar o comportamento dinâmico do sistema dado por (4.14) e (4.15) em torno do ponto de equiĺıbrio não-hiperbólico em questão. No entanto, para estudarmos o comportamento de (4.16) precisamos encontrar uma expressão anaĺıtica para p(., .) e q(., .). O terceiro teorema de Carr (teorema 2.3) mostra como determinar as expressões para variedade central ao menos por uma aproximação polinomial. Considere funções v(., .) e w(., .) definidas em um conjunto aberto U ⊂ R 2 tal que v(0, 0) = w(0, 0) = v′(0, 0) = w′(0, 0) = 0 e o mapa H[v, w] : U → R 2 dado por H[v, w](y3, y4) = D[v, w](y3, y4) [B(y3, y4) + G1[v, w](y3, y4)] − C[v, w](y3, y4) − G2[v, w](y3, y4), (4.17) onde D[v, w](y3, y4) =   vy3 vy4 wy3 wy4   , B(y3, y4) =   0 1 −a6 0     y3 y4   , G1[v, w](y3, y4) =   0 a5a4y3y 2 4 + M(3)(0) 6I y3 4   +   0 −a5(a1v(y3, y4) + a2w(y3, y4))y3 + M(2)(0) 2I y2 4   , Caṕıtulo 4. Estudo de um Caso Particular 57 C[v, w](y3, y4) =   0 1 a1 a2     v(y3, y4) w(y3, y4)   , G2[v, w](y3, y4) =   0 a4a6y 2 3 − a4y 2 4   +      0   a3v(y3, y4) 3 − 1 2 a4M (2)(0)y2 4y3 a1a4a5y 2 3v(y3, y4) + a2a4a5y 2 3w(y3, y4)        . Efetuando os cálculos e separando os termos de grau menor ou igual a 2 obtemos H[v, w](y3, y4) =      −a6y3vy4 + y4vy3 − w   −a6y3wy4 + y4wy3 − a1v −a2w − a4a6y 2 3 + a4y 2 4        + O(3). (4.18) Considere agora os polinômios    p̄(y3, y4) = B1y 2 3 + B2y3y4 + B3y 2 4, q̄(y3, y4) = C1y 2 3 + C2y3y4 + C3y 2 4. (4.19) Substituindo (4.19) em (4.18) temos H[p̄, q̄](y3, y4) =              −a6(B2y3 + 2B3y4) + y1(2B1y3 + B2y4) −C1y 2 3 − C2y3y4 − C3y 2 4        −a6y3(C2y3 + 2C3y4) + y4(2C1y3 + C2y4) −a1(B1y 2 3 + B2y3y4 + B3y 2 4) −a2(C1y 2 3 + C2y3y4 + C3y 2 4) − a4a6y 2 3 + a4y 2 4                 + O(3), =                    (−C1 − a6B2)y 2 3 (−2a6B3 − C2 + 2B1)y3y4 (B2 − C3)y 2 4           (−a6C2 − a1B1 − a4a6 − a2C1)y 2 3 (−a2C2 − 2a6C3 + 2C1 − a1B2)y3y4 (−a2C3 − a1B3 + C2 + a4)y 2 4                    + O(3). (4.20) Caṕıtulo 4. Estudo de um Caso Particular 58 Agora, igualando os coeficientes de y2 3, y3y4, y2 4 a zero obtemos               B1 B2 B3 C1 C2 C3               = R               −a4a6(a1 + 4a6) 4a4a6a2 a4(a1 + 4a6) −4a4a 2 6a2 −4a4a6(a1 + 4a6) 4a4a6a2               (4.21) onde R = 1 (a1 + 4a6)2 + 4a2 2a6 . Portanto, de (4.17) à (4.21) segue que H[p̄, q̄] = O(3), e podemos usar o terceiro teorema de Carr (teorema 2.3) para obter h(y3, y4) =   p(y3, y4) q(y3, y4)   =   p̄(y3, y4) q̄(y3, y4)  + O(3). (4.22) Note que em (4.16) p(y3, y4) e q(y3, y4) são multiplicados por y3 e, em vista disto foi escolhido os polinômios p̄(y3, y4) e q̄(y3, y4) de grau 2, pois termos de ordem maior ou igual a 3 em (4.22) serão multiplicados por y3, resultando em termos de ordem maior ou igual a 4, os quais não influenciam na dinâmica do sistema. Fazendo agora, a mudança de variável y3 = u e y4 = √ a6v temos ẏ3 = u̇, ẏ4 = √ a6v̇. Assim, podemos escrever o sistema (4.16) como   1 0 0 √ a6     u̇ v̇   =   0 1 −a6 0     1 0 0 √ a6     u v   +          0       −a1a5p(u, √ a6)u − a2a5p(u, √ a6) + M (3)(0) 2I a6v 2 + (a6 − 6a6a4a5) 6 u3 +a4a5u √ a6v + M (3)(0) 6I a6 √ a6v 3                +   0 O(4)   . (4.23) Caṕıtulo 4. Estudo de um Caso Particular 59 Usando-se (4.19) e (4.22) escrevemos a variedade central h(y3, y4) nas novas variáveis na forma    p(u, √ a6v) = B1u 2 + √ a6B2uv + a6B3v 2 + O(3), q(u, √ a6v) = C1u 2 + √ a6C2uv + a6C3v 2 + O(3). (4.24) Agora, fazendo N =   1 0 0 √ a6   , então N−1 =   1 0 0 1√ a6   . Portanto, substituindo (4.24) em (4.23) e multiplicando ambos os membros por N−1 obtemos   u̇ v̇   =   0 √ a6 −√ a6 0     u v   +              0           − ( a1a5√ a6 B1 + a2a5√ a6 C1 − √ a6 6 (1 − 6a4a5) ) u3 −(a1a5B2 + a2a5C2)u 2v −(a1a5 √ a6B3 + a2a5C3 − a4a5 √ a6)uv2 + M3(0) 6I a6v 3 + M2(0) 2I √ a6v 2                        +   0 O(4)   . (4.25) Se um sistema bi-dimensional está na forma   ẋ ẏ   =   0 −ω ω 0     x y  +   f(x, y) g(x, y)   , com f(0) = g(0) = 0 e Df(0) = Dg(0) = 0, então o valor do coeficiente k é dado por k = 1 16 [fxxx + fxyy + gxxy + gyyy] + 1 16ω [fxy(fxx + fyy) − gxy(gxx + gyy) − fxxgxx + fyygyy] sendo fxxx = ∂3f ∂x3 , fxyy = ∂3f ∂x∂y2 e assim por diante (ver [8] ou [14]). O valor das derivadas parciais é calculado no ponto de bifurcação. Portanto, de (4.25) obtemos k = 2a2a4a5a 2 6 (a1 + 4a6)2 + 4a2 2a6 + 1 16I M (3)(0)a6, (4.26) Caṕıtulo 4. Resultados Numéricos 60 e se a eq. (4.26) for diferente de zero, temos a bifurcação de Hopf. Além disso, se M (3)(0) = 0 e como a2 < 0, obtemos uma bifurcação de Hopf super-cŕıtica. Usando o resultado anterior, podemos obter algumas informações sobre o caso M(0) 6= 0. Anteriormente fixamos δ0 = −1 e γ0 = 0. Note que (γ0, δ0) = (0,−1) satisfaz (4.3) e (0,−1) /∈ S. Desde que ∂P ∂γ (0,−1) = a2((a1 + a6) 2 + a2 1a6) 6= 0, então usando o teorema da função impĺıcita segue que existe uma função γ(δ), a qual é definida em uma vizinhança de −1, e tal que γ(−1) = 0 e P (γ(δ), δ) = 0. Sobre a continuidade de γ, podemos assumir que (γ(δ), δ) satisfaz (4.3) e (0,−1) /∈ S. Disto e da observação (4.1), obtemos que d e k são funções cont́ınuas de δ e µ. De fato d = d(δ, µ) e k = k(δ, µ) e sobre os cálculos anteriores segue que d(−1, 0) = 1 2 e k(−1, 0) = 2a2a4a5a 2 6 (a1 + 4a6)2 + 4a2 2a6 + a6M (3)(0) 16I . Assim, em uma vizinhança adequada de (−1, 0) temos que d(δ, µ) 6= 0 e k(δ, µ) 6= 0 se a relação (4.26) for diferente de zero. Disto, obtemos uma bifurcação de Hopf. Usando (3.12), temos δ = cos ϕ0 = −(1 − (M(0)/mgr)2)1/2 e a existência da bifurcação de Hopf pode ser escrita na seguinte forma (ver [3]). Se (4.26) for diferente de zero, então para todo M(0) adequadamente pequeno, existe uma bifurcação de Hopf no espaço de fase do sistema dinâmico dado por (3.7) e (3.8). 4.3 Resultados Numéricos 4.3.1 Bifurcação de Hopf para um caso particular Com o intuito de constatar numericamente a bifurcação de Hopf para o caso particular (seção 4.2), simulamos as equações (3.13) fixando para a constante A o valor A = 0, e utilizando vários valores para a constante B. Para os parâmetros f́ısicos do sistema, adotamos os mesmos valores adotados na seção 3.4. Além disso, nesta subseção consideramos que B assume valores negativos. O modelo de torque considerado aqui é o linear já que o torque exponencial para o caso particular A = 0 não é adequado. Inicialmente mostraremos a mudança de sinal na parte real dos autovalores quando B cruza o valor cŕıtico Bc = 0 (tab. 4.1). Caṕıtulo 4. Resultados Numéricos 61 Tabela 4.1: Autovalores para diferentes valores de B. A = 0 e Bc = 0 A = 0 e B = −0.00005 A = 0 e B = 0.00005 autovalores autovalores autovalores 4.385057068i 0.00049019 + 4.38505704i -0.00049019 + 4.38505704i -4.385057068i 0.00049019 - 4.38505704i -0.00049019 - 4.38505704i -0.05 + 0.384057287i -0.05 + 0.38405728i -0.05 + 0.38405728i -0.05 - 0.384057287i -0.05 - 0.38405728i -0.05 - 0.38405728i Para mostrar o comportamento do sistema dinâmico para diferentes valores de B, in- tegramos as equações (3.13) e analisamos as soluções estacionárias. Nas figuras (4.1) e (4.2) mostramos respectivamente os espaços de fases (y1, y2) e (y3, y4) para vários valores do parâmetro B. O diagrama de bifurcação da variável y1 é mostrado na figura (4.3) e (4.4). Figura 4.1: Espaços de fases (y1, y2) para diferentes valores de B Caṕıtulo 4. Resultados Numéricos 62 Figura 4.2: Espaços de fases (y3, y4) para diferentes valores de B Figura 4.3: Diagrama de bifurcação da variável y1 O diagrama de bifurcação é o obtido simulando o sistema dinâmico para vários valores do parâmetro B, e notamos que a bifurcação ocorre quando a inclinação da reta M(y4) = −By4 cruza o valor M ′(y4) = 0. Na figura (4.5) estão mostradas retas caracteŕısticas para diferentes valores de B. Note, pelo diagrama de bifurcação e pela expressão de M ′(y4) que, quando B < 0, temos que M ′(y4) > 0 e, neste caso, as trajetórias para cada valor de B são trajetórias Caṕıtulo 4. Resultados Numéricos 63 Figura 4.4: Ampliação da figura 4.3 em torno da origem. periódicas nos retratos de fase de (y1, y2) e (y3, y4). Para B > 0 temos M ′(y4) < 0, e assim essas trajetórias são assintoticamente estáveis. Figura 4.5: Função torque para o caso particular M(y4) = −By4. As retas 1 (B = 0.00005) e 2 (B = 0.00002) correspondem a M ′(y4) < 0, a reta 3 é o caso cŕıtico M ′(y4) = 0 e as retas 4 (B = −0.00002) e 5 (B = −0.00005) referem-se a M ′(y4) > 0. Caṕıtulo 5 Análise da dinâmica do sistema através do método de Melnikov Neste caṕıtulo, consideramos um tipo de mola com rigidez menor. Em particular, assu- mimos que o coeficiente de não-linearidade da mola d seja estritamente negativo. Em vista disto, surge um par de equiĺıbrio do tipo sela na dinâmica da mola. Um aspecto que será discutido é um posśıvel surgimento de vibrações caóticas nesta estrutura. Para detectar a dinâmica caótica, usaremos o método clássico de Melnikov. A idéia básica desse método, é que se pode provar analiticamente a existência de bifurcações homocĺınicas ou heterocĺınicas em sistemas Hamiltonianos perturbados. Neste sentido, reescreveremos as equações Hamil- tonianas de movimento do vibrador centŕıfugo eqs.(3.28) usando uma função Hamiltoniana com uma perturbação em primeira potência de ε (definiremos ε posteriormente). Além disso, usaremos um método de redução de graus de liberdade dado em [9] e aplicaremos a teoria de Melnikov neste sistema reduzido, usando a função de Melnikov assim como ela é dada em [19]. 5.1 Formulação do Problema num Campo Hamiltoni- ano Perturbado Iniciamos esta seção voltando a atenção para a pequena massa m acoplada à parte gi- ratória do motor. Fisicamente falando, quando o motor entra em rotação, vibrações diferentes na estrutura de sustentação do motor pode ocorrer quando consideramos que o motor gira sem massa ou com a pequena massa m. Desta forma, podemos considera