Soleil-Terre-Lune

1. Introduction

L'objectif est d'écrire et de résoudre les équations différentielles du mouvement du système à trois corps Soleil-Terre-Lune. Il est impossible d'obtenir une solution analytique de ces équations; nous devons donc effectuer une intégration numérique.

Nous considérons tout d'abord le problème plus général du mouvement des corps dans le système solaire.

2. Équations du mouvement

Dans tout ce qui suit, on entend par planète tout corps du système solaire autre que le Soleil, y compris les satellites comme la Lune, les comètes, les astéroïdes.

L'application des lois de la mécanique classique au système solaire suppose l'existence d'un référentiel inertiel à l'échelle du système solaire. Le système solaire (Soleil+planètes) étant soumis à l'influence gravitationnelle des autres corps de la galaxie, le centre de masse C du système solaire est nécessairement fixe dans ce référentiel. On appelle référentiel barycentrique (RB) le référentiel inertiel dans lequel le centre de masse C est fixe. En pratique, ce référentiel est réalisé avec le centre C et deux étoiles assez lointaines pour que leur mouvement soit négligeable sur l'échelle de temps considérée (de l'ordre du siècle). Sachant que l'étoile la plus proche du système solaire est à une distance environ 6000 fois plus grande que celui-ci, le référentiel RB peut être considéré en mouvement de chute libre dans un champ uniforme. Autrement dit, les forces de marée sont négligeables dans ce référentiel et l'on peut ignorer complètement l'influence gravitationnelle des autres corps de la galaxie.

Le référentiel RB est peu pratique car le centre C n'est pas directement observable. On introduit donc le référentiel héliocentrique RH, en mouvement de translation par rapport à RB. RH n'est pas tout à fait galiléen, puisqu'il est accéléré par rapport à RB. Il faudra donc prendre en compte les forces d'inertie d'entraînement dans ce référentiel. Ce réferentiel RH est réalisé par le centre du Soleil et deux étoiles assez lointaines (les trois non alignés pour définir un système rigide de référence).

On introduit aussi le référentiel Terre-Lune (RTL) en translation par rapport à RB mais centré sur le centre de masse du système Terre-Lune, et le référentiel géocentrique (RG) en translation par rapport à RB mais centré sur le centre de masse de la Terre.

Soit M_s la masse du Soleil et m_i la masse des planètes, l'indice i variant de 1 à N, où N est le nombre de planètes considérées. La constante de gravitation est notée G. La position du centre du Soleil est S, celle du centre de la planète i est P_i. Considérons tout d'abord l'accélération du centre du Soleil dans le référentiel inertiel RB. Le Soleil est soumis à l'influence gravitationnelle des planètes donc :

{\vec{a}}_{S / R B} = \sum_{j = 1}^{N} - G m_{j} \frac{\vec{P_{j} S}}{P_{j} S^{3}} (1)

Dans le référentiel non inertiel RH (en translation par rapport à RB), il y a une force d'inertie par unité de masse égale à l'opposée de cette accélération.

En tenant compte de cette force d'inertie et des forces gravitationnelles directes, l'accélération d'une planète dans le référentiel RH s'écrit :

\vec{a_{i}} = - G M_{s} \frac{\vec{S P_{i}}}{S P_{i}^{3}} - \sum_{j \neq i} G m_{j} \frac{\vec{P_{j} P_{i}}}{P_{j} P_{i}^{3}} - \sum_{j} G m_{j} \frac{\vec{S P_{j}}}{S P_{j}^{3}} (2)

En sortant le terme de la planète considérée dans la force d'inertie, on obtient :

\vec{a_{i}} = - G (M_{s} + m_{i}) \frac{\vec{S P_{i}}}{S P_{i}^{3}} - \sum_{j \neq i} G m_{j} (\frac{\vec{P_{j} P_{i}}}{P_{j} P_{i}^{3}} + \frac{\vec{S P_{j}}}{S P_{j}^{3}}) (3)

3. Unités astronomiques internationales

Des unités astronomiques ont été définies :

L'unité de temps est le jour (j), soit 86400 s.
L'unité de masse est la masse du Soleil M_s.
L'unité de longueur astronomique (UA) est le rayon de l'orbite circulaire d'une planète de masse négligeable et non perturbée, dont la vitesse angulaire autour du Soleil est de k radian par jour, où k est la constante de Gauss.

La constante de Gauss a été fixée en 1938 pour que l'UA soit égale au demi-grand axe de l'orbite terrestre à cette date. Sa valeur est :

k = 0, 01720209895 r a d \cdot j^{- 1} (4)

Sur une échelle de temps de plusieurs siècles, le demi-grand axe de l'orbite terrestre évolue très lentement, mais la masse du Soleil reste constante, et donc l'unité astronomique de change pas. La troisième loi de Kepler s'écrit :

G M_{s} = a^{3} k^{2} (5)

Si l'on utilise les unités astronomiques, on a donc :

G M_{s} = k^{2} = 2, 95912208286 \cdot 1 0^{- 4} U A^{3} \cdot j^{- 2} (6)

Pour écrire les équations différentielles, il faut les masses des planètes exprimées en unité de masse solaire. Voici les rapports de la masse du Soleil sur la masse des planètes principales du système Solaire. Pour les planètes comportant des satellites, c'est la masse totale du sytème planète-satellite qui est donnée (à l'exception de la Terre).

Rapports des masses
Planete	Ms/mi
Mercure	6023600.0
Vénus	408523.5
Terre	332946.0
Terre-Lune	328900.5
Lune	27068620.9
Mars	3098710.0
Jupiter	1047.355
Saturne	3498.5
Uranus	22869.0
Neptune	19314.0

L'inverse de ces valeurs donne la masse m_i de chaque planète exprimée en masse solaire.

Voici finalement l'expression de l'accélération d'une planète dans le système d'unités astronomiques :

\vec{a_{i}} = - k^{2} (1 + m_{i}) \frac{\vec{S P_{i}}}{S P_{i}^{3}} - k^{2} \sum_{j \neq i} m_{j} (\frac{\vec{P_{j} P_{i}}}{P_{j} P_{i}^{3}} + \frac{\vec{S P_{j}}}{S P_{j}^{3}}) (7)

4. Équations différentielles

On note (x_i,y_i,z_i) les coordonnées de la planète dans un repère lié au référentiel héliocentrique RH. Les composantes de sa vitesse sont notées (u_i,v_i,w_i). Le système étudié comporte N planètes. Pour chaque planète, voici les 6 équations différentielles :

\begin{matrix} \frac{d x_{i}}{d t} = u_{i} & (8) \\ \frac{d y_{i}}{d t} = v_{i} & (9) \\ \frac{d z_{i}}{d t} = w_{i} & (10) \\ \frac{d u_{i}}{d t} = - k^{2} ((1 + m_{i}) \frac{x_{i}}{r_{i}^{3}} + \sum_{j \neq i} m_{j} (\frac{x_{i} - x_{j}}{r_{i j}^{3}} + \frac{x_{j}}{r_{j}^{3}})) & (11) \\ \frac{d v_{i}}{d t} = - k^{2} ((1 + m_{i}) \frac{y_{i}}{r_{i}^{3}} + \sum_{j \neq i} m_{j} (\frac{y_{i} - y_{j}}{r_{i j}^{3}} + \frac{y_{j}}{r_{j}^{3}})) & (12) \\ \frac{d w_{i}}{d t} = - k^{2} ((1 + m_{i}) \frac{z_{i}}{r_{i}^{3}} + \sum_{j \neq i} m_{j} (\frac{z_{i} - z_{j}}{r_{i j}^{3}} + \frac{z_{j}}{r_{j}^{3}})) & (13) \end{matrix}

où les distances au cube sont :

\begin{matrix} r_{i}^{3} = {(x_{i}^{2} + y_{i}^{2} + z_{i}^{2})}^{3 / 2} & (14) \\ r_{i j}^{3} = {({(x_{i} - x_{j})}^{2} + {(y_{i} - y_{j})}^{2} + {(z_{i} - z_{j})}^{2})}^{3 / 2} & (15) \end{matrix}

5. Conditions initiales

L'institut de mécanique céleste et de calcul des éphémérides (IMCCE) fournit les éphémérides des corps du système solaire.

Nous utiliserons le générateur d'éphémérides Miriade pour obtenir les positions et les vitesses initiales des corps. Voici les paramètres à choisir :

Target : choix d'une planète ou d'un satellite.
Epoch : choix de la date initiale, du nombre de dates, et de l'intervalle de temps.
Reference center : heliocenter.
Advanced parameters : choisir le plan de référence écliptique et les coordonnées rectangulaires.

Obtenir les données sous forme d'un fichier texte. Voici un exemple, avec 5 dates espacées de 1 jour (réglage par défaut) :

# Miriade.ephemcc.results
# Request:
# targetType: Planet, targetNumber: 3, targetName: Earth, Diameter: 12756.27 km, Orbital period: 3.65256000E+02 d, Time_Scale: UTC, Planetary_Theory: INPOP13C, Coordinates: Ecliptic rectangular, frameType: Astrometric J2000, frameCentre: Heliocenter
# Nb rows: 5
Target, Date, X (au), Y (au), Z (au), Xp (au/day), Yp (au/day), Zp (au/day), Observer distance (au)
Earth,   2018-11-14T18:14:07.00,    0.6081696313585,    0.7802823762238,   -0.0000365394517,   -0.0138546789557,    0.0105075516480,    0.0000001487763,    0.9892981797655
Earth,   2018-11-15T18:14:07.00,    0.5942232455120,    0.7906696191502,   -0.0000364055301,   -0.0140373686738,    0.0102664152828,    0.0000001147281,    0.9890701256632
Earth,   2018-11-16T18:14:07.00,    0.5800963459705,    0.8008141872244,   -0.0000363184239,   -0.0142157027089,    0.0100222201737,    0.0000000554127,    0.9888453541314
Earth,   2018-11-17T18:14:07.00,    0.5657932981497,    0.8107130761994,   -0.0000363024983,   -0.0143896628625,    0.0097750751071,   -0.0000000272574,    0.9886241699805
Earth,   2018-11-18T18:14:07.00,    0.5513184813428,    0.8203633898588,   -0.0000363798366,   -0.0145592401045,    0.0095250871543,   -0.0000001305845,    0.9884068800925

Le tableau comporte la date, les coordonnées héliocentriques en unité astronomique (AU), la distance au Soleil, et les composantes héliocentriques de la vitesse, en AU par jour. Voici comment extraire les positions et les vitesses du fichier :

import numpy
data = numpy.loadtxt("Terre.txt",skiprows=5,delimiter=",",usecols=(2,3,4,5,6,7))
Y_terre = data[0] # condition initiale pour la Terre au 20/8/2016  0h00

print(Y_terre)
--> array([  6.08169631e-01,   7.80282376e-01,  -3.65394517e-05,
        -1.38546790e-02,   1.05075516e-02,   1.48776300e-07])

Nous allons intégrer les équations pour le système Soleil-Terre-Lune. Il faut donc les conditions initiales pour la Lune :

# Miriade.ephemcc.results
# Request:
# targetType: Satellite, targetNumber: 10, targetName: Moon, Diameter: 3474.21 km, Orbital period: 2.73217000E+01 d, Time_Scale: UTC, Planetary_Theory: INPOP13C, Coordinates: Ecliptic rectangular, frameType: Astrometric J2000, frameCentre: Heliocenter
# Nb rows: 5
Target, Date, X (au), Y (au), Z (au), Xp (au/day), Yp (au/day), Zp (au/day), Observer distance (au)
Moon,   2018-11-14T18:18:12.00,    0.6099662214182,    0.7783294111015,   -0.0000942317615,   -0.0134485470602,    0.0108859715065,   -0.0000479239391,    0.9888657504113
Moon,   2018-11-15T18:18:12.00,    0.5963807969035,    0.7891374321571,   -0.0001404834598,   -0.0137243922241,    0.0107264591470,   -0.0000442244392,    0.9891450659445
Moon,   2018-11-16T18:18:12.00,    0.5825138099485,    0.7997746062336,   -0.0001819909412,   -0.0140111043378,    0.0105437836588,   -0.0000384573255,    0.9894249808235
Moon,   2018-11-17T18:18:12.00,    0.5683561772417,    0.8102163204757,   -0.0002167602738,   -0.0143050280650,    0.0103350786681,   -0.0000307787965,    0.9896864539636
Moon,   2018-11-18T18:18:12.00,    0.5539027725410,    0.8204352086928,   -0.0002429837565,   -0.0146018666113,    0.0100977233274,   -0.0000214089195,    0.9899102343811

data = numpy.loadtxt("Lune.txt",skiprows=5,delimiter=",",usecols=(2,3,4,5,6,7))
Y_lune = data[0] # condition initiale pour la Lune au 20/8/2016  0h00

print(Y_lune)
--> array([  6.09966221e-01,   7.78329411e-01,  -9.42317615e-05,
        -1.34485471e-02,   1.08859715e-02,  -4.79239391e-05])

On choisit par convention de ranger les variables dans l'ordre suivant : coordonnées de la Terre, vitesses de la Terre, coordonnées de la Lune, vitesse de la Lune. La condition initiale complète est donc obtenue en concaténant les deux précédentes :

Yi = numpy.append(Y_terre,Y_lune)

Dans certains cas, il est nécessaire d'avoir séparément les conditions initiales sur les positions et sur les vitesses :

Xi=[Y_terre[0],Y_terre[1],Y_terre[2],Y_lune[0],Y_lune[1],Y_lune[2]]
Vi=[Y_terre[3],Y_terre[4],Y_terre[5],Y_lune[3],Y_lune[4],Y_lune[5]]

6. Travaux pratiques

On se propose d'intégrer les équations différentielles du mouvement du système à trois corps Soleil-Terre-Lune.

Les quatres parties seront programmées dans quatre scripts distincts, ce qui permettra de comparer les résultats.

6.a. Utilisation de scipy.integrate.ode

(1) Télécharger les fichiers suivants contenant les conditions initiales : Terre.txt, Lune.txt, Jupiter.txt.

(2) Définir les constantes k (constante de Gauss), k², m_T (masse de la Terre) et m_L (masse de la Lune).

(3) Écrire une fonction systeme(Y,t) qui calcule les dérivées des différentes variables. Celles-ci sont rangées dans le tableau Y dans l'ordre suivant :

xT,yT,zT : 0,1,2
uT,vT,wT : 3,4,5
xL,yL,zL : 6,7,8
uL,vL,wL : 9,10,11

La fonction scipy.integrate.odeint s'utilise de la manière suivante :

Y=odeint(systeme,Yi,t,rtol=1e-5,atol=1e-5)

où t est un tableau contenant les instants auxquels les valeurs des variables doivent être mémorisées et renvoyées dans le tableau Y. La méthode numérique implémentée par la fonction ode effectue une adaptation du pas de temps, ce qui permet de contrôler l'erreur locale au cours du calcul. La tolérance relative rtol et la tolérance absolue atol permettent d'ajuster le degré de précision souhaité. On prendra des valeurs égales pour ces deux tolérances. Plus elles sont petites, plus l'erreur locale est faible. Evidemment, une erreur plus faible est obtenue avec plus de pas de calculs et donc un temps d'exécution plus long. La fonction ode est une interface vers le solveur LSODE de la bibliothèque ODEPACK (implémentée en FORTRAN).

Le tableau Y renvoyé par cette fonction contient les valeurs des variables aux instants demandés, rangés en colonne. Par exemple, la variable xT est accessible par Y[:,0] (première colonne).

(3) Effectuer l'intégration numérique avec la fonction ode pour une durée de 365 jours.

(4) Tracer les courbes suivantes :

trajectoire de la Terre, projetée sur le plan XY;
trajectoire de la Lune dans le référentiel géocentrique, projetée sur le plan XY;
distance Terre-Soleil en fonction du temps;
distance Terre-Lune en fonction du temps;

Afin d'avoir la même échelle sur les axes X et Y, utiliser pyplot.axes().set_aspect('equal').

(5) Refaire le calcul avec des tolérances de plus en plus en petites jusqu'à ne plus voir de variations dans les résultats.

(6) D'après ces résultats, le mouvement de la Lune dans le référentiel géocentrique suit-il les lois de Kepler ?

(7) Pour voir l'influence de la Lune sur le mouvement de la Terre, refaire les calculs (à la suite du script) après avoir annulé la masse de la Lune par une valeur très petite. Tracer l'écart entre la position de la Terre avec la Lune et sa position sans la Lune.

6.b. Méthode d'Euler

Cette partie vise à retrouver les résultats précédents avec la méthode d'Euler.

(1) Reprendre les fonctions pas_euler et euler étudiées en cours.

(2) Effectuer l'intégration et reprendre les tracés faits dans la partie précédente, avec un pas de temps h = 0,1.

(3) Abaisser le pas de temps jusqu'à obtenir des résultats semblables à ceux obtenus avec la fonction ode.

6.c. Méthode d'Euler asymétrique

Cette partie vise à retrouver les résultats précédents avec la méthode d'Euler asymétrique. À chaque pas de calcul, les accélérations sont calculées après les nouvelles positions.

(1) Écrire une fonction accel(X) qui calcule les dérivées secondes des différentes variables. Celles-ci sont rangées dans le tableau Y dans l'ordre suivant :

xT,yT,zT : 0,1,2
xL,yL,zL : 3,4,5

(2) Reprendre les fonction pas_eulerA et eulerA étudiées en cours.

(3) Effectuer l'intégration et reprendre les tracés faits dans la partie précédente, avec un pas de temps h = 0,1.

(4) Abaisser le pas de temps jusqu'à obtenir des résultats semblables à ceux obtenus avec la fonction ode.

(5) Comparer à la méthode d'Euler standard.

6.d. Influence de Jupiter

Mettre en place l'intégration des équations différentielles du système Soleil-Terre-Lune-Jupiter, afin de voir si Jupiter a une influence sur le mouvement de la Terre et de la Lune.