Problème à trois corps


Hypothèses :
On considère le Soleil de masse M1, une masse M2 = µ.M1 et une masse m beaucoup plus petite que M1 et M2. On suppose que M1 et M2 sont immobiles. On se place dans le plan Π = M1M2m et on prend l'origine sur M1, M1-M2 donne la direction de l'axe Ox. A l'instant t initial les coordonnées de m sont x0 et y0. Le vecteur vitesse initiale de m est contenu dans Π et ses composantes sont vx0 et vy0. On pose D = OM2, r1 = Om et r2 = M2-m.

Équations du mouvement :
La masse m est soumise aux forces F1 = G.M1.m.r1 / r13 et F2 = G.M2.m.r2 / r23.
Les équations du mouvement de m sont donc :
m.d2x / dt2 = − G.M1.m.x / r13 − G.M2.m.(x − D) / r23.
m.d2y / dt2 = − G.M1.m.y / r13 − G.M2.m.y / r23 .
     d2x / dt2 = − G.M1.x / r13 − G.µM1.(x − D) / r23.
     d2y / dt2 = − G.M1.y / r13 − G.µ.M1.y / r23.
Ces équations admettent une solution analytique mais elle est plutôt pénible à calculer. Ici nous faisons une résolution numérique en intégrant ces équations par la méthode de Runge-Kutta à l'ordre 4.
Les unités SI ne sont pas bien adaptées au calcul et à la représentation graphique. On effectue le changement d'unités suivant :
Comme unité de masse, on prend la masse M du Soleil, comme unité de longueur L on prend l'unité astronomique qui est la distance moyenne du Soleil à la Terre et comme unité de temps T l'année terrestre. La raison de ce changement d'unités est que d'après la troisième loi de Kepler, pour le système solaire on a : G.M.T2 / L3 = 4π2 = K.
Montrer qu'avec ces unités, les équations deviennent :
     d2x / dt2 = − K.[x / r13 −µ.(x − D) / r23].
     d2y / dt2 = − K.[ y / r13 −µ.y / r23].
Le calcul numérique pose des problèmes quand r1 ou r2 deviennent très petits. Pour limiter ceux-ci, on utilise un pas de calcul variable : on prend un pas très petit quand r1 ou r2 sont petits et on augmente le pas quand m s'éloigne de M1 ou de M2.
Pour contrôler la validité des calculs, on utilise le fait que l'énergie mécanique totale de m est constante. A l'instant t elle vaut :
E = ½m.v2 − G.M1.m / r1 − G.M2.m / r2
Montrer qu'avec les unités choisies, l'énergie par unité de masse E / m s'écrit : ε = v2 / 2 − K.[1 / r1 + µ / r2]
A partir des conditions initiales, on calcule ε0. A chaque pas du calcul, on compare la valeur de ε à celle de ε0. Le calcul est arrêté si |(ε − ε0) / ε0| devient trop grand. Cette situation correspond en général à l'impact de m sur M1 ou M2.

Utilisation :
Avec les différents curseurs, régler les valeurs initiales de D, x0, y0, vx0 et vy0. Le programme interdit que la position initiale de m soit trop proche de M1 ou de M2.
Avec la zone de texte choisir la valeur de µ puis cliquer sur le bouton [Départ] pour lancer le calcul.
Le programme trace la trajectoire suivie (en vert), le vecteur vitesse(en rouge) et indique les valeurs des composantes de la position et de la vitesse ainsi que la valeur du pas et celle de l'énergie par unité de masse de m.
Attention : A cause du pas variable la vitesse paraît faible au voisinage des masses (le pas est alors petit) et grande loin des masses (le pas est grand) alors qu'en réalité c'est l'inverse.
Les coordonnées des positions successives sont rangées dans deux tableaux. Quand les tableaux sont pleins, on décale l'origine du tracé de la trajectoire.
Il est possible de prendre une valeur nulle pour µ. On est alors ramené au problème de Kepler et les orbites sont des ellipses ou des branches d'hyperboles.
Le programme affiche dE = 0 si la quantité |ε − ε0 / ε0| est inférieure à 10−3.
Le programme est arrêté si cette quantité dépasse 3.10−2 . Pour les valeurs négatives de l'énergie on a en général un système lié. Il est possible en agissant sur les divers paramètres de choisir une énergie initiale pratiquement nulle. Dans le cas ou µ est nul, les trajectoires sont alors des paraboles.
On peut noter dans certains cas de figure la grande sensibililité du système aux conditions initiales.

Expérimentez en essayant de prévoir la trajectoire en fonction des conditions initiales choisies.