Le problème à trois corps

L'idée de cette page m'est venue en lisant un roman de science-fiction de l'auteur chinois Liu Cixin, qui s'intitule justement "Le problème à trois corps". Dans un autre système stellaire (Alpha Centauri), formé de trois étoiles, une planète abrite une civilisation. Le roman décrit l'influence du comportement chaotique de la trajectoire de leurs étoiles et de leur planète : très dur ! Pour info, Alpha Centauri est vraiment un système stellaire à trois étoiles dont Proxima du Centaure qui est l'étoile la plus proche de notre Soleil (4 années lumière). Ce système stellaire comporterait trois planètes où, a priori, il ne ferait pas bon vivre !
Un comportement chaotique dans un système à trois corps ? Voyons cela...

L'énoncé du problème

La loi de la gravitation universelle

Vous la connaissez tous : lorsque deux corps de masse m1 et m2 sont en interaction gravitationnelle, loin de tout autre corps, alors la force que l'un exerce sur l'autre est égale à \( \vec{F_{12}} = -Gm_1m_2\dfrac{\vec{r_{12}}}{r_{12}^3}\) et \( \vec{F_{12}} = - \vec{F_{21}}\).  G est la constante de gravitation universelle et \( r_{12} \) est la distance entre les corps 1 et 2.

Que se passe-t-il si trois corps sont en interaction gravitationnelle ? La force de gravitation newtonienne est cumulative. En d'autres termes, l'ensemble de l'univers exerce ses forces gravitationnelles sur chacun d'entre vous. Si on se limite à trois corps, et c'est déjà beaucoup, alors chaque corps subit l'attraction des deux autres, ce que je peux exprimer par \( \vec{F_i} = \sum_{i \ne j} \vec{F_{ij}} = -Gm_i  \sum m_j \dfrac{\vec{r_{ij}}}{r_{ij}^3}\).

Pour trois corps, j'ai donc 3 équations vectorielles de ce genre, soit 6 équations scalaires en projetant sur un plan. Jusqu'ici, rien de bien compliqué, rien que l'application des lois de Newton à un système de trois corps en interaction gravitationnelle.

Pour déterminer les trajectoires de ces trois corps, je vais projeter ces 3 équations vectorielles sur un plan et appliquer la seconde loi de Newton, notre PFD préféré, ce qui va me donner un système différentiel composé de 6 équations.

Les équations différentielles d'un système à trois corps

En reprenant mon équation ci-dessus et en projetant, j'obtiens pour chaque corps de masse \( m_i \) et de coordonnées X, l'équation \( \dfrac{d^2 X_i}{dt^2} = -G\sum_{i \ne j} \dfrac{m_j}{r_{ij}^2}\). Ce que je traduis par le système différentiel suivant :

\( \dfrac{d^2 x_1}{dt^2} = -G( \dfrac{m_2 (x_1 - x_2)}{r_{12}^3} + \dfrac{m_3 (x_1 - x_3)}{r_{13}^3})  \)

\( \dfrac{d^2 y_1}{dt^2} = -G( \dfrac{m_2 (y_1 - y_2)}{r_{12}^3} + \dfrac{m_3 (y_1 - y_3)}{r_{13}^3})  \)

\( \dfrac{d^2 x_2}{dt^2} = -G( \dfrac{m_1 (x_2 - x_1)}{r_{12}^3} + \dfrac{m_3 (x_2 - x_3)}{r_{23}^3})  \)

\( \dfrac{d^2 y_2}{dt^2} = -G( \dfrac{m_1 (y_2 - y_1)}{r_{12}^3} + \dfrac{m_3 (y_2 - y_3)}{r_{23}^3})  \)

\( \dfrac{d^2 x_3}{dt^2} = -G( \dfrac{m_1 (x_3 - x_1)}{r_{13}^3} + \dfrac{m_3 (x_3 - x_2)}{r_{23}^3})  \)

\( \dfrac{d^2 y_3}{dt^2} = -G( \dfrac{m_1 (y_3 - y_1)}{r_{13}^3} + \dfrac{m_3 (y_3 - y_2)}{r_{23}^3})  \)

Voilà un joli petit système différentiel. Le problème, c'est qu'il n'est pas intégrable. Poincaré a démontré qu'il n'existait pas de solution générale à ce problème ! Il en existe des solutions analytiques exactes, mais elles souffrent d'un défaut rédhibitoire : il s'agit de séries divergentes qui ne permettent donc pas de prévision à long terme. Nous allons donc intégrer ce système numériquement.

Le script Python

Définition du système différentiel

Je reprends la définition du système ci-dessus, en réarrangeant les signes et je code tout ça dans la fonction TroisCorps :

def TroisCorps(Y,t):

    global m1,m2,m3

    dY = arange(12,dtype=float)

    # extraction des coordonnées

    x1 = Y[0]; y1 = Y[2]

    x2 = Y[4]; y2 = Y[6]

    x3 = Y[8]; y3 = Y[10]

    # calcul des distances r(i,j)

    r12 = ((x1 - x2)**2 + (y1 - y2)**2)**1.5

    r13 = ((x1 - x3)**2 + (y1 - y3)**2)**1.5

    r23 = ((x2 - x3)**2 + (y2 - y3)**2)**1.5

    # calcul de l'accélération

    dY[1] = -m2*(x1 - x2)/r12 - m3*(x1 - x3)/r13

    dY[3] = -m2*(y1 - y2)/r12 - m3*(y1 - y3)/r13

    dY[5] = m1*(x1 - x2)/r12 - m3*(x2 - x3)/r23

    dY[7] = m1*(y1 - y2)/r12 - m3*(y2 - y3)/r23

    dY[9] = m1*(x1 - x3)/r13 + m2*(x2 - x3)/r23 

    dY[11] = m1*(y1 - y3)/r13 + m2*(y2 - y3)/r23

    # retour des dérivées premières

    dY[0] = Y[1];dY[2] = Y[3]

    dY[4] = Y[5];dY[6] = Y[7]

    dY[8] = Y[9];dY[10] = Y[11]

    return dY

Rien de particulier à signaler dans le codage. Le vecteur Y contient donc les valeurs des positions en (x,y) des 3 corps et leur vitesse, qui serviront à odeint() pour l'intégration numérique du système différentiel. A partir du vecteur Y, je pourrai donc tracer les orbites des trois corps. Pour rappel, une orbite est l'ensemble des points x(t) parcourus par t dans le domaine de définition des solutions de l'EDO (l'espace de phase).

Ah si, vous aurez noté que j'ai posé G = 1, ce qui ne nuit pas au modèle...

Les conditions initiales

Dans un premier temps, je vais utiliser des conditions initiales un peu particulières, nous verrons pourquoi plus loin. Les masses des 3 corps sont identiques et égales à 1.

m1 = 1.0

m2 = 1.0

m3 = 1.0

# position initiale des corps - solution exacte

pos0_1 = [0.97000436,-0.24308753]

pos0_2 = [-pos0_1[0],-pos0_1[1]]

pos0_3 = [0.0,0.0]

# vitesse initiale des corps - solution exacte

v0_3 = [-0.93240737,-0.86473146]

v0_2 = [-v0_3[0]/2.0, -v0_3[1]/2.0]

v0_1 = [-v0_3[0]/2.0, -v0_3[1]/2.0]

Vous constatez que les positions et les vitesses initiales ont des valeurs bien précises. Elles correspondent aux conditions initiales nécessaires pour obtenir une solution exacte du problème à trois corps, solution due au physicien Christopher Moore en 1993.

Le script

Le script ne possède aucune difficulté particulière. Nous retrouvons les instructions classiques d'intégration d'un système d'EDO par odeint().

CI = [pos0_1[0],v0_1[0],pos0_1[1],v0_1[1],\

pos0_2[0],v0_2[0],pos0_2[1],v0_2[1],\

pos0_3[0],v0_3[0],pos0_3[1],v0_3[1]]

tmin = 0.0

tmax = 2*pi/3.0

nbpas = 1000

t = linspace(tmin, tmax, nbpas)

Y = odeint(TroisCorps,CI,t)

Je vous laisse découvrir les instructions d'affichage, qui sont tout à fait classiques.

Expériences et analyse

La solution exacte du système à 3 corps

Valeurs des paramètres masses, positions et vitesses initiales

Dans ma première expérience, j'utilise les données de Christopher Moore pour les masses, positions et vitesses initiales.

Orbites du système

Voilà les orbites obtenues :

Trois corps - solution exacte

Ces orbites semblent stables, au moins pendant un intervalle de temps assez long. Mais que devient cette stabilité si l'on modifie les conditions initiales, même faiblement ? C'est ce que nous allons vérifier maintenant.

L'influence des positions initiales

Introduction d'une perturbation sur les positions initiales

Le principe est d'introduire une petite perturbation dans les valeurs des positions initiales. Je vais donc introduire une faible variation paramétrable sur le tuple pos0_1 qui vaut initialement :

pos0_1 = [0.97000436,-0.24308753]

Pour ce faire, j'ajoute à mon code, après avoir défini pos0_1, les instructions :

epsp = 0.0

pos0_1[0] = pos0_1[0]*(1.0 - epsp)

pos0_1[1] = pos0_1[1]*(1.0 - epsp)

Pour que l'influence de la perturbation soit visible, je modifie la durée de la manipulation en choisissant tmax = nbp fois la période initiale \( \dfrac{2 \pi}{3} \).

Orbites obtenues selon l'importance de la perturbation

Je vais commencer par régler mes paramètres avec une perturbation très petite, un millième par exemple, et une durée égale 10^3 périodes. Voici les orbites :

Trajectoire nbp=1000 eps=0.001

Vous constatez que les orbites ne sont plus aussi stables et oscillent plus ou moins autour de leurs valeurs initiales.

J'augmente maintenant la perturbation en la portant à 0.007 par pas de 5*10^-4, avec la même durée d'évolution. Jusqu'à epsp = 0.00065, j'ai un ensemble d'orbites à peu près dense, qui ressemble à ça pour epsp = 0.0065 :

Trajectoire nbp=1000 eps=0.0065

et soudain, pour epsp = 0.007, j'obtiens :

Trajectoire nbp=1000 eps=0.007

Un vrai fouilli d'orbites qui partent dans tous le sens. Pour y voir un peu plus claire, je diminue la durée d'évolution, en la fixant à 300 périodes :

Trajectoire nbp=300 eps=0.007

Les orbites sont très désordonnées mais restent néanmoins assez groupées pour les valeurs de paramètres mentionnées. En modifiant les paramètres, on peut voir le système "exploser" avec parfois l'éjection d'un composant du système.

L'influence de chaque masse composant le système

Modifions la répartition des masses

Conservons les positions initiales de la solution exacte, mais modifions les masses, par exemple en posant m2 = 5. Les orbites, pour une durée d'évolution de 10 périodes sont :

Trajectoire nbp=10 m2 = 5

Même après une durée assez brève, ici 10 périodes, le système est manifestement instable. On arrive même à l'éjection du corps bleu. La modification de la répartition des masses relatives produit un système instable.

Introduction d'une perturbation sur les masses des composants

Revenons maintenant aux conditions initiales de masses de la solution exacte (les 3 masses égales) et affectons une perturbation de masse sur m1, avec la même méthode que pour les positions. Pour une perturbation de 5*10^-4 et une durée de 50 périodes, j'obtiens :

Trajectoire nbp=50 epsm = 0.0005

Pour une perturbation de 10^-3 et une durée d'évolution de 100 périodes, j'obtiens:

Trajectoire nbp=100 epsm = 0.001

Plus l'évolution du système est longue et plus l'orbite de chaque composant devient instable en s'écartant de son orbite initiale. A long terme, le système est susceptible de de disperser.

Notre système est chaotique

En faisant différentes manipulations, par exemple en conservant les positions initiales, sans perturbation, mais en modifiant la valeur relative des masses, on obtient des choses curieuses. Parfois, Python se plante parce que l'intégrateur odeint() rend l'âme. Il faudrait utiliser ode() en employant un autre algo d'intégration plus apte à traiter les problèmes un peu raides. On obtient aussi des situations intéressantes comme celle ci-dessus, obtenue avec m1 = 1 , m2 = m3 = 2.

Trajectoire nbp=10 variation de m

On observe que l'un des corps est carrément expulsé du système après des soubresauts assez incroyables.

Notre système à trois corps présente un comportement que l'on appelle le chaos déterministe. Car vous avez noté qu'il n'y a rien d'aléatoire dans le système différentiel qui gouverne son comportement. Il est simplement très sensible aux conditions initiales. Cette sensibilité est d'autant plus visible lorsque l'on suit l'évolution du système pendant une durée assez longue. Je vous invite à faire plusieurs manips en faisant varier les différents paramètres que j'ai rassemblé en début de code : la répartition des masses, leur position et les perturbations sur ces grandeurs (resp. epsm et epsp). Utilisez aussi plusieurs durées, plus ou moins longues en réglant nbp. Mais attention, il arrive que odeint() se plante. Il utilise l'algo lsoda (Adam implicit et BDF) de la librairie FORTRAN odepack. Avec ode(), on peut utiliser d'autres algo plus appropriés comme dopri5.

Pour la petite histoire, Liu Cixin situait son problème dans Alpha Centauri. Henri Poincaré, lui, se posait la question de la stabilité de notre système solaire. La réponse est qu'à long terme (plusieurs millions ou dizaines de millions d'années) notre système solaire est instable. Comme sur le diagramme ci-dessus, il n'est pas exclu qu'un des corps le constituant soit éjecté ! Le plus beau est que cela ne doit rien au hasard mais bien aux lois de la mécanique...

Le script Python

Le script TroisCorps.py étudié dans cette page est disponible dans le package TroisCorps.zip


Contenu et design par Dominique Lefebvre - www.tangenteX.com mars 2017   Licence Creative Commons    Contact : PhysiqueX ou

Cette œuvre est mise à disposition selon les termes de la Licence Creative Commons Attribution - Pas d’Utilisation Commerciale - Pas de Modification 3.0 France.