Coreografias no Problema de Três Corpos Restrito

Choreographies on the restricted three-body problem

Resumo

O problema de três corpos trata de objetos puntiformes interagindo mutuamente através da força gravitacional de Newton. Ao longo de mais de três séculos, o estudo deste tipo de sistema levou ao desenvolvimento e ao aprimoramento de diversas técnicas matemáticas, tanto analíticas quanto numéricas, para a compreensão de problemas que envolvem sistemas dinâmicos. Este trabalho discute alguns desses resultados aplicados ao problema de três corpos restrito, formulado a partir da segunda lei de Newton e da lei da gravitação universal. Em particular, foi estudado o comportamento e a estabilidade de dois tipos importantes de soluções periódicas dessse problema: o alinhamento em linha reta de L. Euler e a coreografia em forma de oito de C. Moore. O software Mathematica foi utilizado para resolver o sistema dinâmico e gerar as imagens de movimento dos corpos, bem como calcular o conjunto de expoentes de Lyapunov associados a cada solução. Apesar do caráter inerentemente caótico do problema de três corpos observado nos expoentes de Lyapunov, a solução em oito é linearmente estável, como discutido nos trabalhos de C. Simó.

Palavras-chave:
Problema de Três Corpos; Sistemas Dinâmicos; Expoentes de Lyapunov

Abstract

The three-body problem deals with point objects interacting mutually through Newton's gravitational force. For more than three centuries, the study of this type of system has led to the development and improvement of several mathematical techniques, both analytical and numerical, for the understanding of problems involving dynamical systems. This paper discusses some of these results applied to the restricted three-body problem formulated from Newton's second law and the law of universal gravitation. In particular, the behavior and stability of two important types of periodic solutions of this problem were studied: L. Euler's straight-line alignment and C. Moore's figure-eight choreography. Mathematica software was used to solve the dynamical system and to generate the images of movement of the bodies, as well as to calculate the set of Lyapunov exponents associated with each solution. Despite the inherently chaotic character of the three-body problem observed in the Lyapunov exponents, the figure-eight solution is linearly stable, as discussed in C. Simó's works.

Keywords:
Three-body Problem; Dynamical Systems; Lyapunov Exponents

1. Introdução

O problema de três corpos trata, originalmente, de objetos celestes interagindo mutuamente através da força gravitacional e data de 1687, como uma extensão do problema de dois corpos resolvido por Sir Isaac Newton ^[1]. Apesar do nome ser relativamente recente, esse problema foi explorado desde a antiguidade através das observações do Sol, da Terra e da Lua, na busca pela compreensão e previsão de eclipses ^[2]. Ao longo da história, a sua resolução obteve a contribuição de vários físicos e matemáticos famosos e levou a avanços importantes, desde o entendimento da natureza dos corpos celestes até a programação e implementação de métodos numéricos computacionais, os quais também contribuíram com outras áreas do conhecimento.

A formulação do problema parte da lei da Gravitação Universal e da segunda Lei de Newton. O sistema resultante apresenta seis equações por corpo, três para a posição e três para a velocidade, totalizando 18 equações. Apesar de obtidas de modo simples, elas são equações diferenciais não-lineares e acopladas, o que dificulta sua resolução. Este tipo de conjunto de equações é chamado sistema dinâmico [³ -⁵], cuja solução descreve o comportamento de um ponto, a condição inicial, com o passar do tempo. Um modo de simplificar este sistema é restringir o movimento a um plano, o que é conhecido como problema de três corpos restrito [⁶ -⁸] e diminui o número de equações para 12, removendo as seis equações de uma das direções do sistema.

Apesar dessa e de outras simplificações, o problema de três corpos exigiu 80 anos desde a sua formulação até a publicação de uma solução, obtida por Leonhard Euler ^[9], em 1767, na qual os três corpos permanecem dispostos sobre uma linha reta, que rotaciona em torno do centro de massa do sistema. Outra solução foi obtida em 1772 por Joseph Louis Lagrange ^[10], com os três corpos nos vértices de um triângulo equilátero rotacionando. Analiticamente, apenas soluções mais simples como essas foram encontradas inicialmente e o problema continuou sendo estudado em busca de soluções mais gerais.

Entretanto, as limitações da abordagem analítica só ficaram claras bastante tempo depois, após o trabalho de Karl Sundman ^[11], publicado em 1912, o qual propôs uma solução geral construída a partir de uma série infinita nomeada em sua homenagem. Em 1930, David Belorizky ^[12] mostrou que esta série precisaria de cerca de 10⁸⁰⁰⁰⁰ termos para a obtenção da precisão usual nas observações astronômicas para a solução de Lagrange. Ainda assim, o resultado não seria exato, o que mostra a impossibilidade prática de uma solução analítica geral para o problema de três corpos ^[2].

Com as dificuldades que emergiram do trabalho de Sundman, métodos numéricos passaram a ser empregados para a solução do problema. Em 1913, Carl Burrau estudou a chamada solução Pitagórica ^[13], denominada pela configuração inicial do sistema com os corpos inicialmente em repouso nos vértices de um triângulo retângulo. Devido às limitações do cálculo numérico na ocasião, dada a inexistência de computadores, uma característica desta solução permaneceu desconhecida por algumas décadas. Nesta solução e em outras semelhantes, após um determinado tempo um dos corpos é ejetado e um sistema binário se forma. Em outro caso, estudado em 1920 por Ludwig Becker ^[14], um corpo se aproxima de um sistema binário a partir de uma grande distância, causando perturbações sobre ele. Estas duas soluções obtiveram menos destaque por sua falta de periodicidade, como as de Euler e Lagrange. Uma nova solução periódica surgiu apenas em 1993, obtida por Cristopher Moore ^[15], com os três corpos seguindo uma mesma trajetória em formato de oito. Este caso pertence a uma classe de soluções comumente chamadas de coreográficas, as quais têm em particular a característica de que os corpos devem percorrer a mesma trajetória e estarem igualmente espaçados no tempo ^[16].

Outra grande contribuição foi dada por Milovan Šuvakov e Veljko Dmitrašinović ^[17] em 2013, quando obtiveram treze novas soluções periódicas e utilizaram análises algébricas e topológicas para classificá-las em quatro famílias de soluções distintas. Após isso, em 2014, Iasko e Orlov ^[18] obtiveram mais nove soluções para o problema com massas iguais, enquanto que Šuvakov ^[19] encontrou onze soluções na vizinhança da solução em oito. Já em 2015, Hudomal ^[20] apresentou 25 famílias de soluções, incluindo as quatro de Šuvakov e Dmitrašinović. Em 2016, Rose ^[21] encontrou outras 90 soluções periódicas para a configuração isósceles colinear. Mais recentemente, em 2017, Xiaoming Li, Yipeng Jing e Shijun Liao ^[22] obtiveram 1223 novas soluções para o problema com massas diferentes e momento angular nulo. É importante destacar a utilização de computadores na obtenção de todas as soluções mais recentes, o que mostra os avanços nos métodos numéricos utilizados e, em parte, possibilitados pelo problema de três corpos.

O objetivo deste trabalho é utilizar métodos numéricos computacionais para apresentar e discutir duas soluções do problema de três corpos restrito, escolhido por simplicidade e pelo maior número de soluções importantes se comparado ao problema geral. Assim, este trabalho está organizado da seguinte forma. Na seção 2 é apresentado um resumo sobre sistemas dinâmicos e sobre os expoentes de Lyapunov, que auxiliam o estudo proposto. Na seção 3, o problema de três corpos restrito é formulado a partir da segunda lei de Newton ^[1], com base na notação usual de sistemas dinâmicos [³ -⁵], apresentada na seção anterior. Na seção 4, o software Mathematica é utilizado para a obtenção e análise das soluções em linha reta, de Euler ^[9], e no formato de oito, de Moore ^[15]. Por fim, as discussões finais estão na seção 5

2. Sistemas Dinâmicos

Sistemas dinâmicos descrevem o comportamento de todos os pontos de um dado espaço e suas trajetórias com o passar do tempo ^[3]. Quando se trata do espaço de fase na mecânica, para um sistema com $n$ graus de liberdade, existem $N=2n$ dimensões, pois, associadas a cada grau de liberdade, há uma coordenada de posição e outra de velocidade ou de momento. Em geral, na mecânica clássica o espaço de fase é o espaço euclideano ${\mathbb{R}}^N$. Por esse motivo, o vetor geral $\mathit{X}\left(t\right)=(x_1\left(t\right),...,x_N\left(t\right))\in {\mathbb{R}}^N$, chamado de coordenada, guarda toda a informação do sistema em um instante t e as suas componentes $x_i\left(t\right)$ são coordenadas generalizadas e seus respectivos momentos generalizados associados. Além da mecânica, os sistemas dinâmicos são utilizados para modelar problemas em diversas áreas, como na teoria de circuitos ^[23], na meteorologia ^[24] e na cosmologia ^[25]. Desta forma, as componentes do vetor $\mathit{X}\left(t\right)$ não estão estritamente ligadas à mecânica, mas podem representar uma grande variedade de quantidades distintas.

Partindo da suposição de um sistema determinístico, dada uma condição inicial ${\mathit{X}}_{\mathbf{0}}=\mathit{X}\left(t_0\right)$, o estado do sistema pode ser determinado em qualquer instante, anterior ou posterior. Matematicamente, a evolução temporal do sistema pode ser descrita pela equação

(1) $$\dot{\mathit{X}}\left(t\right)=\mathit{C}(\mathit{X}\left(t\right),t;\lambda ),$$

na qual $\mathit{C}$, chamada de campo vetorial, é uma função das coordenadas, do tempo e de um ou mais parâmetros de controle $\lambda$. Os sistemas dinâmicos podem ser classificados em discretos ou contínuos de acordo com a continuidade ou não do parâmetro $t$. Podem ser classificados, ainda, de acordo com a estrutura do campo vetorial. Nesse caso, o sistema é linear ou não-linear seguindo a forma da dependência de $\mathit{C}$ com relação às coordenadas de $\mathit{X}$. Caso $\mathit{C}$ independa explicitamente do tempo, o sistema é dito autônomo e, caso contrário, não-autônomo. Para este segundo caso, é possível incluir uma coordenada extra, $x_{N+1}$, de modo que $x_{N+1}\left(t\right)=t$ para torná-lo um sistema autônomo, o qual é mais simples de se analisar, ao custo de uma dimensão extra no sistema a ser resolvido ^[4]. Por fim, são chamados suaves, caso possam ser derivados uma quantidade arbitrária de vezes.

O parâmetro de controle $\lambda$, na equação (1), indica a presença de uma ou mais constantes que influenciam diretamente no comportamento do sistema. Por exemplo, o sinal das constantes pode determinar a estabilidade das soluções de um dado sistema, bem como tornar seu comportamento caótico, no sentido explicado no decorrer deste trabalho. A análise desses parâmetros indica, em geral, valores críticos ${\lambda }_c$, nos quais ocorrem essas mudanças, chamadas de bifurcações [³ -⁵].

Para tratar dos sistemas dinâmicos suaves, é comum utilizar um mapeamento ${\Phi }_t:{\mathbb{R}}^N\to {\mathbb{R}}^N$ que faz a conexão entre ${\mathit{X}}_0$ e $\mathit{X}\left(t\right)$ através de

(2) $${\text{Φ}}_t\left({\mathit{X}}_{\mathit{0}}\right)=\mathit{X}\left(t\right),$$

chamado de fluxo de fase ^[4]. Utilizando esta relação, a equação (1) se torna

(3) $${\dot{\text{Φ}}}_t\left({\mathit{X}}_{\mathit{0}}\right)=\mathit{C}\left({\text{Φ}}_t\left({\mathit{X}}_{\mathit{0}}\right)\right),$$

e o conjunto $\left\{{\mathit{}\text{Φ}}_t\left({\mathit{X}}_0\right)\right|t\in \mathbb{R}\}$ representa uma solução do sistema, cujo estado inicial é ${\mathit{X}}_0$.

Um resultado importante ^[4] é mostrado a partir da equação

(4) $$\Lambda \left({\mathit{}\text{Φ}}_t\left({\mathit{X}}_0\right)\right)=\nabla ·\mathit{C},$$

a qual relaciona a taxa de mudança de um volume arbitrário no espaço de fase, $\Lambda$, em um estado ${\mathit{}\text{Φ}}_t\left({\mathit{X}}_0\right)$, com a divergência do campo vetorial. Esta equação mostra que um sistema conservativo, quando $\nabla ·\mathit{C}=0$, é caracterizado por um fluxo de fase que preserva os volumes, enquanto que um sistema dissipativo, quando $\nabla ·\mathit{C}<0$, leva à contração destes volumes. No entanto, esta característica é local e, por esse motivo, uma análise global da dinâmica do sistema exige uma média dos valores de $\Lambda$ ao longo de toda a trajetória.

Em sistemas dissipativos, os volumes podem até mesmo tender a zero com o fluxo de fase. Outra possibilidade é a contração em apenas algumas direções, o que torna um volume $N$-dimensional em algum objeto geométrico de menor dimensão. Estes são dois exemplos dos chamados atratores ^[4], os quais podem ser entendidos como subespaços do espaço de fase para os quais as trajetórias, dentro de um certo limite, tendem com o fluxo de fase.

Os atratores podem ser de cinco tipos: pontos fixos, ciclos limite, quasi-periódicos, caóticos e estranhos [⁴,²⁶]. Pontos fixos são soluções estacionárias ${\mathit{X}}_0$ que obedecem

(5) $${\mathit{}\text{Φ}}_t\left({\mathit{X}}_0\right)={\mathit{X}}_0,\forall t,$$

de modo que soluções vizinhas podem se afastar dele com o fluxo, o que o caracteriza como instável, permanecer a uma mesma distância, quado ele é estável, ou mesmo tender para ele, quando é assintoticamente estável. É a estabilidade assintótica que auxilia na definição de um atrator e a da sua chamada bacia de atração, o subconjunto do espaço de fase no qual todas as trajetórias tendem à ele.

Os ciclos limites também apresentam esta propriedade, mas são soluções periódicas de modo que

(6) $${\mathit{}\text{Φ}}_t\left(\mathit{X}\right)={\mathit{}\text{Φ}}_{t+T}\left(\mathit{X}\right),\forall T>0,$$

onde $T$ representa o período da solução. Um atrator quasi-periódico, é similar a um ciclo limite e seu caso mais básico pode ser entendido como uma solução periódica em duas direções distintas, como no movimento em um toro. Um movimento quasi-periódico ocorre quando as duas frequências, neste caso de circulação no toro, ${\omega }_1$ e ${\omega }_2$, formam uma razão irracional, fazendo com que sua trajetória nunca retorne ao mesmo ponto, apesar de completar movimentos em ambas as direções.

Por fim, um atrator caótico é caracterizado por possuir uma grande sensibilidade às condições iniciais e, quando pode ser descrito por uma geometria fractal, é chamado também de atrator estranho. Em geral, estas duas propriedades se manifestam juntas, mas existem exemplos de casos separados ^[27].

Esta grande sensibilidade às condições iniciais define o caos em um sistema dinâmico. Isto implica na imprevisibilidade, mesmo nos sistemas determinísticos propostos neste artigo, e faz com que medições com diferenças ínfimas, da ordem de ${10}^{-8}$ ou menores, levem a resultados completamente diferentes. Devido ao erro intrínseco nas medições que podem ser feitas,os efeitos do caos em um determinado sistema são sempre consideráveis. Isto é o que justifica previsões meteorológicas possuírem boa precisão apenas poucos dias no futuro, mesmo conhecendo bastante sobre o comportamento atmosférico.

Os tipos de atratores em um sistema podem ser determinados através de um conjunto de expoentes estudados por Aleksandr Lyapunov. Nomeados em sua homenagem, seus sinais indicam qualitativamente a dinâmica de um sistema arbitrário e seus valores, quantitativamente, a média logarítmica da taxa de expansão ou contração da distância entre duas trajetórias infinitesimalmente próximas no instante inicial [⁴,²⁶]. Para o estudo da estabilidade de uma trajetória que parte de ${\mathit{X}}_0$, é utilizada uma trajetória em sua vizinhança, partindo de ${\mathit{X}}_0+{\mathit{}\xi }_0$, e analisada a propagação da perturbação ${\mathit{}\xi }_0$ com o fluxo. Matematicamente, após um tempo $t$,

(7) $${\mathit{}\xi }_t={\mathit{}\text{Φ}}_t({\mathit{X}}_0+{\mathit{}\xi }_0)-{\mathit{}\text{Φ}}_t\left({\mathit{X}}_0\right).$$

Nesta equação, a expansão de ${\mathit{}\text{Φ}}_t({\mathit{X}}_0+{\mathit{}\xi }_0)$ em série de Taylor em torno do ponto ${\mathit{X}}_0$ leva a

(8) $${\mathit{}\text{Φ}}_t({\mathit{X}}_0+{\mathit{}\xi }_0)={\mathit{}\text{Φ}}_t\left({\mathit{X}}_0\right)+D_{{\mathit{X}}_0}{\mathit{}\text{Φ}}_t\left({\mathit{X}}_0\right)·{\mathit{}\xi }_0,$$

onde $D_{{\mathit{X}}_0}{\mathit{}\text{Φ}}_t\left({\mathit{X}}_0\right)$ representa a derivada com relação à ${\mathit{X}}_0$ do fluxo de fase de ${\mathit{X}}_0$. Substituir esta equação na (7), resulta em

(9) $${\mathit{}\xi }_t=D_{{\mathit{X}}_0}{\mathit{}\text{Φ}}_t\left({\mathit{X}}_0\right)·{\mathit{}\xi }_0.$$

Com isso, o maior expoente de Lyapunov de ordem 1 é definido através de

(10) $${\sigma }_{{\mathit{X}}_0,{\mathit{}\xi }_0}=\underset{t\to \infty }{lim}\frac{1}{t}ln|D_{{\mathit{X}}_0}{\mathit{}\text{Φ}}_t\left({\mathit{X}}_0\right)·{\mathit{}\xi }_0|\equiv {\sigma }_1.$$

O limite nesta equação garante que o resultado seja o maior expoente para um dado sistema, para o qual foi explicitada a dependência com relação a ${\mathit{X}}_0$ e a ${\mathit{}\xi }_0$. No entanto, dada uma base linearmente independente de $N$ vetores ${\mathit{v}}_i$ para o espaço de fase, existem $N$ expoentes do tipo

(11) $${\sigma }_i={\sigma }_{{\mathit{X}}_0,{\mathit{v}}_i},$$

com os índices organizados de modo que

(12) $${\sigma }_1\ge {\sigma }_2\ge \cdots \ge {\sigma }_N.$$

Este conjunto completo de expoentes de Lyapunov exige um método mais complicado para ser obtido ^[26], o qual envolve repetidas iterações de equações do tipo

(13) $${\sigma }_p=\underset{k\to \infty }{lim}\frac{1}{kT}\sum _{i=1}^{k}ln\left(\right|{\mathit{u}}_p^i\left|\right),$$

onde ${\mathit{u}}_p^i$ indica o $p$-ésimo vetor linearmente independente da base ${\mathit{v}}_p$ ortogonalizado durante o processo de ortonormalização de Gram-Schmidt ^[26], durante a $i$-ésima iteração.

Através da implementação numérica da equação (13), é possível obter o conjunto completo de expoentes de Lyapunov para um determinado sistema dinâmico. Da combinação de sinais pode-se descobrir o tipo de atrator presente no espaço de fase do sistema [⁴,²⁶]. Por exemplo, em $N$ dimensões, o conjunto de expoentes $({\sigma }_1,\cdots ,{\sigma }_N)$ indica

(14) $$\left\{\begin{array}{c}(-,\cdots ,-),\text{ponto}\text{fixo,}\hfill \\ (0,-,\cdots ,-),\text{c.}\text{limite,}\hfill \\ (0,\cdots ,0,-,\cdots ,-),\text{toro}{T}^{N-1}\text{,}\hfill \\ (+,\cdots ,+,0,\cdots ,0,-,\cdots ,-),\text{caos}{C}^{N-2}\text{.}\hfill \end{array}\right.$$

Na seção seguinte é resumida a formulação do problema de três corpos, descrito através de um sistema dinâmico suave, autônomo e não-linear.

3. O Problema de Três Corpos Restrito

Para a formulação do problema de três corpos restrito, são utilizadas a segunda lei de Newton e a lei da Gravitação Universal ^[1]. A segunda pode ser definida da seguinte maneira ^[6]: se dois corpos têm massa $m_1$ e $m_2$, e se estão separados por uma distância $r$, o corpo 2 irá atuar no centro de massa do corpo 1 com uma força direcionada ao centro de massa do corpo 2 de módulo

(15) $$F_{12}=\frac{G{m}_1{m}_2}{r^2},$$

onde $F_{12}$ denota a força sobre o corpo 1 exercida pelo corpo 2 e G é a constante da gravitação universal e equivale a, aproximadamente, $6,674\times {10}^{-11}{m}^3·k{g}^{-1}·s^{-2}$, no Sistema Internacional de Unidades (SI). Nesse caso, os corpos são considerados puntiformes, com toda a massa concentrada em seus centros de massa. Apesar de restritivo, isso é geralmente válido para os corpos celestes, já que seus tamanhos são muito pequenos em relação distâncias envolvidas em seus movimentos.

Dado um sistema de coordenadas inercial, vetorialmente a posição dos corpos 1 e 2 são representadas pelos vetores ${\mathit{r}}_1$ e ${\mathit{r}}_2$, da origem ao centro de massa do respectivo corpo. Com isso, a distância entre os corpos é o módulo do vetor ${\mathit{r}}_{12}={\mathit{r}}_1-{\mathit{r}}_2$, o qual também determina a direção ${\widehat{\mathit{r}}}_{12}={\mathit{r}}_{12}/\left|{\mathit{r}}_{12}\right|$ entre os centros de massa. Assim, a lei da gravitação pode ser dada por

(16) $${\mathit{F}}_{\mathit{12}}=-\frac{G{m}_1{m}_2}{|{\mathit{r}}_{12}{|}^2}\frac{{\mathit{r}}_{12}}{|{\mathit{r}}_{12}|},$$

onde o sinal negativo indica que a força sobre o corpo 1 está no sentido de 1 para 2. Em um sistema com mais corpos, a força gravitacional resultante sobre o corpo 1 é a soma vetorial das forças de interação com os demais corpos.

Assim, com três corpos de massas $m_1$, $m_2$ e $m_3$, a força resultante sobre o corpo 1 será

(17) $${\mathit{F}}_1=-G{m}_1{m}_2\frac{{\mathit{r}}_{12}}{|{\mathit{r}}_{12}{|}^3}-G{m}_1{m}_3\frac{{\mathit{r}}_{13}}{|{\mathit{r}}_{\mathit{13}}{|}^3}.$$

Seguindo este princípio, pode-se generalizar a equação (17) para $n$ corpos. Assim,

(18) $${\mathit{F}}_i=-\sum _{j\ne i}^{n}G{m}_i{m}_j\frac{{\mathit{r}}_{ij}}{|{\mathit{r}}_{ij}{|}^3}.$$

Em coordenadas cartesianas, a posição de cada corpo tem a forma

(19) $${\mathit{r}}_i=\left[\begin{array}{c}{x}_i\\ y_i\\ z_i\end{array}\right].$$

Devido à maior complexidade de movimentos possíveis em três dimensões, inicialmente tratou-se do problema de três corpos restrito para a análise de suas características básicas. Esta simplificação não deixa de ser fisicamente válida, pois os planetas do sistema solar orbitam o sol aproximadamente em um mesmo plano, por exemplo.

Desta forma, no problema restrito, uma das coordenadas deve ser constante para todos os corpos e, preferencialmente, nula por simplicidade. Então, a restrição ao plano $x-y$ implica que $z_i=0$ e as posições podem ser descritas simplesmente por

(20) $${\mathit{r}}_i=\left[\begin{array}{c}{x}_i\\ y_i\end{array}\right].$$

Com isso, os vetores de distância entre dois dos corpos são

(21) $${\mathit{r}}_{ij}=\left[\begin{array}{c}{x}_i\\ y_i\end{array}\right]-\left[\begin{array}{c}{x}_j\\ y_j\end{array}\right]=\left[\begin{array}{c}{x}_i-x_j\\ y_i-y_j\end{array}\right],$$

e os seus módulos são da forma

(22) $$|{\mathit{r}}_{ij}|={[{(x_i-x_j)}^2+{(y_i-y_j)}^2]}^{\frac{1}{2}}.$$

A segunda lei de Newton é

(23) $${\mathit{F}}_r=m\mathit{}\ddot{\mathit{r}},$$

onde ${\mathit{F}}_r$ é a força resultante sobre um corpo, $m$ é a sua massa e $\mathit{r}$ a sua posição, em coordenadas cartesianas no plano $x-y$. Assim,

(24) $$\left[\begin{array}{c}{F}_x\\ F_y\end{array}\right]=m\left[\begin{array}{c}\ddot{x}\\ \ddot{y}\end{array}\right],$$

onde $F_x$ e $F_y$ são as componentes da força e $x$ e $y$, da posição.

Para o corpo $i$, então, utilizando as equações (18) e (24), as equações de movimento são

(25) $$m_i\left[\begin{array}{c}{\ddot{x}}_i\\ {\ddot{y}}_i\end{array}\right]=-\sum _{j\ne i}^{n}G{m}_i{m}_j\frac{1}{{|r_{ij}|}^{\frac{3}{2}}}·\left[\begin{array}{c}{x}_i-x_j\\ y_i-y_j\end{array}\right].$$

Isto leva ao sistema de equações

(26) $$\left\{\begin{array}{c}{\ddot{x}}_i=-\sum _{i\ne j}^{n}G{m}_j{\displaystyle \frac{x_i-x_j}{{|r_{ij}|}^{\frac{3}{2}}}},\\ {\ddot{y}}_i=-\sum _{i\ne j}^{n}G{m}_j{\displaystyle \frac{y_i-y_j}{{|r_{ij}|}^{\frac{3}{2}}}}.\end{array}\right.$$

Introduzindo, então, a relação ${\dot{x}}_i=v_{x_i}$, esse sistema de equações diferenciais passa a ser de primeira ordem, dado por

(27) $$\left\{\begin{array}{c}{\dot{x}}_i=v_{x_i},\hfill \\ {\dot{y}}_i=v_{y_i},\hfill \\ {\dot{v}}_{x_i}=-\sum _{j\ne i}^{n}G{m}_j{\displaystyle \frac{x_i-x_j}{{|r_{ij}|}^{\frac{3}{2}}}},\hfill \\ {\dot{v}}_{y_i}=-\sum _{j\ne i}^{n}G{m}_j{\displaystyle \frac{y_i-y_j}{{|r_{ij}|}^{\frac{3}{2}}}}.\hfill \end{array}\right.$$

Esse sistema mostra todas as coordenadas do $i$-ésimo corpo sob efeito da gravidade de $n-1$ outros corpos. Para $n$ corpos é necessário resolver um conjunto de $4n$ equações diferenciais ordinárias não-lineares e acopladas, o que é impraticável analiticamente mesmo para $n=3$, como discutido por Belorizky ^[12].

Este sistema, para três corpos, pode ser representado de forma reduzida pelo sistema dinâmico

(28) $$\dot{\mathit{X}}=\mathit{F},$$

(29) $$\begin{array}{ccc}& & \mathit{X}=\left[\begin{array}{c}{x}_i\\ y_i\\ v_{x_i}\\ v_{y_i}\end{array}\right],\hfill \\ & & \text{e}\hfill \\ & & \mathit{F}=\left[\begin{array}{c}{v}_{x_i}\\ v_{y_i}\\ -\sum _{j\ne i}^{n}G{m}_j{\displaystyle \frac{x_i-x_j}{{|r_{ij}|}^{\frac{3}{2}}}}\\ -\sum _{j\ne i}^{n}G{m}_j{\displaystyle \frac{y_i-y_j}{{|r_{ij}|}^{\frac{3}{2}}}}\end{array}\right],\hfill \end{array}$$

com o índice $i$ variando de 1 a 3. Para resolver numericamente esse sistema de 12 equações, é necessário definir as constantes $G$, $m_1$, $m_2$ e $m_3$ e as condições iniciais, o que determina a solução a ser obtida. Apesar de quaisquer soluções podem ser obtidas, nem todas são de interesse, como as que terminam na colisão entre os corpos. Na seção seguinte são discutidas duas soluções importantes para este problema, obtidas numericamente neste trabalho com o software Mathematica.

4. Soluções do Problema de Três Corpos Restrito

Partindo do que foi discutido sobre sistemas dinâmicos e o problema de três corpos, nesta seção são apresentadas duas soluções obtidas para este problema ao longo da história. Uma é a solução em linha reta obtida em 1763, por Leonhard Euler ^[9]. A outra é a solução em oito, obtida numericamente com o auxílio computacional por Cristopher Moore, em 1993 ^[15], e comprovada pouco tempo depois pelos matemáticos Alain Chenciner e Richard Montgomery ^[28].

A solução em linha reta, a primeira obtida após a formulação do problema por Newton, originalmente foi obtida de modo analítico. Nela, os corpos ficam dispostos sobre uma mesma reta, a qual rotaciona em torno do centro de massa do sistema, de modo que a força resultante sobre cada corpo seja nula e as distâncias relativas entre eles sejam constantes.

Numericamente, esta solução pode ser obtida com valores específicos para as condições iniciais e as constantes envolvidas. Para isso, é necessária uma precisão considerável, mostrada abaixo ^[29], para que os corpos se comportem como desejado.

(30) $${\mathit{X}}_0=\left[\begin{array}{c}{x}_{1_0}\\ y_{1_0}\\ x_{2_0}\\ y_{2_0}\\ x_{3_0}\\ y_{3_0}\\ v_{x_{1_0}}\\ v_{y_{1_0}}\\ v_{x_{2_0}}\\ v_{y_{2_0}}\\ v_{x_{3_0}}\\ v_{y_{3_0}}\end{array}\right]=\left[\begin{array}{c}-1\hfill \\ 0\hfill \\ 1.254953728\hfill \\ 0\hfill \\ 2.745046272\hfill \\ 0\hfill \\ 0\hfill \\ -0.3660350371\hfill \\ 0\hfill \\ 0.4593570344\hfill \\ 0\hfill \\ 1.004783114\hfill \end{array}\right].$$

Para estas condições iniciais, são utilizados $G=1$, $m_1=2$ e $m_2=m_3=0.5$, para resolver o sistema (29) numericamente através do comando “NDSolve” do software Mathematica. Estes valores para a constante da gravitação universal e para as massas são requisitos para que a solução em linha reta seja obtida com as condições inicias apresentadas. Ao definir $G=1$, no entanto, as unidades de comprimento, tempo e massa $uc$, $ut$ e $um$, respectivamente) devem obedecer

(31) $$\frac{u{c}^3}{um·u{t}^2}=\frac{1}{6,674\times {10}^{-11}}\frac{m^3}{kg·s^2}.$$

Uma característica do problema de três corpos, é a abrangência originada da equação (31). Definidas as unidades, existem sistemas de comportamento similar em diferentes escalas [⁶,⁷]. Por exemplo, eles podem representar três pequenos satélites a altas velocidades, tomando $uc={10}^{4}m$, $ut=1s$ e $um=66,74kg$ ou três objetos maiores com grandezas da ordem de $uc={10}^{9}m$, $ut={10}^{6}s$ e $um=6,674·{10}^{4}kg$. Isto se relaciona também com as unidades naturais ^[30], como as unidades de Planck, nas quais constantes universais, como $G$, são fixadas como 1 e as demais grandezas são definidas a partir delas.

As figuras a seguir, produzidas com o Mathematica basicamente através dos comandos “ParametricPlot” e “Show”, mostram quaisquer destes casos, de acordo com as unidades utilizadas. Nestas imagens, os corpos foram representados por círculos proporcionais às suas massas, mas é importante ressaltar que os corpos são infinitesimais de acordo com definições do problema e esta é apenas uma opção didática para facilitar a observação da posição dos corpos em um dado instante.

As figuras a seguir representam os corpos em um mesmo plano $x-y$. Elas mostram as trajetórias que os corpos seguem de fato em seus movimentos. A Figura 1 mostra os instantes iniciais da solução.

Figura 1
Instantes iniciais da solução em linha reta, com as condições iniciais (30). O parâmetro $t$ foi utilizado de 0 a $1.5$ para mostrar o início do movimento dos três corpos. Está representada, também, a linha reta sobre a qual os corpos se encontram.

O período desta solução foi estimado em $T=17.1655$, o que está representado na Figura 2.

Figura 2
Um período da solução em linha reta, com as condições iniciais (30). O parâmetro $t$ foi utilizado de 0 a $17.1655$, o que evidencia a trajetória circular e o alinhamento dos corpos.

Como é de se esperar pela complexidade do problema de três corpos, esta solução é instável. Apesar disto não ser o caso para todas as soluções, esta característica está ligada ao caos que surge, com frequência, em sistemas não-lineares com muitas dimensões e até em sistemas com uma dimensão em casos particulares ^[4]. Nestes casos, o sistema é bastante sensível às condições iniciais e sua precisão, o que ocorre nesta solução como pode ser visto na Figura 3, simplesmente adotando um intervalo maior para o parâmetro tempo.

Figura 3
Com as mesmas condições iniciais (30), o parâmetro t foi utilizado de 0 a 100 para mostrar a instabilidade da solução.

Esta instabilidade é uma característica intrínseca desta solução e é evidenciada devido à imprecisão inevitável das condições iniciais (30). Esta situação ilustra um dos desafios de se lidar com sistemas caóticos, pois sempre que são feitas medições, elas estão sujeitas uma precisão limitada, o que dificulta previsões precisas sobre o comportamento do sistema. Por exemplo, seja $\delta$ uma perturbação numérica a ser adicionada em uma componente das condições iniciais (30), o comportamento dos corpos se altera completamente mesmo para valores da ordem de ${10}^{-8}$, como é visível na Figura 4.

Figura 4
O parâmetro t foi utilizado de 0 a 100 com uma perturbação $\delta ={10}^{-8}$ somada em $v_{y_{3_0}}$ nas condições iniciais (30).

Isto implica, também, que em um sistema caótico duas medidas levemente diferentes para as condições iniciais levam a resultados bastante distintos. Esta análise pode ser feita através dos expoentes de Lyapunov, dados pelas equações (13), e obtidos pelo algoritmo de Marco Sandri ^[26]. Para este sistema de 12 dimensões, há 12 equações desse tipo e, para obter os expoentes neste trabalho, elas passaram por $k=1000$ iterações em intervalos $T=0.05$, nas quais o sistema foi resolvido com passos $dt=0.001$. Com estes parâmetros, o algoritmo levou à

(32) $$\left\{\begin{array}{c}{\sigma }_1=0.0936397,\hfill \\ {\sigma }_2=0.0886711,\hfill \\ {\sigma }_3=0.0109841,\hfill \\ {\sigma }_4=-0.0011041,\hfill \\ {\sigma }_5=-0.00770106,\hfill \\ {\sigma }_6=-0.00837444,\hfill \\ {\sigma }_7=-0.00905421,\hfill \\ {\sigma }_8=-0.00949284,\hfill \\ {\sigma }_9=-0.0130357,\hfill \\ {\sigma }_{10}=-0.0175909,\hfill \\ {\sigma }_{11}=-0.0388666,\hfill \\ {\sigma }_{12}=-0.088075.\hfill \end{array}\right.$$

A obtenção de três expoentes positivos indica um comportamento caótico nesta solução, o que está de acordo com a análise do movimento dos corpos no plano $x-y$. É possível analisar a distância entre duas soluções com o passar do tempo, para estudar seu comportamento. Assim, seja a distância entre a solução original $\mathit{X}\left(t\right)$ e uma solução perturbada ${\mathit{X}}_{\mathit{p}}\left(t\right)$, dada por $D\left(t\right)=|\mathit{X}\left(t\right)-{\mathit{X}}_{\mathit{p}}\left(t\right)|$ ou, explicitamente, por

(33) $$\begin{array}{c}D\left(t\right)=[{(x_1-x_{p_1})}^2+{(x_2-x_{p_2})}^2+{(x_3-x_{p_3})}^2+\hfill \\ {(y_1-y_{p_1})}^2+{(y_2-y_{p_2})}^2+{(y_3-y_{p_3})}^2+\\ {(v_{x_1}-v_{x_{p_1}})}^2+{(v_{x_2}-v_{x_{p_2}})}^2+{(v_{x_3}-v_{x_{p_3}})}^2+\\ \hfill {(v_{y_1}-v_{y_{p_1}})}^2+{(v_{y_2}-v_{y_{p_2}})}^2+{(v_{y_3}-v_{y_{p_3}})}^2{]}^{\frac{1}{2}}.\end{array}$$

Considerando $\mathit{X}\left(t\right)$ a solução com condições iniciais (30) e ${\mathit{X}}_{\mathit{p}}\left(t\right)$, a solução na Figura 4, com uma perturbação $\delta ={10}^{-8}$ somada em $v_{y_{3_0}}$ nas mesmas condições (30). Para $t\in [0,50000]$, o comportamento de $D\left(t\right)$ com o tempo é mostrado na Figura 5.

Figura 5
A figura mostra a distância entre a solução em linha reta e uma solução perturbada na sua vizinhança. O parâmetro $t$ foi utilizado de 0 a 50000 para a análise da distância entre as soluções, em função do tempo. Neste caso, as soluções evoluem de modo distinto, principalmente para tempos grandes.

Isto mostra que, apesar de distintas, as soluções ainda se mantém, inicialmente, em uma região similar para este caso. Apesar disso, elas se distanciam muito em alguns instantes, dadas as proporções da solução, principalmente para $t>10000$.

A solução em oito, consiste nos três corpos seguindo uma mesma trajetória, de formato característico. Esta foi a primeira das soluções deste tipo, as quais passaram a ser denominadas coreografias [²⁸,³¹] e foram identificadas, também, para sistemas com mais corpos [¹⁶,³²].

Para este caso os três corpos devem ter massas iguais, obedecendo as condições iniciais ^[19]

(34) $${\mathit{X}}_0=\left[\begin{array}{c}{x}_{1_0}\\ y_{1_0}\\ x_{2_0}\\ y_{2_0}\\ x_{3_0}\\ y_{3_0}\\ v_{x_{1_0}}\\ v_{y_{1_0}}\\ v_{x_{2_0}}\\ v_{y_{2_0}}\\ v_{x_{3_0}}\\ v_{y_{3_0}}\end{array}\right]=\left[\begin{array}{c}-1\hfill \\ 0\hfill \\ 1\hfill \\ 0\hfill \\ 0\hfill \\ 0\hfill \\ 0.3471128135672417\hfill \\ 0.532726851767674\hfill \\ 0.3471128135672417\hfill \\ 0.532726851767674\hfill \\ -2\times 0.3471128135672417\hfill \\ -2\times 0.532726851767674\hfill \end{array}\right].$$

Com essas condições, devem ser utilizados $G=m_1=m_2=m_3=1$ para resolver numericamente o sistema. Desse modo, os três corpos percorrem a trajetória desejada e as figuras seguintes os mostram representados no plano $x-y$. Como no caso anterior, elas representam as trajetórias percorridas no espaço e a Figura 6 mostra os instantes iniciais da solução em oito.

Figura 6
Instantes iniciais da solução em oito, com as condições iniciais (34). O parâmetro $t$ foi utilizado de 0 a $1.5$, mostrando a formação da trajetória em oito.

Aumentando o intervalo do parâmetro $t$, pode-se observar a trajetória fechada, após um período $T=6.3250$ ^[19], na Figura 7.

Figura 7
Um período da solução em oito. O parâmetro $t$ foi utilizado de 0 a $6.3250$ e a figura mostra que todos os corpos seguem a mesma trajetória.

Esta solução periódica foi a primeira estável para o problema de três corpos ^[31]. Neste caso, são necessárias perturbações consideráveis para promover, apenas visivelmente, a rotação e/ou translação da trajetória em oito e maiores ainda para que o sistema sistema se desfaça. Por exemplo, o efeito de uma perturbação da ordem de ${10}^{-2}$ em apenas uma direção, $y_{2_0}$, neste sistema é mostrado na Figura 8.

Figura 8
O parâmetro $t$ foi utilizado de 0 a 100 com uma perturbação $\delta ={10}^{-2}$, para maior visibilidade dos seus efeitos, somada em $y_{2_0}$ nas condições iniciais (34).

Neste caso, a figura em oito rotaciona no sentido anti-horário. Analisando a distância $D\left(t\right)$, entre a solução original e esta, é notável que após algum tempo a solução perturbada completa sua rotação e retorna a uma região similar à figura em oito. Isto está ilustrado na Figura 9.

Figura 9
A figura mostra a distância $D\left(t\right)$ entre as soluções. O parâmetro $t$ foi utilizado de 0 a 1000. É possível perceber a aproximação da solução perturbada com relação à original. As regiões mais espessas do gráfico correspondem a oscilações mais acentuadas de $D\left(t\right)$.

Uma perturbação $\delta ={10}^{-2}$, agora em $v_{x_{3_0}}$, causa o efeito visível na Figura 10.

Figura 10
O parâmetro $t$ foi utilizado de 0 a 100 com uma perturbação $\delta ={10}^{-2}$ somada à $v_{x_{3_0}}$ nas condições iniciais (34).

Neste caso a solução se desloca no sentido positivo do eixo $x$. A Figura 11 mostra o comportamento da distância $D\left(t\right)$ entre esta solução e a original.

Figura 11
O parâmetro $t$ foi utilizado de 0 a 1000. A solução perturbada, em geral, se distancia da original. As regiões mais espessas do gráfico, como antes, correspondem à oscilações acentuadas de $D\left(t\right)$.

É lógico que há um limite para que estas soluções não colapsem completamente, como pode ser observado na Figura 12, na qual a perturbação foi $\delta =0.5$.

Figura 12
O parâmetro $t$ foi utilizado de 0 a 10 com uma perturbação $\delta =0.5$ somada à $x_{3_0}$ nas condições iniciais (34). Neste caso os corpos deixam a região na qual se mantinham para perturbações menores.

Esta solução perturbada ejeta um dos corpos em uma direção e um sistema binário na direção oposta, como mostrado na Figura 13.

Figura 13
O parâmetro $t$ foi utilizado de 0 a 20 com uma perturbação $\delta =0.5$ somada à $x_{3_0}$ nas condições iniciais (34). A ejeção dos corpo 2 do sistema pode ser observada na parte inferior e o sistema binário resultante, na parte superior.

Este resultado é importante pois sua ocorrência é bastante típica em outras configurações iniciais para o problema [⁶,¹³]. Isto indica, por exemplo, que sistemas binários de estrelas podem ter se originado de um sistema com três corpos após a ejeção de um deles, o que auxilia na compreensão daquilo que é possível observar a partir da Terra.

Neste caso, a distância $\mathit{D}\left(t\right)$, como definida anteriormente, se comporta de forma aproximadamente linear após a ejeção, já que a solução original se mantém próxima à origem e a perturbada afasta os corpos com velocidade constante, em decorrência da conservação de momento. Isto é visível na Figura 14.

Figura 14
O parâmetro $t$ foi utilizado de 0 a 100, o que permite a observação do comportamento de $D\left(t\right)$ após a ejeção.

Para a obtenção dos expoentes de Lyapunov, novamente são feitas $k=1000$ iterações das 12 equações (13), em intervalos $T=0.05$, resolvendo o sistema com passos $dt=0.001$ durante as iterações. Com isso, o algoritmo de Marco Sandri resulta em

(35) $$\left\{\begin{array}{c}{\sigma }_1=0.114005,\hfill \\ {\sigma }_2=0.107316,\hfill \\ {\sigma }_3=0.0250382,\hfill \\ {\sigma }_4=0.0109859,\hfill \\ {\sigma }_5=0.00651626,\hfill \\ {\sigma }_6=0.00420007,\hfill \\ {\sigma }_7=-0.0101587,\hfill \\ {\sigma }_8=-0.0103498,\hfill \\ {\sigma }_9=-0.0223216,\hfill \\ {\sigma }_{10}=-0.0344755,\hfill \\ {\sigma }_{11}=-0.0853039,\hfill \\ {\sigma }_{12}=-0.105452.\hfill \end{array}\right.$$

Com seis expoentes positivos indicando um sistema caótico, espera-se que a solução em oito, de referência, e uma solução perturbada, em sua vizinhança, se afastem exponencialmente no espaço de fase com o passar do tempo. Vale destacar que, para uma solução caótica, como a discutida anteriormente, mesmo imprecisões nas condições iniciais ou perturbações muito pequenas acabariam com a sua periodicidade, principalmente para tempos grandes, o que não ocorre neste caso.

O problema está nas limitações do método numérico utilizado neste trabalho, o qual, para tempos grandes apresenta erros consideráveis, bem como no hardware disponível para a execução dos cálculos. A solução em oito é um caso notável, pois contraintuitivamente, dada a complexidade do sistema, trabalhos como o de Carles Simó ^[31] mostram que a estabilidade observada nas figuras apresentadas existe. Vale destacar que os resultados do algoritmo ^[26] utilizado descreveram satisfatoriamente soluções de problemas mais simples e para tempos menores nos testes realizados.

5. Conclusão

A solução em linha reta é bastante sensível às condições iniciais e a perturbações, o que indica o caos inerente a ela. Com as condições iniciais utilizadas, com tempos maiores que dois períodos, a trajetória se desfaz, como mostrado na Figura 3. Com uma perturbação da ordem de ${10}^{-8}$, os corpos descrevem uma trajetória bastante diferente da original, o que está na Figura 4. Por esse motivo, a distância entre essas soluções cresce bastante, principalmente para tempos grandes, de acordo com a Figura 5.

Por outro lado, a solução em oito é estável ^[31]. As soluções perturbadas permanecem na vizinhança da trajetória original, na Figura 6, de modo que as perturbações causam translações e rotações da figura em oito, como nas Figuras 8 e 10, respectivamente. Nesses casos, a distância aumenta pouco para tempos grandes, como visto nas Figuras 9 e 11, além de apresentar certa periodicidade. O comportamento periódico se mantém, a menos que haja uma perturbação da ordem de ${10}^{-1}$, para a qual a solução se desfaz, a exemplo da Figura 13. A Figura 14 mostra que a distância entre essa solução e a original cresce de forma aproximadamente linear, devido a ejeção de um dos corpos do sistema.

Para o primeiro caso foi obtido um conjunto de três expoentes de Lyapunov positivos e nove negativos. Isto indica o caos e a instabilidade das soluções, como parece intuitivo por conta da complexidade do problema de três corpo e das perturbações analisadas. No entanto, é notável no segundo caso que apesar de seis expoentes positivos, até certo limite para as perturbações, as soluções se comportam estavelmente. Este fato está mais ligado ao método numérico utilizado neste trabalho, que foi limitado pelo hardware disponível, o que impediu cálculos mais precisos e para tempos maiores. Vale destacar que, para a solução em linha reta de Euler, por exemplo, a combinação de sinais para os expoentes representa bem o sistema.

Agradecimentos

À Carles Simó, pelos esclarecimentos acerca da estabilidade da solução em oito, e, por fim, ao PIBIC e ao CNPq, pelo apoio financeiro.

Referências

• ^[1]
  I. Newton, Philosophiæ Naturalis Principia Mathematica (Royal Society e Typis Streater, Londres, 1687), p. 510, 3ª ed.
• ^[2]
  M. Valtonen, J. Anosova, K. Kholshevnikov, A. Mylläri, V. Orlov e K. Tanikawa, The Three-body Problem from Pythagoras to Hawking (Springer, New York, 2016), p. 188, 1ª ed.
• ^[3]
  M.W. Hirsch, S. Smale e R.L. Devaney, Differential Equations, Dynamical Systems & An Introduction to Chaos (Academic Press, San Diego, 2003), p. 425, 2ª ed.
• ^[4]
  W. Greiner, Classical Mechanics: Systems of Particles and Hamiltonian Dynamics (Springer, New York, 2010), p. 580, 2ª ed.
• ^[5]
  J.A. Villate, Introdução aos Sistemas Dinâmicos: Uma Abordagem Prática com Maxima, disponível em http://www.villate.org/doc/sistemasdinamicos/sistdina m-1_2.pdf
  » http://www.villate.org/doc/sistemasdinamicos/sistdina m-1_2.pdf
• ^[6]
  M. Valtonen e H. Karttunen, The Three-Body Problem (Cambridge University Press, New York, 2006), p. 356.
• ^[7]
  C. Marchal, The Three-Body Problem (Elsevier Science, New York, 1990), p. 592.
• ^[8]
  F.R. Moulton, An Introduction to Celestial Mechanics (The MacMillan Company, New York, 1914), p. 478, 2ª ed.
• ^[9]
  L. Euler, Novi Comm. Acad. Sci. Imp. Petrop. 11, 144 (1767).
• ^[10]
  J.L. Lagrange, Essai sur le Problème des Trois Corps, disponível em http://sites.mathdoc.fr/cgi-bin/oeitem?id=OE_LAGRANGE__6_229_0
  » http://sites.mathdoc.fr/cgi-bin/oeitem?id=OE_LAGRANGE__6_229_0
• ^[11]
  K.F. Sundman, Acta Mathematica 36, 105 (1912).
• ^[12]
  D. Beloriszky, Bulletin Astronomique 6, 417 (1930).
• ^[13]
  C. Burrau, Astronomical Notes 13, 8 (1913).
• ^[14]
  L. Becker, Monthly Notices of the Royal Astronomical Society 80, 590 (1920).
• ^[15]
  C. Moore, Physical Review Letters 70, 3675 (1993).
• ^[16]
  G.I. Depetri, Coreografias no Problema de N corpos Dissertação de Mestrado, Universidade de São Paulo, São Paulo (2011).
• ^[17]
  M. Šuvakov e V. Dimitrašinović, Physical Review Letters 110, 1 (2013).
• ^[18]
  P.P. Iasko e V.V. Orlov, Astronomy Reports 58, 869 (2014).
• ^[19]
  M. Šuvakov, Celestial Mechanics and Dynamical Astronomy 119, 369 (2014).
• ^[20]
  A. Hudomal New periodic solutions to the three-body problem and gravitational waves Dissertação de Mestrado, Universidade de Belgrado, Belgrado (2015).
• ^[21]
  D.B. Rose, Geometric phase and periodic orbits of the equal-mass, planar three-body problem with vanishing angular momentum Tese de Doutorado, Universidade de Sydney, Sydney (2016).
• ^[22]
  L. Xiaoming, J. Yipeng e L. Shijun The 1223 new periodic orbits of planar three-body problem with unequal mass and zero angular momentum, disponível em: https://arxiv.org/abs/1709.04775, acesso em 09/10/2017.
  » https://arxiv.org/abs/1709.04775
• ^[23]
  T.S. Parker e L.O. Chua, Practical Numerical Algorithms for Chaotic Systems (Springer-Verlag, New York, 1989), p 348.
• ^[24]
  E.N. Lorenz, Journal of the Atmospheric Sciences 20, 130 (1963).
• ^[25]
  C.G. Böhmer e N. Chan, in: Dynamical and Complex Systems editado por S. Bullett, T. Fearn e F. Smith (World Scientifc, Londres, 2017), p. 121.
• ^[26]
  M. Sandri, The Mathematica Journal 6, 78 (1996).
• ^[27]
  C. Grebogi, E. Ott, S. Pelikan e J.A. Yorke, Physica D: Nonlinear Phenomena 13, 261 (1984).
• ^[28]
  A. Chenciner e R. Montgomery, Annals of Mathematics 152, 881 (2000).
• ^[29]
  Y. Nakamura Applications by Dr. Yasuyuki Nakamura, disponível em: http://www.maplesoft.com/applications/Author.aspx?mid=927, acesso em 01/05/2017.
  » http://www.maplesoft.com/applications/Author.aspx?mid=927
• ^[30]
  D. Trancanelli, Rev. Bras. Ensino Fís. 38, e2505 (2016).
• ^[31]
  C. Simó, in: Proceedings of an International Conference on Celestial Mechanics, Evanston, 1999, (American Mathematical Society, Providence, 2002), p. 209.
• ^[32]
  C. Simó, Progress in Mathematics 201, 101 (2001).

Datas de Publicação

• Publicação nesta coleção
  2019

Histórico

• Recebido
  21 Dez 2017
• Revisado
  30 Jul 2018
• Aceito
  12 Set 2018