Investigando a Dinâmica do Rotor Duplo Pulsado: um laboratório dinâmico para sistemas caóticos discretos com espaço de fase 4D

Investigating the Dynamics of the Kicked Double Rotor: a dynamical laboratory for discrete chaotic systems with 4D phase space

Resumos

O estudo de sistemas dinâmicos não lineares possibilita compreender, aplicar e prever fenômenos em diversas áreas do conhecimento científico, indo desde a meteorologia e a dinâmica de corpos celestes até o crescimento populacional e as sinapses cerebrais e muitas outras. Tendo em vista sua grande aplicabilidade, é importante conhecer ferramentas básicas para análise de sistemas dinâmicos. No entanto, embora muitos sistemas físicos e de engenharia sejam descritos por modelos matemáticos multidimensionais, grande parte da literatura introdutória sobre o assunto se limita ao estudo de sistemas com espaço de fase uni ou bidimensional. Portanto, esta contribuição objetiva apresentar, de forma didática e empregando uma abordagem computacional, ferramentas básicas para análise de sistemas com espaço de fase com dimensão mais elevada, explorando a obtenção e a interpretação de diagramas de bifurcação, cálculo de pontos de equilíbrio, análise de estabilidade e cálculo de expoentes de Lyapunov. Para isso utilizaremos um sistema de tempo discreto com espaço de fase 4-dimensional conhecido como Rotor Duplo Pulsado como laboratório dinâmico.

Palavras-chave:
Mapas caóticos; Análise de estabilidade; Expoentes de Lyapunov

The study of nonlinear dynamical systems makes it possible to understand, apply and predict phenomena in several areas of scientific knowledge, ranging from meteorology and the dynamics of celestial bodies to population growth and brain synapses, anad many more. In view of its wide applicability, it is important to know basic tools used for analyzing dynamical systems. However, although many physical and engineering systems are described by multidimensional mathematical models, much of the introductory literature on the subject is limited to the study of systems with one or two-dimensional phase space. Therefore, this contribution aims to present, in a didactic way and employing a computational approach, basic tools to analyse systems with higher dimensional phase space, exploring how to obtain and interpret bifurcation diagrams, compute equilibrium points and analyse their stability and compute Lyapunov exponents. To do so, we will use a discrete time system with 4-dimensional phase space known as Pulsed Double Rotor as a dynamical laboratory.

Keywords
Chaotic maps; Stability analysis; Lyapunov exponents.

1. Introdução

O conceito de Sistema Dinâmico pode ser entendido analisando a etimologia de cada palavra separadamente. O termo “sistema” remonta a um conjunto de corpos ou objetos que possuem uma conexão devido a uma certa interdependência de modo que exista uma relação de causa e efeito entre os mesmos, isto é, os acontecimentos, ou fenômenos, que envolvem um elemento do conjunto geram consequências nos demais elementos. Já o termo “dinâmico” foi introduzido por G. W. Leibniz em seus estudos e se refere a dependência temporal das grandezas que compõem o sistema [¹]. Dessa forma, a definição de Sistemas Dinâmicos refere-se a um conjunto de possíveis estados de um sistema que são determinados por uma regra matemática que define o estado atual em função do estado passado [²].

Exemplos clássicos de Sistemas Dinâmicos incluem modelos para o movimento de corpos celestes interagindo por meio da força gravitacional que depende da massa dos corpos envolvidos e da posição relativa entre eles que, por sua vez, depende do tempo. Mais recentemente, a abordagem de Sistemas Dinâmicos tem sido utilizada para investigar a atividade neural no cérebro humano, as variações no sistema econômico global e a dinâmica das populações [³].

A possibilidade de conhecer o estado presente de um sistema apenas conhecendo suas condição inicial e as equações que o regem, isto é, aplicando uma regra determinística, motivou os estudos sobre a Teoria dos Sistemas Dinâmicos durante anos. Com o passar do tempo, avanços importantes foram obtidos como, por exemplo, o surgimento da teoria do Caos e da teoria de controle moderno, entre muitas outras aplicações advindas dos estudos de muitos pesquisadores que se dedicaram a essa área, dentre os quais se destacam cientistas, matemáticos e engenheiros de renome, indo desde o próprio Isaac Newton até Henri Poincaré e Edward Lorenz [¹, ²].

Com isso, torna-se cada vez mais importante e relevante conhecer conceitos e ferramentas para estudar Sistemas Dinâmicos de modo a identificar e interpretar a resposta de modelos matemáticos de uma vasta gama de sistemas físicos e de engenharia cujas soluções podem ser diversas, incluindo pontos fixos, órbitas periódicas e até mesmo trajetórias caóticas.

Neste artigo são apresentadas ferramentas para o estudo de sistemas dinâmicos e seus algoritmos. Especificamente, discutem-se métodos para se construir diagramas de bifurcação, cálculo de pontos fixos, análise de estabilidade e cálculo de expoentes de Lyapunov para caracterização de dinâmica caótica. Esses métodos são, então, aplicados a um sistema dinâmico não linear com 2 graus de liberdade modelado por um sistema de quatro equações de tempo discreto, conhecido como Rotor Duplo Pulsado (RDP).

Uma vez que a maior parte da bibliografia disponível se limita a tratar sistemas com espaço de fase uni ou bidimensionais, o objetivo principal desta contribuição é expor, de forma didática, o uso de ferramentas aplicadas ao estudo do RDP. Assim, ao apresentar e aplicar algumas das principais ferramentas de análise de Sistemas Dinâmicos a um sistema com espaço de estados 4-dimensional, espera-se contribuir para o ensino e para estudo de outros modelos de sistemas dinâmicos com mais graus de liberdade e auxiliar futuros estudos na área. Adicionalmente, objetiva-se apresentar novos resultados que aprofundam a compreensão sobre a dinâmica do Rotor Duplo Pulsado, visto que este sistema possui uma muito dinâmica rica e ainda é pouco estudado.

O artigo está organizado da seguinte forma. Na Seção ² apresenta-se o sistema conhecido como Rotor Duplo Pulsado, incluindo um breve histórico, o modelo físico e as equações que o regem. Na seção ³ apresentam-se as principais ferramentas de análise de sistemas dinâmicos utilizados neste trabalho, bem como sua aplicação ao caso específico do Rotor Duplo Pulsado como forma de exemplo didático. Por fim, na seção ⁴, são feitas as considerações finais.

2. O Rotor Duplo Pulsado

O sistema dinâmico não linear rotativo conhecido como Rotor Duplo Pulsado (RDP), introduzido originalmente por Romeiras et al. [³], é um modelo matemático de tempo discreto que descreve a evolução temporal de um sistema mecânico composto por duas hastes finas de massa desprezível que giram em torno de um pivô fixo. Uma das hastes está sujeita a um forçante impulsivo aplicado periodicamente como mostrado na Fig. 1.

Figura 1
Representação esquemática do Rotor Duplo Pulsado.

Mais especificamente, a primeira haste tem comprimento l₁ e está anexada em uma das suas extremidades ao pivô fixo no ponto P₁ e na outra extremidade ao ponto P₂ que se localiza no ponto médio da segunda haste que tem comprimento 2l₂, conforme o diagrama esquemático da Fig. 1. As duas hastes são livres para rotacionar horizontalmente e suas posições angulares são dadas pelos ângulos θ₁(t) e θ₂(t) que representam, respectivamente, a orientação angular da primeira e da segunda haste no instante de tempo t. Uma massa m₁ é fixada no ponto de contato entre as hastes, P₂, e duas massas m₂/2 são fixadas em cada extremidade da segunda haste nos pontos P₃ e P₄. Assim, a energia cinética do sistema é dada por

(1) $$K(t)=\frac{1}{2}(I_1{\stackrel{.}{\theta }}_1^2+I_2{\stackrel{.}{\theta }}_2^2),$$

onde ${\stackrel{.}{\theta }}_1$ e ${\stackrel{.}{\theta }}_2$ são as velocidades angulares e I₁ e I₂ são os momentos de inércia de cada haste, respectivamente, com $I_1=(m_1+m_2)l_1^2$ e $I_2=m_2{l}_2^2$.

No ponto P₃ são aplicados impulsos constantes de intensidade ρ sempre na mesma direção nos instantes t = T,2T,3T,etc. Isso equivale ao rotor estar na presença de um campo homogêneo ligado periodicamente de forma instantânea e desligado em seguida. Assim, a energia potencial do sistema pode ser escrita como

(2) $$U(t)=(l_1\text{cos}{\theta }_1+l_2\text{cos}{\theta }_2)f(t),$$

onde f(t) = ∑_kρδ(t−kT), k = 1,2,3,…, representa matematicamente os impulsos aplicados periodicamente, com δ(t−kT) denotando a função Delta de Dirac.

Na ausência de dissipação a Lagraneana do sistema é

(3) $$\begin{array}{l}L=K-U\\ L=\frac{1}{2}(I_1{\dot{\theta }}_1^2+I_2{\dot{\theta }}_2^2)-(l_{1}cos{\theta }_1+l_{2}cos{\theta }_2)f(t),\end{array}$$

e as equações de movimento são dadas por

$$\frac{d}{dt}\left(\frac{\partial L}{\partial \stackrel{.}{{\theta }_i}}\right)-\frac{\partial L}{\partial {\theta }_i}=0,i=1,2,$$

ou, explicitamente,

(4) $$I_1{\ddot{\theta }}_1-l_{1}f(t)\text{sen}{\theta }_1=0$$ $$I_2{\ddot{\theta }}_2-l_{2}f(t)\text{sen}{\theta }_2=0.$$

Pode-se introduzir uma dissipação de energia no sistema considerando o atrito nos pontos de pivô P₁ e P₂. Nesse caso, o atrito é modelado como uma dissipação viscosa com a função de dissipação de Rayleigh dada por

(5) $$ℱ=\frac{1}{2}{\nu }_1{I}_1{\stackrel{.}{\theta }}_1^2+\frac{1}{2}{\nu }_2{I}_2{({\stackrel{.}{\theta }}_2-{\stackrel{.}{\theta }}_1)}^2,$$

onde ν₁I₁ e ν₂I₂ representam coeficientes de atrito viscoso em P₁ e P₂, respectivamente, com ν₁ e ν₂ dados por constantes reais positivas. Agora, as equações de movimento são dadas por

$$\frac{d}{dt}\left(\frac{\partial L}{\partial \stackrel{.}{{\theta }_i}}\right)-\frac{\partial L}{\partial {\theta }_i}=-\frac{\partial ℱ}{\partial {\stackrel{.}{\theta }}_i},i=1,2,$$

ou, explicitamente,

(6) $$\begin{array}{r}{I}_1{\ddot{\theta }}_1+({\nu }_1{I}_1+{\nu }_2{I}_2){\dot{\theta }}_1-{\nu }_2{I}_2{\dot{\theta }}_2-l_{1}f(t)\text{sen}{\theta }_1=0\\ I_2{\ddot{\theta }}_2-{\nu }_2{I}_2{\dot{\theta }}_1+{\nu }_2{I}_2{\dot{\theta }}_2-l_{2}f(t)\text{sen}{\theta }_2=0.\end{array}$$

Uma vez que f(t) = 0 para t≠kT,k = 1,2,…, as equações (6) são lineares para todo t tal que (k−1)T < t < kT e podem ser resolvidas pelos métodos usuais para resolução de equações diferenciais ordinárias com coeficientes constantes. No entanto, neste trabalho estamos interessados nos efeitos não lineares introduzidos pela aplicação periódica dos impulsos.

Notamos que, em t = kT, o forçante externo altera de forma instantânea a velocidade angular das hastes, mas não suas posições, de forma que a velocidade angular de cada haste é descontínua em kT. Matematicamente, podemos escrever

(7) $$\begin{array}{l}{\theta }_i(k{T}^{+})={\theta }_i(k{T}^{-})={\theta }_i(kT)\\ {\dot{\theta }}_i(k{T}^{+})-{\dot{\theta }}_i(k{T}^{-})=\frac{\rho l_i}{I_i}sen({\theta }_i(kT)),\end{array}$$

onde θ_i(kT⁻) e $\stackrel{.}{{\theta }_i}(k{T}^{-})$, com i = 1,2, representam as posições e as velocidades imediatamente antes da aplicação do forçante, enquanto que θ_i(kT⁺) e $\stackrel{.}{{\theta }_i}(k{T}^{+})$ representam as posições e as velocidades imediatamente após a aplicação do forçante.

Essas considerações permitem a modelagem do sistema em termos de 4 equações de diferenças que descrevem a evolução temporal do modelo mecânico em intervalos de tempo discreto imediatamente após a aplicação do forçante periódico [³, ⁴]. Assim, obtém-se o sistema de equações (8) que envolve as posições angulares θ₁ e θ₂ e as velocidades angulares ${\stackrel{.}{\theta }}_1$ e ${\stackrel{.}{\theta }}_2$ de cada haste:

(8) $$\begin{array}{l}{\theta }_1^{(k+1)}={\theta }_1^{(k)}+M_{11}{\dot{\theta }}_1^{(k)}+M_{12}{\dot{\theta }}_2^{(k)}\\ {\theta }_2^{(k+1)}={\theta }_2^{(k)}+M_{21}{\dot{\theta }}_1^{(k)}+M_{22}{\dot{\theta }}_2^{(k)}\\ {\dot{\theta }}_1^{(k+1)}=\frac{\rho l_1}{I_1}sen{\theta }_1^{(k+1)}+L_{11}{\dot{\theta }}_1^{(k)}+L_{12}{\dot{\theta }}_2^{(k)}\\ {\dot{\theta }}_2^{(k+1)}=\frac{\rho l_2}{I_2}sen{\theta }_2^{(k+1)}+L_{21}{\dot{\theta }}_1^{(k)}+L_{22}{\dot{\theta }}_2^{(k)}.\end{array}$$

As variáveis θ₁ e θ₂ são moduladas no intervalo [0,2π] e M_ij e L_ij são termos de matrizes constantes definidas por:

(9) $$L={\displaystyle \sum _{i=1}^2}{W}_i{e}^{{\lambda }_{i}T},$$ $$M={\displaystyle \sum _{i=1}^2}{W}_i{\displaystyle \frac{e^{{\lambda }_{i}T}-1}{{\lambda }_i}},$$ $$W_1=\left(\begin{array}{cc}\hfill a\hfill & \hfill b\hfill \\ \hfill b\hfill & \hfill d\hfill \end{array}\right),W_2=\left(\begin{array}{cc}\hfill d\hfill & \hfill -b\hfill \\ \hfill -b\hfill & \hfill a\hfill \end{array}\right),$$ $$a={\displaystyle \frac{1}{2}}\left(1+{\displaystyle \frac{{\nu }_1}{\mathrm{\Delta }}}\right),b=-{\displaystyle \frac{{\nu }_2}{\mathrm{\Delta }}},$$ $$d={\displaystyle \frac{1}{2}}\left(1-{\displaystyle \frac{{\nu }_1}{\mathrm{\Delta }}}\right),$$ $${\lambda }_{1,2}=-{\displaystyle \frac{1}{2}}({\nu }_1+2{\nu }_2\pm \mathrm{\Delta }),$$ $$\mathrm{\Delta }=\sqrt{{\nu }_1^2+4{\nu }_2^2}.$$

Há diversos parâmetros envolvidos na modelagem do RDP, o que torna possível estudar diversas condições de operação do modelo físico. Os valores dos parâmetros utilizados neste trabalho serão definidos na próxima seção.

Vale destacar que, o sistema de equações (8) é uma extensão quadridimensional do rotor simples, modelado pelo Mapa Padrão, um sistema bidimensional amplamente estudado e descrito na literatura. Mais detalhes sobre a obtenção do mapa do RDP podem ser encontrados nas Refs. [³, ⁵]. Adicionalmente, a Ref. [⁶] apresenta uma dedução semelhante para o rotor simples.

3. Ferramentas básicas para análise de sistemas dinâmicos não lineares aplicadas ao RDP

Nesta seção são apresentadas algumas ferramentas fundamentais para a investigação de sistemas dinâmicos não lineares. As ferramentas descritas são, então, aplicadas ao RDP de modo a exemplificar a aplicação das mesmas e facilitar o entendimento.

Notadamente, o estudo de sistemas dinâmicos envolve, em grande parte, a utilização de ferramentas computacionais que podem ser implementadas em diversas linguagens de programação como Fortran, C++, Python ou MatLab. A seguir, junto com a descrição de cada ferramenta, são apresentados algoritmos básicos para implementação das mesmas. Neste trabalho as implementações foram realizadas em Python e/ou Fortran.

3.1. Diagramas de Bifurcação

Diagramas de bifurcação são ferramentas extremamente úteis na análise de sistemas dinâmicos pois descrevem de maneira gráfica a evolução de uma determinada variável dinâmica em função de um parâmetro p do sistema. Enquanto a variável é expressa no eixo vertical, p evolui no eixo horizontal e todos os demais parâmetros do sistema são mantidos constantes [²].

A obtenção destes diagramas deve levar em conta as seguintes etapas gerais:

1. Escolher os parâmetros e as variáveis dinâmicas a analisar.

2. Definir o intervalo de variação do parâmetro p e determinar o valor do incremento que será utilizado.

3. Definir quantas condições iniciais serão utilizadas para investigar cada valor p do parâmetro.

4. Realizar a construção numérica do diagrama propriamente dita: para cada p, evolui-se o sistema a partir do conjunto de condições iniciais por um certo número de iterações predefinido, desprezando o transiente adequado, e guardam-se os pontos finais.

5. Plotar em um gráfico os pontos salvos em arquivo.

Ao definir o valor inicial, o valor final e o incremento de p, deve-se levar em conta quantos valores do parâmetro serão analisados de forma a cobrir adequadamente o intervalo objetivando um gráfico sem “buracos” que acarretariam uma representação visual pobre do sistema.

Suponha que se deseja construir um diagrama que descreva a evolução de um sistema considerando p ∈ [0,3]. Um incremento Δ = 0.01 implica que 300 valores de p serão analisados, enquanto que um incremento Δ = 0.03 proporciona uma resolução gráfica três vezes menor. Vale ressaltar que o custo computacional e de armazenamento do primeiro caso também é maior, de forma que o custo benefício de utilizar um dado incremento deve ser avaliado caso a caso.

Ao construir um diagrama é possível evoluir uma ou mais condições iniciais para cada valor do parâmetro investigado. Em ambos os casos pode-se seguir a evolução de um conjunto de condições iniciais, isto é, utilizar o último ponto da evolução temporal de uma dada condição inicial para um determinado p como a nova condição inicial para p + Δ, ou pode-se ainda reiniciar as condições iniciais ao conjunto original a cada novo incremento do parâmetro variável. Vale notar que, diferentemente do que ocorre em casos unidimensionais clássicos, como é o caso do Mapa Logístico [²], em dimensão mais alta é comum a coexistência de diversos estados assintóticos para um mesmo p, sendo que cada um deles têm a sua própria bacia de atração. Assim, são necessárias condições iniciais nas diferentes bacias de atração para visualizar de forma adequada a evolução do sistema. A escolha do número de condições iniciais deve ser criteriosa pois quanto maior o número de condições iniciais escolhidas, maior será o tempo computacional consumido. No Apêndice, o Algoritmo 1 escrito em pseudocódigo apresenta o procedimento básico para obter um diagrama de bifurcação utilizando várias condições iniciais que são reinicializadas a cada novo incremento do parâmetro variável.

A definição do número de iterações também depende do sistema e de sua complexidade. Para sistemas unidimensionais, como o Mapa Logístico [²], 100 iterações são suficientes, enquanto que para sistemas bidimensionais como o Mapa Padrão pode ser necessário utilizar valores da ordem de 10³ iterações, especialmente logo após o acontecimento de bifurcações globais. A resposta transiente do sistema deve ser desprezada, salvando a variável dinâmica que deseja-se estudar e o valor do parâmetro para uma ou mais iterações após alcançar o regime permanente. Novamente, o número de pontos salvos após o transiente depende de alguns fatores. Por exemplo, suponha que para um dado parâmetro p um atrator caótico coexista com algumas soluções periódicas de período elevado e que, entre todas as condições iniciais escolhidas apenas uma está na bacia de atração do atrator caótico. Nesse casso, se um número pequeno de pontos forem salvos após desprezar o transiente, pode-se interpretar incorretamente os pontos obtidos como sendo parte de uma solução periódica, por exemplo, deixando, assim, de identificar a solução caótica.

Construção de diagramas de bifurcação do RDP. A seguir, exemplifica-se a construção de diagramas de bifurcação utilizando o RDP com os seguintes parâmetros fixados: ν₁ = ν₂ = 1, T = 1, I₁ = I₂ = 1, l₂ = 1, $l_1=\frac{1}{\sqrt{2}}$ e m₁ = m₂ = 1. A escolha dos valores é baseada nas referências [⁴] e [³] e o forçante ρ é o parâmetro variável escolhido.

O intervalo de analise do parâmetro foi definido como 5≤ρ≤7.5, no qual foram analisados 700 valores de ρ, ou seja, Δρ = 2.5/700. Além disso, para cada ρ foram escolhidas 625 condições iniciais diferentes, considerando pares equidistantes de θ_1,2 ∈ [0,2π), com Δθ_1,2 = 2π/25, e ${\stackrel{.}{\theta }}_{1,2}=0$.

É possível escolher qualquer uma das 4 variáveis dinâmicas para construir o diagrama de bifurcação. Inicialmente, considera-se a evolução da variável θ₂.

A definição do transiente para o RDP possui certa particularidade, pois este sistema possui transientes consideravelmente elevados após transformações globais sofridas por atratores caóticos, o que aumenta o tempo computacional gasto na construção do diagrama. A Figura 2 mostra os diagramas de bifurcação obtidos considerando as condições inicias descritas anteriormente e variando o transiente desprezado. No primeiro gráfico, Figura 2(a), foram desprezadas 6×10³ iterações, enquanto que no segundo, Figura 2(b), o valor do transiente considerado foi de 10⁶. Por fim, na Figura 2(c) foram utilizados dois valores de transientes diferentes: 10⁸ iterações para ρ ∈ [7.0142,7.25] e 10⁶ iterações no restante do intervalo. Em todos os casos salvou-se apenas 1 iteração após o transiente para cada condição inicial.

Figura 2
Diagramas de bifurcação do RDP considerando diferentes transientes.

Pode-se observar que, apesar das condições iniciais serem as mesmas, os diagramas são visualmente distintos no intervalo de 6.8≤ρ≤7.5. Em todos os casos percebe-se que o sistema evolui a partir de uma solução periódica que se inicia a partir da perda de estabilidade de um ponto fixo na região anterior à explorada, em ρ≈4.2771. Essa solução periódica perdura coexistindo com pequenos atratores caóticos a medida que ρ aumenta, como pode ser visto, por exemplo, em ρ≈5.25,6.0 e 6.25. Em ρ≈6.4 a órbita de período 2 bifurca novamente e desenvolve uma cascata de duplicação de período até o caos.

A principal diferença nos diagramas de bifurcação ocorre para ρ > 6.8. Nessa região pode-se notar a importância da escolha correta do transiente. Na Figura 2(a), observa-se que a região a partir de ρ = 6.8 aparenta ter comportamento caótico. Ao aumentar o transiente para 10⁶ iterações, conforme visto na Figura 2(b), percebe-se que a região caótica vista anteriormente desaparece a partir de ρ≈7.0142, dando lugar a um par de soluções periódicas que sofrem cascatas de duplicação de período, originando dois atratores caóticos distintos. Na região entre 7.0142≤ρ≤7.25 da Figura 2(b) ainda podem ser percebidos pontos que correspondem à resposta transiente do sistema. Com efeito, ao aumentar transiente para 10⁸ iterações na Figura 2(c), o diagrama de bifurcação fica mais limpo nessa região, exibindo apenas o par de órbitas periódicas que coexistem no espaço de soluções.

3.2. Cálculo de Pontos Fixos

Seja uma função f em ℝ^m. Um ponto fixo de um sistema dinâmico é dado por x* tal que f(x*) = x*, onde f(x) é a regra matemática que define o sistema e x ∈ ℝ^m é o vetor de variáveis dinâmicas do sistema. De maneira análoga, órbitas periódicas são soluções de um sistema dinâmico que se repetem após um certo período k, isto é, ${\mathbf{\text{x}}}_{n+1}^{*}=\mathbf{\text{f}}({\mathbf{\text{x}}}_n^{*})$ para n = 1,…,k−1 e ${\mathbf{\text{x}}}_1^{*}=\mathbf{\text{f}}({\mathbf{\text{x}}}_k^{*})$ [⁷].

Pontos fixos do RDP. Os pontos fixos do RDP são calculados a partir da equação (8), resolvendo o sistema de equações

(10) $$\begin{array}{l}{\theta }_1^{*}={\theta }_1^{*}+M_{11}{\dot{\theta }}_1^{*}+M_{12}{\dot{\theta }}_2^{*}\\ {\theta }_2^{*}={\theta }_2^{*}+M_{21}{\dot{\theta }}_1^{*}+M_{22}{\dot{\theta }}_2^{*}\\ {\dot{\theta }}_1^{*}=\frac{\rho l_1}{I_1}sen{\theta }_1^{*}+L_{11}{\dot{\theta }}_1^{*}+L_{12}{\dot{\theta }}_2^{*}\\ {\dot{\theta }}_2^{*}=\frac{\rho l_2}{I_2}sen{\theta }_2^{*}+L_{21}{\dot{\theta }}_1^{*}+L_{22}{\dot{\theta }}_2^{*}.\end{array}$$

Nesta seção, serão considerados novamente os mesmos valores de parâmetros fixos descritos na seção ^3.1, variando apenas o parâmetro ρ.

Realizando as devidas manipulações algébricas de forma a determinar a expressão analítica de cada variável dinâmica tem-se que os pontos fixos do RDP são dados por:

(11) $${\stackrel{.}{\theta }}_1^{*}={\displaystyle \frac{2\pi n_1{M}_{22}-2\pi n_2{M}_{12}}{M_{11}{M}_{22}-M_{12}{M}_{21}}},$$

(12) $$\stackrel{.}{{\theta }_2^{*}}={\displaystyle \frac{2\pi n_2{M}_{11}-2\pi n_1{M}_{2}1}{M_{11}{M}_{22}-M_{12}{M}_{21}}},$$

(13) $${\theta }_1^{*}={\text{sen}}^{-1}\left({\displaystyle \frac{{\stackrel{.}{\theta }}_1^{*}(1-L_{11})-L_{12}{\stackrel{.}{\theta }}_2^{*}}{J_1\rho }}\right),$$

(14) $${\theta }_2^{*}={\text{sen}}^{-1}\left({\displaystyle \frac{\stackrel{.}{{\theta }_2^{*}}(1-L_{22})-L_{21}{\stackrel{.}{\theta }}_1^{*}}{J_2\rho }}\right),$$

onde,

(15) $$n_1=\pm 1,n_2=\pm 2,J_{1,2}=\frac{l_{1,2}}{I_{1,2}},$$

sendo n_1,2 valores inteiros relacionados a rotação do sistema, uma vez que são utilizadas variáveis angulares moduladas entre 0 e 2π. As combinações de n_1,2 utilizadas neste trabalho foram [n₁,n₂] = [1,2] e [n₁,n₂] = [−1,−2], uma vez que valores inteiros diferentes dos considerados não fornecem soluções reais no intervalo estudado do parâmetro ρ.

Pode-se observar que ${\stackrel{.}{\theta }}_1^{*}$ e $\stackrel{.}{{\theta }_2^{*}}$, dados pela equações (11) e (12), dependem de parâmetros mantidos constantes, variando apenas em função de n_1,2. Por outro lado, as variáveis ${\theta }_1^{*}$ e ${\theta }_2^{*}$ dependem explicitamente de ρ.

Outro aspecto particular deste sistema se deve a função sen⁻¹(x) presente nas equações (13) e (14). Essa função possui domínio no intervalo D = [−1,1] e imagem em $I=[-\frac{\pi }{2},\frac{\pi }{2}]$, como pode ser observado na Figura 3.

Figura 3
Gráfico da função sen⁻¹ (x).

Devido a imagem da função ser definida no intervalo de $[-\frac{\pi }{2},\frac{\pi }{2}]$ e os valores de θ_1,2 no sistema estarem modulados entre 0 e 2π, tem-se os dois casos a seguir, conforme ilustrado na Figura 4.

Figura 4
Representação do círculo trigonométrico unitário.

Caso 1:

$$0\le x\le 1\{\begin{array}{cc}{\text{sen}}^{-1}(x)=\theta \hfill & \\ {\text{sen}}^{-1}(x)=\pi -\theta \hfill & \end{array}$$

Caso 2:

$$-1\le x<0\{\begin{array}{cc}{\text{sen}}^{-1}(x)=\theta +2\pi \hfill & \\ {\text{sen}}^{-1}(x)=\pi -\theta \hfill & \end{array}$$

No Caso 1, quando 0≤x≤1, há duas soluções possíveis. A primeira corresponde ao valor de θ no primeiro quadrante, isto é, $\theta \in [0,\frac{\pi }{2}]$, enquanto que a segunda, corresponde ao valor de θ no segundo quadrante, conforme visto na Figura 4(a). Considerando a modulação de θ entre 0 e 2π, essas duas soluções possíveis correspondem a dois pontos fixos com coordenadas θ_1,2 ∈ [0,π].

Já no Caso 2, quando −1≤x < 0, há outras duas soluções possíveis: um valor de θ no quarto quadrante e um segundo valor de θ no terceiro quadrante, conforme ilustrado na Figura 4(b). Agora, considerando a modulação de θ entre 0 e 2π, essas duas soluções possíveis correspondem a dois pontos fixos com coordenadas θ_1,2 ∈ [π,2π].

Dessa forma, obtém-se os valores possíveis de θ em todos os quadrantes. Portanto, para cada valor de ρ e um dado par (n₁,n₂), existe um total de 4 pontos fixos.

Analisando o argumento da função sin⁻¹⁡(x) é possível determinar para quais valores de ρ o sistema possui ou não pontos fixos, uma vez que valores de x fora do domínio da função sin⁻¹ (x), isto é, fora do intervalo D = [−1,1], não fornecem soluções reais para as equações (13) e (14).

Para que a a função sin⁻¹⁡(x) possua solução real temos as seguintes relações: primeiramente, partindo da equação (13), tem-se a equação (16) que está relacionada aos valores de ${\theta }_1^{*}$:

(16) $$1\ge \left|{\displaystyle \frac{{\stackrel{.}{\theta }}_1^{*}(1-L_{11})-L_{12}{\stackrel{.}{\theta }}_2^{*}}{J_1\rho }}\right|$$ $${\rho }_{{\theta }_1^{*}}\ge \left|{\displaystyle \frac{{\stackrel{.}{\theta }}_1^{*}(1-L_{11})-L_{12}{\stackrel{.}{\theta }}_2^{*}}{J_1}}\right|.$$

Por outro lado, utilizando a equação (14) obtém-se os valores relacionados a ${\theta }_2^{*}$:

(17) $$1\ge \left|{\displaystyle \frac{{\stackrel{.}{\theta }}_2^{*}(1-L_{22})-L_{21}{\stackrel{.}{\theta }}_1^{*}}{J_2\rho }}\right|$$ $${\rho }_{{\theta }_2^{*}}\ge \left|{\displaystyle \frac{{\stackrel{.}{\theta }}_2^{*}(1-L_{22})-L_{21}{\stackrel{.}{\theta }}_1^{*}}{J_2}}\right|.$$

Para os parâmetros utilizados neste trabalho, os valores de ${\stackrel{.}{\theta }}_{1,2}^{*}$ dos pontos fixos para os pares de números de rotação escolhidos são mostrados na Tabela 1. Utilizando esses valores nas equações (16) e (17), encontram-se os valores mínimos de ρ para os quais as equações (13) e (14) têm solução real. No caso, tem-se ${\rho }_{{\theta }_1^{*}}\ge 1.25607\times {10}^{-15}$ e ${\rho }_{{\theta }_2^{*}}\ge 6.283185$. Uma vez que ${\rho }_{{\theta }_2^{*}}>{\rho }_{{\theta }_1^{*}}$, concluí-se que a função sin⁻¹⁡(x) relacionada a variável θ₂ implica que apenas é possível encontrar analiticamente pontos fixos para o sistema quando ρ≥6.28315.

TABLE 1
Coordenadas constantes dos pontos fixos do RDP em função de n_1,2.

Caso	n 1	n 2	θ 1	${\stackrel{.}{\theta }}_1$	${\stackrel{.}{\theta }}_2$
1	1	2	0	5.81964	16.19396
2	1	2	π	5.81964	16.19396
3	−1	−2	2π	−5.81964	−16.19396
4	−1	−2	π	−5.81964	−16.19396

Além dos valores de ${\stackrel{.}{\theta }}_{1,2}^{*}$ dos pontos fixos serem independentes de ρ para ρ≥6.28315, devido às simetrias do modelo, os valores de ${\theta }_1^{*}$ também independem do forçante, como indicado na Tabela 1. Assim, apenas a coordenada ${\theta }_2^{*}$ dos pontos fixos varia com ρ.

Os valores de ${\theta }_2^{*}$ são calculados utilizando a equação (14). A Figura 5 mostra os valores de ${\theta }_2^{*}$ dos quatro pontos fixos possíveis dos sistema, calculados analiticamente, sobrepostos a um digrama de bifurcação para ilustrar a coexistência dos pontos fixos com outras estados assintóticos.

Figura 5
Pontos fixos do RDP calculados analiticamente e sobrepostos a um diagrama de bifurcação do sistema mostrando a variável θ₂ em função de ρ. As cores referem-se aos casos listados na Tabela 1, sendo o amarelo referente ao caso 1, o azul referente ao caso 2, o vermelho ao caso 3 e o verde ao caso 4.

Vale destacar que o diagrama de bifurcação mostrado na Figura 5 foi gerado iterando condições iniciais com velocidades angulares iniciais nulas, isto é, ${\stackrel{.}{\theta }}_{1,2}=0$ e que essas condições iniciais não fazem parte das bacias de atração dos pontos fixos calculados analiticamente.

3.3. Análise de Estabilidade de Pontos Fixos

A análise de estabilidade consiste em determinar se as soluções assintóticas de um sistema são estáveis ou instáveis. De forma geral, a estabilidade assintótica de pontos fixos está relaciona à situação na qual sucessivas iterações do mapa fazem com que trajetórias na vizinhança do ponto fixo se aproximem dele. Por outro lado, pontos fixos assintoticamente instáveis são aqueles para os quais, após sucessivas iterações do mapa, as trajetórias na vizinhança tendem a se afastar do ponto. Essa noção simplificada não abrange todas as possíveis situações dinâmicas, uma vez que o estudo da estabilidade de soluções é uma questão muito mais complexa com grande importância para aplicações práticas [⁸].

No caso de mapas unidimensionais dados genericamente por x_{n + 1} = f(x_n), a linearização do sistema em torno de um ponto fixo x^∗ fornece η_{n + 1} = f′(x^∗)η_n, onde η representa uma perturbação em torno de x^∗ e f′(x^∗) é a derivada de f(x) com respeito a x avaliada em x^∗. Assim, se o autovalor ou multiplicador λ = f′(x^∗) for tal que |λ| < 1, o ponto fixo é dito estável. Por outro lado, se |λ| > 1 o ponto fixo é instável, enquanto que para o caso marginal |λ| = 1 o teste é inconclusivo, sendo necessário recorrer a análises de ordens superiores [²].

Ao considerar mapas em ℝ^m, com m > 1, a análise de estabilidade linear dos pontos fixos do sistema exige a avaliação dos autovalores da matriz Jacobiana J do sistema. Assim, seja um mapa f em ℝ^m com um ponto fixo p*. A matriz Jacobiana J do sistema é a matriz de derivadas parciais de primeira ordem de f. Se a magnitude de cada autovalor de J(p*) é menor que 1, então p* é estável. Por outro lado, se a magnitude de cada autovalor de J(p*) é maior que 1, p* é instável. Quando m > 1, podem ocorrer simultaneamente autovalores de J(p*) com magnitude menor que 1 e autovalores de J(p*) com magnitude maior que 1. Nesse caso, haverá, simultaneamente, direções estáveis e direções instáveis. Se nenhum dos autovalores de J(p*) tiver magnitude igual à 1, então p* é um ponto fixo hiperbólico. Se p* é hiperbólico e se pelo menos um autovalor de J(p*) tem magnitude maior que 1 e pelo menos um valor próprio tem magnitude menor que 1, então p* é chamado de ponto de sela e é instável. Quando o ponto fixo não é hiperbólico é necessário recorrer a análises de ordens superiores [²].

Análise de estabilidade dos pontos fixos do RDP. Para realizar a análise de estabilidade dos pontos fixos do RDP é necessário determinar a matriz Jacobiana do sistema dada por:

$$\mathbf{\text{J}}=\left[\begin{array}{cccc}\hfill \frac{\delta F_1}{\delta {\theta }_1}\hfill & \hfill \frac{\delta F_1}{\delta {\theta }_2}\hfill & \hfill \frac{\delta F_1}{\delta {\stackrel{\_}{\theta }}_1}\hfill & \hfill \frac{\delta F_1}{\delta {\stackrel{\_}{\theta }}_2}\hfill \\ \hfill \frac{\delta F_2}{\delta {\theta }_1}\hfill & \hfill \frac{\delta F_2}{\delta {\theta }_2}\hfill & \hfill \frac{\delta F_2}{\delta {\stackrel{\_}{\theta }}_1}\hfill & \hfill \frac{\delta F_2}{\delta {\stackrel{\_}{\theta }}_2}\hfill \\ \hfill \frac{\delta F_3}{\delta {\theta }_1}\hfill & \hfill \frac{\delta F_3}{\delta {\theta }_2}\hfill & \hfill \frac{\delta F_3}{\delta {\stackrel{\_}{\theta }}_1}\hfill & \hfill \frac{\delta F_3}{\delta {\stackrel{\_}{\theta }}_2}\hfill \\ \hfill \frac{\delta F_4}{\delta {\theta }_1}\hfill & \hfill \frac{\delta F_4}{\delta {\theta }_2}\hfill & \hfill \frac{\delta F_4}{\delta {\stackrel{\_}{\theta }}_1}\hfill & \hfill \frac{\delta F_4}{\delta {\stackrel{\_}{\theta }}_2}\hfill \end{array}\right],$$

onde F_i, com i = 1,2,3,4 são as quatro equações que descrevem o RDP dadas pelo sistema de equações (8). Assim,

$$\mathbf{\text{J}}=\left[\begin{array}{cccc}\hfill 1\hfill & \hfill 0\hfill & \hfill M_{11}\hfill & \hfill M_{12}\hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill M_{21}\hfill & \hfill M_{22}\hfill \\ \hfill C_{1}cos({\theta }_1)\hfill & \hfill 0\hfill & \hfill L_{11}+M_{11}{C}_{1}cos({\theta }_1)\hfill & \hfill L_{22}+M_{12}{C}_{1}cos({\theta }_1)\hfill \\ \hfill 0\hfill & \hfill C_{2}cos({\theta }_2)\hfill & \hfill L_{21}+M_{21}{C}_{2}cos({\theta }_2)\hfill & \hfill L_{22}+M_{22}{C}_{2}cos({\theta }_2)\hfill \end{array}\right]$$

com C₁ = ρL₁/I₁ e C₂ = ρL₂/I₂.

Após determinar a matriz Jacobiana para cada ponto fixo, são calculados os autovalores da matriz em cada caso, determinando a estabilidade das soluções. A Figura 6(a) mostra as magnitudes dos autovalores para o ponto fixo identificado como Caso 2 na Tabela 1. Na figura, as curvas vermelha e amarela representam o par de autovalores complexos conjugados, enquanto que as curvas em cor azul e verde mostram o par de autovalores reais.

Figura 6
(a) Magnitudes dos autovalores da matriz Jacobiana calculados para o ponto fixo identificado como Caso 2 na Tabela 1, com as curvas vermelha e amarela representando o par de autovalores complexos conjugados e as curvas azul e verde representando o par de autovalores reais. (b) Estabilidade dos pontos fixos listados na Tabela 1, com a cor verde correspondendo a estabilidade e a cor vermelho correspondendo a instabilidade.

Adicionalmente, a Figura 6(b) apresenta o resultado da análise de estabilidade dos pontos fixos listados na Tabela 1. Na figura observa-se que os pontos fixos denominados Casos 1 e 3, mostrados em vermelho, são sempre estáveis desde o seu surgimento em ρ = 2π até ρ = 7.5 quando a análise é interrompida. Por outro lado, os pontos fixos denominados Casos 2 e 4 são soluções estáveis para ρ ∈ [2π,7.01418] e nesse intervalo do espaço de parâmetro coexistem com soluções caóticas. As bacias de atração dos pontos fixos são menos prevalentes no espaço de estados que o atrator caótico, de forma que a dinâmica do sistema é predominantemente caótica. Quando os dois pontos fixos se tornam instáveis em ρ≈7.01418, os únicos estados assintóticos do sistema são duas órbitas periódicas de período 2 que coexistem no espaço de estados. A medida que ρ aumenta, as órbitas sofrem uma cascata de de duplicações de período, desenvolvendo uma conhecida rota para o caos.

3.4. Expoentes de Lyapunov

De forma simplificada, um sistema caótico é definido apresenta dependência sensível às condições iniciais, o que pode levar a uma completa perda de previsibilidade do sistema [²]. Isso significa que, em um sistema que apresenta comportamento caótico, duas condições iniciais que diferem apenas por uma quantidade infinitesimal podem apresentar um comportamento totalmente divergente a medida que o sistema evolui no tempo. Por isso, é importante conhecer um método que permita caracterizar quantitativamente a ocorrência de fenômenos caóticos.

No final do século XIX, o matemático e físico russo Aleksandr Mikhailovich Lyapunov, ao defender sua tese “O problema geral da estabilidade do movimento ” na Universidade de Moscou, desenvolveu conceitos de estabilidade que resultaram posteriormente no chamado método dos Expoentes de Lyapunov para a caracterização de regimes caóticos [⁹].

Esse método consiste em medir o quanto duas trajetórias inicialmente próximas divergem com o passar do tempo [²]. Para isso, considerando sistemas unidimensionais, utiliza-se uma função δ(t) que mede a distância entre as trajetórias a cada iteração do sistema, de forma que δ(0) seja a distância inicial, δ(1) a distância após a primeira iteração e δ(n) a distância após a n-ésima iteração, conforme ilustrado na Figura 7.

Figura 7
Representação gráfica do conceito dos Expoentes de Lyapunov.

Para caracterizar a presença de caos em um sistema, a função δ(0) deve possuir um valor infinitesimal, da ordem de 10⁻¹⁵, e a função δ(n) deve assumir um valor alto quando n→∞. Caso contrário, δ(n) permanece aproximadamente igual a δ(0) [¹⁰]. Assim, vale a a seguinte relação:

(18) $$||\delta (n)||\approx ||\delta (0)||e^{\lambda n},$$

onde λ é chamado expoente de Lyapunov.

Analisando a equação (18) percebe-se que, quando n→∞, há três possibilidades em função de λ.

• Se λ > 0, e^λn→∞, portanto ||δ(n)||→∞. Assim, fica caracterizada uma divergência exponencial entre as trajetórias e, consequentemente, a presença de caos.

• Se λ < 0, e^λn→1/∞ que por sua vez tende a zero, portanto ||δ(n)||→0. Assim, não fica caracterizada uma dependência sensível às condições iniciais, de forma que a trajetória não é caótica.

• Se λ = 0, ocorre uma indeterminação (∞×0) o que pode indicar a iminência de uma bifurcação ou o início de um regime caótico.

Reescrevendo a equação (18) para isolar o expoente de Lyapunov no limite n→∞, tem-se:

(19) $$\lambda =\underset{n\to \mathrm{\infty }}{lim}\frac{1}{n}\ln \left|\frac{\delta (n)}{\delta (0)}\right|.$$

Sabendo que δ(n) = fⁿ(x₀ + δ(0))−fⁿ(x₀), ou seja, que a função δ(n) pode ser escrita como a diferença entre a n-ésima iteração do sistema fⁿ e a condição inicial x₀ e conhecendo a condição inicial infinitesimalmente próxima x0 + δ(0), tem-se:

(20) $$\lambda =\underset{n\to \mathrm{\infty }}{lim}\frac{1}{n}\ln \left|\frac{f^n(x_0+\delta (0))-f^n(x_0)}{\delta (0)}\right|.$$

A equação 20 remete à definição de derivadas, de forma que

(21) $$\lambda =\underset{n\to \mathrm{\infty }}{lim}\frac{1}{n}\ln \left|{(f^n)}^{\prime }(x_0)\right|.$$

Expandindo o termo derivável e utilizando a regra da cadeia, o expoente de Lyapunov é dado formalmente por:

(22) $${(f^n)}^{\prime }(x_0)={\displaystyle \prod _{i=0}^{n-1}}{f}^{\prime }(x_i)$$

(23) $$\lambda =\underset{n\to \mathrm{\infty }}{lim}{\displaystyle \frac{1}{n}}\ln \left|{\displaystyle \prod _{i=0}^{n-1}}{f}^{\prime }(x_i)\right|$$

(24) $$\lambda =\underset{n\to \mathrm{\infty }}{lim}\left\{{\displaystyle \frac{1}{n}}\ln \left|{\displaystyle \sum _{i=0}^{n-1}}{f}^{\prime }(x_i)\right|\right\}.$$

Para sistemas multidimensionais o cálculo dos expoentes de Lyapunov segue o mesmo princípio, mas é necessário observar que para cada trajetória obtém-se um espectro de m valores, com m sendo a dimensão do sistema. Nesse caso, para determinar se uma trajetória tem comportamento caótico, observa-se o valor do maior expoente de Lyapunov, λ_max, já que este ocorre na direção de maior afastamento entre a trajetória de referência e a trajetória teste. Assim, de forma análoga ao caso unidimensional, para λ_max > 0 verifica-se a presença de caos no sistema [¹⁰].

De forma prática, o procedimento básico descrito a seguir pode ser utilizado para obter numericamente o expoente de Lyapunov máximo [¹⁰, ¹¹].

1. Escolher uma condição inicial na bacia de atração de um dos atratores do sistema.

2. Iterar o mapa até que a órbita esteja no atrator. Se a condição inicial for escolhida no atrator, este passo pode ser omitido.

3. Definir uma condição inicial $x_0^a$ no atrator para a trajetória de referência e uma segunda condição inicial $x_0^b$ a uma distância infinitesimal d₀ de $x_0^a$ para a trajetória teste.

4. Iterar o sistema a partir dos estados $x_0^a$ e $x_0^b$ para obter novos estados $x_1^a$ e $x_1^b$, respectivamente.

5. Calcular a distância d₁ entre $x_1^a$ e $x_1^b$.

6. Calcular o logaritmo natural de d₁/d₀.

7. Reajustar a trajetória teste de forma que sua distância à trajetória de referência seja d₀, na mesma direção de d₁.

8. Repetir o procedimento dos itens 4 a 7 diversas vezes e calcular a média da soma dos valores encontrados no item 6.

A escolha do número de iterações transiente depende de algum conhecimento prévio sobre o sistema em estudo. Para a maioria dos sistemas, algumas centenas de iterações são suficientes para atingir um estado assintótico, mas pode ser necessário aumentar consideravelmente este valor próximo a bifurcações ou transformações globais no espaço de fases.

Após desprezar o transiente, a ultima iteração do mapa é tomada como $x_0^a$ e define-se $x_0^b$ a uma distância d₀ de $x_0^a$. Essa distância depende da precisão e da arquitetura do sistema computacional empregado. Este último passo é executado até que o valor de λ convirja para uma dada precisão desejada.

Em geral, a convergência é lenta, mas no caso de mapeamentos de dimensão baixa alguns milhares de iterações são suficientes para obter uma estimativa com precisão de cerca de dois algarismos significativos. Para melhorar o tempo de convergência, ao calcular a média, pode-se descartar os primeiros valores obtidos para λ, de forma a assegurar que as trajetórias estão de fato orientadas ao longo da direção de máxima expansão.

É interessante introduzir um critério de parada adicional baseado no número de repetições do cálculo da média. Finalmente, para garantir a validade do resultado obtido para um dado estado assintótico é importante verificar se este independe das condições iniciais escolhidas na bacia, do valor de d₀ e do número de iterações incluídas na média.

O procedimento descrito pode ser traduzido em pseudocódigo conforme mostrado no Algoritmo 2 contido no Apêndice. Vale ressaltar que o expoente de Lyapunov não tem significado algum para órbitas ilimitadas devido aos erros numéricos que são acumulados no cálculo nesse caso.

Expoentes de Lyapunov do RDP. A fim de exemplificar o uso dos expoentes de Lyapunov considerou-se a dinâmica do RDP com ρ ∈ [5,7.5]. Conforme discutido anteriormente, no caso do RDP podem ser identificados diferentes regimes em termos do tempo de transiente, especialmente após bifurcações ou crises de atratores caóticos. Além disso, atratores que coexistem para um mesmo par de parâmetros do sistema podem apresentar transientes característicos muito distintos. Assim, optou-se por descartar em torno de 10⁵ iterações transientes para garantir que o cálculo do expoente de Lyapunov fosse iniciado a partir de um estado assintótico e não durante um longo transiente caótico.

Escolheu-se d₀ da ordem de 10⁻⁸ em apenas uma das variáveis dinâmicas, sendo as distâncias calculadas a partir da distância euclidiana no espaço 4-dimensional,

$$d_1=\sqrt{\sum _i^4{(x_{1_i}^a-x_{1_i}^b)}^2},$$

onde o subíndice i denota cada uma variáveis dinâmicas do sistema. Neste caso, x₁ = θ₁, x₂ = θ₂, $x_3={\stackrel{.}{\theta }}_1$ e $x_4={\stackrel{.}{\theta }}_2$.

A cada iteração do mapa, a órbita b foi reajustada de modo que sua distância da órbita a fosse d₀ ao longo da mesma direção de d₁. Assim, a cada nova iteração as variáveis de estado da órbita b são reinicializadas como:

$$x_{0_i}^b=x_{1_i}^a+(d_0/d_1)(x_{1_i}^b-x_{1_i}^a).$$

No caso do RDP o cálculo da média foi repetido por aproximadamente 50 mil vezes, descartando-se 5 mil valores iniciais, até obter uma convergência adequada. Os máximos expoentes de Lyapunov obtidos ao seguir uma condição inicial desde ρ = 5.0 é mostrado na Figura 8.

Figura 8
Gráfico com os maiores expoentes de Lyapunov, em laranja, sobreposto ao diagrama de bifurcação parcial do sistema.

No gráfico mostrado, a coordenada final da trajetória para um dado valor de ρ foi utilizada como condição inicial para o cálculo do expoente para o valor seguinte de ρ. Assim, seguimos a dinâmica de uma único estado assintótico que inicialmente é periódico e que evolui para o caos a medida que ρ aumenta. Vale ressaltar que para 6.8⪅ρ⪅7.02 observa-se um salto nos valores positivos encontrados para os expoentes de Lyapunov. Esses valores, estão associados à ocorrência de um longo transiente caótico após a crise do atrator caótico que desaparece em ρ≈6.8, pois o transiente de 10⁵ iterações não é suficiente para atingir os estados assintóticos periódicos que existem nessa região do espaço de parâmetro. Assim, nesse caso, o uso do expoente de Lyapunov como forma de caracterizar a presença de caos no sistema deve ser avaliado com base no conhecimento de aspectos adicionais da dinâmica.

4. Considerações finais

Neste trabalho foram abordadas diversas ferramentas utilizadas no estudo de sistemas dinâmicos não lineares caóticos, apresentando aplicações destas ferramentas ao Rotor Duplo Pulsado. Em particular, foram obtidos diagramas de bifurcação e avaliados os transientes para diferentes valores do parâmetro forçante sistema. Além disso, foram calculados pontos de equilíbrio, analisando sua estabilidade. Por fim, foi realizado o cálculo de expoentes de Lyapunov, discutindo sua aplicabilidade.

Espera-se que este estudo possa contribuir para outros trabalhos envolvendo sistemas dinâmicos representados por mapas com mais de duas dimensões, uma vez que a maior parte da literatura costuma abordar estas ferramentas somente até sistemas uni e bidimensionais.

Agradecimentos

O primeiro autor agradece o apoio financeiro da FAPESP (processo 2018/25001-3) e a segunda autora agradece o apoio financeiro do CNPq (processo 422282/ 2018-9) e da FAPESP (processo 2018/00059-9). Os autores também registram seu agradecimento à CAPES por manter o Portal de Periódicos e à UNESP por proporcionar acesso a essa plataforma e pelo apoio institucional.

Apêndice

Algoritmos

REFERÊNCIAS

• ^1.
  L.H.A. Monteiro, Sistemas Dinâmicos (Editora Livraria da Física, São Paulo, 2011).
• ^2.
  K. Alligood, T. Sauer e J. Yorke, Chaos: An Introduction to Dynamical Systems (Springer, New York, 1996).
• ^3.
  F.J. Romeiras, C. Grebogi, E. Ott e W.P. Dayawansa, Physica D: Nonlinear Phenomena 58, 165 (1992).
• ^4.
  U. Feudel, C. Grebogi, L. Poon e J. Yorke, Chaos, Solitons and Fractals 9, 171 (1998).
• ^5.
  C. Grebogi, E. Kostelich, E. Ott e J. Yorke, Physica D: Nonlinear Phenomena 25, 347 (1987).
• ^6.
  E. Ott, Chaos in dynamical systems (Cambridge University Press, Cambridge, 2002).
• ^7.
  T. Parker e L. Chua, Practical Numerical Algorithms for Chaotic Systems (Springer, New York, 1989).
• ^8.
  P. Holmes e E. Shea-Brown, Scholarpedia 1, 1838 (2006).
• ^9.
  J.J. O’Connor e E. Robertson, “Aleksandr Mikhailovich Lyapunov”, disponível em http://www-history.mcs.st-andrews.ac.uk/Biographies/Lyapunov.html, acessado em 08/08/2019.
  » http://www-history.mcs.st-andrews.ac.uk/Biographies/Lyapunov.html
• ^10.
  S. Strogatz, Nonlinear Dynamics And Chaos (Perseus Books, Reading, 1994).
• ^11.
  J. Sprott, “Numerical Calculation of Largest Lyapunov Exponent”, disponível em http://sprott.physics.wisc.edu/chaos/lyapexp.htm, acessado em 25/06/2020.
  » http://sprott.physics.wisc.edu/chaos/lyapexp.htm

Datas de Publicação

• Publicação nesta coleção
  31 Mar 2021
• Data do Fascículo
  2021

Histórico

• Recebido
  06 Jan 2021
• Revisado
  03 Mar 2021
• Aceito
  08 Mar 2021