mécanique céleste avec scilab

Le transit de Vénus de juin dernier  m'a redonné envie de m'intéresser à la mécanique céleste. Le sujet est mathématiquement difficile car il est en général impossible de résoudre les systèmes d'équations qu'on y rencontre, mis à part le cas extrêmement simple du problème à deux corps, quoique même dans ce cas la solution n'est pas aussi simple que cela à exprimer analytiquement. Par contre les problèmes de mécanique céleste sont très amusants à résoudre numériquement, surtout si l'on dispose d'un logiciel de calcul efficace comme Scilab! Dans la pratique une fois les équations posées il suffit d'utiliser la fonction ode pour résoudre l'équation et afficher la trajectoire aussi complexe soit-elle, voici un exemple de problème à trois corps résolu numériquement :

1 modélisation du problème

 Si on veut simuler l'évolution gravitationnel d'un système composé de N corps (numérotés de i=1 jusqu'à N) le but est de calculer leurs positions $X_i(t)$ au cours du temps (t). Le principe fondamental de la dynamique permet d'écrire les équations que vérifient les positions, l'accélération ${d^2\over dt^2}X_i $ subie par le $i^\text{ième}$ corps est la somme des forces gravitationnelles exercées par les N-1 autres corps :

$$ m_i{d^2\over dt^2}X_i(t) = \sum_{j\neq i} -{{\cal G} m_i m_j\over \vert X_i(t)-X_j(t)\vert^3}(X_i(t)-X_j(t))=m_i F_{G,i}(X_1,\dots,X_N)~~~~i=1,\dots, N$$

il est évident que les forces gravitationnelles exercées par les corps ayant la plus grosse masse $m_j$ vont avoir l'effet le plus important, et que dans certains cas on peut négliger l'effet des corps les plus légers sur les corps les plus massifs pour simplifier les équations. Dans le cas du système solaire  il y a

• un corps principal le soleil dont la masse est très supérieure aux autres corps
• des corps de masse intermédiaire les planètes
• des corps négligeables (astéroïdes, comètes, satellites, ...) 

du fait de sa masse, le soleil  est peu influencé par les autres corps, on peut donc le considérer comme fixe (à l'origine des coordonnées). Pour les planètes, elles subissent l’attraction du soleil mais on peut négliger l'action d'autres corps plus petits. Par contre pour calculer la trajectoire des comètes, il faudra tenir compte de l'action du soleil et des différentes planètes dans la force gravitationnelle.

 Pour résoudre numériquement un tel problème il faut ensuite le mettre sous la forme d'une équation différentielle d'ordre 1:
$$ {d\over dt}U =F(U,t)$$
 C'est possible en faisant intervenir la vitesse ${d\over dt} X_i$ des différents objets comme une variable à part entière :
$$U=\begin{pmatrix}X_1\\ {d\over dt}X_1\\ \vdots \\ X_N\\ {d\over dt}X_N \end{pmatrix}
~~et~~
F(U,t)=
\begin{pmatrix}U_2\\ F_{G,1}(U_1,\dots,U_{2N-1})\\ \vdots \\ U_{2N}\\ F_{G,N}(U_1,\dots,U_{2N-1})\ \end{pmatrix}$$

car pour $i=1,\dots, N$ on a
$$ {d\over dt}U_{2i-1}= {d\over dt}X_i=U_{2i}$$
 et
$${d\over dt}U_{2i}= {d^2\over dt^2}X_i=F_{G,i}(X_1,\dots,X_N)$$

2 résolution numérique

Pour résoudre numériquement une équation différentielle d'ordre un  écrite sous la forme $ {d\over dt}U =F(U,t)$ il existe dans la plupart des langages informatiques une fonction "ode" qui calcule automatiquement la liste des valeurs prise par la solution U(t) à certaines dates t ! Si vous utilisez Scilab la fonction "ode" s'utilise comme ceci  U=ode(U0,t0,t,F)  avec

• U= liste des valeurs de la solution 
• t0= date du début de la simulation
• U0=valeur de U(t) à t=t0 
• t= listes des dates auxquelles calculer la valeur de U(t) 
• F= la fonction F(U,t)

pour résoudre le problème il faut donc commencer par programmer la fonction calculant les forces gravitationnelles de chaque corps :

  
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

je mesuis placé dans le cadre d'un mouvement plan (u est à deux dimensions) mais on pourrait facilement rajouter une troisième dimension. Il reste à programmer la fonction F :

 function [du]=force(t,u,masse0,masse1)
         u1=[u(1);u(2)]
        du1=[u(3);u(4)]
         u2=[u(5);u(6)]
        du2=[u(7);u(8)]
       ddu1=force_g(t,u1,masse0)
       ddu2=force_g(t,u2,masse0)+force_g(t,u2-u1,masse1)
         du=[du1;ddu1;du2;ddu2]
endfunction

ici je me suis restreint au cas simple d'une planète géante (disons Jupiter) de masse m1 et d'une comète (de masse négligeable) le tout en orbite autour du soleil (de masse m0). Pour obtenir une animation il faut maintenant fixer les valeurs numériques des différents objets (positions, vitesses, masses) et récupérer les résultats de calculs d'ode puis les afficher. C'est avec ce script que j'ai obtenu l'animation au début de ce billet :

  
//constantes 
G=0.04;//constante gravitationnelle en unités astronomique*
m0=1000;//masse soleil
m1=1;//masse Jupiter
dt=0.05;// intervalle entre les différentes dates
T=500;//durée de la simulation
dx=0.5;dy=0.5;// taille du dessin des planètes
//coordonnées de la planète au début 
x1=5;y1=0;vx1=0;vy1=2.5;
//coordonnées de la comète  au début 
x2=6;y2=6;vx2=-2;vy2=-0.5;
// valeur initiale de U 
u0=[x1; y1; vx1; vy1; x2; y2; vx2; vy2];
//dates de calcul de U
t=[0:dt:T];
//calcul de la solution 
u=ode(u0,0,t,list(force,m0,m1));
//récupération des coordonnées correspondant à la position (X,Y) des deux objets
X=[u(1,:)',u(5,:)'];Y=[u(2,:)',u(6,:)'];

//animation graphique
plot2d(X,Y,[5 2],rect=[-10,-10,10,10],frameflag=4)
rect=[min(X),min(Y),max(X),max(Y)]
N=length(t);
plot2d(0,0,-9,rect=rect,frameflag=3)
f=gcf()
winnum=winsid() ;
toolbar(winnum(1),'off');
f.pixmap="on"
for k=2:N
   drawlater()
   clf()
   plot2d(0,0,-9,rect=rect,frameflag=3)
   plot2d(X(1:k,:),Y(1:k,:),[5 2],rect=rect,frameflag=0)
   P1=[X(k,1)-dx/2;Y(k,1)+dy/2;dx;dy;0;360*64];//la planète en rouge
   P2=[X(k,2)-dx/2;Y(k,2)+dy/2;dx;dy;0;360*64];//la comète en bleu
   xfarcs([P1,P2],[5,2])
   xstring(-4,-6,'t='+string(t(k)))
   drawnow()
end

pour les valeurs numériques il est plus simple de travailler en unités "astronomiques" à savoir :

• unité de distance 1ua=150 millions de kilomètres
• unité de temps 1 an=365 jours
• unité de masse masse de jupiter = $1.9 \times 10^{27}$ kg

ce qui donne pour la constante gravitationnelle une valeur d'environ G=0.04 .
3 Quelques exemples 

Après on peut s'amuser à modifier les paramètres des planètes ou le  nombre de planètes pour obtenir d'autres animation. Par exemple si au départ  la comète a une position et une vitesse initiale proche de la planète on va obtenir qu'elle se satellise autour de la planète principale. Cela permet de visualiser les perturbations auquelles sont soumises les orbites des satellites comme la lune autour de la terre :

autre exemple en ajoutant une dimension et en initialisant les vitesses de la planètes et de la comète pour qu'elles ne soient pas dans le même plan , on peut même réaliser des animations en 3D :

j'ai amélioré la visualisation de la trajectoire en 3D en faisant apparaitre le plan orbital de la planète et en faisant tourner le point d'observation.