Desenvolvimento de um aplicativo android para a análise do circuito de Chua-Matsumoto

Development of an android application for the analysis of the Chua-Matsumoto circuit

Um aplicativo Android para a análise do circuito de Chua-Matsumoto foi desenvolvido. Este aplicativo fornece a trajetória de solução para o sistema dinâmico do circuito com condição inicial (x₀,y₀,z₀). O sistema de equações diferenciais não-lineares é resolvido pelo método de Runge-Kutta-Fehlberg. O aplicativo possui três telas de entrada de dados: uma para as condições iniciais (x₀,y₀,z₀), uma para os parâmetros do sistema (α, β,a,b), e uma para os parâmetros do método (t_max,h_max,h_min,tol). Há uma tela utilizada para traçar as curvas das coordenadasx,y e z como funções do tempo. Também, há uma tela utilizada para traçar as projeções da trajetória nos planosxy, xz e yz e para exportar todos os dados para um arquivo txt. Além disso, há uma tela utilizada para exibir os expoentes de Lyapunov calculados pelo aplicativo.

Palavras-chave:
circuito de Chua-Matsumoto; sistemas dinâmicos; atratores caóticos; aplicativos Android

Abstract

An Android application for the analysis of the circuit of Chua-Matsumoto was developed. This application provides the solution trajectory to the dynamical system of the circuit with the initial condition (x₀,y₀,z₀). The system of nonlinear differential equations is solved by Runge-Kutta-Fehlberg method. The application has three data entry screens: one for the initial conditions (x₀,y₀,z₀), one for the parameters of the system (α,β,a,b), and one for the parameters of the method (t_max,h_max,h_min,tol). There is a screen used to plot the curves of the coordinatesx,y and z as functions of the time. Also, there is a screen used to plot the projections of the trajectory into thexy, xz and yz planes and to export all data to a txt file. Furthermore, there is a screen used to display the Lyapunov exponents calculated by the application.

Keywords:
circuit of Chua-Matsumoto; dynamical systems; chaotic attractors; Android applications

1. Introdução

A ideia de caos relaciona-se à existência de respostas imprevisíveis para um sistema determinístico. Todavia, não existe uma definição única para o caos. Pode-se afirmar que um sistema que apresenta o mesmo comportamento dinâmico complexo do mapa de ferradura de Smale (apresentar uma infinidade contável de órbitas periódicas de período arbitrariamente grande, uma infinidade incontável de órbitas não-periódicas, órbita densa, e sensibilidade às condições iniciais) é caótico.

Durante muitos anos não se soube se o caos representava apenas um conceito matemático, materializado na forma de erros de aproximação nos cálculos numéricos das equações que descrevem sistemas físicos. Na década de 1980, o circuito de Chua-Matsumoto possibilitou a visualização de um atrator caótico na tela do osciloscópio, corroborando experimentalmente que o caos representa o comportamento dinâmico de sistemas físicos reais ^[1]. A importância de estudar o circuito de Chua-Matsumoto está relacionada à maior variedade de fenômenos que ele apresenta, os quais têm sido verificados experimentalmente em laboratório, em comparação com outros sistemas comumente estudados, como Rössler ou Lorenz ^[2]. A criptografia baseada em caos é uma possível aplicação prática do circuito de Chua-Matsumoto ^[3].

O uso de dispositivos que utilizam a plataforma Android encontra-se disseminado entre os estudantes de todos os níveis, trazendo ao professor o desafio de incorporar tais dispositivos ao ensino, em vez de simplesmente condená-los. O presente trabalho foi desenvolvido nesse sentido. Nele apresenta-se o desenvolvimento de um aplicativo em plataforma Android (executável em um Smartphone) para simulações computacionais do circuito de Chua-Matsumoto. O aplicativo é gratuito e está disponível no endereço eletrônicohttp://www.amazon.com.au/Circuit-of -Chua-Matsumoto-and-Chaos.

O aplicativo soluciona o circuito de Chua-Matsumoto, para uma dada condição inicial (x₀,y₀,z₀), através do método Runge-Kutta-Fehlberg. É possível escolher os valores das condições iniciais (x₀,y₀,z₀), dos parâmetros do circuito (α,β,a e b), assim como, os parâmetros do método (t_max, h_max, h_min, tol). É possível plotar os gráficos das coordenadas x, y ez em função do tempo. Também há a possibilidade de visualização das projeções da trajetória nos planos xy, xzeyz. É possível realizar o cálculo dos expoentes de Lyapunov. Além disso, pode-se exportar os dados para arquivos no formato txt, os quais são armazenados na pasta Documents.

Na seção 2, o objeto de estudo do aplicativo, o circuito de Chua-Matsumoto, é descrito. A seção 3 contém alguns conceitos da teoria de sistemas dinâmicos, os quais são necessários para que o usuário seja capaz de entender os resultados gerados pelo aplicativo. A seção 4 descreve o aplicativo e apresenta resultados obtidos com o aplicativo em um Smartphone, os quais são comparados com simulações computacionais executadas em um computador de mesa (desktop) rodando o software MATLAB. A seção 5 é a conclusão do trabalho.

2. Modelo do circuito de Chua-Matsumoto

O circuito mostrado na Fig. 1 é conhecido na literatura como circuito de Chua-Matsumoto, o qual é constituído de um indutor de indutância L, um resistor de condutânciaG, dois capacitores de capacitânciasC₁ e C₂, e um resistor não-linear, N_R, conhecido como diodo de Chua. Os parâmetros dos componentes de circuito,C₁,C₂,G, L, são números reais positivos.

Figura 1
Circuito de Chua-Matsumoto.

O diodo de Chua é um dispositivo não-linear, formado por três regiões lineares (vN_R< −B_p, − B_p ≤ vN_R ≤ B_P e vN_R> BP) de resistências negativas, cuja curva característica corrente-voltagem é dada por

em que B_P é a tensão de limiar. Se o valor absoluto da tensão nos terminais deN_R,vN_R, excede a tensão de limiar, a inclinação da curva característica é G_b; caso contrário, a inclinação da curva característica éG_a (ver Fig. 2).

Figura 2
Curva característica corrente-voltagem do diodo de Chua.

2.1. Modelo em coordenadas adimensionais

O circuito de Chua-Matsumoto é modelado matematicamente, em coordenadas adimensionais, através do sistema de equações diferenciais ordinárias tridimensional

Em que a função f (x) é dada por

As correntes iC₁ eiC₂fluem, respectivamente, através dos capacitoresC₁ eC₂. As tensões nos terminais desses capacitores são vC₁evC₂, como ilustrado na Fig. 1. Além disso, v^C1 é igual a v^NR, pois C₁eN_R possuem os mesmos terminais. A capacitância mede a capacidade de armazenamento de energia elétrica, assim como, a oposição à variação de tensão (inércia de tensão) nos terminais de um capacitor, quando atravessado por uma corrente elétrica; o que é expresso matematicamente como

A corrente i_G flui através do resistor de condutânciaG. Pela primeira Lei de Ohm, tem-se

Aplicando a lei de Kirchhoff das correntes (conservação de carga) respectivamente aos nós l₁ e l₂, obtém-se

Aplicando a Lei de Kirchhoff das tensões (conservação de energia) à malha constituída pelo indutor L e pelo capacitorC₂, obtém-se a equação

As variáveis vN_R,iN_R,vC₂,i_L,t, presentes nas Eqs. (6), (7) e (8) podem ser trocadas por variáveis adimensionais x, y, z, τ, utilizando as fórmulas

Realizando essas mudanças de coordenadas, obtém-se o sistema tridimensional mostrado na Eq. (2), ondeα e β são dados por

A seu turno, a Eq. (3) pode ser obtida substituindo a Eq. (9)naEq. (1), onde

2.2. Implementação do diodo de Chua

O diodo de Chua não é vendido comercialmente, mas pode ser emulado com componentes convencionais ^[4-8]. Na Ref. ^[5], Matsumoto implementou umdiodo de Chua com um amplificador operacional (UA741), um par de diodos, sete resistores e uma fonte de alimentação DC de ±9 V, de acordo com o esquema de circuito mostrado na Fig. 3, resultando em um componente N_Rcom a curva característica da Eq. (1), com os seguintes valores: Ga ≈ −0, 8 mA/V,Gb ≈ −0, 5 mA/V e ±Bp = ±1 V.

Figura 3
Esquema do circuito do diodo de Chua implementado por Matsumoto^[5].

O diodo de Chua também pode ser obtido através de componentes nanoeletrônicos denominados de memoristores ^[9-11]. Embora o conceito de memoristor tenha sido desenvolvido por Leon Chua no início da década de 1970, o primeiro memoristor foi construido por pesquisadores da HP Labs apenas em 2008, e ficou conhecido como memoristor da HP ^[13]. Em 2012, um circuito caótico baseado em memoristores HP foi publicado ^[14].

3. Sistemas dinâmicos

As variáveis x(t),y(t) ez(t), cujas derivadas temporais aparecem no lado esquerdo das equações diferenciais naEq. (2), são as variáveis de estado do sistema. Como as variáveis de estado variam com o tempo, o sistema é designado como sistema dinâmico. As variáveis de estado formam o vetor de estado dado por

As variáveis de um sistema dinâmico que não representam variáveis de estado; comoα, β, a e b; recebem a denominação deparâmetros do sistema. Se esses parâmetros não variam com o tempo t, o sistema é dito a parâmetros fixos.

Pode-se pensar no vetor de estado r(t) como ovetor posição de uma partícula, em que as variáveis x, y e z representam as coordenadas da posição. Uma trajetória é descrita pela partícula com o passar do tempo t.

O espaço de fases (ou espaço de estados) é o conjunto de todas as posições r(t) com t ∈ ℝ. Para a Eq. (2), o espaço de fases é o ℝ³. A fim de solucionar o sistema, deve-se fornecer, em um certo instante de tempo t₀, o vetor de estador₀ =r(t₀). Uma solução é uma trajetória no espaço de fases, iniciada em r₀. Esse problema é conhecido como problema de valor inicial com condição inicialr₀. Para garantir a existência e unicidade da solução de um problema de valor inicial, as funções h_i e ∂h_i / ∂x_i(i = 1, 2, 3) devem ser contínuas.

Um sistema de equações diferenciais é dito autônomo, se a variável independente t não aparece explicitamente nas equações. Portanto, o sistema (2) é autônomo. Um sistema não-autônomo de ordem n, ṙ (t) =h(r(t),t), pode ser transformado em um sistema autônomo de ordemn + 1, ṙ' (t) =h(r'(t)), substituindo o vetor de estado r(t) pelo vetor de estador'(t) = (r(t),t).

Para x ≥ 1, |x−1| = x−1 e |x+1| = x+1; para x ≤ −1, |x − 1| = −x + 1 e |x + 1| = −x − 1; para |x| ≤ 1, |x−1| = −x+1 e |x+1| =x+1; portanto, a Eq. (3) pode ser reescrita como

evidenciando a existência de três diferentes domínios para as soluções da Eq. (2), a saber, D₋₁ = {(x, y, z) ∈ ℝ³| x ≤ −1},D₁ = {(x, y, z) ∈ ℝ³ | x ≥ 1} eD₀ = {(x, y, z) ∈ ℝ³ | − 1 ≥ x ≤ 1}. O planoU₋₁ = {(x, y, z) ∈ ℝ³ | x = −1} é a fronteira entreD_-1 e D₀, e o planoU₁ = {(x, y, z) ∈ ℝ³ | x = 1} é a fronteira entreD₁ e D₀.

No domínio D_-1, f(x) =bx - a + b, entãoẋ = α(y−x − f(x)) = (αb−α)x+αy+α(a-b), portanto, o sistema de equações pode ser escrito na forma matricial como

onde p = α (a −b) e v = b. Definindo a matriz B dados por e o vetor coluna

é possível escrever a equação matricial (15) como

em que o vetor B é chamado de entrada do sistema, o vetorré chamado de saída do sistema e a matriz é chamada de matriz do sistema.

No domínio D₁, f(x) =bx + a - b, entãoẋ = α(y−x − f(x)) = (-αb−α)x+αy+α(-a+b), portanto, o sistema de equações pode ser escrito na forma matricial através da Eq. (15) com p = −a + b e ν=b.

No domínio D₀, f(x) =ax, então ẋ=α(y − x−f(x)) = (−αa−α)x + αy, portanto, o sistema de equações pode ser escrito na forma matricial através da Eq. (15) com p =0 eν = a.

O sistema dinâmico (2) é contínuo, haja vista o vetor de estado r(t) variar continuamente com o tempo t. Um sistema é discreto quando há um conjunto de instantes de tempo {t₀,t₁,t₂, ...,t_j, ...}, comj ∈ ℕ, fora dos quais nunca ocorrem variações no vetor de estado.

O princípio da superposição estabelece que, ser₁(t) é a saída produzida pela entradaB₁(t) er₂(t) é a saída produzida pela entradaB₂(t), então

onde k₁ e k₂ são constantes. Um sistema é dito linear quando obedece ao princípio da superposição, caso contrário, é dito não-linear. Malgrado o sistema (2) para o domínio ℝ³ seja não-linear, ele assume a forma linear (17) para os domíniosD_-1, D₀eD₁, o que pode ser demonstrado utilizando a expressão (18).

Diz-se que um sistema é dissipativo, quando o volume que ele ocupa no espaço de fase contrai-se com o passar do tempo, o que pode ser expresso através da equação

Portanto, −α(b +1)− 1 < 0 e −α(a +1)− 1 < 0 são as condições para que o circuito de Chua-Matsumoto seja dissipativo emD₋₁ ∪D₁ eD₀, respectivamente.

3.1 Pontos de equilíbrio e estabilidade no sentido de Lyapunov

A análise do sistema nos seus pontos de equilíbrio(oupontos fixos) é importante para estudar a estabilidade do sistema. Os pontos de equilíbrio são aqueles em que a partícula está parada, ou seja, se a partícula encontra-se no ponto de equilíbrio r* em um certo instante de tempo t^*, então a partícula encontrar-se-á no mesmo ponto r^* para qualquert>t^*. Matematicamente, pode-se dizer que pontos de equilíbrio r* são aqueles para os quais

Para sistemas de equações diferenciais que descrevem um circuito elétrico, um ponto de equilíbrio corresponde a uma solução DC (direct current) do circuito.

Um ponto regular ou ordinário é todo ponto do espaço de fases que não é ponto de equilíbrio.

Um ponto de equilíbrio r* é assintoticamente estável se todas as trajetórias com condições iniciais contidas numa esfera de raioδ centrada em r* tendem ar* com o passar do tempo. Quando δ é finito o ponto de equilbrio é localmente assintoticamente estável. Quando δ é infinito o ponto de equilíbrio é globalmente assintoticamente estável. Ambos os pontos são atratores. O conjunto das condições iniciais que tendem a um mesmo atrator formam a bacia de atração do atrator.

Um ponto de equilíbrio r* é marginalmente estável se todas as trajetórias com condições iniciais contidas numa esfera de raioδ centrada em r* não tendem ar* com o passar do tempo, mas continuam dentro de uma esfera de raio ∈ centrada em r*.

Um ponto de equilíbrio r* é instável se pelo menos uma trajetória com condição inicial contida numa esfera de raio δcentrada emr* escapa de uma esfera de raio ∈ centrada emr* em um tempo finito.

Os pontos de equilíbrio do sistema da Eq. (2), obtidos quandoẋ =ẏ =ż = 0, são:P₋ = (−θ, 0,θ), emx < −1;P₀ =O = (0, 0, 0), em −1 <x < 1; eP₊ = (θ, 0,−θ) em x > 1; em queθ = (b−a)/(b + 1) > 1 e, portanto, a sentença [(a < −1)∧(b > −1)]∨[(a < −1)∨(b > −1)] deve ser verdadeira, para que existam os pontos de equilíbrioP_- eP₊.

3.2. Linearização

Dado um ponto de equilíbrio P cuja posição ér_P, as expansões em séries de Taylor das componentes do vetor h = (ẋ,ẏ,ż) da Eq. (2) são

onde i =1, 2, 3. Comoh₁(r) =ẋ,h₁(r_P) =ẋ_P,h₂(r) =ẏ,h₂(r_P) = ż_P,h₃(r) = ż eh₃(r_P) =ż_P, as três equações obtidas podem ser expressas na forma matricial como

A Eq. (22) é uma versão linearizada do sistema (2), válida para pontos infinitesimalmente próximos do ponto P. A matriz do sistema linearizado é denominada de matriz jacobiana do sistema, a qual para o circuito de Chua-Matsumoto pode ser obtida a partir da Eq. (2), sendo expressa como

onde ν = b para |x|≥ 1, eν =a para |x|≤ 1. A Eq. (23) coincide com a Eq. (15), pois o circuito de ChuaMatsumoto é um sistema não-linear constituído por três regiões lineares. AEq. (22) pode ser reescrita como

onde r' = r − r_P.

Sabe-se quer(t)=ke^λt é solução da equação diferencial de primeira ordemṙ(t) =λr(t). Por analogia, pode-se pensar em uma solução r'(t) =ve^λtpara a Eq. (24), resultando em λve^λt= ve^λt, logo

onde v é o autovetor associado a λ. Como o circuito de Chua-Matsumoto é tridimensional,₁, λ₂ e λ₃, e três autovetores v₁,v₂ev₃, respectivamente associados. Esses autovalores são obtidos como soluções da equação seculardet(λ é a matriz identidade de ordem 3. Nessa equação, λ é um autovalor de − ,e possui três autovalores λ) = 0, a qual resulta em

onde ν = b para a regiãoD₋₁ ∪ D₁, que contém os pontos de equilíbrio P_-eP₊;e ν =apara a região D₀, que contém o ponto de equilíbrioP₀ = O.

Se o conjunto solução da Eq. (26)não é formado por três autovalores reais, então esse conjunto possui um autovalor real λ₁ = γ e dois autovalores complexos conjugados λ₂ =σ + iω e λ₃ = σ −iω. Nesse caso, se as soluções da Eq. (24) são linearmente independentes, a solução geral é a combinação linear delas, resultando na expressão

em que os autovetores v₂η_R+iη_Iev₃η_R−iη_Iconstituem um par complexo conjugado, com parte realη_R e parte imagináriaη_I. A matriz cujas colunas são os autovetoresv₁v₂v₃) é denominada de matriz modal.= (

Definindo o vetor w(t), dado por

é possível reescrever a Eq. (27)como

Derivando a Eq. (28), tem-se

Como w(0) = [C1C2C3]^T, fazendo t = 0 na Eq. (29), obtém-se a expressão

A partir da Eq. (31), são obtidos o coeficiente real C₁ e o par complexo conjugadoC₂cCe^ΦC eC₃cCe^−ΦC. Reescrevendo a Eq. (27) em termos decC, -Φ_C,η_R eη_I , resulta em

A trajetória r'(t), t ∈ ℝ, pode ser decomposta em duas componentesr'_Rer'_C, dadas por

A trajetória descrita por r'_R é uma linha reta (pois o versor de r'_R não varia com o tempo). O conjunto dos pontos dessa trajetória forma oautoespaço correspondente ao autovalor real γ no ponto fixo P, representado porE_R(P). Por outro lado,r'_C é um vetor tangente ao plano gerado pelos vetores de base η_Reη_R. O conjunto dos pontos desse plano forma o autoespaço correspondente aos autovalores complexos σ ± iω no ponto fixo P, representado porE_C(P). Nota-se que os vetores de base são mantidos quando realizadas as derivações da Eq. (33) ou da Eq. (34) com relação ao tempot. Como consequência disso, se r'∈E^C(P) (ou r'∈E^R(P)) em um instantet = 0, então r'∈E^C (P) (ou r'∈E^R(P)) para qualquert > 0; isto é, uma trajetória iniciada em um autoespaço, evolui nele.

A partir da Eq. (33) conclui-se que, se γ < 0, entãor'_R → 0 quandot → ∞, ou seja, as trajetórias convergem assintoticamente para o ponto fixo PemE_R(P). Todavia, seγ > 0, as trajetórias afastam-se deP.

Da Eq. (34), vê-se que, quandoσ < 0 e ω ≠ 0, as trajetórias emE_C(P) descrevem uma espiral que converge assintoticamente para o ponto fixo P (aEq. (34) é uma elipse multiplicada por uma exponencial que decresce com o tempo), que nesse caso é umfoco espiral atrativo (estável). Seσ > 0 eω ≠ 0, a trajeté uma espiral que se afasta do ponto fixoP, o qual representa um foco espiral repulsivo (instável).

Os autovetores v₁que geramE_R(P) são obtidos a partir da equação (γv₁= 0. Fazendo v₁ = (1,, − ) chega-se ao autovetor expresso como)

Portanto, a equação da reta E_R(P) é dada por

Utilizando [(σ ± iω)v_2,3 = 0 com v_2,3 = (1η_R = (v₂ +v₃)/2eη_I = (v₂−v₃)/(2i), encontra-se, −]), e calculando os autovetores ,,)=() = (,,

Da Eq. (26) conclui-se que

Substituindo a Eq. (38) na Eq. (37), pode-se obter, a partir do produto vetorial η_R ×η_I, um vetor normal ao planoE_C(P). Após algumas manipulações algébricas, com a utilização de Eq. (26) com λ =γ, chega-se à equação deE_C(P dada por

3.3. Forma de Jordan real

Da Eq. (29), tem-se

resultando na variável real w₁ e no par complexo conjugado w₂w_R + iw_Iew₃w_R − iw_I. A matriz é dada por

Da Eq. (30), tem-se que

A matriz do sistema (42) é aforma de Jordan real do circuito de Chua-Matsumoto daEq. (2), cuja validade restringe-se a um dos domínios lineares (as variáveisw_R(t),w_I(t),w₁(t) obtidas emD₀ não são as mesmas obtidas emD₁).

Utilizando a Eq. (40), o vetor de base v₁ do subespaçoE_R(O) transforma-se emv_∼1= (w_R,w_I,w₁) = (0, 0, 1), e os vetores de baseη_R e η_I do subespaço E_C(O) transformam-se emη_∼R = (w_R,w_I, w₁) = (1/2, 0, 0) e η_∼I = (w_R, w_I, w₁) = (0, −1/2, 0). Portanto, o sistema (42) é obtido a partir de (2), através de uma transformação de coordenadas que coloca E_R(O) sobre um dos eixos, no caso w₁,eE_C(O) sobre o plano perpendicular a esse eixo, no caso w₁ = 0.

3.4. Ponto fixo do tipo sela-foco e órbita homoclínica

A trajetória que passa por um ponto de equilíbrio P é referida como órbita homoclínica se tende a Pquandot → ∞ e t → − ∞. O pontoP é denominado de ponto homoclínico. Para facilitar a compreensão, utilizar-se-á o exemplo ilustrado na Fig. 4.

Figura 4
Órbita homoclínica no ponto fixo O = (0,0,0), para os parâmetrosα = 11.0917459, β= 14.2857143, a = -1.1428571 e b = -0.7142857. A trajetória foi obtida através da Eq. (27).

Se σ γ < 0 para um ponto de equilíbrioO com autovalores λ₁ = γ, λ₂ =σ + iω, λ₃ = σ−iω; então esse ponto é do tipo sela, pois dados os subespaços E_R(O) eE_C(O), se O é assintoticamente estável em um deles, será instável no outro. Além disso, seω ≠ 0, então as trajetórias emE_C(O) são focos espirais. Um ponto fixo do tipo sela-foco é um ponto fixo do tipo sela comω≠ 0.

Para o ponto fixo O, a posição dos pontos da trajetória ér(t) =r'(t)+O =r'(t). Caso o ponto fixo fosseP_+/-, a posição dos pontos da trajetória seriar(t) =r'(t)+P_+/−.

Uma órbita homoclínica é obtida no ponto fixo O = (0,0,0), para os parâmetrosα = 11.0917459, β = 14.2857143,a = -1.1428571 e b = -0.7142857. Os autovalores de O são λ₁ =γ, λ₂ = σ+ iω, λ₃ = σ−iω comγ = 2.835,σ = −1.125 e ω=2.593. Comoσ γ < 0 e ω = 0,O é um ponto fixo do tipo sela-foco. A trajetória mostrada na Fig. 4 foi obtida atrav´ es daEq. (27).

Na Fig. 4, a trajetória no subespaço instável E_R(OM) parte deO para M₁. O pontoM₁ é a intersecção E^R(O) ∩ U1 (obtido fazendo x = 1 na Eq. (36)). Na regiãoD₁ a trajetória é uma curva que parte deM₁ para M2. O pontoM2 pertence à reta formada pela intersecção dos planosE^C(O)∩U₁. Finalmente, partindo de M₂, a trajetória descreve uma espiral que converge assintoticamente para O. Em síntese, a trajetória tende ao ponto fixo O quando t→ ∞ et → − ∞, ou seja, a trajetória é uma órbita homoclínica.

Em 1986, Chua, Komuro e Matsumoto apresentaram uma prova rigorosa de que é possível encontrar órbitas homoclínicas para o circuito de Chua-Matsumoto através da escolha adequada dos parâmetros α, β,a e b^[15].

3.5. Seção de Poincaré e mapa de primeiro retorno

Dado um sistema tridimensional, pode-se analisar os pontos de intersecção entre as trajetórias e uma seção transversal às mesmas (denominada de seção de Poincaré), em vez das próprias trajetórias. Vários pontos de interseção entre uma trajetória e a seção de Poincaré são obtidos ao longo do tempo, dando origem a uma versão bidimensional discreta do sistema contínuo tridimensional. Um mapa de primeiro retorno (ou mapa de Poincaré) transforma os pontos de uma seção de Poincaré nos pontos obtidos quando a trajetória atravessa novamente a seção no mesmo sentido. Uma análise bastante completa do circuito de Chua-Matsumoto utilizando mapas de primeiro retorno é apresentada na Ref. ^[15].

Em Fig. 5, é mostrada uma seção de Poincaré, em x = 1, para α = 11.0917459,β = 14.2857143, a = - 1.1428571 eb = -0.7142857. No tempo t₀, há 400 pontos distanciados de 0.0001, em torno do pontoM₁. Em um tempot₁, posterior, os pontos evoluem para pontos próximos aM₂. Nota-se que nesse caso, em que as trajetórias não passam pelo ponto homoclínico, a trajetória não apresenta sensibilidade às condições iniciais.

Figura 5
seção de Poincaré, em x = 1, para α = 11.0917459, β = 14.2857143, a = -1.1428571 e b = -0.7142857. No tempo t₀, há 400 pontos distanciados de 0.0001, em torno do ponto M₁. Em um tempo t₁, posterior, os pontos evoluem para pontos pr´oximos a M₂.

Na Fig. 6, é mostrada uma seção de Poincaré, em x = 1, para α = 11.0917459,β = 14.2857143, a = -1.1428571 eb = -0.7142857. A seção mostra dois instantes consecutivos de um mapa de primeiro retorno. Em t₀, há 400 pontos (20×20) distanciados entre si de 0.0001, em torno do pontoM₁. Em t₁, posterior, o sistema evoluiu para pontos distantes entre si. Nesse caso, em que a trajetória passa nas proximidades do ponto homoclínico O, a trajetória apresenta sensibilidade às condições iniciais.

Figura 6
Seção de Poincaré, em x = 1, para α = 11.0917459, β = 14.2857143, a = -1.1428571 e b = -0.7142857. A seção mostra dois instantes consecutivos de um mapa de primeiro retorno. Em t₀, há 400 pontos (20 × 20) distanciados entre si de 0.0001, em torno do ponto M₁. Em t₁, posterior, o sistema evoluiu para pontos distantes entre si.

3.6. Ferradura de Smale

O quadrado unitário D na Fig. 7 é transformado por f em um objeto em forma de ferradura (ferradura de Smale). A transformaçãof realiza a expansão, contração e dobra em três direções distintas, produzidas pela equação de evolução temporal dada por

Figura 7
Expansão, contração e dobra do quadrado unitário realizadas pela transformação f.

A transformação inversa f_-1 é mostrada na Fig. 8.

Figura 8
Expansão, contração e dobra do quadrado unitário realizadas pela transformação inversa f^-1.

A medida que ocorrem sucessivas iterações, a área de intersecção da ferradura comD reduz-se e divide-se em um maior número de regiões desconexas, como ilustrado na Fig. 9. Dá-se o mesmo quando sucessivas iterações são retrocedidas no tempo, como ilustrado naFig. 10.

Figura 9
Aplicação de f sobref(D), ou seja,f²(D).

Figura 10
Aplicação de f^-1sobref^-1(D), ou seja,f^-2(D).

Se todas as transformações f forem realizadas sobreD, obter-se-á infinitos segmentos de reta unitários verticais. Para transformações inversas f-1, obter-se-á infinitos segmentos de reta unitários horizontais. O conjunto das órbitas que permanecem em D após todas as transformações é obtido pelos pontos de intersecção entre os segmentos unitários horizontais e verticais, e é dado por

Em vez de analisar a evolução das órbitas, é possível analisar a evolução dos itinerários que informam em qual região deD(H₀ ouH₁) a órbita encontra-se a cada iteração. Esses itinerários podem ser representados por sequências bi-infinitas s = {...s₋₂s₋₁.s₀s₁s₂...}, com s_k = 0, se a órbita passa porH₀,e s_k = 1, se órbita passa por H₁.

Pode-se definir a transformação de desvio σ(s)= {...s₀.s₁s₂s₃...} em Σ (conjunto de todas as sequências bi-infinitas). A transformaçãoσ faz cada termo s_n da sequência avançar paras_n+1, representando a lei de evolução temporal em Σ, que corresponde a f em Λ. É possível provar que as conclusões obtidas para a dinâmica deσ em Σ também são válidas paraf em Λ.

Uma sequência bi-infinita s = {...s₋₁.s₀s₂...} possui uma órbita periódica, se existe T∈ ℕ tal que σ^T(s)= s, onde T é o período. Por exemplo, s = {...101101.101101...} = {101.101} é uma órbita periódica de período 3. Para qualquer número natural escrito na forma binária χ, é possível obter uma órbita periódica s = {χ∙χ}. Portanto, há uma infinidade contável de órbitas periódicas de período arbitrariamente grande.

A cada sequência bi-infinita s = {...s₋₁.s₀s₂...}corresponde uma sequência infinitas={0.s₀s₋₁s₂s₋₂...}. Se s é não-periódica, 0.s₀s₋₁s₂s₋₂... é um número irracional; portanto, há uma infinidade incontável de órbitas não-periódicas.

A distância entre duas sequências bi-infinitas s = {...s₋₂s₋₁.s₀s₁s₂...} e s' = {...s'₋₂s'₋₁.s'₀s₁s'₂...} é dada por d(s, s') = δ_k/2^|k|, ondeδ_k = 0, se s_k=s'_k , ou δ_k = 1, ses_k ≠ s'_k. Afirmar que s e s' estão próximos (dado um ε > 0, d(s, s') < ε), equivale a dizer: dado um N ∈ ℕ,s_k =s'_kpara |k| < N, ondeN = N (ε).

Uma sequência s possui órbita densa quando, dados s' ∈ Σ e ε > 0, existe um número de iterações T tal qued(σ^T(s),s') < ε. A sequênciaT de iterações, tem-se s_k=s'_k para |k| <N, onde s=σ^T (N=N(ε) ∈ ℕ. = {...{11}{01}{1}.{0}{00}{10}...} (formada por todas as sequências finitas com 0's e 1's) possui órbita densa, pois, após um certo número ) e

A transformação σ em Σ possui sensibilidade às condições iniciais. Dados s, s ∈ Σ e ε > 0, tais que d(s,s) < ε, ou seja,s_k = s_kpara |k| < N; apósNiterações, não será mais possível afirmar qued(σ^N(s),σ^N(s)) < ε, pois |k| >N, e não será possível garantir ques_k = s_k.

Como a dinâmica de f em Λ pode ser descrita em termos da dinâmica de σ em Σ, a transformação da ferradura de Smale possui:

• Uma infinidade contável de órbitas periódicas de período arbitrariamente grande;

• Uma infinidade incontável de órbitas não periódicas;

• Órbita densa;

• Sensibilidade às condições iniciais.

3.7. Teorema de Shilnikov

O teorema de Shilnikov foi publicado em 1965 ^[16]. Segundo esse teorema ^[15], dado um sistema tridimensional ṙ = h(r) contínuo e linear por partes que possui um ponto fixo do tipo sela-foco na origem com um autovalor real γ e um par complexo conjugado de autovaloresσ ± iω, com σ γ < 0 eω = 0; se |σ| < |γ| (condição de Shilnikov) e existe uma órbita homoclínica através da origem, então existe uma ferradura de Smale próxima a essa órbita homoclínica.

É possível encontrar um sistema de coordenadas em que uma órbita homoclínica, como aquela mos trada na Fig. 4, tenha o subespa¸ co instável ER(0) sobre o eixo z e o subespaço estável E_C(0) sobre o planoxy , como ilustrado na Fig. 11. Isso pode ser feito utilizando o conjunto dos autovetores como base, como em (42).

Figura 11
Ferradura de Smale próxima a uma trajetória homoclínica.

Um conjunto de pontos de uma seção de Poincaré tende a contrair-se para um único ponto, quando a seção é tomada cada vez mais próxima a um ponto fixo assintoticamente estável. Se em vez disso, o ponto fixo for instável, os pontos da seção afastarse-ão. Se o ponto fixo for homoclínico, a seção será contraída na direção do subespaço estável e será expandida na direção do subespaço instável. Além disso, se o o sistema for dissipativo, a seção terá de ser dobrada para que caiba no espaço de fases. Portanto, o retângulo (que é um conjunto de condições iniciais de várias trajetórias) na Fig. 4 transforma-se com as mesmas propriedades da ferradura de Smale. Nesse caso, o circuito de Chua-Matsumoto possuirá: uma infinidade contável de órbitas periódicas de período arbitrariamente grande, uma infinidade incontável de órbitas não-periódicas, órbita densa, e sensibilidade às condições iniciais.

3.8. Expoentes de Lyapunov

Pode-se pensar em um cubo centrado em um ponto de uma órbita homoclínica, que seja formado por pontos que representam condições iniciais. A medida que o tempo passa, os pontos desse cubo evoluem de uma forma que, a dimensão do cubo na direção da órbita deve ser mantida constante, pois caso diminuísse, convergiria para um ponto, e caso aumentasse, os pontos afastar-se-iam uns dos outros e da órbita, com o passar do tempo.

Considerando que as dimensões do cubo sãod_j(t)=d_j(t₀)e^Λj(t−t0), onde j = 1, 2, 3 e Λ₁,Λ₂ eΛ₃são os expoentes de Lyapunov. Por exemplo, sed₂(t) é a dimensão na direção de uma órbita homoclínica, então Λ₂ = 0, ed₂(t) =d₂(t₀) não varia com o tempo (o expoente na direção da trajetória deve ser nulo, pois as trajetórias nunca se afastam do atrator caótico, nem tendem a um ponto fixo). Sed₁(t) é a dimensão na direção do subespaço instável, então Λ₁ > 0, ed₁(t) =d₁(t₀)e^Λ1(t−t0)cresce exponencialmente com o tempo. Sed₃(t) é a dimensão na direção do subespaço estável, então Λ₃ < 0, ed₃(t) =d₃(t₀)e^Λ3(t−t0)diminui exponencialmente com o tempo. O volume do cubo é

Um sistema caótico é dissipativo, pois as trajetórias não se afastam do atrator caótico, portanto,

4. Aplicativo Circuit of Chua and Chaos

Nessa seção serão descritos os métodos numéricos implementados no aplicativo. Além disso, serão apresentados alguns resultados obtidos através de simulações com o programa Circuit of Chua and Chaos, realizadas em umSmartphone. Esses resultados serão comparados aos obtidos com simulações utilizando o programa MATDS (para MATLAB),¹ realizadas em um computador de mesa. As liguagens utilizadas no desenvolvimento do aplicativo foram XMLeJava. Todas as telas do aplicativo foram desenvolvidas emXML, e toda a lógica de programação foi implementada emJavainclusive os algoritmos do método Runge-Kutta-Fehlberg e da técnica empregada no cálculo dos expoentes de Lyapunov.

É possível fazer simulações sem realizar a entrada de dados, pois os valoresα = 15.0, β; = 25.58,a= -1.14, b = -0.71x₀ =1.6,y₀ = 0.0,z₀ = −1.5 et_max = 45.0,h_max = 0.01,h_min = 1.0E − 16 etol = 1.0E − 15 estão previamente selecionados. A escolha de outros valores podem ser realizadas nas telas de entradas de dados. Na Fig. 12 é mostrada a tela de entrada de dados para os parâmetros do circuito de Chua. Há outras telas para escolher a condição inicial e os parâmetros do método RungeKutta-Fehlberg.

Figura 12
Tela para entrada de dados em que os parâmetros do circuito de Chua-Matsumoto (α,β,a e b) podem ser escolhidos.

4.1. Simulação computacional

O método de Runge-Kutta-Fehlberg foi realizado a partir da generalização para o caso tridimensional de um algoritmo (para o caso unidimensional) apre sentado na Ref. ^[17](algoritmo 5.3 da p. 275). O cálculo dos expoentes de Lyapunov foi realizado com base no algoritmo apresentado na Ref. ^[18]. Ambos os algoritmos tiveram de ser adaptados ao caso do circuito de Chua-Matsumoto.

O exemplo de simulação apresentado aqui utiliza os valores previamente selecionados no aplicativo, isto é

Os parâmetros do método Runge-Kutta-Fehlbergt_max,h_max,h_minetol representam, respectivamente, o tempo máximo de simulação, o passo temporal máximo, o passo temporal mínimo e a tolerância do método.

A tela utilizada para plotar os gráfico da projeção da trajetória no planoxy está mostrada na Fig. 13, para os parâmetros mostrados na Eq. (47). Pode-se escolher também o planoyz ou xz. A tela utilizada para plotar os gráficos das curvas x(t) está mostrada naFig. 14. Pode-se plotar tambémy(t) ez(t) em função do tempo.

Figura 13
Tela de plotagem das projeções da trajetória. Foi plotada a projeção da trajetória no plano xy, para os parâmetros mostrados na Eq. (47).

Figura 14
Tela de plotagem das curvas em função de t. Foi plotado o gráfico da curva x(t), para os parâmetros mostrados na Eq. (47).

Na Fig. 15 é mostrada a tela de cálculo dos expoentes de Lyapunov, a qual apresenta os valores Λ₁ = 0.36, Λ₂ = 0.00, e Λ₃ = −4.15, para os parâmetros mostrados na Eq. (47). A soma dos expoentes é negativa, indicando que o sistema é dissipativo. Como visto anteriormente, esses valores correspondem a uma solução caótica. Utilizando o aplicativo MATDS, foram obtidos os valores Λ₁ = 0.39, Λ₂ = 0.00, e Λ₃ = −4.12. Vê-se que esses resultados estão de acordo.

Figura 15
Tela de cálculo dos expoentes de Lyapunov, para os parâmetros mostrados na Eq. 47.

Os pontos de equilíbrio, para os parâmetros da Eq. (47), são P₊ = (1.48, 0, −1.48), P₀ = (0, 0, 0) eP₋ = (−1.48, 0, 1.48). Os autovalores deP₀ são λ₁ = −1.1259 + 3.8418i, λ₂ = −1.1259 − 3.8418ie λ₃ =3.3517. Os autovalores deP₊ eP_- são λ₁ =0.31038 + 4.30581i, λ₂ =0.31038 − 4.30581i e λ₃ = −5.9708. Nota-se que todos os pontos de equilíbrio são do tipo sela-foco e obedecem à condição de Shilnikov.

Na Fig. 16 é mostrada a trajetória do atrator de rolo duplo (double-scroll) para o circuito de Chua-Matsumoto com os parâmetros mostrados na Eq. (47). Os resultados obtidos com o aplicativoCircuit of Chua and Chaos foram comparados aos obtidos com o programa MATDS. Os resultados são coincidentes.

Figura 16
Atrator de rolo duplo (double-scroll) para o circuito de Chua-Matsumoto com os parâmetros mostrados na Eq. (47). Os resultados obtidos com o aplicativo Circuit of Chua and Chaos são comparados aos obtidos com o programa MATDS. Os resultados são coincidentes.

Nas Figs. 17, 18 e 19 são mostrados, respectivamente, os gráficos das curvas dex(t),y(t)ez(t) em função de t. Todos os resultados estão em concordância.

Figura 17
Curvas de x em função de t obtidos com o aplicativo Circuit of Chua and Chaos (linhas) comparadas às curvas obtidas com o programaMATDS(pontos na forma de x), os parâmetros mostrados naEq. (47). Os resultados são coincidentes.

Figura 18
Curvas de y em função de t obtidos com o aplicativo Circuit of Chua and Chaos (linhas) comparadas às curvas obtidas com o programaMATDS(pontos na forma de x), os parâmetros mostrados naEq. (47). Os resultados são coincidentes.

Figura 19
Curvas de z em função detobtidos com o aplicativo Circuit of Chua and Chaos(linhas) comparadas às curvas obtidas com o programaMATDS (pontos na forma dex), os parâmetros mostrados na Eq. (47). Os resultados são coincidentes.

5. Conclusão

Realizar simulações com o circuito de Chua-Matsumoto é uma forma eficiente de tornar mais atrativo o estudo de sistemas dinâmicos (uma disciplina que em geral é ofertada tanto para os cursos de engenharia quanto para os cursos de física e matemática), auxiliando no aprendizado de conceitos relacionados à disciplina. O uso da plataforma Android possui algumas vantagens em termos da realização de simulações com relação às outras tecnologias, tais como: o número de usuários ativos do Android ultrapassou a quantidade de um bilhão de usuários em 2014; ² o custo dos dispositivos é menor; os dispositivos possuem tela sensível ao toque; os dispositivos apresentam maior comodidade de manuseio (ocupam pouco espaço, não apresentam superaquecimento, apresentam um menor consumo de energia elétrica).

O aplicativo Circuit of Chua and Chaos foi desenvolvido para o estudo de um sistema que apresenta um certo grau de complexidade, uma vez que as soluções podem ser caóticas. O aplicativo mostrouse eficiente na obtenção de expoentes de Lyapunov; projeções das trajetórias nos planosxy,xz e yz; e gráficos dex(t),y(t)ez(t). Os resultados do aplicativo foram exportados para aquivos txt e comparados com os resultados obtidos com o programa MATDS, utilizado por muitos dos pesquisadores que investigam o circuito de Chua-Matsumoto, inclusive o próprio Leon Chua, que dá nome ao circuito, como pode ser visto na Ref.^[9]. A partir da compara¸ cão, verificou-se que os resultados obtidos com o aplicativo coadunam com os obtidos com o MATDS.

Pretende-se desenvolver outros programas de simulação para plataformaAndroid, de forma a cobrir o conteúdo da disciplina de sistemas dinâmicos que normalmente é ministrado em um curso regular de graduação.

Referências

• ^[1]
  Luiz Henrique Alves Monteiro, Sistemas Dinâmicos(Ed. Livraria da Física, São Paulo, 2011), 3a ed.
• ^[2]
  Stephen Lynch, Dynamical Systems with Applications using MATLAB(r) (Birkhäuser, Basileia, 2014).
• ^[3]
  Ahmed G. Radwan and Mohammed E. Fouda , On the Mathematical Modeling of Memristor, Memcapacitor, and Meminductor(Springer International Publishing, Suiça 2015).
• ^[4]
  T. Matsumoto, IEEE Trans. on Circ. & Syst. 31,1055 (1984).
• ^[5]
  T. Matsumoto, L.O. Chua and M. Komuro, IEEE Trans. on Circ. & Syst. 32, 797 (1985).
• ^[6]
  T. Matsumoto, L.O. Chua and M. Komuro, Journal Circuit Theory and Appl. 14, 117 (1986).
• ^[7]
  M.P. Kennedy, Frequenz 46, 66 (1992).
• ^[8]
  J.M. Cruz and L.O. Chua, IEEE Trans. on Circ. & Syst.46, 985 (1992).
• ^[9]
  M. Itoh and L.O. Chua, Int. J. Bif. and Chaos 18,3183 (2008).
• ^[10]
  B. Muthuswamy, Int. J. Bif. and Chaos 453, 1335 (2010).
• ^[11]
  B. Muthuswamy and L.O. Chua, Int. J. Bif. and Chaos20,1567 (2009).
• ^[12]
  L.O. Chua, IEEE Trans. on Circ. Theory 5, 507 (1971).
• ^[13]
  D.B. Strukov, G.S. Snider, D.R. Stewart and S.R. Williams, Nature453, 80 (2008).
• ^[14]
  A. Buscarino, L. Fortuna, M. Frasca and L.V. Gambuzza, Chaos: Interdiscip. J. Nonlinear Sci. 22, 023136 (2012).
• ^[15]
  L.O. Chua, M. Komuro and T. Matsumoto, IEEE Trans. on Circ. and Systems 33, 1072 (1986).
• ^[16]
  L.P. Shil'nikov, Dokl. Akad. Nauk SSSR 160, 558 (1965).
• ^[17]
  Richard L. Burden e J. Douglas Faires, Análise Numérica (CENGAGE Learning, São Paulo, 2008) 8a ed.
• ^[18]
  A. Wolf, J.B. Swift, H.L. Swinney and J.A. Vastano, Physica D16, 285 (1985).

• 1
  ver http://www.math.rsu.ru/mexmat/kvm/matds/.
• 2
  ver http://www.techspot.com/news/

Datas de Publicação

• Publicação nesta coleção
  Mar 2016

Histórico

• Recebido
  27 Jul 2015
• Aceito
  28 Nov 2015