une solution exacte du problème à trois corps

Il y a quelques jours sur google+  je suis tombé sur une image gif animée (comme celle ci-dessous) montrant une solution exacte du problème à 3 corps (pour les forces gravitationnelles) découverte par Cris Moore :

 une question m'est alors venu à l'esprit : "Peut il vraiment exister, quelque part dans l'univers, 3 étoiles dans une telle configuration ?" . Pour un mathématicien la réponse à cette question ne relève pas du hasard mais d'une notion mathématique générale dans l'étude des systèmes dynamiques : "la stabilité des solutions". La question étant très complexe un logiciel de calcul numérique peut être très utile pour se faire une idée de la réponse.

1 modélisation du problème 

Comme je l'ai déjà expliqué dans un précédent billet mécanique céleste avec scilab le principe fondamental de la dynamique donne facilement l'équation différentielle qui gouverne le système à 3 corps :
$$ {d^2\over dt^2}X_i =\sum_{j\neq i}-{{\cal G}m_j(X_i-X_j)\over |X_i-X_j|^3}=\sum_{j\neq i} F_g(X_i-X_j,m_j) ~~~~i=1,2,3$$
avec 6 conditions initiales : $X_i(0), {d\over dt}X_i(0) $ données pour $i=1,2,3$.  Pour résoudre ce système numériquement sous Scilab il suffit de l'écrire sous forme d'un système d'ordre 1 :
$$ U=\begin{pmatrix}X_1\\ {d\over dt}X_1\\ X_2\\ {d\over dt}X_2\\ X_3\\{d\over dt}X_3 \end{pmatrix}
~~~~{d\over dt}U =F(U)=\begin{pmatrix}
U_2\\
F_g(U_1-U_3,m_2)+F_g(U_1-U_5,m_3)\\
U_4\\
F_g(U_3-U_1,m_1)+F_g(U_3-U_5,m_3)\\
U_6\\
F_g(U_5-U_3,m_2)+F_g(U_5-U_1,m_1)\\
 \end{pmatrix} $$
et d'utiliser la commande ode. Si on se place dans le cas d'un mouvement plan les termes $U_i$ sont des vecteurs de ${\mathbb R}^2$ contenant les potions $x_i,y_i$ et les vitesses  $vx_i,vy_i$  des 3 corps :
$$ U_1=\begin{pmatrix}x_1\\y_1\end{pmatrix},~~U_3=\begin{pmatrix}x_2\\y_2\end{pmatrix},~~U_5=\begin{pmatrix}x_3\\y_3\end{pmatrix}~~~~
U_2=\begin{pmatrix}vx_1\\vy_1\end{pmatrix},~~U_4=\begin{pmatrix}vx_2\\vy_2\end{pmatrix},~~U_6=\begin{pmatrix}vx_3\\vy_3\end{pmatrix} $$

 Les données initiales permettant d'obtenir la trajectoire découverte par Cris Moore sont données dans l'article "A remarkable periodic solution of the three-body problem in the case of equal masses"A. Chenciner, & R. Montgomery


$$x_1=−x_2=0.97000436~~~y_1=-y_2=−0.24308753,~~~x3=y_3=0; $$
et
$$vx_3=−2vx_1=−2vx_2=−0.93240737~~~vy_3=−2vy_1=−2vy_2−0.86473146$$

avec un choix d'unités telles que $m_1=m_2=m_3=1$ et ${\cal G}=1$. Mon script, initialisé avec ces données, donne bien la trajectoire attendue avec une période $T\approx 6.3$ bonne approximation de la valeur théorique.  Voici le script Scilab utilisé pour obtenir l'animation au début de ce billet :

// ********************************
/// définition des fonctions
// ********************************
//forces gravitationnelles
function [u2]=force_g(t, u, masse)
    module=-G*masse*((u(1)^2+u(2)^2)^(-3/2))
    u2=[module*u(1); module*u(2)]
endfunction
//fonction décrivant le système
function [du]=force(t, u, masse1, masse2, masse3)
        u1=[u(1);u(2)]
        du1=[u(3);u(4)]
        u2=[u(5);u(6)]
        du2=[u(7);u(8)]
        u3=[u(9);u(10)]
        du3=[u(11);u(12)]
        ddu1=force_g(t,u1-u2,masse2)+force_g(t,u1-u3,masse3)
        ddu2=force_g(t,u2-u1,masse1)+force_g(t,u2-u3,masse3)
        ddu3=force_g(t,u3-u1,masse1)+force_g(t,u3-u2,masse2)
        du=[du1;ddu1;du2;ddu2;du3;ddu3]
endfunction
// ********************************
/// données initiales
// ********************************
//constantes
G=1;//constante gravitationnelle
m0=1;;m=[m0,m0,m0];//masse des 3 corps
dt=0.02;//pas de temps
T=40;//durée totale de la simulation
t=[0:dt:T];//intervalle de temps
N=length(t);//nombre de pas de temps
dx=0.07;dy=0.07;// taille affichage des 3 corps
//coordonnées des 3 corps à t=0
x3=0; y3=0; vx3=-0.93240737; vy3=-0.86473146;
x1=0.97000436;y1=-0.24308753;vx1=-vx3/2;vy1=-vy3/2;
x2=-x1;y2=-y1;vx2=-vx3/2;vy2=-vy3/2;
u0=[x1; y1; vx1; vy1; x2; y2; vx2; vy2;x3; y3; vx3; vy3];
// ********************************
/// calcul de la solution
// ********************************
u=ode(u0,0,t,list(force,m(1),m(2),m(3)));// calcul par ode
X=[u(1,:)',u(5,:)',u(9,:)'];Y=[u(2,:)',u(6,:)',u(10,:)'];//récupération des trajectoires
// ********************************
/// l'animation 
// ********************************
rect=[min(X),min(Y),max(X),max(Y)]// taille fenêtre
clf();plot2d(0,0,-1,rect=rect,frameflag=3)//initialisation fenêtre graphique
tail=6/dt;//longueur de la trajectoire tracée
for k=2:N
   drawlater()
   //affichage trajectoires
   clf(); plot2d(0,0,0,rect=rect,axesflag=0,frameflag=3)
   plot2d(X(max(1,k-tail):k,:),Y(max(1,k-tail):k,:),[5 2 3],rect=rect,axesflag=0,frameflag=0)
   //affichage des 3 corps
   P1=[X(k,1)-dx/2;Y(k,1)+dy/2;dx;dy;0;360*64];
   P2=[X(k,2)-dx/2;Y(k,2)+dy/2;dx;dy;0;360*64];
   P3=[X(k,3)-dx/2;Y(k,3)+dy/2;dx;dy;0;360*64];
   xfarcs([P1,P2,P3],[5,2,3])
   //trajectoires en gras
   A=gca(); A.children(2).children([1:3]).thickness=3;
   //affichage du temps
   xinfo('t='+string(t(k)))
   drawnow()
end
//trajectoire globale 
clf()
plot2d(X,Y,[5 2 3],rect=rect,frameflag=3,axesflag=0)
 A=gca();A.children(1).children([1:3]).thickness=2;

 2 étude de la stabilité

Maintenant que nous pouvons simuler les trajectoire des trois corps on peut se faire une idée de la stabilité du système en modifiant légèrement les conditions initiales pour voir comme évoluerai le système.  pour celà il suffit de rajouter une partie "aléatoire"  aux conditions initiales du problème juste avant de lancer ode.
 La première idée est de regarder l'influence des positions initiales  des 3 corps sur le mouvement :    

//perturbations des données initiales 
pert=8;// amplitude de la perturbation en %
u0=u0+(pert/100)*rand(u0);//position/vitesses

en faisant varier le paramètre pert de 0 à 10%  on va voir la trajectoire idéale se déformer :

   

ces quelques essais montre que bien que  la trajectoire ne soit plus périodique pour une perturbation même faible (<1%) l'aspect général de la trajectoire reste comparable à la trajectoire idéale .... enfin jusqu'à un certain niveau de perturbation (10%) où là le comportement du système devient pour le moins chaotique. La deuxième idée qui m'est venue est de regarder l'influence des masses des 3 corps sur le système. j'ai donc simulé le mouvement pour des masses s'écartant légèrement les unes des autres  avec :

//perturbations des données initiales 
pert=8;// amplitude de la perturbation en %
u0=u0+(pert/100)*rand(u0);//position/vitesses
m=(1+ pert*rand(1:3)/100)*m0// masses

 
voilà une des animations obtenues avec ce type de perturbations :

  je trouve le résultat étonnant : la trajectoire globale est assez chaotique  mais pourtant le système semble rester  mieux confiné que si les 3 masses étaient égales!  Par contre on remarque qu'à de très nombreuses reprises 2 des 3 corps s'approchent très près l'un de l'autre, à tel point qu'on peut se demander si une collision est fortement probable pour un tel système? L'exemple ci-dessous n'a rien de particulier, on obtient le même genre de phénomènes à quasiment tous les essais, mais dans certains cas un  des 3 corps est éjecté du système :

Conclusion, à mon avis, il est difficile de croire qu'un tel système existe dans l'univers, une perturbation minime des paramètres de départ  suffisant à le déstabiliser.