Analyse de la trajectoire d'une comete

Analyse de l'orbite d'une comète (ou d'un satellite)

Description du problème

Le système différentiel à intégrer

Considérons le système simplifié constitué du Soleil et d'une comète orbitant autour de celui-ci. Nous négligerons dans la suite toutes les forces autres que l'attraction gravitationnelle du Soleil sur la comète. Cette problématique est bien sur très semblable à l'étude de l'orbite d'un satellite autour de la Terre à une altitude suffisamment élevée.
Nous savons que la comète subit de la part du Soleil une force qui s'exprime par (les vecteurs sont en gras, comme d'habitude):

F = -G*M*m*r/r³

Notez le signe -, qui traduit le caractère attractif de la force. Les paramètres et variables sont très classiques : M la masse du Soleil, m celle de la comète et r le rayon vecteur. On ne présente plus G.
Appliquons maintenant la seconde loi de Newton, ce qui me donne finalement:

g = -G*M*r/r³

Vous aurez sans doute remarqué que j'ai fait une approximation, tout à fait justifiée, en ramenant la masse réduite du système Soleil+comète à la masse du Soleil!
Nous obtenons donc, en projetant sur les axes Ox et Oy d'un référentiel héliocentrique supposé galiléen pour notre problème, le système différentiel suivant:

d²x/dt² = -K*x/(x² + y²)^3/2

d²y/dt² = -K*y/(x² + y²)^3/2

avec K = G*M, que nous allons bien sur adimensionner.

Vu la situation, nous allons utiliser des unités astronomiques (UA). L'unité de longueur sera la distance moyenne Terre-Soleil, soit environ 1.496*10^-11 m, notée L. L'unité de temps sera l'année terrestre, c'est à dire le temps nécessaire à la Terre pour parcourir son orbite autour du Soleil, orbite supposée circulaire, soit 2*π*L. On peut donc lier l'unité de temps à l'unité de longueur. Avec ces unités caractéristiques, le coefficient K prend la valeur K = 4*π² UA, valeur que vous retrouverez dans le programme.
Nous procéderons à l'intégration numérique de ce système, en précisant les constantes d'intégration à l'aide des valeurs initiales des données (x,y) et (vx,vy).

La conservation de l'énergie

Nous savons que la force F(r) est conservatrice et qu'elle dérive d'un potentiel U(r), tel que F(r) = -dU(r)/dr * r/r. L'énergie totale de la comète est donc E = mv²/2 - Km/r dans notre système d'unités.
Nous calculerons donc continument l'énergie potentielle et l'énergie cinétique de la comète pour vérifier la conservation de l'énergie, avec:

Ec = 0.5*m*(vx*vx + vy*vy)
Ep = -K*m/sqrt(x*x + y*y)

Par convention, de manière toute à fait arbitraire, nous fixerons m = 1.

Quelques rappels sur les ellipses

Une ellipse est décrite par:

son demi-grand axe a,
son demi-petit axe b,

Quelques petites formules :

l'excentricité e = sqrt(1 - b²/a²)
l'aphélie = (1 + e)*a
le périhélie = (1 - e)*a
la vitesse de la comète sur son orbite elliptique v = sqrt(K*(2/r - 1/a))
la période T de la comète répond à T² = a³ dans notre système d'unités.

Nous voilà prêt pour écrire notre programme simulation.

Simulation de la trajectoire d'une comète

Le programme

Le programme complet en C, qui utilise la librairie graphique Dislin pour le tracé graphique, présente quelques particularités. Son code source Comete.cpp est téléchargeable.
Les grandes composantes du programme sont communes à tous mes programmes. J'utilise mon schéma RK4 pour l'intégration du système différentiel, qui est décrit comme d'habitude dans la routine Derivee:

//**********************************************

// Description du système différentiel pour RK4

//**********************************************

void Derivee(double y[], double dydt[])

{

double d;

d = pow((y[0]*y[0]+y[1]*y[1]),1.5);

dydt[0] = y[2]; // dx/dt

dydt[1] = y[3]; // dy/dt

dydt[2] = -K*y[0]/d; // d²x/dt² = -K*x/(x²+y²)^3/2

dydt[3] = -K*y[1]/d; // d²y/dt² = -K*y/(x²+y²)^3/2

}

Il comprend toutefois quelques nouveautés dans l'utilisation de la librairie dislin, comme le tracé de deux graphiques sur la même page et le tracé de la légende des 3 courbes d'évolution des énergies.
Les fonctions Min et Max pourront également être réutilisées, en changeant ou pas le type des variables.

Résultats de simulation et commentaires

Voici un premier exemple de simulation. Il s'agit de la simulation d'orbite de la comète Encke dont les caractéristiques principales sont:

périhélie : 0.339 UA
aphélie: 4.09 UA
période: 3.3 UA
excentricité: 0.847

Le programme demande 3 valeurs : l'aphélie, la vitesse à l'aphélie et la période. L'aphélie et la période sont des données que l'on trouve facilement dans les tables. Ce n'est pas le cas de la vitesse à l'aphélie que l'on devra calculer au préalable par la formule trouvée ci dessus.
Je lance donc le programme comete avec les données:

et j'obtiens les graphiques suivants:

Le référentiel est centré sur le Soleil.
Le mouvement commence à t=0 en x = r0, soit à l'aphélie de la comète, c'est à dire à sa plus grande distance du Soleil. Le calcul de l'orbite se fait dans le sens trigonométrique. Sans surprise, l'orbite obtenue est une ellipse. Les axes sont gradués en UA.
Le diagramme d'évolution de l'énergie est gradué en abscisse en % de la période et en ordonnées en UA d'énergie.
Notez les données issues de la simulation (imprimées dans la fenêtre de commandes): la valeur de a, le demi-grand axe, b le demi-petit axe et l'excentricité e de l'orbite. Pour cette dernière, la valeur issue de la simulation est de 0.846. Les tables astronomiques donnent une valeur comprise selon les sources entre 0.847 et 0.849. Le programme n'est pas mauvais !
Vous pourrez également calculer que la relation entre la période et le demi-grand axe (T² = a³) est très correctement vérifiée. Regardons maintenant de plus près les tracés de l'orbite et des énergies.
Nous remarquons que l'orbite est fermée, ce qui dénote d'une bonne stabilité du calcul et du schéma numérique. Pour vous en convaincre, modifiez le problème en utilisant la méthode d'Euler au lieu du RK4. Vous noterez une instabilité des résultats qui se traduiront par une orbite ouverte et un diagramme d'énergie incohérent. Vous pouvez aussi essayer de changer le pas de calcul deltat, en le diminuant, pour vérifier l'influence sur le résultat.
Concernant l'évolution de l'énergie de la comète, remarquons que:

l'énergie totale de la comète est constante et négative. Elle est constante parce que la force en 1/r² est conservatrice et aussi parce que le schéma est stable! Elle est négative car l'orbite est fermée. C'est une caractéristique de ce type d'orbite!
l'énergie cinétique présente un maximum au périhélie et un minimum à l'aphélie. C'est bien sur l'inverse pour l'énergie potentielle, ce qui est assez compréhensible....

Vous pouvez à votre tour tracer les orbites d'autres comètes. Il suffit de trouver sur le net leurs caractéristiques...