UNIVERSIDADE ESTADUAL PAULISTA Instituto de Geociências e Ciências Exatas Campus de Rio Claro O USO DA ANÁLISE DE FOURIER, DE WAVELETS E DOS EXPOENTES DE LYAPUNOV NO ESTUDO DE UM SISTEMA DINÂMICO NÃO-IDEAL COM ATRITO SECO E EXCITAÇÃO EXTERNA Natale Chierice Júnior Orientador: Prof. Dr. José Roberto Campanha Dissertação de Mestrado elaborada junto ao Programa de Pós-Graduação em Física - Área de Concentração em Física Aplicada, para obtenção do Título de Mestre em Física. Rio Claro (SP) 2007 621 Chierice Júnior, Natale C533u O uso da análise de Fourier, de wavelets e dos expoentes de Lyapunov no estudo de um sistema dinâmico não-ideal com atrito seco e excitação externa / Natale Chierice Júnior. – Rio Claro : [s.n.], 2007 94 f. : il. Dissertação (mestrado) – Universidade Estadual Paulista, Instituto de Geociências e Ciências Exatas Orientador: José Roberto Campanha 1. Física aplicada. 2. Caos. 3. Atrator. 4. Bifurcação. 5. Seção de Poincaré. 6. adere-desliza. I. Título. Ficha catalográfica elaborada pela STATI – Biblioteca da UNESP Campus de Rio Claro/SP Comissão Examinadora Prof. Dr. José Roberto Campanha Instituição: IGCE/RC Prof. Dr. José Manoel Balthazar Instituição: DEMAC/RC Prof. Dr. Reyolando M. L. R. F. Brasil Instituição: ESCOLA POLITÉCNICA-USP/SP Natale Chierice Júnior Rio Claro, 19 de março de 2007. Resultado: Aprovado i Dedico este trabalho à minha Família. ii Agradecimentos Agradeço às seguintes pessoas e entidades: Ao Prof. Dr. José Roberto Campanha, pela preciosa orientação deste trabalho, amizade e incentivo e pelos ensinamentos transmitidos. À Comissão Permanente do Magistério da Aeronáutica – COPEMA, da Academia da Força Aérea, pelo apoio e pela dispensa semanal concedida para que pudesse concluir este trabalho. Ao Departamento de Física da UNESP de Rio Claro, pelo apoio e facilidades proporcionadas. Aos professores do curso de pós-graduação em Física, área de concentração Física Aplicada, pela amizade e ensinamentos. Aos colegas do curso de Pós – Graduação: Luciano Ferro, Luiz Roberto Salomão e Sidney Jorge Schinaider. A todos que colaboraram de alguma forma na realização deste trabalho. Obrigado a todos. Obrigado a Deus. iii Sumário Índice iv Resumo vi Abstract vii Lista de figuras viii Nomenclatura x I – Introdução 1 II – Caos e métodos matemáticos no estudo de sistemas dinâmicos 7 III – Equações que descrevem o sistema 36 IV – Simulações numéricas do sistema 50 V – Conclusão 75 iv Índice Capítulo 1 – Introdução 1 1.1 – Movimento adere-desliza 1 1.2 – Sistema proposto 3 1.3 – Sistema ideal e não-ideal 4 1.4 – Objetivo 5 Capítulo 2 – Caos e métodos matemáticos no estudo de sistemas dinâmicos 7 2.1 – Caos 7 2.2 – Breve histórico do caos 8 2.3 – Métodos matemáticos no estudo do sistema dinâmicos 11 2.3.1 – Série de Fourier 11 2.3.2 – Transformada de Fourier 13 2.3.3 – Transformada discreta de Fourier 15 2.3.4 – Wavelets 17 2.3.5 – Transformada wavelet contínua 27 2.3.6 – Transformada wavelet discreta 28 2.3.7 – Obtenção da transformada wavelet 30 2.3.8 – Expoentes de Lyapunov 33 Capítulo 3 – Equações que descrevem o sistema 36 3.1 – Introdução 36 3.2 – Forças que agem no bloco 37 3.2.1 – Força de atrito 37 3.2.2 – Força mola 41 3.2.3 – Força externa 42 3.3 – Equação do movimento do bloco 42 3.4 – Equação do movimento do motor CC 43 3.5 – Equações que descrevem o sistema 47 v Capítulo 4 – Simulações numéricas do sistema 50 4.1 – Introdução 50 4.2 – Atribuição de valores aos coeficientes do sistema 51 4.3 – Integração numérica do sistema 52 4.4 – Resultados das simulações 53 4.5 – Seção de Poincaré 61 4.6 – Bifurcação 65 4.7 – Análise da interação do movimento do bloco e da fonte de energia 66 4.7.1 – Freqüência da força externa igual a 14 Hz 67 4.7.2 – Freqüência da força externa igual a 12,6 Hz 70 4.7.3 – Comparações dos resultados 74 Capítulo 5 – Conclusão 75 5.1 – Comentários sobre os métodos utilizados 75 5.2 – Sugestões para trabalhos futuros 76 Referências Bibliográficas 77 Apêndice 81 Apêndice A 81 Apêndice B 85 Apêndice C 88 Apêndice D 90 Apêndice E 92 vi Resumo As oscilações mecânicas quando interferem no comportamento de um sistema mecânico estão relacionadas à transferência de energia devido ao atrito. A dinâmica desses sistemas com atrito pode ser prejudicada com o surgimento de movimentos caóticos. O estudo do comportamento dinâmico dessas oscilações mecânicas é o objetivo deste trabalho e para isto propomos um sistema não-ideal que descreve um modelo físico que trata do movimento de um bloco e de um motor elétrico de corrente contínua. O bloco preso a um extremo de uma mola com o outro extremo preso a um suporte fixo está apoiado em uma correia movimentada pelo motor elétrico. Sofrendo influências da força de atrito, da força da mola e de uma força externa que age harmonicamente, o bloco muitas vezes interfere na velocidade angular do motor, causando comportamentos caóticos no sistema. Com simulações numéricas estudamos o sistema, usando a transformada rápida de Fourier, transformada wavelet, expoentes de Lyapunov, diagrama de bifurcação, seção de Poincaré, trajetórias de plano de fase e gráficos da posição do bloco em função do tempo, em busca das freqüências que fazem o bloco oscilar em movimentos periódicos e caóticos. A importância desse estudo está em mostrar que métodos distintos conduzem a um mesmo resultado. Palavras-chave: sistema não-ideal, atrito, movimentos caóticos, diagrama de bifurcação. vii Abstract The mechanical oscillations when they interfere in the behavior of a mechanical system are related to the transfer of energy due to the friction. The dynamics of such systems with friction can be harmed by the appearance of chaotic movements. The study of the dynamic behavior of those mechanical oscillations is the objective of this work and for this we proposed a non-ideal system that describes a physical model that treats the movement of a block and a direct current motor. The block locked to the end of a spring with the other end locked to a fixed support is rested in a belt moved by a direct current motor. Suffering influences of the friction force, the spring force and the external force that act harmoniously, the block many times interferes in the angular speed of the motor, causing chaotic behaviors in the system. With numeric simulations, we studied the system using the fast Fourier transform; wavelet transform, Lyapunov exponents, bifurcation diagram, Poincaré section, phase plane trajectories and graphs of the block position in time function, looking of the frequencies that make the block to oscillate in periodic and chaotic movements. The importance of such study is to show that different methods lead to a same result. Keyword: non-ideal system, friction, chaotic movement, bifurcation diagram. viii Lista de Figuras Figura 1-1 – Modelo de um sistema dinâmico não-ideal com atrito seco e excitação externa 04 Figura 2-1 – Wavelet de Haar 20 Figura 2-2 – Wavelet chapéu mexicano 20 Figura 2-3 – Wavelet de Morlet 21 Figura 2-4 – Compressão da wavelet de Haar 23 Figura 2-5 – Dilatação da wavelet de Morlet 23 Figura 2-6 – Translação de wavelets 24 Figura 2-7 – Wavelets de Haar geradas a partir da wavelet “mãe” 25 Figura 2-8 – Escalogramas de funções com freqüências previamente conhecidas com escalas em Hz e tempo em segundo 26 Figura 2-9 – Análise wavelet em uma seção do sinal 31 Figura 2-10 – Análise wavelet em nova seção do sinal 31 Figura 2-11 – Análise wavelet com uma wavelet em nova escala 32 Figura 2-12 – Correspondência entre escalas e freqüências na análise wavelet 32 Figura 2-13 – Evolução de um elemento de volume esférico de raio 0(x0) em torno de um ponto inicial x0 34 Figura 3-1 – Modelo de um sistema dinâmico não-ideal com atrito seco e excitação externa 37 Figura 3-2 – Força de atrito com s e k diferentes 38 Figura 3-3 – Gráfico da função h(u) com s e k iguais 40 Figura 3-4 – Gráfico da função h*(u) = tanh(30u) com s e k iguais 41 Figura 3-5 – Esquema genérico de um motor elétrico CC com carga inercial I 43 Figura 3-6 – Diagrama do torque de um motor CC versus velocidade angular 46 Figura 4-1 – Gráficos da simulação do sistema dinâmico proposto para = 12 Hz 54 Figura 4-2 – Escalograma da transformada wavelet para = 12 Hz 54 Figura 4-3 – Gráficos da simulação do sistema dinâmico proposto para = 12,5 Hz 55 Figura 4-4 – Escalograma da transformada wavelet para = 12,5 Hz 55 Figura 4-5 – Gráficos da simulação do sistema dinâmico proposto para = 13 Hz 56 Figura 4-6 – Escalograma da transformada wavelet para = 13 Hz 56 ix Figura 4-7 – Gráficos da simulação do sistema dinâmico proposto para = 13,5 Hz 57 Figura 4-8 – Escalograma da transformada wavelet para = 13,5 Hz 57 Figura 4-9 – Gráficos da simulação do sistema dinâmico proposto para = 14 Hz 58 Figura 4-10 – Escalograma da transformada wavelet para = 14 Hz 58 Figura 4-11 – Gráficos da simulação do sistema dinâmico proposto para = 14,5 Hz 59 Figura 4-12 – Escalograma da transformada wavelet para = 14,5 Hz 59 Figura 4-13 – Gráficos da simulação do sistema dinâmico proposto para = 15 Hz 60 Figura 4-14 – Escalograma da transformada wavelet para = 15 Hz 60 Figura 4-15 – Passagem do movimento periódico para o caótico 62 Figura 4-16 – Movimento caótico do sistema dinâmico proposto 63 Figura 4-17 – Passagem do movimento caótico para o periódico 64 Figura 4-18 – Diagrama de bifurcação 65 Figura 4-19 – Gráficos da simulação do sistema dinâmico proposto para = 14 Hz e 0 = 18m2 kg/s2 68 Figura 4-20 – Gráficos da simulação do sistema dinâmico proposto para = 14 Hz e 0 = 12m2 kg/s2 69 Figura 4-21 – Gráficos da simulação do sistema dinâmico proposto para = 14 Hz e 0 = 2,2m2 kg/s2 70 Figura 4-22 – Gráficos da simulação do sistema dinâmico proposto para = 12,6 Hz e 0 = 18m2 kg/s2 71 Figura 4-23 – Gráficos da simulação do sistema dinâmico proposto para = 12,6 Hz e 0 = 12m2 kg/s2 72 Figura 4-24 – Gráficos da simulação do sistema dinâmico proposto para = 12,6 Hz e 0 = 2,2m2 kg/s2 73 x Nomenclatura CC – corrente contínua C(a,b) – transformada wavelet fa – força de atrito fk – força de atrito cinético fmo – força da mola fs – força de atrito estático fe – força externa f – amplitude da força externa Hz – hertz i – corrente do circuito do motor elétrico I – carga inercial do motor elétrico k – constante de elasticidade da mola ka – coeficiente de atrito viscoso do motor elétrico L – indutância do motor elétrico m – massa do bloco r – raio das rodas envolvidas pela correia R – resistência do circuito do motor elétrico u – velocidade relativa de deslizamento V – voltagem aplicada no motor elétrico Vemf – força contra eletromotriz x – posição do bloco x – velocidade do bloco x – aceleração do bloco – expoente de Lyapunov – freqüência angular da força externa 0 – velocidade angular máxima do motor elétrico – freqüência natural da oscilação do bloco – deslocamento angular do motor elétrico – velocidade angular do motor elétrico – parâmetro da tangente hiperbólica xi – coeficiente de atrito k – coeficiente de atrito cinético s – coeficiente de atrito estático – torque do motor elétrico 0 – torque de estol do motor elétrico CAPÍTULO 1 Introdução 1.1 Movimento adere-desliza Muitos fenômenos físicos que notamos, intuitivamente, ocorrem quando existe transferência de energia mecânica de uma forma para outra. Exemplo disto está, dentre muitos, no movimento oscilante da suspensão de um carro, no movimento do solo pela passagem de um ônibus nas proximidades e na propagação do som de uma música associada às vibrações acústicas. Essas oscilações (vibrações) presentes nesses fenômenos estão relacionadas à transferência de energia associada a dois efeitos devido ao atrito: a dissipação de energia e a auto-excitação. A dissipação de energia reduz e a auto-excitação aumenta a quantidade de energia de um sistema. Na dinâmica dos sistemas com atrito ocorrem os ciclos limites adere-desliza e movimentos caóticos que podem prejudicar o funcionamento das máquinas (PONTES JÚNIOR, 2003). Nos exemplos citados, temos as seguintes formas de energia: 2 1 a mola da suspensão do carro está relacionada à energia potencial elástica de deformação do material e o amortecedor é o dissipador (transformador) de energia mecânica para energia térmica, que é verificada pelo aquecimento do óleo do amortecedor; 2 as rodas do ônibus em movimento sobre o pavimento transmitem energia potencial elástica para o solo; 3 a música é a transmissão de ondas de pressão no ar (pressão acústica) que são percebidas por meio de nosso sistema auditivo. Em alguns problemas de engenharia, as oscilações (vibrações) mecânicas são indesejadas e o controle delas traz benefícios ao sistema. Nos sistemas com o fenômeno das oscilações induzidas por atrito pode ocorrer geração de ruído e/ou oscilações auto-excitadas que prejudicam o funcionamento desses sistemas. A dinâmica de uma superfície de contato com atrito produz uma alternância de estados do movimento, que são conhecidas como modo de aderência (stick em inglês) e modo de deslizamento (slip em inglês). Tal fenômeno pode ser facilmente observado no comportamento dinâmico entre o pneu e o pavimento (estrada) durante a frenagem brusca de um veículo. As oscilações induzidas por atrito que apresentam uma alternância de modos de oscilar são denominadas de oscilações adere-desliza (stick-slip em inglês). O movimento adere-desliza é indesejável, na maioria dos casos, por causa do seu efeito prejudicial sobre a operação e desempenho de máquinas. Sua eliminação ou controle é um dos principais interesses dos engenheiros envolvidos no projeto e operação de máquinas com componentes deslizantes. Podemos citar como exemplo de um sistema de oscilações adere-desliza, com necessidade de controle, o sistema de perfuratrizes de poços de petróleo, onde uma broca muito longa é movida para perfuração de rochas. O atrito com o solo, ao longo do comprimento da broca, sobrecarrega o sistema gerando oscilações devido à elasticidade da broca. Em razão deste fenômeno torna-se necessário a utilização de um sistema de controle do motor da perfuratriz que suprima as oscilações torcionais adere-desliza da broca. Um outro exemplo de fenômeno onde estão presentes as oscilações adere-desliza é nos pneus 3 do trem de pouso de uma aeronave quando tocam o solo. Em um primeiro momento os pneus deslizam e depois aderem ao solo. Em máquinas com transmissão de movimento ou de potência que utilizam componentes flexíveis, denominadas correias ou esteiras, pode aparecer um fenômeno de interação auto-excitadora de oscilações adere-desliza devido ao atrito. A correia ou a esteira transmite o movimento entre rodas ou interage com outros mecanismos através do contato entre suas superfícies. A análise da dinâmica desses sistemas é fundamental no controle das vibrações buscando funcionalidade e conforto ambiental devido à emissão de ruídos. Os efeitos do atrito em sistemas dinâmicos são de grande importância no mundo tecnológico. Atrito é a principal origem do amortecimento em palhetas de turbinas, em pinos flexíveis de juntas de estruturas espaciais como também em robôs com andadores e agarradores (FEENY, 1994). O atrito também está presente no ranger de portas, no som produzido por instrumentos musicais como no violino, no sistema de freios, em máquinas- ferramentas da indústria de manufaturas, em sistemas articulados de robótica, em rodagens de trens em curvas gerando ruídos e muitos outros fenômenos. Por muitos anos, vem sendo estudado analiticamente, numericamente e experimentalmente. Daí a importância deste trabalho no estudo do comportamento de um sistema dinâmico com oscilações adere- desliza e atrito seco. 1.2 Sistema proposto O sistema dinâmico que analisamos consiste de um bloco sobre uma correia inextensível envolvendo duas rodas que giram pela ação de um motor elétrico não-ideal. O bloco preso a um extremo de uma mola linear e o outro extremo preso a um suporte fixo, oscila sobre a correia em movimento. Age sobre ele a força da mola, a de atrito e a externa que consiste de uma força periódica. Este modelo está ilustrado na figura 1-1 e essa situação é o que geralmente ocorre na prática em ciência da engenharia (BALTHAZAR, 1999). 4 Um modelo semelhante a este já foi estudado por NARANAYAN, S. e JAYARAMAN, K., 1991 e neste estudo os autores consideraram constante a velocidade da correia e não-lineares a força da mola e de amortecimento. Em busca de um sistema dinâmico mais simples com oscilações adere-desliza e atrito seco, desprezamos a força de amortecimento do bloco, consideramos variável a velocidade da correia e linear a força da mola. Isto facilitou na compreensão do fenômeno já que o número de parâmetros a serem analisados foi reduzido. Este trabalho é uma continuidade da pesquisa: análise dinâmica por wavelets em um sistema com fricção seca e amortecimento (PEREIRA, 2002), realizada por DANILO CARLOS PEREIRA em outubro de 2002 sob orientação do prof. Dr. JOSÉ ROBERTO CAMPANHA. 1.3 Sistema ideal e não-ideal Um aspecto que consideramos neste trabalho é a interação da dinâmica do sistema mecânico e a sua fonte de energia. Em problemas de engenharia, em geral, não Figura 1-1 Modelo de um sistema dinâmico não-ideal com atrito seco e excitação externa. x dt d m rI fcos( t) k 0 5 é considerada a influência do movimento do sistema oscilante sobre sua fonte de energia. Entretanto, em muitos problemas práticos de engenharia foram observados que a fonte de energia do sistema é influenciada pela resposta do sistema. Esta constatação invalida a formulação tradicional da teoria das oscilações o que torna necessária uma formulação, mais realista, que considere a interação das variáveis de estado da fonte de energia e das variáveis de estado do sistema mecânico (KONONENKO, 1969), (BALTHAZAR et al., 2002). O modelo tradicional, onde a existência do fenômeno de interação da dinâmica do sistema mecânico e a fonte de energia não são consideradas, é denominado de modelo de sistema dinâmico com fonte de energia ideal (máquina ideal). O modelo que considera a interação dinâmica entre o sistema mecânico e a fonte de energia é denominado de modelo de sistema dinâmico com fonte de energia não-ideal (máquina não-ideal) (PONTES JÚNIOR, 2003). Essa interação foi, pioneiramente, detectada por SOMMERFELD, 1903 e recebeu o nome de efeito Sommerfeld. Neste trabalho usamos o modelo de sistema dinâmico com fonte de energia não-ideal. 1.4 Objetivo O objetivo deste trabalho é analisar o comportamento do modelo proposto por meio das séries temporais geradas quando integramos as equações do movimento pelo método numérico de Runge-Kutta. Para essa análise usamos a transformada de Fourier (Apêndice A), a transformada wavelet (Apêndice B), os expoentes de Lyapunov (Apêndice C), os diagramas de bifurcação (Apêndice D) e as seções de Poincaré (Apêndice E). Também usamos os gráficos que dão a posição do bloco em função do tempo (Apêndice A) e as trajetórias de plano de fase (Apêndice A) definidas pela posição versus velocidade do bloco. Calculamos os expoentes de Lyapunov pelo algoritmo LET (SIU, 1998), que é um método já existente, usado no cálculo de expoentes de Lyapunov e disponível para ser usado no MATLAB 7.0. 6 Observamos o movimento de oscilação do bloco sobre a correia, o movimento rotacional do motor, a influência do movimento adere-desliza do bloco sobre o motor e a resposta do motor a esses movimentos. Neste modelo, a força de atrito depende da velocidade relativa de deslizamento descrita pela diferença entre a velocidade da correia e a velocidade do bloco. Para atingir esse objetivo, no capítulo 2, fizemos um breve histórico sobre caos e descrevemos os métodos matemáticos: transformada de Fourier, transformada wavelet e os expoentes de Lyapunov. No capítulo 3, descrevemos as equações diferenciais que tratam do deslocamento angular e velocidade do motor elétrico CC e do deslocamento e velocidade do bloco sobre a correia. No capítulo 4, atribuímos valores aos coeficientes do sistema, fizemos simulações numéricas, usamos para estudo dos casos limítrofes a seção de Poincaré e o diagrama de bifurcação e no capítulo 5, sugerimos a continuidade do trabalho com simulações numéricas do sistema para valores diferentes de massa do bloco e para valores diferentes de amplitude da força externa. Sugerimos também, como continuidade lógica deste trabalho, a construção experimental do modelo aqui proposto. CAPÍTULO 2 Caos e métodos matemáticos no estudo de sistemas dinâmicos 2.1 Caos Uma maneira de estudar o comportamento de certos fenômenos está na formulação de modelos matemáticos (sistemas) desses fenômenos que analisem, no tempo, a variação e evolução dinâmica desses sistemas. Na formulação desses modelos matemáticos as equações diferenciais têm ampla aplicação. Elas podem ser usadas em todo tipo de fenômeno físico, químico, biológico ou social que envolva taxas de variação nos parâmetros do sistema. As soluções obtidas em sistemas de equações diferenciais por meio de métodos numéricos e algoritmos revelam detalhes do comportamento do sistema e, consequentemente, do comportamento do fenômeno que esse modelo representa. Quando por meio de um processo numérico, integramos as equações diferenciais que definem o fenômeno e, para valores muito próximos da condição inicial, não obtemos os mesmos resultados, dizemos que o sistema é caótico e com grande sensibilidade a pequenas mudanças na condição inicial. A principal característica nesses sistemas é que comportamentos passados não se repetem. 8 Citamos como exemplo, as variações climáticas, os fluidos aquecidos, a dinâmica de populações, os osciladores mecânicos, etc. Existem fenômenos, como a queda de um objeto, o som, o movimento dos astros e muitos outros que não são caóticos. 2.2 Breve histórico do caos Apresentaremos a seguir um breve histórico do caos e detalhes sobre o assunto podem ser encontrados em, PEITGEN, 1992, STEWART, 1991, RUELLE, 1993. O estudo sobre sistemas caóticos se desenvolveu com maior intensidade, salvo poucas exceções, a partir de 1975. O norte-americano Edward Lorenz, pesquisador de meteorologia do Instituto de Tecnologia de Massachusetts, no início da década de 1960, formulou um sistema de equações diferenciais, com vista à aplicação em meteorologia, descartando algumas variáveis do modelo proposto por B. Saltzman. Por serem equações de difícil resolução analítica, procurou resolve-las, numericamente, em um computador, com certas substituições intuitivas. Notou que para valores próximos atribuídos à condição inicial, os resultados obtidos se repetiam por um curto espaço de tempo e depois disso o comportamento desses resultados se tornava irregular e imprevisível. Essa irregularidade e imprevisibilidade na evolução temporal de alguns sistemas (não lineares) são denominadas “caos”. Em 1971, os matemáticos David Ruellle, francês, e Floris Takens, holandês, apresentaram um mecanismo que, por intermédio do qual, poderiam aparecer soluções turbulentas da equação usada para descrever o comportamento de um líquido – a equação de Navier-Stokes – quando um dos parâmetros sofria alteração. A turbulência, embora sendo um fenômeno facilmente produzido quando, por exemplo, abrimos o tampão de uma banheira cheia de água, permanece misteriosa e controversa. No entanto, Henri Poincaré (1854-1912) um matemático proeminente e astrônomo teórico que estudou sistemas dinâmicos foi quem reconheceu a impossibilidade da predição, considerando que o conhecimento do estado inicial de um sistema é cercado de incerteza. Ele descreveu como segue: “para pequenas diferenças nas condições iniciais pode ocorrer grandes 9 diferenças no final do fenômeno. Um pequeno erro na formação pode produzir um enorme erro no final. Predições tornam-se impossíveis e nós temos um fenômeno fortuito”. No início da década de 1970, Robert May, físico australiano do Instituto de Estudos Avançados em Princeton, fez uso, na biologia, de um modelo matemático usado para a compreensão da dinâmica das populações. A equação do segundo grau y = kx(1 x), (2.1) concebida pelo matemático belga Pierre François Verhulst em 1845, foi usada por Robert May para simular comportamento de epidemias. Não se sabe por que esta equação ficou conhecida pelo nome de equação logística, mas é um fato notável que ela tenha sido muito usada no estudo de populações no século 20. Robert May verificou que os resultados obtidos no estudo dessa equação coincidiam com os dados reais de uma epidemia, na qual o número de infectados, em um determinado momento, aumentava abruptamente e num momento posterior diminuía drasticamente. Se admitirmos k o valor da razão de crescimento da população em estudo e, no instante n, xn a porcentagem de infectados e 1 – xn a porcentagem de não infectados dessa população, então a população no instante n + 1, pode ser obtida pela equação de recorrência: xn+1 = kxn (1 xn). (2.2) conhecida pelo nome de mapa logístico. Considerando uma população inicial de 10% de infectados, isto é, x0 = 0,1, os valores x1, x2, x3, serão encontradas por meio do mapa logístico (2.2), e representarão, respectivamente, as porcentagens da população infectada, no fim do 10 primeiro período, no fim do segundo período, no fim do terceiro período e assim por diante. Robert May notou que, para valores crescentes de k, a população, em função desse crescimento, num primeiro momento apresentava tamanho constante, depois oscilava, regularmente, entre valores altos e baixos, e mais adiante, oscilava subitamente duas vezes mais depressa até o momento em que esses valores passavam a mudar irregularmente sem qualquer previsibilidade. Esse resultado foi surpreendente porque enfraquecia um dos dogmas fundamentais da ciência: equações matemáticas eram consideradas a forma mais elevada para expressar princípios da natureza, e as soluções de equações matemáticas que descrevem sistemas naturais eram tidas como repetíveis, não importando quem fizesse o cálculo ou quantas vezes fossem eles repetidos. Esta é a aplicação básica da matemática à ciência — fazer previsões precisas e repetíveis. Robert May demonstrou que equações escritas para descrever processos naturais podem, sob certas circunstâncias, dar resultados imprevisíveis. Desde a descoberta de Robert May, comportamentos caóticos foram encontrados em muitas outras áreas, tais como epidemias, ritmos cardíacos, ciclos econômicos e fluxo de fluidos. No mapa logístico, para certos valores do parâmetro de controle, o sistema mostra um comportamento regular, mas ao atingirmos certo valor crítico deste parâmetro, o sistema passa a exibir bruscamente o comportamento caótico. Nessa região caótica, qualquer incerteza inicial na especificação de x0 crescerá exponencialmente com o crescimento do número n de iterações. A descoberta do caos em sistemas físicos levou a um novo entendimento das leis da natureza. O caos é inevitável, mesmo para sistemas físicos muito simples. Por um lado, existe uma ordem não esperada dentro do caos, em virtude das simetrias no movimento regular que o suporta. Sistemas não caóticos são raros, embora sirvam quase sempre de base para nossa compreensão física da natureza. Por outro lado, para a classe dominante de sistemas caóticos, erros iniciais de observação em geral crescem exponencialmente e o determinismo torna-se sem sentido numa curta escala de tempo. Assim, se a precisão infinita deve ser abandonada, é necessário que atentemos para o fato de que as equações determinísticas, não garantem a capacidade de previsibilidade, em virtude das incertezas nas condições iniciais dos sistemas a que se aplicam. 11 2.3 Métodos matemáticos no estudo de sistemas dinâmicos As simulações numéricas no sistema dinâmico proposto geram séries temporais que, analisamos pelos seguintes métodos matemáticos: a) Transformada rápida de Fourier; b) Transformada wavelet; c) Expoentes de Lyapunov. Na análise periódica do diagrama do plano de fase usamos o mapa de Poincaré e na análise da súbita mudança qualitativa no comportamento dinâmico do sistema quando fazemos uma pequena mudança no valor de um parâmetro usamos o diagrama de bifurcação. Esses dois métodos serão explorados no capítulo 4. 2.3.1 Série de Fourier A evolução no tempo de um sistema dinâmico pode ser representada por uma função f(t) dependente do tempo t ou por uma série temporal, com a condição de que t seja escolhido, em intervalos regulares de tempo. Dependendo do tipo da função f(t), sua representação pode ser feita de dois modos diferentes. Se f(t) é periódica, então o espectro deve ser expresso como uma combinação linear de senos e cossenos, cujas freqüências são múltiplos inteiros de uma freqüência básica. Essa combinação linear é denominada série de Fourier. Porém, se f(t) não é periódica, o espectro deve ser expresso pela transformada de Fourier de f(t) (BAKER, 1996). Essa representação é útil em dinâmicas caóticas, pois a transformada de Fourier é, em geral, uma função de valores complexos e, muitas vezes, é preferível transforma-la em uma função de valor real, usando para isso a raiz quadrada do módulo da transformada. Essa função real é chamada de espectro de potência de f(t). 12 O matemático francês, Joseph Fourier, em 1807 demonstrou que toda função periódica f(t), de período 2L, definida no intervalo [ L, L] com t R pode ser expressa como uma somatória das funções trigonométricas seno e cosseno. Desse modo, a equação matemática que define f(t) assume a forma: 0n nnnn )]t(senb)tcos(a[)t(f (2.3) onde an e bn são chamados coeficientes de Fourier, cujos valores são deduzidos do fato de que as funções seno e cosseno formam uma base ortogonal e são dados por: L L 0 dtf(t) L 1 a ; L L nn 0ndt,t)cos(f(t) L 1 a ; (2.4) L L nn 0ndt,t)sen(f(t) L 1 b . Se o intervalo [ L, L] é dividido em n partes iguais, então cada subintervalo resultante dessa divisão, terá período igual a: n L2 Tn . (2.5) Consequentemente, a freqüência angular desse subintervalo é: n nn T 2 f2 (2.6) Substituindo (2.5) em (2.6) obtemos: 13 L n n L2 2 n (radianos por unidade de tempo) (2.7) O período Tn, relacionado com a freqüência angular é obtido da equação (2.6) e dado por: n n 2 T (unidades de tempo) (2.8) Substituindo (2.7) em (2.3), a série de Fourier assume a forma: 0n nn )]t L n (senb)t L n cos(a[)t(f (2.9) Dizemos que a equação (2.9) é uma representação espectral de f(t). O objetivo de usar a série de Fourier na análise de uma função periódica é obter a amplitude, em função da freqüência para cada onda senoidal que compõe o sinal analisado e com isso, construir uma série discreta de freqüências a partir de uma série temporal. Utilizando a fórmula de Euler, podemos reescrever a série de Fourier na forma complexa: n ti n necf(t) (2.10) em que os coeficientes cn são obtidos por: 3,2,1,0n,dte)t(f L2 1 c L L ti n n (2.11) 2.3.2 Transformada de Fourier 14 A série de Fourier é usada para funções periódicas, enquanto que, a transformada de Fourier é para funções não periódicas. Se a função é periódica tal, que f(t) = f(t + nT) com n sendo inteiro positivo ou negativo e T sendo a periodicidade básica, então as freqüências das várias componentes espectrais são todas inteiras, múltiplas da freqüência básica 0 = 2 /T. Assim, a série de Fourier, na forma complexa, que representa a função f(t) pode ser escrita da seguinte forma: n n tin 0ec)tf( (2.12) onde Cn são as amplitudes das componentes de freqüências n 0. dtef(t) 2 c 0 0 00 tin n . (2.13) A transformada de Fourier é uma extensão da série de Fourier em que a periodicidade básica T de f(t) pode ficar infinitamente grande. Assim, quando T tende para infinito, o espaço entre as componentes das freqüências torna-se infinitesimal e o espectro discreto das componentes da freqüência torna-se contínuo. Por essa razão n 0 tende para que é uma variável contínua e cn tende para c( )d onde d é um pequeno intervalo de freqüência e c( ) é a amplitude dependente da freqüência ou transformada de Fourier (BAKER, 1996). Nessas condições a equação (2.12) é então escrita na forma: de)(c)t(f ti (2.14) e a equação dos coeficientes (2.13) torna-se: 15 dtde)t(f 2 1 d)(c ti (2.15) de onde deduzimos o valor de c( ), que é igual a: tde)t(f 2 1 )(c ti (2.16) As equações (2.14) e (2.16) são, respectivamente, a função original f(t) e a transformada de Fourier. Com a equação (2.16) temos condição de converter um conjunto de sinais digitais, no domínio do tempo, em um conjunto de pontos no domínio das freqüências como também, reconstruir os sinais gerados pela função original, multiplicando os coeficientes de Fourier por senóides e cossenóides de freqüências apropriadas. 2.3.3 Transformada discreta de Fourier Uma série temporal não corresponde a uma função contínua e infinita, e sim, a um conjunto discreto de valores que muitas vezes obtemos por meio da simulação computacional de um modelo matemático. O método de Runge-Kutta, aplicado na resolução do problema proposto neste trabalho, gera uma série temporal discreta com um número finito, n, de pontos, para intervalos constantes de tempo. Consequentemente, para o cálculo dos coeficientes de Fourier, as integrais apresentadas em (2.4) são calculadas usando um método numérico de integração que será explicado a seguir. Os pontos = t0, t1, t2, tn = dividem o intervalo [ , ] em n subintervalos iguais, cujo comprimento é definido por: 16 n 2 t (2.17) Designando os valores da função f(t) nos pontos t0, t1, t2, tn, respectivamente, por f(t0), f(t1), f(t2), f(tn) e utilizando, a fórmula dos retângulos (PISKOUNOV, 1997) determinamos os coeficientes de Fourier, substituindo as integrais apresentadas em (2.4) por somatórios. Assim, n 1k k )f(tt 1 a 0 (2.18) Substituindo (2.17) em (2.18) temos: n 1k k )f(t n 2 a 0 (2.19) Analogamente, calculamos os coeficientes an e bn: )cos(nt)f(t n 2 a k 1k kn n (2.20) )sen(nt)f(t n 2 b k 1k kn n (2.21) Para o caso discreto, a transformada de Fourier dada pela equação (2.16) tem sua integral substituída por um somatório e obtemos: ti 1k k ne)f(t n 1 nc n (2.22) A transformada de Fourier inversa é dada por 17 ti n kecf(t) 1k k (2.23) Mesmo usando computadores modernos, certos cálculos numéricos tomam muito tempo. Com o objetivo de reduzir esse tempo, Tukey e Cooley, em 1965, desenvolveram um algoritmo chamado Fast Fourier Transform (FFT) para o cálculo da transformada discreta de Fourier permitindo, de forma rápida, encontrar o espectro de freqüência de um sinal. Esse algoritmo reduz o número de operações computacionais. Para uma quantidade de N pontos ele transforma 2N2 operações em uma quantidade de 2Nlog2N. A partir do espectro do sinal, podemos caracterizar se o comportamento do sistema é periódico ou caótico. No entanto, deve ser aqui salientado que quando passamos do domínio do tempo para o domínio da freqüência a informação temporal é perdida. 2.3.4 Wavelets Apresentamos a seguir um breve histórico das wavelets e detalhes sobre o assunto podem ser encontrados em MORETTIN, 1999, DAUBECHIES, 1992. A palavra wavelet derivada do vocábulo em francês ondelette, que denota o diminutivo de onda, surgiu em meados dos anos 80, a partir dos trabalhos de um grupo de pesquisadores franceses. O termo ondelette foi introduzido pelo engenheiro geofísico Jean Morlet (1980), sendo a base matemática de suas idéias formalizada pelo físico teórico de mecânica quântica Alex Grossmann, que o ajudou a formalizar a transformada wavelet em sua forma contínua. Os dados estudados por Morlet exibiam conteúdos de freqüências que mudavam rapidamente ao longo do tempo. Nesse caso a transformada de Fourier não era adequada 18 como método de análise. Na realidade eles redescobriram e deram uma interpretação ligeiramente diferente do trabalho de Alberto Caldéron sobre análise harmônica de 1946. Foi Yves Meyer, um matemático francês, que em 1984 ressaltou uma semelhança no trabalho de Morlet e Caldéron. Isto o levou, em 1985, a construir funções básicas (wavelets ortonormais) que, quando usadas na transformada wavelet, geravam um domínio de tempo e freqüências satisfatórias na solução de problemas. Mas Meyer verificou que, cerca de cinco anos antes, J. O. Strömberg já tinha descoberto as mesmas wavelets. Em 1909, num apêndice de sua tese de doutorado, o matemático alemão, Alfred Haar já propunha um conjunto de funções base de wavelets ortonormais. Suas wavelets eram de pequeno uso prático e de pobre domínio de freqüência. Em 1930, Paul Levey utilizando o trabalho de Haar desenvolveu funções básicas ortonormais para estudar os sinais aleatórios do movimento browniano. Na década de 90, Ingrid Daubechies, uma estudante de graduação da Universidade de Bruxelas, desenvolveu sistemas para discretização de parâmetros de tempo e escala da transformada wavelet. Estes sistemas apresentavam uma maior liberdade na escolha das funções básicas. Daubechies, com Stephane Mallat desenvolveram a transição da análise de sinais contínuos para discretos. Em 1986, Mallat e Meyer desenvolveram a idéia de análise de multiresolução (MRA) para transformada wavelet discreta (DWT). Daubechies utilizando os trabalhos de Mallat formalizou a teoria moderna de wavelet desenvolvendo as bases ortonormais de wavelets suaves com suportes compactos. Nos últimos anos, tem se verificado muita pesquisa para outras funções básicas wavelets com diferentes propriedades e modificações no algoritmo de MRA. Em 1992, Daubechies com Albert Cohen, Jean Feauveau construíram as wavelets biortogonais com suporte compacto que são utilizados por muitos pesquisadores sobre as funções de base ortonormais. Desde então as wavelets tem-se desenvolvido com o aparecimento de novas funções-base de wavelets, novas aplicações e novos métodos, constituindo um novo e sólido campo de pesquisas. 19 Definição de wavelet Wavelet é uma “pequena onda” que, matematicamente, expressamos por uma função (t), que pode ser suave ou não, simétrica ou não e pode ter expressão matemática simples ou não. Ela deve integrar para zerar, “ondeando” acima e abaixo do eixo t e, principalmente, assegurar rapidez e facilidade no cálculo da transformada wavelet direta e indireta. Citaremos a seguir alguns exemplos de funções wavelets, para melhor ilustrar a definição dada. Wavelet de Haar A função (t) que define a wavelet de Haar é: contráriocaso0, 1 t 2 1 se1, 2 1 t0se1, (t) (2.24) e seu gráfico está representado na figura 2-1. 20 Wavelet chapéu mexicano Essa wavelet recebe esse nome pelo fato de seu gráfico ser parecido ao de um chapéu mexicano e é definida pela função: 22te)t(1(t) 2 (2.25) e seu gráfico está representado na figura 2-2. (t) 0,5 –1 0 1 1 t Figura 2-1 Wavelet de Haar Figura 2-2 Wavelet chapéu mexicano 21 Wavelet de Morlet Essa wavelet também denominada gaussiana modulada é uma função complexa, para 0 fixo e é definida por: 22 0 tti ee(t) (2.26) A equação (2.26), pela fórmula de Euler, pode ser decomposta em uma parte real e outra imaginária. 2/2tet)isent(cos(t) 00 (2.27) onde a parte real da wavelet de Morlet é dada pela equação: t)cos(e(t) 0 2/2t (2.28) Em (2.28), se atribuímos à freqüência 0 um valor, como por exemplo 0 = 5, temos a parte real da wavelet de Morlet, cujo gráfico está apresentado na figura 2-3. Figura 2-3 Wavelet de Morlet 22 Escala e translação Na análise de Fourier, o objetivo básico está em aproximar uma função f(t) por uma combinação linear de componentes senoidais, cada uma com uma freqüência dada. O conjunto { n(t) = eint, n = 0, 1, 2, } de funções ortogonais, de período 2 , forma a base para a análise de Fourier. Na realidade, esse conjunto é gerado por dilatações de uma única função (t) = eit, ou seja, n(t) = (nt) para qualquer n inteiro. O fato básico é que toda função periódica, de período 2 , de quadrado integrável, é gerada por uma superposição de dilatações inteiras da função (t). No entanto a análise de Fourier é apropriada para analisar os processos estacionários. Para processos não estacionários e processos com características especiais, outros sistemas ortogonais podem ser úteis, como as wavelets (MORETTIN, 1979). Enquanto a análise de Fourier consiste em decompor um sinal em senos e cossenos de várias freqüências a análise wavelet é a decomposição por meio de escalas (dilatação e compressão) e de translações. Para interpretarmos a dilatação (ou compressão) e uma wavelet, consideremos (t) a wavelet de Haar definida em (2.24). Calculando (2t) e (4t) temos: contráriocaso0, 2 1 t 4 1 se1, 4 1 t0se1, t)2( contráriocaso0, 4 1 t 8 1 se1, 8 1 t0se1, t)4( (2.29) 23 Os gráficos das wavelets de Haar (t), (2t) e (4t) estão apresentados na figura 2-4: Observando os domínios das wavelets de Haar (t), (2t) e (4t) na figura 2-4, notamos que cada domínio é obtido multiplicando o domínio da wavelet mãe por, respectivamente, 1, 1/2 e 1/4. Estes valores são denominados fatores de escala e são representados, genericamente, por j2 1 a com j = 0, 1. 2... . O valor a é denominado fator de escala. Analogamente, o fator de escala é também aplicado em outras wavelets. Assim, fatores de escala decrescentes fazem a compressão da wavelet e fatores de escala crescentes fazem a sua dilatação. Isso pode ser observado quando aplicamos essas mesmas escalas na wavelet de Morlet conforme figura 2-5. (4t) (2t) (t) Figura 2-5 Dilatação da wavelet de Morlet. Wavelets de Morlet nas escalas (a) 1/4, (b) 1/2 e (c) 1. (a) (b) (c) 0 10,5 0 1 0 10,5 –1 0 1 –2 (t) (4t) (a) (c) a = 1 a = 1/4 0 10,5 –1 0 -2 (2t) (b) a = 1/2 Figura 2-4 Compressão da wavelet de Haar. Wavelets de Haar nas escalas (a) 1, (b) 1/2 e (c) 1/4. -1 -2 2 1 2 2 24 Se quisermos deslocar uma wavelet para um ponto qualquer b, afastado da origem do eixo, basta calcularmos (t – b), conforme mostra a figura 2-6. Na dilatação, (ou compressão) de uma wavelet (t) é preciso usar uma escala a que satisfaça a condição de a < 1 para contrair e de a > 1 para dilatar a wavelet, enquanto que, na translação temos que usar o parâmetro b para mudar a origem do eixo t para t – b. Alterando escala e posição da wavelet (t), denominada wavelet “mãe”, geramos uma família de wavelets, denominadas wavelets “filhas”, que são, matematicamente, definidas por: a bt a(t) 2/1 ba, a, b R (2.30) onde, usualmente, os valores atribuídos à a e b são: a = 2 j e b = k 2 j, j, k Z. Usando a wavelet de Haar em (2.29) e atribuindo valores a j e k, obtemos algumas wavelets “filhas”. Por exemplo: 1. 0,0 = 20/2 (20 t 0) = (t); 2. 1,0 = 21/2 (21 t 0) = 2 (2t); 3. 1,1 = 21/2 (21 t 1) = 2 (2t 1); 4. 2,0 = 22/2 (22 t 0) = 2 (4t); 5. 2,1 = 22/2 (22 t 1) = 2 (4t 1); 6. 2,2 = 22/2 (22 t 2) = 2 (4t 2); 7. 2,3 = 22/2 (22 t 3) = 2 (4t 3). (t) (t – 4 ) (t – 8 ) Figura 2-6 Translação de wavelets. 25 O gráfico da wavelet “mãe” e de algumas wavelets “filhas” estão apresentados na figura 2-7. Nos gráficos da figura 2-7, à medida que os valores atribuídos a j crescem, notamos que as wavelets “filhas” tornam-se mais estreitas e mais alongadas, isto é, com amplitudes maiores. Os valores crescentes atribuídos a k, fazem as wavelets “filhas” se movimentarem para a direita, dentro do intervalo de tempo t. Dessa forma, quando j e k 0,0(t) 0 10,5 –1 0 1 –2 Figura 2-7 Wavelets de Haar geradas a partir da wavelet “mãe”. 1,0(t) 2 1 0 –1 –2 0 10,5 1,1(t) 2 1 0 –1 –2 0 10,5 2,3(t) 2 1 0 –1 –2 0 10,5 2,0(t) 2 1 0 –1 –2 0 10,5 2,1(t) 2 1 0 –1 –2 0 10,5 2,2(t) 2 1 0 –1 –2 0 10,5 2 26 Figura 2-8 Escalogramas de funções com freqüências previamente conhecidas com escalas em Hz e tempo em segundo. Em (a) a função cos(8 t) e em (b), a função sen ( t). (a) (b) assumem uma grande quantidade de valores inteiros, o sinal a ser analisado que estiver dentro do intervalo de tempo t será varrido pelas wavelets “filhas”. A transformada wavelet gera os coeficientes wavelets. Cada um dos coeficientes dj, k multiplicados por apropriadas wavelets, dependentes de escala e posição, obtidas de uma wavelet “mãe” (t) reconstroem o sinal original ou função. Quanto maior o coeficiente maior é a semelhança entre a wavelet e o sinal, e quanto menor o coeficiente menor a semelhança. Assim, reconstruímos a função f(t) usando a relação que segue: 2,12,12,02,01,11,11,01,00,00,0 dddddf(t) (2.31) Um sinal pode também ser representado por meio de um diagrama tempo-escala conhecido como escalograma. Nele, os coeficientes wavelets são representados por intensidades de cores, onde, usualmente, as cores frias estão associadas aos coeficientes menores e as cores quentes aos maiores. Os gráficos da figura 2-8 com escala em Hz e tempo em segundos (zero até 85s) mostram os coeficientes wavelets calculados para as funções (a) cos (8 t) e (b) sen ( t) por meio das cores. As cores quentes se concentram nas proximidades da escala 4 Hz no escalograma (a) e 0,5 Hz no escalograma (b) indicando, respectivamente, as freqüências das funções cos(8 t) e sen( t). 27 2.3.5 Transformada wavelet contínua A transformada wavelet apareceu em sua forma contínua com os trabalhos de dois pesquisadores franceses, J. Morlet, um geofísico, e A. Grossmann, um físico teórico (MORETINN, 1999). Na análise de Fourier é feita uma decomposição do sinal em funções senos e cossenos com freqüências diversas e na análise usando a transformada wavelet a decomposição do sinal é feita por meio da função wavelet que será modificada com dilatações e compressões (uso da escala) e será transladada ao longo de todo sinal. Com isto, teremos indicação das freqüências presentes no sinal, em que tempo elas ocorrem (informação não fornecida com a transformada de Fourier) e, também, evita a possibilidade de escolha errada da janela na transformada de Gabor. Segundo DAUBECHIES, 1992, LIMA, 2002, MORETTIN, 1979, podemos definir a transformada wavelet como: td)(tf(t))b,a(C ba, (2.32) ) a bt ( a 1 (t)ba, (2.33) onde (t) é a função (wavelet) escolhida para se fazer a análise do sinal f(t). O valor a é o parâmetro de escala e b o parâmetro de localização com a e b R e a 0. Mudando o valor de a tem-se o efeito de dilatação (a > 1) ou de contração (a < 1), enquanto que mudanças no parâmetro b têm o efeito de analisar a função f(t) em torno deste ponto. A função (t) deve possuir as seguintes características (DAUBECHIES, 1992): a) satisfazer as condições de normalização; 28 1dt(t) 2 (2.34) b) decair suficientemente rápido para se obter uma boa localização; dt(t) (2.35) c) satisfazer a condição de admissibilidade; dC 2)( (2.36) onde ( ) é a transformada de Fourier de (t). Isso é o que garante a existência da transformada wavelet inversa: 2a dadb ) a bt ( a 1 b)C(a, C 1 f(t) (2.37) d) a média ser igual a zero, (t) comporta-se como uma onda (daí o nome wavelet). 0dt(t) (2.38) 2.3.6 Transformada wavelet discreta As propriedades de wavelets, descritas anteriormente, de maneira geral, também, são verdadeiras na análise de wavelet discreta. 29 Consideremos o conjunto de valores T 1 0 X X X X (2.39) que representam amostras de um sinal observado ou gerados por uma função f(t). Supondo T = 2M , M > 0, inteiro, definimos a transformada wavelet discreta de X, com respeito a wavelet “mãe” (t), como: )T/t(Xd k,j 1T 0t tk,j (2.40) que denotamos, simplesmente k,jd . Essa transformada é calculada para j = 0,1, , M – 1 e k = 0, 1, , 2j – 1 (MORETTIN, 1999). Assim, na transformada wavelet discreta, os sinais são representados por séries do tipo: j k kj,kj, (t)df(t) (2.41) em que as funções bases dadas por: k)t(2(t) j k,j j, k Z, (2.42) são chamadas wavelets “filhas” e são geradas a partir de dilatações e translações de uma única função (t) também denominada wavelet “mãe”. 30 As funções wavelets (t)k,j são obtidas de (t) por uma dilatação e uma translação. As funções { j,k(t), j, k Z} formam uma base que não precisa ser necessariamente ortogonal. Uma das vantagens de se trabalhar com bases ortogonais é que elas permitem a reconstrução perfeita do sinal original a partir dos coeficientes da transformada de forma mais econômica. Consideremos numa base ortogonal gerada por )t( onde os coeficientes de wavelets são dados por dtk)t2f(t)2d j2j/ kj, (2.43) A transformada wavelet é a aplicação que associa um sinal aos seus coeficientes wavelets. kj,df(t) . (2.44) Os índices j e k representam, respectivamente, escala e translação. Assim, cada um desses coeficientes dj,k multiplicados por apropriadas wavelets, dependentes de escala-posição, obtidas de uma wavelet “mãe” (t), reconstroem o sinal original. A expansão de uma função ou sinal por funções base wavelets não é única. Existem várias funções wavelets diferentes que podem ser usadas como funções base. 2.3.7 Obtenção da transformada wavelet 31 Figura 2-9 Análise wavelet em uma seção do sinal. Wavelet Toolbox-(MISITI, 1996). Figura 2-10 Análise wavelet em nova seção do sinal. Wavelet Toolbox-(MISITI, 1996). Para obtermos a transformada wavelet de um sinal, devemos proceder conforme segue: passo 1: tomar uma wavelet mãe e compará-la com uma seção no começo do sinal original; passo 2: calcular o coeficiente C a partir da transformada wavelet. Este coeficiente representa a comparação da wavelet mãe com esta seção do sinal. Quanto maior o valor de C maior é a semelhança entre a wavelet e o sinal da seção tomada. passo 3: tomar a próxima seção do sinal e novamente repetir os passos 1 e 2 até cobrir, de seção em seção, todo o sinal; passo 4: ajustar uma nova escala para a wavelet (dilatar ou comprimir) e repetir os passos de 1 até 3; 32 passo 5: repetir os passos de 1 a 4 para todas as escalas da wavelet (tantas vezes quanto for possível ou desejado). Quando terminamos todos os passos descritos acima, produzimos coeficientes wavelets em diferentes escalas e seções do sinal. Os coeficientes wavelets das escalas maiores estão associados a wavelets mais dilatadas e, escalas menores, a wavelets mais comprimidas. Note que quanto mais dilatada for a wavelet, maior será a seção do sinal com o qual ela estará sendo comparada e, deste modo, as características mais visíveis estão sendo medidas pelos coeficientes wavelets. Assim, existe uma correspondência entre as escalas wavelets e a freqüência revelada pela análise wavelet, conforme vemos na figura 2-12. Figura 2-11 Análise wavelet com uma wavelet em nova escala. Wavelet Toolbox-(MISITI, 1996). escalas pequenas wavelets mais comprimidas detalhes rapidamente variáveis alta freqüência no sinal escalas grandes wavelets mais dilatadas características mais visíveis mudando lentamente baixa freqüência no sinal Figura 2-12 Correspondência entre escalas e freqüências na análise wavelet. Wavelet Toolbox-(MISITI, 1996). 33 2.3.8 Expoentes de Lyapunov Os expoentes de Lyapunov, nome em homenagem ao matemático russo, A. M. Lyapunov, (BAKER, 1996), (FERRARA, 1994) desempenham um papel crucial na descrição do comportamento de sistemas dinâmicos. Eles medem a taxa média de divergência ou convergência das trajetórias do espaço de fase a partir de pontos iniciais próximos. Por essa razão, eles podem ser usados para analisar a estabilidade dos ciclos limites e para checar a sensível dependência às condições iniciais indicando a presença de atratores caóticos (SANDRI, 1996). Consideremos, inicialmente, sistemas contínuos com n equações diferenciais ordinárias (FERRARA, 1994) e imaginemos uma pequena hiper-esfera de condições iniciais no espaço de fase para escalas de tempo suficientemente pequenas. O efeito da dinâmica do sistema será distorcer este conjunto para um hiper-elipisóide esticando ao longo de algumas direções e contraindo ao longo de outras. A taxa assintótica de expansão do eixo maior que corresponde à direção mais instável do fluxo é medida pelo maior expoente de Lyapunov 1. Em geral, se ordenarmos os eixos e os expoentes de Lyapunov em ordem decrescente pela magnitude, 1 ... n e 1 ... n, cada i quantifica a taxa exponencial média de expansão ou contração para o i-ésimo eixo i. Mais formalmente, consideremos dois pontos iniciais próximos x0 e y0 em um espaço de fase, conforme figura 2-13 onde y0 é uma pequena hiper-esfera de raio 0(x0), representada matematicamente na relação (2.45), cuja finalidade é de teste aos estados iniciais vizinhos em torno do ponto x0 de uma linha de fluxo. )(xxy 0000 (2.45) Com o passar do tempo, o fluxo deforma a hiper-esfera que se transforma em um hiper-elipsóide com eixos principais i(t), para i = 1, 2,..., n, e que está representado na figura 2-13. 34 Os expoentes de Lyapunov medem o crescimento exponencial dos eixos principais i(t) e são definidos por n,2,1,i )( (t) ln t 1 limlim 0000 x0x i )(t i (2.46) Em geral os i dependem do estado inicial x0, mas em muitos casos eles são constantes ao longo de uma significativa região do espaço de fase (FERRARA, 1994). Da equação 2.46 obtemos: t 00i ie)(x(t) . (2.47) Concluímos que (FERRARA, 1994): a) a existência de um ou mais expoentes de Lyapunov positivos define uma instabilidade orbital nas direções associadas; 0(x0) 1(t) 2(t) x0 Figura 2-13 Evolução de um elemento de volume esférico de raio 0(x0) em torno de um ponto inicial x0. Depois de um tempo t a esfera torna-se um elipsóide com eixos principais 1(t) e 2(t). Na figura está representado o caso bidimensional. 35 b) para uma solução caótica, associada a um atrator estranho, a dependência sensitiva as condições iniciais implica na existência de pelo menos um expoente de Lyapunov i > 0; c) para uma solução periódica ou quasi-periódica podemos esperar que deslocamentos na direção perpendicular ao movimento diminuam com o tempo, enquanto que ao longo da trajetória eles não devem se alterar, correspondendo a um simples deslocamento do ponto inicial. Segue, portanto, de (2.47) que no caso de solução periódica ou quasi-periódica i < 0 nas direções perpendiculares ao movimento e i = 0 ao longo da trajetória. É possível identificar um atrator pelos sinais dos expoentes de Lyapunov. Assim, a dinâmica de um atrator para um sistema contínuo com quatro equações diferenciais será de um movimento periódico se o espectro de Lyapunov apresentar um expoente nulo e os demais negativo, caótico, se apresentar um positivo, um nulo e os outros dois negativos e ponto fixo se todos forem negativos (SANDRI, 1996). Esta classificação está descrita abaixo: Dinâmica dos atratores Espectros de Lyapunov Ponto fixo , , , Movimento periódico 0, , , Movimento caótico +, 0, , . Assim, um sistema contínuo com quatro equações diferenciais ordinárias será caótico, com base no espectro de Lyapunov, se tiver um expoente positivo, um nulo e os outros dois negativos e periódicos se o espectro tiver um nulo e os demais negativos. CAPÍTULO 3 Equações que descrevem o sistema 3.1 Introdução Quando um sistema recebe energia de uma fonte, geralmente, admitimos que a fonte não reaja ao sistema. Assim, se consideramos a fonte de energia como sendo um motor girando com velocidade angular constante, usualmente, admitimos que esta velocidade permaneça constante e que o sistema não interfere no movimento rotacional do motor (KONONENKO, 1969). No entanto, em problemas práticos da ciência de engenharia, geralmente, a excitação é dependente da resposta estrutural do sistema dinâmico. Isto é, a excitação manifesta influência no sistema e o sistema age na fonte de excitação. Neste caso, o sistema dinâmico é dito não-ideal e a excitação depende da resposta do sistema. Este fenômeno é conhecido na literatura científica como efeito Sommerfeld, em homenagem ao primeiro pesquisador do assunto. No modelo proposto, conforme figura 3-1, nós admitimos constante a voltagem que alimenta o motor elétrico CC (motor de corrente contínua), a qual implica em velocidade angular constante. No entanto, essa velocidade nem sempre é constante, 37 podendo oscilar devido à influência do movimento de oscilação adere-desliza do bloco no motor elétrico. No estudo deste modelo consideramos uma correia inextensível e sobre ela um bloco preso a um extremo de uma mola linear com o outro extremo preso a um suporte fixo. Esse bloco oscila sobre a correia movimentada pelo motor elétrico. Influenciando na oscilação do bloco temos a força de atrito, a força da mola e a força externa, que serão detalhadas a seguir. 3.2 Forças que agem no bloco 3.2.1 Força de atrito Apesar de difícil a tarefa de modelar o atrito seco em um oscilador, visto que a física do atrito seco não é completamente entendida (FENNY, 1994), dos diferentes modelos existentes para reproduzir esse fenômeno, usamos o modelo que envolve a lei de Coulomb que trata do atrito entre dois sólidos deslizantes, publicada por Augustin de Coulomb e que diz: x dt d m r k I fcos( t) 0 Figura 3-1 Modelo de um sistema dinâmico não-ideal com atrito seco e excitação externa. 38 1 – a força de atrito é independente da área de contato aparente entre as superfícies deslizantes; 2 – a força de atrito é proporcional à força normal N, isto é, à força perpendicular às superfícies deslizantes que pressiona os dois sólidos juntos. O fator de proporcionalidade é chamado coeficiente de atrito; 3 – o atrito cinético fk, que é a força para manter um corpo deslizante em uma velocidade constante, não depende da velocidade de deslizamento. Ele é menor ou igual ao atrito estático fs que é a força para iniciar o deslizamento. Temos então que: fs = sN, fk = kN, com k s. (3.1) No sistema dinâmico proposto modelamos a força de atrito fa, causada pela interação do bloco e da correia em movimento, em função da velocidade relativa de deslizamento, da força normal N e dos coeficientes de atrito estático s e cinético k (HECKL, 1996). A velocidade relativa de deslizamento que é a diferença da velocidade do bloco é definida por xru onde representa a velocidade angular do motor e x representa a velocidade e r o raio da roda envolta pela correia. Em um modelo tradicional, a força de atrito fa, figura 3-2, é proporcional à força normal N e é definida pela equação: fa = N (3.2) Figura 3-2 Força de atrito com s e k diferentes. fa u 39 No modelo em estudo, segundo a lei de Coulomb, o coeficiente de atrito varia em função da variação da força normal na área em contato. Isso acontece porque existe um coeficiente de atrito estático s e um coeficiente de atrito cinético k (FENNY, 1994). Quando o bloco tem a mesma velocidade da correia, a velocidade relativa de deslizamento é nula. Nessa condição, o atrito estático s toma o valor necessário para equilíbrio das forças que agem no bloco. Consequentemente, para manter essa condição de equilíbrio a força de atrito passa a assumir muitos valores já que o atrito estático está variando e é definida pela equação (FENNY, 1994): fa(u)= h(u) s N (3.3) onde 1 0 e h(u) = 1 para u < 0. Quando admitimos (LYANG, 1965) s = k = e consideramos que N = mg, a força de atrito fa(u), que é descontínua e com muitos valores quando u = 0, assume a forma: fa(u) = h(u) mg (3.5) com h(u) definida pela equação: 40 0u1 Rye0u1y1 0u1 )u(h (3.6) e representada pelo gráfico da figura 3-3: Para evitar a descontinuidade da lei de atrito de Coulomb e poder usar a integração numérica de funções contínuas, substituímos a equação (3.6) pela equação (LIANG, 1995): u)tanh((u)h* (3.7) onde é um parâmetro que determina a proximidade entre a descontinuidade do modelo de Coulomb e essa aproximação contínua. Utilizando esse modelo, a força de atrito convencional passa a ser definida pela equação: u)mgtanh(µfa (3.8) Atribuindo a (unidade em s/m) o valor de 30 a equação (3.8) passa, então, a ser escrita na seguinte forma: Figura 3-3 Gráfico da função h(u) com s e k iguais. h(u) 1 1 u 41 )utanh(30gmµfa (3.9) A figura 3-4, obtida com recursos do MATLAB 7.0 para u variando entre 1,3 e 1,3, mostra a função h*(u) = tanh(30u), substituta da função h(u). 3.2.2 Força da mola A força da mola também contribui para o movimento oscilatório do bloco e obedece a lei de Hooke que é enunciada assim: quando um sólido é deformado, ele resiste a isto com uma força proporcional à deformação experimentada, desde que esta não seja demasiado grande. Em se tratando de uma deformação unidimensional, essa lei pode ser escrita como: F = – k x (3.10) Figura 3-4 Gráfico da função h*(u) = tanh (30u) com s e k iguais. 42 onde x é a compressão ou alongamento da mola a partir do seu estado indeformado e k a constante de proporcionalidade. O sinal menos (–) indica que a força se opõe, ou resiste, à deformação. Ao tomarmos a constante elástica k, denominada constante elástica linear, a força da mola (força restauradora) que age no bloco, consequentemente, será linear e assume a forma: fmo = – k x (3.11) 3.2.3 Força externa O bloco sofre ainda influência de uma força externa de excitação harmônica que definimos pela equação: fe = fcos( t) (3.12) onde f e são, respectivamente, a amplitude e freqüência desta excitação harmônica. 3.3 Equação do movimento do bloco A equação que descreve o movimento do bloco obedece à segunda lei de Newton e tem a forma: emoa fffxm (3.13) onde m é a massa do bloco que oscila sobre a correia, x seu deslocamento em relação a posição de equilíbrio, fa a força de atrito, fmo a força da mola, e fe a força externa. 43 carga inercial I Substituindo (3.9), (3.11) e (3.12) em (3.13), obtemos a equação que descreve o movimento do bloco, definida por: )tcos(fkx)u30tanh(mgxm (3.14) 3.4 Equação do movimento do motor CC O motor CC (motor de corrente contínua) é um transdutor eletromecânico bidirecional que pode converter, em qualquer direção, a energia elétrica em mecânica e energia mecânica em elétrica. O esquema desse motor está descrito na figura 3-5 e na dinâmica idealizada para ele consideramos o campo magnético constante e denotamos por R a resistência do circuito e por L a indutância do rotor. A seguir, desenvolvemos as equações diferenciais que descrevem a dinâmica desse sistema eletromagnético. A relação entre o potencial elétrico e a força mecânica no movimento do rotor, através de um campo magnético, é definida pela lei de indução de Faraday e pela lei de Ampère. i + vel. angular torque L + motor CC R atrito Figura 3-5 Esquema genérico de um motor elétrico CC com carga inercial I. 44 A aplicação da voltagem V gera a corrente i no circuito do motor provocando, assim, um torque mecânico no rotor. Em contrapartida, é criada no circuito certa resistência R a essa corrente. O trabalho principal de todo motor com imã permanente é adquirir um torque constante pela manipulação da bobina e da corrente. O torque é diretamente proporcional a corrente e é definido pela equação: ik (3.15) onde k é geralmente conhecido como torque constante. A força contra-eletromotriz Vemf é diretamente proporcional à velocidade angular do motor elétrico e é definida pela equação: Vemf = kemf (3.16) onde kemf é geralmente conhecida como voltagem constante e depende de certas propriedades físicas do motor elétrico. Em motores de corrente contínua o torque e a voltagem constantes tem o mesmo valor (k = kemf = k) e por essa razão os subscritos emf e são usualmente desprezados. A voltagem total pode então ser escrita como: V = iR + L dt di + k (3.17) Admitindo que a corrente do circuito mude muito pouco, a derivada da corrente pode ser desprezada e, consequentemente, igualada a zero. Este procedimento é representado, conforme segue, pela equação: 45 dt di = 0 (3.18) Substituindo a equação (3.18) na equação (3.17) obtemos a equação: R kV i (3.19) Substituindo a equação (3.19) na equação (3.15) obtemos a equação: R k R kV 2 (3.20) Na velocidade zero o torque do motor é máximo. Este torque é denominado torque de estol e é definido pela equação: R kV 0 (3.21) Com a substituição da equação (3.21) na equação (3.20) obtemos a equação do torque do motor: R k2 0 (3.22) Substituindo na equação (3.22) o torque por zero, isto é fazendo = 0, obtemos a velocidade máxima do motor 0, definida pela equação: 20 k R 0 (3.23) Substituindo a equação (3.23) na equação (3.22) obtemos a equação geral do torque, definida por: 46 )(1 0 0 (3.24) Esta equação pode ser graficamente exibida e é comumente conhecida como diagrama do torque versus velocidade angular do motor CC. O gráfico da figura 3-6 mostra o torque gerado pelo motor elétrico em função da velocidade. Para voltagens constantes um segmento de reta une o torque de estol 0 à velocidade máxima 0. O declive dessa linha é definido pela razão 0 0 . O diagrama do torque de um motor CC versus sua velocidade angular é uma representação útil, porque podemos obter o desejado torque de carga nesse diagrama e verificar se o motor está habilitado em produzir esse torque. Obtemos a equação da parte mecânica do motor usando a segunda lei de Newton que diz que a carga inercial I vezes a derivada da velocidade angular é igual à soma de todos os torques sobre o eixo do motor. O resultado é a equação: k dt d I a (3.25) onde ka é o coeficiente de atrito viscoso. 0 0 0 0 Figura 3-6 Diagrama do torque de um motor CC versus velocidade angular. 47 Substituindo a equação (3.24) em (3.25) obtemos: k)(1 dt d I a 0 0 (3.26) Substituindo as relações dt d e em (3.26) obtemos: k) 1 (1I a 0 0 (3.27) onde (3.27) é a equação diferencial que descreve o movimento de rotação do motor. 3.5 Equações que descrevem o sistema Em muitos problemas práticos de engenharia foram observados que a fonte de energia do sistema é influenciada pela resposta do sistema. Esta constatação mostra a necessidade da construção de modelos que considerem a interação das variáveis de estado da fonte de energia e das variáveis de estado do sistema mecânico (KONONENKO, 1969), (BALTHAZAR et al., 2002). O modelo utilizado neste trabalho trata de um sistema dinâmico não-ideal com atrito seco e excitação externa onde interagem o movimento de oscilação adere- desliza do bloco e o movimento de rotação do motor elétrico. Quando atrelamos o motor CC à roda que movimenta a correia, surge um torque no eixo do motor CC devido à força de atrito existente entre o bloco e a correia. Essa interação da força de atrito e do movimento rotacional do motor dá origem à máquina não-ideal, objeto de estudo deste trabalho. No modelo aqui proposto, o movimento do bloco interfere no movimento rotacional do motor elétrico o qual responde interferindo no movimento de oscilação do bloco. 48 A ação da força de atrito no motor, devido ao torque, é definida pela equação: )utanh(mgµrfam (3.28) Sob ação da força fam no motor elétrico, obtemos a equação diferencial que descreve o movimento rotacional do motor e é definida por: u)µmgtanh(30rk) 1 (1I a 0 0 (3.29) As equações diferenciais (3.27) e (3.29) descrevem o modelo dinâmico proposto que é definido pelo sistema: )u(tanh I mgµr I k II t)(cos m f x m k )utanh(gµx a 0 00 (3.30) onde a primeira equação diferencial do sistema descreve o movimento do bloco sobre a correia e a segunda o movimento rotacional do motor. Na primeira equação diferencial do sistema, efaf )(cos m f x m k )u(tanhgx mf t (3.31) fe é a força externa que excita harmonicamente o bloco, fm a força da mola que age no bloco e fa a força de atrito dependente da velocidade relativa de deslizamento. A diferença da velocidade da correia e da velocidade do bloco é denominada velocidade relativa de deslizamento e é definida pela equação: 49 xru (3.32) A força de atrito que usamos neste trabalho é definida por: )utanh(mgµfa (3.33) que é uma aproximação suave da lei de atrito de Coulomb. Ao parâmetro , que determina uma proximidade entre a descontinuidade do modelo de Coulomb e essa aproximação suave (LIANG, 1995), atribuímos o valor = 30s/m. A força que movimenta harmonicamente o bloco é definida pela equação: )tcos( m f fe (3.34) onde é a freqüência da força externa que variaremos para analisar o comportamento do sistema dinâmico proposto. Na segunda equação do sistema, amf )utanh(mgµ I r m I k II )( a 0 00 (3.35) )m( é a função que define as características do motor elétrico CC e que admitimos aqui ser uma função linear e fam é a força que age no movimento rotacional do motor e que sofre interferência das oscilações adere-desliza que ocorrem no bloco. CAPÍTULO 4 Simulações numéricas do sistema 4.1 Introdução As séries temporais obtidas por integração numérica, traduzem o comportamento do sistema dinâmico em estudo. Com elas analisamos o movimento oscilatório do bloco sobre a correia, a interferência dessas oscilações no movimento rotacional do motor e o sistema dinâmico como um todo. Fazemos essa análise com os gráficos que dão a posição e velocidade do bloco em função do tempo, com as trajetórias do plano de fase, com o espectro de freqüência (análise de Fourier – FFT), com os coeficientes wavelets (análise wavelet) e com a taxa de divergência ou convergência das trajetórias do plano de fase que partem de pontos iniciais próximos (expoentes de Lyapunov). O plano de fase é um caso particular do espaço de fase que é definido como todo espaço matemático, com coordenadas ortogonais que representam as direções de cada variável, necessárias para especificar o estado instantâneo do sistema. Cada ponto do espaço de fase representa um estado possível para o sistema e por ele passa só uma 51 trajetória. No caso em que o espaço de fase é um plano de fase, o estado de uma partícula se movendo em uma dimensão é especificado pela sua posição x e velocidade v (BAKER, 1996). As trajetórias do plano de fase que usamos para análise deste trabalho são traçadas usando a posição e velocidade do bloco. A análise de Fourier e a análise wavelet, ambas com o objetivo de estudar as freqüências do sinal, serão feitas, respectivamente, pela transformada rápida de Fourier (FFT) e pela transformada wavelet. Com recursos do MATLAB 7.0 usamos o programa do Apêndice A para ilustrar os gráficos: posição do bloco em função do tempo, as trajetórias do plano de fase e o espectro de freqüências. Com o programa do apêndice B calculamos a transformada wavelet e com o do apêndice C, os expoentes de Lyapunov. 4.2 Atribuição de valores aos coeficientes do sistema O sistema de equações diferenciais que usamos para descrever o modelo proposto neste trabalho é, conforme as deduções feitas no Capítulo 3, definido por: )]x(r[tanh I mgrµ I k II t)(cos m f x m k )]x(r[tanhµ gx a 0 00 (4.1) e na equação do movimento do bloco o valor de m k define a freqüência natural de oscilação do bloco. Ao sistema descrito em (4.1) atribuímos os valores numéricos que seguem: 52 e obtemos, assim, o seguinte sistema: )]x0702,0(30tanh[137,68017,662,140 t)cos(3,47x5177,4)]x0702,0(30tanh[06,1x (4.3) que será integrado pelo método numérico de Runge-Kutta, residente no software MATLAB 7.0. Nos dados apresentados em (4.2) notamos que a freqüência natural do bloco é definida por Hz32,13 . 4.3 Integração numérica do sistema Transformando o sistema de equações diferenciais ordinárias de segunda ordem (4.3) em um sistema de equações diferenciais ordinárias de primeira ordem pela mudança das variáveis x , x , e por outras, definidas conforme segue, 1xx 2xx 3x 4x obtemos o sistema: 622,0µ o = 2,5 m2 kg/s2 2s/m8,9g I = 0,017778 m2 kg m0702,0r 0 = 78,3 rad/s (4.2) k = 1015kg/ s2 ka = 0,281 m2kg/s m = 5,72 kg u = (0,0702 x )m/s f = 19,86 kgm/s2 Hz32,13 m k m/s30 53 )]xx0702,0(30tanh[137,68x017,6140,62x xx )tcos(3,47x5177,4)]xx0702,0(30 tanh[06,1x xx 2444 43 1242 21 (4.4) Com o método numérico de Runge-Kutta de quarta ordem, instalado no software MATLAB 7.0, integramos o sistema (4.4) no intervalo de zero a 180 segundos, com intervalos de integração de 1/256 segundos e condições iniciais: x1 = 0, x2 = 0, x3 = 0 e x4 = 0. Desprezamos a parte transiente, tomamos os últimos 8192 pontos para obtermos quatro séries temporais que tratam, separadamente, da posição do bloco e sua velocidade, da posição rotacional do motor e sua velocidade. Analisamos o comportamento dinâmico do sistema por meio dessas séries usando os métodos matemáticos propostos neste trabalho. 4.4 Resultados das simulações Nas simulações, variamos a freqüência da força externa no intervalo de 11 Hz a 15 Hz, com passos de 0,1 Hz. Dentre essas simulações, escolhemos uma seqüência crescente de freqüências, que faz o bloco passar de um movimento periódico para um caótico e depois de um movimento caótico para um periódico. A seqüência das freqüências escolhidas é {12; 12,5; 13; 13,5; 14; 14,5; 15}. Exibimos a seguir, nas figuras de 4-1 até 4-14, em (a) os gráficos da posição do bloco em função do tempo, em (b) as trajetórias da velocidade do bloco em função da sua posição, em (c) os espectros de freqüências da posição do bloco, em (d) os diagramas dos expoentes de Lyapunov e em (e) os escalogramas da transformada wavelet. Calculamos os expoentes de Lyapunov (SIU, 1998) com t variando de zero até 1000 segundos e usamos a wavelet de Morlet no cálculo da transformada wavelet (COMPO, 1995). No escalograma, as cores quentes indicam os maiores coeficientes wavelets e as cores frias os menores coeficientes. 54 Na figura 4-1, para = 12 Hz, a característica do gráfico em (a) é de um movimento periódico e em (b), também, onde as trajetórias permanecem confinadas numa fina faixa, conforme ilustrado. Em (c), a freqüência fundamental e suas freqüências múltiplas, também caracterizam movimento periódico do bloco. Em (d) com 1 = 0,00283715, 2 = 1, 39047, 3 = 2, 88453 e 4 = 12, 4264 temos a confirmação do movimento periódico com um expoente de Lyapunov nulo e os demais negativos (SANDRI, 1996). Em (e), figura 4-2, o gráfico da transformada wavelet também caracteriza movimentos periódicos. (a) (c) (d) (b) Figura 4-1 Gráficos da simulação do sistema dinâmico proposto para = 12 Hz. (a) posição do bloco em função do tempo, (b) trajetórias do plano de fase, (c) espectro de freqüência e (d) expoentes de Lyapunov. Figura 4-2 Escalograma da transformada wavelet para = 12 Hz. (e) 55 Na figura 4-3, para = 12,5 Hz, a característica do gráfico em (a) é de um movimento periódico e em (b), também, onde as trajetórias permanecem confinadas numa fina faixa, conforme ilustrado. Em (c), a freqüência fundamental e suas freqüências múltiplas, também caracterizam movimento periódico do bloco. Em (d) com 1 = 0,00272898, 2 = 1, 57962, 3 = 1, 58342 e 4 = 13, 6154 temos a confirmação do movimento periódico com um expoente de Lyapunov nulo e os demais negativos (SANDRI, 1996). Em (e), na figura 4-4, o gráfico da transformada wavelet também caracteriza movimentos periódicos. (a) (c) Figura 4-3 Gráficos da simulação do sistema dinâmico proposto para = 12,5 Hz. (a) posição do bloco em função do tempo, (b) trajetórias do plano de fase, (c) espectro de freqüência e (d) expoentes de Lyapunov. (b) (d) Figura 4-4 Escalograma da transformada wavelet para = 12,5 Hz. (e) 56 Na figura 4-5, para =13 Hz, em (a), nada podemos afirmar quanto à periodicidade do movimento do bloco. Em (b) e (c) os gráficos dão indícios de movimento caótico do bloco. Em (d) com 1 = 0, 0047211, 2 = 0, 119831, 3 = 1, 3453 e 4 = 13, 3847 temos a confirmação do movimento periódico com um expoente de Lyapunov nulo e os demais negativos (SANDRI, 1996). Em (e), na figura 4-6, o gráfico da transformada wavelet também dá indícios de movimento periódico. Figura 4-6 Escalograma da transformada wavelet para = 13 Hz. (a) (b) (c) Figura 4-5 Gráficos da simulação do sistema dinâmico proposto para = 13 Hz. (a) posição do bloco em função do tempo, (b) trajetórias do plano de fase, (c) espectro de freqüência e (d) expoentes de Lyapunov. (d) (e) 57 Na figura 4-7, para =13,5 Hz, em (a) o movimento do bloco não tem características periódicas. Em (b) e (c) os gráficos dão indícios de movimento caótico do bloco. Em (d) com 1 = 0, 404333, 2 = 0, 00276072, 3 = 2, 1918 e 4 = 12, 50767 temos a confirmação do movimento caótico com um expoente de Lyapunov positivo, um nulo e os demais negativos (SANDRI, 1996). Em (e), na figura 4-8, o gráfico da transformada wavelet também dá indícios de movimento caótico. (a) (b) (c) Figura 4-7 Gráficos da simulação do sistema dinâmico proposto para = 13,5 Hz. (a) posição do bloco em função do tempo, (b) trajetórias do plano de fase, (c) espectro de freqüência e (d) expoentes de Lyapunov. (d) Figura 4-8 Escalograma da transformada wavelet para = 13,5 Hz. (e) 58 Na figura 4-9, para =14 Hz, em (a) o movimento do bloco não tem características periódicas. Em (b) e (c) o gráfico dá indícios de movimento caótico do bloco. Em (d) com 1 = 0, 339166, 2 = 0, 0011096, 3 = 1, 91544 e 4 = 12, 1101 temos a confirmação do movimento caótico com um expoente de Lyapunov positivo, um nulo e os demais negativos (SANDRI, 1996). Em (e), na figura 4-10, o gráfico da transformada wavelet também dá indícios de movimento caótico. (a) (b) (c) Figura 4-9 Gráficos da simulação do sistema dinâmico proposto para = 14 Hz. (a) posição do bloco em função do tempo, (b) trajetórias do plano de fase, (c) espectro de freqüência e (d) expoentes de Lyapunov. Figura 4-10 Escalograma da transformada wavelet para = 14 Hz. (d) (e) 59 Na figura 4-11, para =14,5 Hz, a característica do gráfico em (a) é de movimento periódico no bloco e em (b) também, onde as trajetórias permanecem confinadas em três faixas, conforme ilustrado. Em (c) o gráfico dá indícios de movimento caótico do bloco. Em (d) com 1 = 0, 0025521, 2 = 1, 8951, 3 = 5, 9865 e 4 = 8, 3223 temos a confirmação do movimento periódico com um expoente de Lyapunov nulo e os demais negativos (SANDRI, 1996). Em (e), na figura 4-12, o gráfico da transformada wavelet também dá indícios de movimento periódico. (b)(a) (c) Figura 4-11 Gráficos da simulação do sistema dinâmico proposto para = 14,5 Hz. (a) posição do bloco em função do tempo, (b) trajetórias do plano de fase, (c) espectro de freqüência e (d) expoentes de Lyapunov. Figura 4-12 Escalograma da transformada wavelet para = 14,5 Hz. (d) (e) 60 Na figura 4-13, para = 15 Hz., a característica do gráfico em (a) é de um movimento periódico do bloco e em (b), também, onde as trajetórias permanecem confinadas numa fina faixa, conforme ilustrado. Em (c), a freqüência fundamental e suas freqüências múltiplas, também caracterizam movimento periódico do bloco. Em (d) com 1 = 0,00323494, 2 = 1, 91622, 3 = 2, 63337 e 4 = 12,0558 temos a confirmação do movimento periódico com um expoente nulo e os demais negativos (SANDRI, 1996). Em (e), na figura 4-14, o gráfico da transformada wavelet também caracteriza movimentos periódicos. (a) (b) (c) Figura 4-14 Escalograma da transformada wavelet para = 15 Hz. Figura 4-13 Gráficos da simulação do sistema dinâmico proposto para = 15 Hz. (a) posição do bloco em função do tempo, (b) trajetórias do plano de fase, (c) espectro de freqüência e (d) expoentes de Lyapunov. (d) (e) 61 4.5 Seção de Poincaré Uma maneira de observar o comportamento de sistemas dinâmicos é por meio da seção de Poincaré, artifício inventado por Henri Poincaré como um meio de simplificação do diagrama do espaço de fase desses sistemas dinâmicos. Ele é construído pela visão do diagrama do espaço de fase de forma estroboscópica, de tal modo que o movimento é observado periodicamente. Para o sistema dinâmico proposto o período estroboscópico é o período da força externa que age no bloco (BAKER, 1996). O método de Poincaré consiste em cortar ou secionar o atrator espiral em intervalos regulares de tempo, observando nessas seções o plano definido pela posição e velocidade do bloco. Se esta seção é feita em intervalos correspondentes ao movimento da força externa então, os desenhos estroboscópicos são todos mostrados em forma de ponto. A posição e velocidade do bloco são as mesmas quando a seção é feita num período T = 2 / ( = 2 f f = /2 T = 2 / ) (BAKER, 1996). Para um comportamento periódico do sistema dinâmico, a seção de Poincaré apresenta um conjunto de pontos coincidentes, que se repetem e estão concentrados em pequenas regiões do mapa enquanto que para um comportamento caótico os pontos não são coincidentes, não se repetem e são distribuídos no mapa de forma dispersa. Com, aproximadamente, 1250 pontos resultantes do movimento do motor de 300 a 900 segundos em passos de 0,001s, gerados pelo programa do Apêndice E, ilustramos a passagem do movimento periódico para o caótico em valores vizinhos a 13 Hz (figura 4-15), alguns movimentos caóticos em 13,5 Hz, 13,7 Hz e 14 Hz (figura 4-16) e a passagem do movimento caótico para o periódico em valores vizinhos à 14,5 Hz (figura 4-17). 62 Figura 4-15 Passagem do movimento periódico para o caótico. Em (a) os pontos se concentram em duas regiões da seção de Poincaré; em (b) os pontos começam a se dispersar e em (c) estão totalmente dispersos tomando uma grande região da seção de Poincaré. (a) (b) (c) 63 Figura 4-16 Movimento caótico do sistema dinâmico proposto. Em (a), (b) e (c) os pontos estão totalmente dispersos tomando uma grande região da seção de Poincaré. (c) (a) (b) 64 Figura 4-17 Passagem do movimento caótico para o periódico. Em (a) os pontos estão totalmente dispersos tomando uma grande região da seção de Poincaré, em (b) os pontos começam a se concentrar em duas regiões e em (c) estão concentrados em duas regiões da seção de Poincaré. (a) (b) (c) 65 4.6 Bifurcação Uma bifurcação ocorre quando uma pequena mudança feita no valor de um parâmetro (parâmetro de bifurcação) de um sistema dinâmico causa uma súbita mudança qualitativa no comportamento dinâmico desse sistema (STROGATZ, 1994). O comportamento do sistema dinâmico bloco-correia-motor pode ser globalmente observado quando verificamos a posição do bloco no início de cada ciclo para um determinado intervalo de freqüências da força externa. Para isto integramos o sistema numericamente calculando para várias freqüências desse intervalo a posição do bloco. O gráfico gerado pela posição do bloco no começo de cada ciclo versus freqüência é denominado diagrama de bifurcação e está representado na figura 4-18 para as freqüências da força externas que variam de 11,8 Hz até 15 Hz em passos de 0,1 Hz. Na integração do sistema consideramos o tempo variando de 300 a 900 segundos, em passos de 0,001s. Quando fazemos variar a freqüência da força externa é possível notar as posições que o bloco ocupa no começo de cada ciclo, examinando o gráfico da figura 4-18. Figura 4-18 Diagrama de bifurcação. Os valores da posição do bloco no começo de cada ciclo dependem da freqüência da força externa que varia de 11,8 Hz até 15 Hz. 66 Notamos no diagrama de bifurcação (figura 4-18) que, quando varia de 11,8 Hz até as proximidades de 13 Hz e nas proximidades de 14,5 Hz até 15 Hz, o bloco ocupa duas posições no início de cada ciclo, caracterizando um comportamento periódico do sistema dinâmico. Para os valores restantes de isto é, entre 13 Hz e 14,5 Hz, o bloco parece ocupar um número maior de posições no início de cada ciclo, dando indícios do comportamento caótico do sistema dinâmico proposto. Nas regiões limítrofes, onde o sistema passa de um comportamento periódico para caótico e vice-versa, o mapa de Poincaré dá maiores detalhes dessa região. 4.7 Análise da interação do movimento do bloco e da fonte de energia Devido ao atrito, as oscilações que ocorrem no bloco sobre a correia podem influenciar no movimento rotacional do motor elétrico como, também, o movimento rotacional do motor elétrico pode influenciar nas oscilações do bloco. Com o objetivo de analisar o comportamento do sistema dinâmico proposto, em situações onde a interferência do movimento de oscilação do bloco no movimento rotacional do motor elétrico e as respostas do motor elétrico a essa interferência são visíveis, diminuímos a potência do motor elétrico alterando seu torque de estol e conservamos todos os outros parâmetros conforme foram definidos em (4.2). Escolhemos duas freqüências da força externa sendo uma delas igual a 14 Hz, onde o comportamento do sistema é caótico e a outra igual a 12,6 Hz onde o comportamento do sistema é periódico. Os comportamentos do sistema nessas freqüências são conhecidos, pois já foram estudados neste capítulo em 4.4. Integramos o sistema, numericamente, nas freqüências da força externa iguais a 14 Hz e 12,6 Hz e conservamos todos os outros parâmetros conforme já foram definidos em (4.2), com exceção dos torques de estol do motor elétrico, aos quais foram dados valores decrescentes com a finalidade de diminuir a potência do motor elétrico. Os torques escolhidos estão definidos a seguir: 67 1) 0 = 18 m2 kg/s2; 2) 0 = 12 m2 kg/s2; (4.5) 3) 0 = 2,2 m2 kg/s2. Analisamos os resultados dessa integração numérica e ilustramos da figura 4-19 até a figura-4-24, em (a) o gráfico das trajetórias do plano de fase (posição e velocidade do bloco), em (b) o espectro de freqüência em (c) o diagrama de potência do motor elétrico e em (d) o diagrama da velocidade angular do motor elétrico. 4.7.1 Freqüência da força externa igual a 14 Hz Analisamos neste item o comportamento dinâmico do sistema bloco- correia-motor na freqüência da força externa igual a 14 Hz e nos torques que estão propostos em (4.5). a.1) 0 = 18 m2 kg/s2 Neste torque a potência máxima do motor elétrico é igual a 1409,4 watts e sua velocidade angular máxima é igual a 78,3 rad/s. Nessas condições, o motor elétrico é capaz de movimentar a correia sem sofrer grandes influências do atrito. Isto é verificado quando observarmos os gráficos da figura 4-19 onde em (c) a potência do motor elétrico de, aproximadamente, 352 watts e em (d) sua velocidade angular de, aproximadamente, 40 rad/s se mantêm inalteradas. Isto caracteriza que o movimento do bloco sobre a correia não influencia no movimento rotacional do motor elétrico como, também, o motor elétrico não influencia no movimento de oscilação do bloco. Os gráficos da figura 4-19 que ilustram, em (a) as trajetórias do plano de fase confinadas numa fina faixa, em (b) apenas uma freqüência e em (d) velocidade 68 angular constante do motor elétrico revelam, que o sistema bloco-correia-motor tem comportamento periódico. a.2) 0 = 12 m2 kg/s2 Neste torque a velocidade angular máxima do motor elétrico continua sendo igual a 78, 3 rad/s e sua potência máxima diminui para 939,6 watts. Nessas condições, o motor elétrico já começa sofrer alguma influência do atrito. Isto é verificado ao observarmos os gráficos da figura 4-20 onde em (c) a potência do motor elétrico oscila entre, aproximadamente, 200 watts e 230 watts, em (d) sua velocidade angular também oscila entre, aproximadamente, 24 rad/s e 33 rad/s. Isto é uma característica da interferência do movimento de oscilação do bloco no motor elétrico e das respostas do motor elétrico a essa interferência. Figura 4-19 Gráficos da simulação do sistema dinâmico proposto para = 14 Hz e 0 = 18m2 kg/s2. (a) plano de fase, (b) espectro de freqüência, (c) potência do motor elétrico e (d) velocidade angular do motor elétrico. (a) (b) (d) (c) 69 Os gráficos da figura 4-20 que ilustram, em (a) as trajetórias do plano de fase confinadas numa fina faixa, em (b) a freqüência fundamental e suas freqüências múltiplas e em (d) a velocidade angular do motor elétrico oscilando, periodicamente, revelam que o sistema bloco-correia-motor tem comportamento periódico. a.3) 0 = 2,2 m2 kg/s2 Diminuindo ainda mais o torque de estol do motor elétrico para 0 = 2,2 m2 kg/s2 e conservando a mesma velocidade angular máxima em 78, 3 rad/s, temos um motor elétrico com potência máxima igual a 172,26 watts. Com essas características, o motor elétrico já começa a sofrer grande influência do atrito. Notamos isto quando observarmos os gráficos da figura 4-21 onde em (c) a potência do motor elétrico oscila entre, aproximadamente, 2 watts e 27 watts e em (d) sua velocidade angular também oscila entre, aproximadamente, 1 rad/s e 15 rad/s. A interferência do bloco no movimento Figura 4-20 Gráficos da simulação do sistema dinâmico proposto para = 14 Hz e 0 = 12m2 kg/s2. (a) plano de fase, (b) espectro de freqüência, (c) potência do motor elétrico e (d) velocidade angular do motor elétrico. (a) (c) (d) (b) 70 rotacional do motor elétrico e a resposta que dá o motor elétrico a esse tipo de interferência estão presentes neste caso, pois as oscilações de potência e velocidade do motor elétrico são claramente observadas na figura 4-21 (c) e (d). Os gráficos da figura 4-21 que ilustram em (a) as trajetórias do plano de fase com características caóticas, em (b) um número grande de freqüências diferentes, e em (d) velocidade angular oscilante do motor elétrico com características caóticas, revelam que o sistema bloco-correia-motor tem comportamento caótico. 4.7.2 Freqüência da força externa igual a 12,6 Hz Neste item analisamos o sistema bloco-correia-motor com a freqüência da força externa igual a 12,6 Hz, repetindo os mesmos procedimentos usados em 4.7.1. (c) (a) (b) (c) Figura 4-21 Gráficos da simulação do sistema dinâmico proposto para = 14 Hz e 0 = 2,2m2 kg/s2. (a) plano de fase, (b) espectro de freqüência, (c) potência do motor elétrico e (d) velocidade angular do motor elétrico. 71 b.1) 0 = 18 m2 kg/s2 Neste torque e na velocidade angular máxima igual a 78, 3 rad/s o motor elétrico atinge potência máxima igual ao caso a.1. Essa potência dá condições ao motor elétrico de movimentar a correia sem sofrer grandes influências do atrito. Notamos isto quando observamos os gráficos da figura 4-22 onde em (c) a potência do motor elétrico igual a, aproximadamente, 352 watts e em (d) sua velocidade angular, aproximadamente, igual a 40 rad/s se mantêm inalteradas. Isto mostra que o movimento de oscilação do bloco não interfere no movimento rotacional do motor elétrico como também o motor elétrico não interfere no movimento de oscilação do bloco. Os gráficos da figura 4-22 que ilustram, em (a) as trajetórias do plano de fase confinadas numa faixa, em (b) apenas uma freqüência e em (d) velocidade angular constante do motor e