samedi 13 octobre 2012

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.


 

11 commentaires:

  1. Bonjour,
    j'ai le problème suivant, merci de m'aider.
    Soit l'équation $$x'(t) = y(t) - x^3(t) +x(t),y'(t) = -x(t).$$ avec pour tout $t_0 \in \mathbb{R},$ on a $(x(t_0),y(t_0)) = (x_0,y_0).$

    -Soit $(x(t),y(t))$ une trajectoire de l'équation. Comment montrer que $(-x(t),-y(t))$ est aussi une trajectoire de l'équation?

    RépondreSupprimer
    Réponses
    1. pose $u(t)=-x(t)$ et $v(t)=-y(t)$ exprime ensuite $u'(t)$ et $v'(t)$ en fonction de $u(t)$ et $v(t)$ (bien sûr il faudra passer par $x(t)$ et $y(t)$ ... ) et tu trouveras évidement que :
      $$u'(t) = v(t) - u^3(t) +u(t),v'(t) = -u(t)$$
      ce qui veut bien dire que $(u(t),v(t))$ est une trajectoire de l'équation.

      Supprimer
  2. Bonjour,
    Est-ce que la solution trouvée a une relation avec la lemniscate de Bernoulli avec a=1,15 ?
    Le post est très intéressant.
    Merci.

    RépondreSupprimer
    Réponses
    1. Effectivement c'est une très bonne remarque! Même si la courbe me faisait penser à une lemniscate je n'avais essayé de superposer les deux courbes. J'arrive à une superposition parfaite avec l'orbite en traçant avec scilab :
      t=[-10:0.01:10];
      x=1.081*(t+t.^3)./(1+t.^4);
      y=(t-t.^3)./(1+t.^4);
      plot(x,y,'-c')
      j'ai pris le paramètre a=1.081... par ce que c'est le maximum de la coordonnée X sur l'orbite (faire max(X) et aussi min(X)). En cherchant un peu je suis tombé sur un papier
      "Choreographic Three Bodies on the Lemniscate" de Toshiaki Fujiwara, Hiroshi Fukuda, Hiroshi Ozaki que tu trouveras à l'adresse http://arxiv.org/abs/math-ph/0210044 . A première vue la constante devrait être a=1 mais j'ai l'impression qu'en fait les résultats de la page 4 mènent à a=(4/(2+sqrt(3)))^2=1.1487483 ... ce qui serait très proche de ce que tu proposes.

      Supprimer
  3. Très intéressant! Mais on dirait que plus personne sur internet, même d'un niveau doctorat à l'université n'est capable d'écrire la langue française correctement. MI-NA-BLE!

    RépondreSupprimer
    Réponses
    1. Merci pour cette remarque constructive, courageux "anonyme" ! Hé oui il y a certainement beaucoup fautes dans mes textes, faute d'être relus et par ce qu'ils sont écrits trop rapidement sur mon temps libre. L'orthographe n'a jamais été mon fort , j'ai rarement eu plus de 0/20 en dicté au collège, mais ça ne m'a pas empêché d'avoir mon bac C mention très bien avec 17/20 au bac de Français , comme quoi c'est pas si important l'orthographe :-)

      J'ai eu la malchance d'être élève à une époque où on ne savait qualifier que de "nul" ou "fainéant" quelqu'un qui ne réussissait pas. Et dans mon cas on ne s'est jamais posé la question de la dyslexie même si dans une dictée j'écrivais plusieurs fois "ent" pour marquer le pluriel d'un nom et "s" pour un verbe alors que j'avais de très bon résultats en grammaire et conjugaison .... crois moi ce texte "minable" est déjà un progrès, et je continuerai a essayer de m'améliorer, avec mes pauvres moyens, et absolument pas grâce aux gens comme toi.

      Supprimer
  4. Bonjour,
    Merci une autre fois par la rapidité de la réponse et de m’avoir consacré un peu de temps libre. Je suis très content de votre blog. Je le lit souvent vos articles de mathématiques et d’astronomie que vous écrivez, dans mon petit temps libre que je dispose. Grace à ce blog j’ai appris beaucoup de concept de mécanique céleste, Scilab, et mathématique. L’approximation que j’ai trouvée par tâtonnement est :
    t=linspace(0,10,100)
    r=sqrt(1.15*cos(2*t))
    y=r.*sin(t)
    x=r.*cos(t);
    plot(x,y);

    Je lirai avec tranquillité le pdf.
    Ne faites pas attention au commentaire policier, sans aucune direction. Tout le monde est libre de dire ce qu’il veut, mais dans ce cas, comme on disait quand on était au collège, c’est un hors sujet, qui était vraiment très pénalisé. J’espère que ce type de commentaire n’arrête pas la diffusion de ce type de blog.
    Pardon par mon Français que je ne pratique pas depuis 23 ans quand j’ai fini mon Bac C.

    RépondreSupprimer
  5. Bonjour, votre site m'intéresse beaucoup je cherchais à savoir comment je pouvais faire des perturbations non pas avec scilab mais avec python pour le problème à trois corps, sauriez-vous comment faire? Merci d'avance.

    RépondreSupprimer
    Réponses
    1. On peut faire avec python les mêmes choses qu'avec scilab, il faut regarder du coté des bibliothèques scipy
      -l'ensemble des opérations de base ,calculs arithmétiques manipulations de listes/matrices/tableaux, se traduisent assez facilement normalement (bibliothèque "Numpy")
      - pour les tracés avec la commande "plot" de la librairie python "matplotlib" s'effectuent d'une manière très proche de ceux effectués en scilab avec plot2d (d'ailleurs on peut aussi utiliser "plot" avec scilab )
      - il reste à regarder le solveur d'équation différentielles disponible en python, il s'appelle aussi ode et marche d'une manière similaire à celui de scilab : http://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.ode.html




      Supprimer
  6. Bonjour,

    qu'en est-il pour le problème à deux corps, par ex. le soleil et jupiter ? Comment résoud-on la partie theorique pour déterminer les conditions initiales pour les positions et les vitesses en n'ayant que la distance entre le soleil et jupiter et leur masse comme données ?
    Merci !

    RépondreSupprimer
    Réponses
    1. je ne suis pas sûr de comprendre la question. Pour deux corps de masse donnée connaître la distance qui les sépare ne permet pas de déterminer la trajectoire, il faut en plus connaître la vitesse relative des 2 corps. En fait il y a une infinité de trajectoire qui passe par ce point suivant la vitesse ... La position et la vitesse forment les conditions initiales du problème. Tu trouveras peut être des réponses dans les autres billets de ce blog :
      - un modèle simplifié pour le problème à 2 corps
      - résolution numérique d'un problème à 3 corps

      Supprimer

Pour écrire des formules mathématiques vous pouvez utiliser la syntaxe latex en mettant vos formules entre des "dollars" $ \$....\$ $ par exemple :
- $\sum_{n=1}^\infty {1\over n^2}={\pi^2\over 6}$ s'obtient avec \sum_{n=1}^\infty {1\over n^2}={\pi^2\over 6}
- $\mathbb R$ s'obtient avec {\mathbb R} et $\mathcal D$ s'obtient avec {\mathcal D}
- pour les crochets $\langle .,. \rangle$ dans les commentaires utilisez \langle .,. \rangle
vous pouvez écrire du html dans les commentaires :
- italique <i> ... </i> gras <b> ... </b>
- lien <a href="http://adresse "> .... </a>