Caos no Sistema de Três Corpos

Caos no Sistema de Três Corpos Geferson Lucatelli Universidade Federal do Rio Grande - FURG Instituto de Matemática, Estatística e Física 26 de março de 2015

1

Introdução

O estudo de sistemas não-lineares e o Caos vem cada vez mais sendo estudados. Uma área que despertou muito interesse para pesquisadores. A partir desses estudos, notou-se que muito dos fenômenos que nos rodeiam, se comportam de modo não linar e também caótico, e muito poucos sistemas são não-lineares. Um simples ocilador, que tanto é estudado, é um sistema nãolinear, e podemos torná-lo caótico, e mais, nós não resolvemos a equação que o governa. O que fazemos é uma aproximação para sin θ ≈ θ, para resolvê-lo de forma linear. O estudo desses sistemas só foi possível hoje em dia por computadores. Neste presente trabalho, o estudo a ser feito é o sistema de três corpos, sendo constituído por duas massas e uma partícula teste, onde essa ficará evoluindo sob interação das anteriores. Será feito o cálculo do expoente de Lyapunov, que descreve o quão as órbitas vão se afastando, devido a uma muito pequena modificação nas condições iniciais da partícula teste.

2

Física do Sistema

O sistema descrito aqui, são de três corpos. Dois corpos de massa m1 e m2 que são considerados como não pontuais portanto possuem um raio, respectivamente, R1 e R2 . Esses dois corpos estão fixos em dois pontos num espaço de duas dimensões, digamos, xy. Com o sistema sendo caótico, o movimento ficará sempre confinado nesse plano, pois a componente z do momento será conservada. Mas o momento angular não se conserva nas componentes x e y. O outro corpo, é a partícula teste, que podemos considerar ter uma massa muito pequena, em efeito as outras, e designo por m. 1

A equação que rege o movimento da interação entre corpos com massa m é dada pela Lei de Newton da Gravitação, F=−

Gm1 m2 ˆr |r12 |3

(1)

onde o vetor ˆr dá a direção em que a força está atuando. Para um sistema de três corpos, precisamos decompor está última equação para a partícula teste interagindo com cada uma das massas. Essa decomposição é feita em relação a uma massa, seja m1 , e depois em relação à outra, m2 . Primeiro vamos fazer com a massa m1 . Concidere um ângulo θ1 que m faz com m1 em relação a uma linha que liga à elas com qualquer eixo paralelo ao eixo y. Então a força Fy atuando em m é Fy = F cos θ1 , e na direção x, é Fx = F sin θ1 , onde F é o módulo da força dada pela equação (1). Vamos definir algumas variáveis, x e y são as coordenadas da massa m teste em relação á origem, e suas respectivas velocidades nas componentes são vx e vy . A posição do corpo de massa m1 é fixada em uma coordenada, x1 e y1 , no qual, ficará constante. As duas massas não estão orbitando entre si neste problema. Bom, aplicando o teorema de Pitágoras, temos a distância p de m até m1 , que é r1 = (x − x1 )2 + (y − y1 )2 . Note que temos que fazer a diferença entre a posição de m1 e de m com a origem, por isso (x − x1 ) e (y − y1 ). Vamos também reescrever o seno e o cosseno em coordenadas cartesianas, pois, mesmo não parecendo, é melhor de se trabalhar. Temos então x − x1 y − y1 sin θ1 = (2) cos θ1 = r1 r1 A partir da equação (1), escrevemos o sistema de equações, nas componentes Gmm1 Gmm1 x − x1 d2 x cos θ1 = − (3) m 2 =− 2 dt r1 r12 r1 m

d2 y Gmm1 Gmm1 y − y1 = − sin θ = − 1 2 dt2 r1 r12 r1

(4)

Esse sistema é somente a relação entre m1 e m. Para fazermos isso com m2 o processo é exatamente o mesmo, de forma a obtermos m

d2 x Gmm2 Gmm2 x − x2 =− cos θ2 = − 2 2 dt r2 r22 r2

(5)

d2 y Gmm2 Gmm2 y − y2 =− sin θ2 = − 2 2 dt r2 r22 r2

(6)

m

2

onde θ2 é o ângulopentre a linha que liga m à m2 com qualquer eixo paralelo ao eixo y, e r2 = (x − x2 )2 + (y − y2 )2 é a distância entre elas. A soma total das forças que está interagindo com m é a soma das componentes em x, equações (3) com (5), e em y, (4) com (6). O fator comum que aparece em todas as equações é m e pode ser eliminado. Então, a partir disso, chegamos as equações que descrevem o problema de três corpos,     Gm2 d2 x Gm1 (7) = − 3 (x − x1 ) + − 3 (x − x2 ) dt2 r1 r2 e

d2 y = dt2

    Gm2 Gm1 − 3 (y − y1 ) + − 3 (y − y2 ) r1 r2

(8)

Antes de prosseguirmos, não devemos esquecer, de que estamos considerando m1 e m2 como não pontuais, e de fato, obrigatoriamente isso deve ser feito, pois senão, quando a massa teste se aproximar de um dos corpos, a distância seria muito pequena, pois seria pontual, e a força iria crescer muito. Devemos fazer uma modificação em r1 e r2 , que é exatamente adicionar o raio R1 de m1 e R2 de m2 . Isso é simples, p p r2 = (x − x2 )2 + (y − y2 )2 + R2 r1 = (x − x1 )2 + (y − y1 )2 + R1 É interesante discutir as equações (7) e (8), para saber em qual direção a força estará atuando para uma dada posição x. Podemos escolher apenas um termo dessas equações. Pode ser o primeiro da primeira equação. Se a distância x for maior que x1 , que implica m estar a direita de m1 num eixo paralelo à x, o termo (x − x1 ) é positivo, mas como existe um sinal negativo na frente da expressão, a força aponta no sentido −x, que é exatamente o que esperamos que acontece-se. A partir disso, todas as componentes estão corretas para um dado y e x, portanto as equações são válidas.

3

Integrador Numérico

O integrador numérico que foi utilizado para resolver este problema, é o método de Leap-Frog, que foi introduzido em código Python. Este método consiste em calcular posições e velocidade em tempos intercalados. Este método foi escolhido, pois reduz o número de cálculos. Implementar esse código em Runge-Kutta de quarta ordem não seria viável, iria consumir muita memória computacional, devido ao grande número de equações, e visto que não traria muitos benefícios em relação ao erro do integrador.

3

As equações do método de Leap-Frog são, de uma certa forma 1 xi+1 = xi + vxi δt + axi δt2 2

(9)

para a posição, num instante i + 1 e 1 (10) vxi+1 = vxi + (axi + axi+1 )δt 2 e para a velocidade num instante i + 1. Note as acelerações intercaladas na equação da velocidade. A aceleração axi+1 é calculada a partir da posição naquele instante, e a velocidade é então calculada com as duas acelerações. Essas acelerações são os termos do lado direito das equações (7) e (8). Esses cálculos para a componente y são executados da mesma forma que as duas equações anteriores. As equações até agora foram escritas nas unidades do SI. Mas podemos simplificar isso, para diminuir o esforço computacional com os cálculos. As massas m1 e m2 podem ser definidas como sendo menor que 1, e da mesma forma para a distância entre elas e a massa m, como sendo algumas unidades de distância. A escala temporal não necessariamente precisa ter unidades, na verdade nenhuma dessas quantidades precisa. Isso não é muito relevante, desde que é usado a proporção certa entre essas quantidades. Por exemplo para a velocidade com que m translada, os outro corpos pode ser escrita como sendo 2π, que vem de uma órbita circular de raio R unitário 1, e leva 1 unidade tempo para completar uma volta em torno de uma massa M unitária 1, então v = 2πR/T = 2π. Dessa forma, para encontrar G igualamos a força de atração, equação (1) com a resultante centrípeta F = mv 2 /R, √ mv 2 GM m = −→ v = G 2 R R

(11)

como v = 2π, G = 4π 2 .

4

O Sistema de dois Corpos

Vamos agora por em prática o que foi proposto nas equações acima. É interessante em primeiro caso analizar uma condição dessa interação que não seja caótica, para mostrar como isso é limitado e é pouco comum. Esses tipos de sistemas estão a deriva de entrarem em caos. Consideremos primeiro para um teste m1 = 0.0 e m2 = 0.009. E também o raio de m2 como sendo R2 = 0.06. Façamos isso para analizar uma órbita estável. Para as coordenadas de m2 , temos x2 = −1.5 e y2 = 1.5. A massa m, terá suas 4

condições iniciais como sendo x = 16.0, y = 0.0, vx = 0.0 e vy = 0.1. A Figura (1) mostra a trajetória da massa teste durante 500 intervalos de tempo, com um passo do integrador de δt = 0.001. O ponto no interior da elipse é o corpo central. Esta órbita é eliptica, fechada, estável, com um

Figura 1: Órbita estável, fechada para as condições iniciais especificadas. período de translação de ≈ 415.0 unidades de tempo. Não tem movimento de presseção. Esta é uma órbita, sob uma boa aproximação, típica para o geral e a maioria dos corpos que orbitam o Sol e qualquer estrela. Disse sob boa aproximação, porque estamos desconsiderando a interação com outros corpos, e também mesmo que se esses corpos existisem, eles não têm quase nenhuma influência sob os demais, devido as grandes distâncias que os separam. Talvez no passado o Sistema Solar tenha sido caótico por conta de existir uma mais forte interação entre os corpos. Vamos ver agora uma órbita fechada, estável e com dois períodos, um de translação e um perído de presseção, que é a variação do ângulo com que a órbita gira depois de uma volta tranladada por m. Para isso, mantemos todas as condições iniciais iguais que anteriormente, exceto x. Definimos x = 6.0. Veja a figura (2) abaixo onde mostra esse movimento.

5

Figura 2: Órbita estável, aberta, movimento de presseção. Período duplo, Tθ , e Tr Nessa figura vemos que a órbita não se repete depois de uma translação completa de m. A órbita está tendo um movimento de presseção, e ela tem portanto dois períodos, o de translação de m e o de presseção. Essa órbita com período duplo voltará a se repetir se os períodos são comensuráveis, ou seja, se a razão Tθ /Tr pode ser expressa como uma razão de números inteiros. Uma outra maneira de dizer se a órbita é fechada é analizar a razão δθ a = (12) 2π b onde a e b são números inteiros, de maneira que a razão entre a variação de θ sobre um ciclo de 2π sempre possa ser expressa por uma razão de números inteiros. A variação do ângulo depois de uma volta completada é de aproximadamente δθ ≈ 0.14 radianos.

5

O Sistema de Três Corpos

Vamos agora partir a caminho do chaos, e observar a partícula m evoluindo em torno de m1 e m2 . Inicialmente, vamos encontrar um regime que não seja caótico, e depois, ajustar os parâmetros, para que o sistema fique caótico. Vamos manter as posições de m1 e m2 com estavam antes. O que vai se 6

modificado é a massa de m1 = 0.009, idêntica à m2 . As condições iniciais para m agora é, x = 6.0, y = 0.0, vx = 0.0 e vy = 0.4. Para essas condições, a órbita é estável, veja a figura (3) .

Figura 3: Órbita estável para o sistema de três corpos, o tempo de simulação é de 1000 intervalos de tempo. Notamos que esta órbita também tem um movimento de presseção, contudo ela é estável a não entra em caos. Mas até onde isso fica estável? O que poderia acontecer neste sistema para deixa-lo caótico. Podemos dar uma resposta rápida - pouca coisa é nescessária. Conforme vamos aumentando a complexibilidade do problema, mais interações sobre m, o mesmo começa a adiquirir uma cara de ser instável. Isso pode ser explicado pelo fato de que a equação (1), a equação Newtoniana que governa o sistema, não ser linear. Em condições especiais, onde há pouca interação entre os corpos, ela se comporta de maneira linear, da mesma forma que fazemos para um pêndulo, sin θ ≈ θ, onde para pequenas amplitudes, sin θ é linear. Conforme veremos, um ponto forte do caos, é ser sensível as condições iniciais. E agora que estamos mudando parâmetros, vamos ver como ele começa a surgir. Para que o sistema se torne caótico, vamos fazer uma pequena modificação em vy e em x, sendo agora iguais à 0.3 e 5.0 respectivamente. Mas, surge uma pergunta, o porque dessas condições. Veja, como foi afirmado, para pequenas interaçãoes, a equação da força é linear, mas veja que a posição inicial de m está mais próxima de uma das massas. E também, não podemos colocar m 7

muito perto de m1 ou m2 pois senão ela seria jogada para longe. Mas claro, para escolher essas condições, foram feitas várias simulações. Veja a Figura (4) com o resultado dessa simulação.

Figura 4: Órbita instável, caótica para condições iniciais vy = 0.3 e x = 5.0, com 1000 unidades de tempo de simulação Veja, os parâmetros definidos fizeram o sistema entrar em caos (não é possível afirmar corretamente isso, ainda). Isso mostra que esse tipo de problema, é comumente caótico. A maioria dos problemas físicos envolvem equações não-lineares. As equações lineares são casos extremamete partículares. Mas cuidado, para o caos surgir de uma equação não-linear, não é simplesmente sob quaqlquer parâmetro. É preciso determinar esses parâmetros para que surja caos. Então pode-se afirmar com toda a certeza, essas equações, sempre tem um caminho a se seguir que nos levam ao caos. Vamos usar outras condições, que levam a um sistema caótico, e deixá-lo evouluir mais tempo. Seja elas então m1 = m2 = 0.005, R1 = R2 = 0.1, x = 5.23, y = 0.0, vx = 0.0 e vy = 0.2. Vamos plotar o diagrama de faze do sistema, e também a amplitude, pode ser em x, em função do tempo. O diagrama de faze de x versus vx é mostrado na Figura (5), e de y versus vy é mostrado na Figura (6). A amplitude x em função de t é mostrado na Figura (7).

8

Figura 5: Diagrama de faze do sistema de x por vx de para as condições caóticas.

Figura 6: Diagrama de faze do sistema de y por vy para as condições caóticas.

9

Figura 7: Amplitude pelo tempo, do mesmo sistema das duas figuras anteriores. É muito interessante analizar o diagrama de faze de um sistema caótico, pois a partícula de massa m nunca vai passar na mesma posição com a mesma velocidade, e este plano vai ficar totalmente preenchido. Mas notamos que ele mantem uma forma, e isso é fundamental, pois nos diz que o sistema está conservando a energia. No gráfico do deslocamento em x, observa-se que a função que o descreve não é harmônica e não tem períodos pois é caótica. De modo grosseiro, podemos dizer que existe períodos sim, mas um número muito grande, talvez infinitos períodos, e portanto não há maneira de escrever essa função numa forma de uma soma de funções elementares finitas. Vejamos em seguida a característica mais fundamental do caos.

6

Sensibilidade as Condições Iniciais

Quão longe duas trajetórias podem ficar separadas por uma pequena distância, |δx(t)| devido a uma pequena variação |δx| na condição inicial? Em outras palavras, temos certeza se em um tempo t, ou depois de um número n de interações, as trajetórias estarão ainda próximas uma da outra?, ou serão zero? No caos, um sistema é super instável, e qualquer mudança, mesmo sendo pequena, mas muito pequena (até menor que 10−15 ), traz condições 10

extremamente diferentes entre duas órbitas que saíram nessa diferença. Vamos agora criar dois sistemas idênticos, ou melhor, duas partículas, m01 e m02 idênticas. As suas coordenadas serão representadas por x01 , y10 , x02 e y20 , 0 0 0 0 , respectivamente. Para as condições, e vy2 , vy1 , vy1 e suas velocidades por vx1 vamos usar as mesmas que anteriormente, m01 terá para a posição (5.23, 0) e velocidade (0, 0.2), e m02 , posição (5.23, 0), velocidade (0, 0.2). Porém, vamos fazer uma modificação na coordenada inicial de m02 , um incremento na posição x02 de um fator δx = 0.00001, e deixar o sistema evoluir. O que podemos esperar disso? Como havíamos comentado, quando um sistema está a beira do caos, qualquer mudança gera consequências totalmente diferentes. Vamos ao ponto, a figura (8) mostra as duas órbitas. A Figura (9), mostra a amplitude x em função de t, das duas massas, para verificarmos a diferença de caminho de cada uma.

Figura 8: Órbitas de m01 (em verde) e m02 (em vermelho), com uma diferença inicial na posição x de δx = 0.00001 Note que as órbitas, depois de um tempo t, são totalmente diferentes uma da outra, Figura (9), apesar de ter apenas uma mísera diferença em apenas uma condição inicial, isso é incrível. Essa análize é uma ótima ferramenta para descobrir se um sistema é caótico. Sim, podemos afirmar com toda a certeza que o nosso sistema é caótico. Uma característica, talvez a mais fote, do caos, é exatamente isso, mudanças inicias que parecem insignificantes, geram trajetórias totalmente diferentes. 11

Figura 9: Amplitude de m01 (em verde) e m02 (em vermelho), com uma diferença inicial na posição x de δx = 0.00001. Note que a variação brusca ocorre no instante t≈ 482, e parece voltar a ser igual, mas diverge totalmente em t ≈ 680.

12

Mas vamos além, queremos "testar"mais o nosso sistema para ver se ele é caótico. Utilizando as mesmas condições iniciais, com um δx menor ainda, δx = 0.0000000001 = 1−10 , e vejamos o que ocorre, Figura (10) Note que as órbitas demoram um pouco mais para ficarem diferentes, mas mesmo assim ficam e muito, isso é caos. Concluimos que as trajetórias se separam de maneira exponencial, na próxima seção será apresentado isso.

Figura 10: O desvio brusco da órbita ocorre próximo a t ≈ 1065 . Por curiosidade, podemos fazer essa análize, da diferença das condições iniciais mara o primeiro problema, da Seção (4), onde teríamos chegado a conclusão de que o sistema não era caótico, e agora podemos testar isso. Utilizemos δx = 0.00001, e comparamos as duas trajetórias. Veja abaixo, na Figura (11). Vimos o que esperávamos, as órbitas não separam, pois o sistema não é caótico. Vamos descrever agora, de que maneira as duas órbitas se separam.

13

Figura 11: Amplitude na direção x em função do tempo das duas massas, sendo um diferenciando do outro por um δx. Como as tragetórias não divergem, o sistema é não-caótico.

6.1

Expoente de Lyapunov

Uma maneira de quantificar isso, que é o caos determinístico, é estudar o expoente de Lyapunov. A definição disso diz que num sistema caótico, as trajetórias se afastam de forma exponencial. Esse afastamento obedece a uma expressão da forma |δx(t)| = δxeλt (13) onde δx é a diferença entre a condição inicial dos dois sistemas, e λ é o chamado de expoente de Lyapunov. Se dividirmos ambos os lados da equação (13) e aplicarmos o logaritmo natural da mesma forma, temos ln |

1 δx(t) δx(t) |= ln(eλt ) −→ λ = ln | | δx t δx

(14)

A diferença δx(t) é a diferença entre as órbitas, x01 −x02 no instante t. Podemos reescrever a equação (14) como uma soma de todas as diferenças, divididas pela diferença inicial δx, portanto t

1X δxi ln | | t→∞ t δx i=0

λ = lim

14

(15)

Se λ for negativo, as duas órbitas não divergem, ou seja , não é caótico. Se λ é maior que zero, a diferença cresce exponencialmente, como está na equação (13). A Figura (12) mostra o gráfico de λ em função de t.

Figura 12: Veja que inicialmente, λ = ln |x01 (t) − x02 (t)|, é negativo, ou seja as órbitas estão próximas. Mas com o tempo, λ cresce, e fica positivo, onde o sistema diverge. O gráfico nos mostra o que queríamos, a diferença entre as órbitas cresce de maneira exponencial. Isso é mostrado no intervalo de t = 0 até t ≈ 1000, onde λ cresce de forma linear, pois é o logarítmo natural de uma exponencial, que é uma reta. Depois disso, a diferença se estabiliza, pois como vimos em figuras anteriores, o sistema fica confinado de modo caótico em torno dos dois corpos, m1 e m2 , e não há como crescer mais. Veja uma última figura, onde é mostrado λ em função do tempo para o problema de três corpos que não é caótico. Esperamos que λ seja muito negativo, o que realmente é, veja a Figura (13). Nota-se que o pequeno aumento da diferença entre as órbitas se deve ao erro do integrador numérico. Analizar a diferença entre as órbitas em sistemas é muito importante, pois é possível, desta maneira, do expoente de Lyapunov, quantificar o que se está estudando.

7

Conclusões

Fortes conclusões que podemos tirara partir de um breve estudo sobre o caos, é que ele pode estar presente em muitos fenômenos físicos. O estudo sobre 15

Figura 13: A diferença, λ = ln |x01 (t) − x02 (t)| é mínima e estabiliza, pois o sistema não é caótico e as órbitas não se separam. o caos despertou muito interesser em físicos, matemáticos e biólogos a partir da década de 70, pois o domínio computacional entrou em área, resolvendo problemas que antes eram muito difíceis, ou senão impossíveis de se resolverem. O caos reafirma a seguinte situação, de que não podemos determinar nada com exatidão, pois uma pequena mudança em qualquer parâmtro muda completamente a situação do que estamos lidando, em comparação com outro parâmetro. Depois que um sistema caótico evolui, nunca podemos saber de onde eele partiu, e também o quanto tempo ele está evoluindo. Ressaltando a ideia de um pêndulo caótico. Imagine que queiramos construir um pêndulo onde podemos fazer um estudo sobre caos a partir dele. Mas, em paralelo, resolvemos a sua equação computacionalmente, para comaparar a prática e a teoria. Podemos dizer que é muito difícil de se obter os mesmos resultados, pois como vimos , diferenças de escalas de nano unidades de medida, alteram totalmente a evolução de um sitema, como na diferença das posições das duas massas, δx = 1×−10 . Então, no momento de ajustar experimentalmente os parâmetros do pêndulo construído, não se tem possibilidade de certezas de suas condições iniciais, ou seja, se são exatamente as mesmas que a teoria no computador.

16

8

Referências

[1] Acheson D., From Calculus to Chaos: An Introduction to Dynamics, Editora Oxford University Press, 1997. [2] Taylor R. J., Mecânica Clássica, 2 ed., Editora Bookman, 2013. [3] Thornton S.T. et al, Classical Dynamics of Particles and Systems, 5 ed, Editora Thomson. [4] Yorke J.A. et al, CHAOS: An Introduction to Dynamical Systems, Editora Springer, 2000.

17