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...
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.
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.
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...
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 v 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, téléchargeable dans la bibliothèque des codes de TangenteX, 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.
Dans ma première expérience, j'utilise les données de Christopher Moore pour les masses, positions et vitesses initiales.
Voilà les orbites obtenues :
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.
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} \).
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 103 périodes. Voici les orbites :
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 :
et soudain, pour epsp = 0.007, j'obtiens :
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 :
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.
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 :
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.
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 :
Pour une perturbation de 10-3 et une durée d'évolution de 100 périodes, j'obtiens:
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.
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.
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...