jeudi 30 août 2012

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.

9 commentaires:

  1. Pas mal! Bravo. J'avais un satellite à lancer justement ;)

    RépondreSupprimer
  2. désolé, moi je ne lance que des comètes :-)

    RépondreSupprimer
  3. Bonjour,

    Tout d'abord, un grand merci pour avoir réalisé ce travail car il m'a été utile pour développer un petit programme. Je me suis permis de modifier un petit peu le code pour suivre l’évolution d’un satellite qui s’insérerait en orbite géostationnaire. J’ai donc pris en compte les différentes mises à feu à des instants donnés et j’obtiens un résultat satisfaisant pour le novice de Scilab que je suis. Sans votre code comme support, je n’y serais jamais parvenu.
    Toutefois, j’aimerais savoir comment vous vous y prenez pour passer de la 2D à la 3D.
    Dans l’attente de votre réponse,

    Un lecteur intéressé par vos travaux.

    RépondreSupprimer
  4. Merci pour le commentaire diluxe, j'espère pouvoir jeter un oeil à ton code ! Le passage de la 2D à la 3D se fait assez naturellement en passant de 2 coordonnées d'espaces x,y à 3 x,y,z. Par exemple les données de départ du système :
    u0=[x1; y1; vx1; vy1; x2; y2; vx2; vy2];
    vont devenir
    u0=[x1; y1; z1; vx1; vy1; vz1; x2; y2; z2; vx2; vy2; vz2];
    ce qui va modifier la manière d'écrire la fonction [du]=force(t,u,masse0,masse1)
    u1=[u(1),u(2),u(3)]
    du1=[u(4),u(5),u(6)]
    ....
    tu peux trouver un exemple plus travaillé de simulation 3D dans mon livre sur scilab , je les remis sur mon blog dans le billet :
    http://rouxph.blogspot.fr/2013/07/le-livre-scilab-de-la-theorie-la.html
    si tu veux le code complet de l'exemple est accessible en ligne librement dans la première partie de mon livre :
    http://www.d-booker.net/HTML/scilab/base/intro/index.html
    si tu veux plus d'explications c'est décrit en détails dans la partie 4 de mon livre sur la réalisation d'animations.

    Bonne lecture !

    RépondreSupprimer
  5. Bonjour,
    Je vais vous faire parvenir mon code que n'est jamais que du "trifouillage" de matrices. Son unique but était de voir comment s'insère un satellite en orbite géostationnaire mais il mériterait une amélioration de son ergonomie car il manque cruellement de flexibilité. Dans l'idéal, il faudrait que je parvienne à passer le tout en 3D et que le satellite puisse fournir une poussée en continu ( comme c'est le cas pour les satellites à propultioon électrique). Pour l'heure je n'ai réussi qu'à imposer une série de mises à feu à chaque fois que le satellite passe par son apogée.
    Quoi qu'il en soit je vous remercie vivement pour vos conseils!


    PS: le code arrive dans un autre post

    RépondreSupprimer
  6. Voilà le code!
    Bonne chance pour comprendre quelque chose dans ce fouilli où je me perds moi-même
    Si ça vous intéresse et que vous ne comprenez pas ce que j'ai voulu faire, n'hésitez pas à me contacter ou à me donner des conseils pour l'améliorer.

    Merci pour tout!

    le lien : http://dl.free.fr/vnvTm2TXl

    RépondreSupprimer
  7. Merci pour ce post. Je l'ai trouvé très intéressant.

    RépondreSupprimer
  8. Bonjour, je vous remercie énormément pour ce billet plus que très intéressant !
    Je voudrais néanmoins savoir si c'était possible que vous précisiez quelles ont été les modifications apportées aux conditions initiales qui vous ont permises d'obtenir un comportement de la comète similaire à un satellite tel qu'on peut le voir dans l'image animée à la fin de votre article.
    Merci d'avance :-)

    RépondreSupprimer
  9. désolé pour la réponse tardive! J'ai retrouvé les paramètres de cette simulation :

    //constantes
    G=0.04;
    m0=1000;
    m1=1;
    dt=0.002;
    T=20;
    //coordonnées de planète
    x1=5;y1=0;vx1=0;vy1=2.5;
    //coordonnées du satelltite
    x2=5.05;y2=0;vx2=0;vy2=3.5;


    les paramètres pour la planète principale sont choisis pour avoir une orbite circulaire (3ème loi de Kepler). Pour le satellite on choisit une position légèrement décalée vers l'extérieur, avec une vitesse parallèle à celle de la planète et un peu supérieure de sorte qu'elle double la planète comme ça doit se passer pour le système Terre-Lune et effectivement les lois de la physique font que ça se passe ainsi ..

    RépondreSupprimer

Pour écrire des formules mathématiques vous pouvez utiliser la syntaxe latex en mettant vos formules entre des "dollars" $ \$....\$ $ par exemple :
- $\sum_{n=1}^\infty {1\over n^2}={\pi^2\over 6}$ s'obtient avec \sum_{n=1}^\infty {1\over n^2}={\pi^2\over 6}
- $\mathbb R$ s'obtient avec {\mathbb R} et $\mathcal D$ s'obtient avec {\mathcal D}
- pour les crochets $\langle .,. \rangle$ dans les commentaires utilisez \langle .,. \rangle
vous pouvez écrire du html dans les commentaires :
- italique <i> ... </i> gras <b> ... </b>
- lien <a href="http://adresse "> .... </a>