Hypothèses :
On considère le Soleil de masse M₁, une masse M₂ = µ.M₁ et une masse m beaucoup plus petite que M₁ et M₂. On suppose que M₁ et M₂ sont immobiles. On se place dans le plan Π = M₁M₂m et on prend l'origine sur M₁, M₁-M₂ donne la direction de l'axe Ox. A l'instant t initial les coordonnées de m sont x₀ et y₀. Le vecteur vitesse initiale de m est contenu dans Π et ses composantes sont v_x0 et v_y0. On pose D = OM₂, r₁ = Om et r₂ = M₂-m.

Équations du mouvement :
La masse m est soumise aux forces F₁ = G.M₁.m.r₁ / r₁³ et F₂ = G.M₂.m.r₂ / r₂³.
Les équations du mouvement de m sont donc :
m.d²x / dt² = − G.M₁.m.x / r₁³ − G.M₂.m.(x − D) / r₂³.
m.d²y / dt² = − G.M₁.m.y / r₁³ − G.M₂.m.y / r₂³ .
     d²x / dt² = − G.M₁.x / r₁³ − G.µM₁.(x − D) / r₂³.
     d²y / dt² = − G.M₁.y / r₁³ − G.µ.M₁.y / r₂³.
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.T² / L³ = 4π² = K.
Montrer qu'avec ces unités, les équations deviennent :
     d²x / dt² = − K.[x / r₁³ −µ.(x − D) / r₂³].
     d²y / dt² = − K.[ y / r₁³ −µ.y / r₂³].
Le calcul numérique pose des problèmes quand r₁ ou r₂ deviennent très petits. Pour limiter ceux-ci, on utilise un pas de calcul variable : on prend un pas très petit quand r₁ ou r₂ sont petits et on augmente le pas quand m s'éloigne de M₁ ou de M₂.
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.v² − G.M₁.m / r₁ − G.M₂.m / r₂
Montrer qu'avec les unités choisies, l'énergie par unité de masse E / m s'écrit : ε = v² / 2 − K.[1 / r₁ + µ / r₂]
A partir des conditions initiales, on calcule ε₀. A chaque pas du calcul, on compare la valeur de ε à celle de ε₀. Le calcul est arrêté si |(ε − ε₀) / ε₀| devient trop grand. Cette situation correspond en général à l'impact de m sur M₁ ou M₂.

Utilisation :
Avec les différents curseurs, régler les valeurs initiales de D, x₀, y₀, v_x0 et v_y0. Le programme interdit que la position initiale de m soit trop proche de M₁ ou de M₂.
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é |ε − ε₀ / ε₀| est inférieure à 10⁻³.
Le programme est arrêté si cette quantité dépasse 3.10⁻² . 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.