Projet: Evolution d'une galaxie

On veut simuler l'évolution d'une galaxie, en particulier sa structure, au cours du temps. Il s'agit en fait d'un problème à N corps, N étant grand (quelques centaines pour commencer). Les notions de C++ vues en cours seront (ré)introduites et utilisées au fur et à mesure. On a d'abord scindé le problème en parties, chacune étant traitée par un binôme.

Structure générale du code et définition des classes de base

Les classes de base à considérer sont :

• la classe "UEtoile" définissant une étoile. Pour le moment, on y inclut la masse, la position et la vitesse de l'étoile. Un premier squelette de classe peut être trouvé ici : [UEtoile.h] et [UEtoile.C]. Utilisez-le (voir son utilisation ci-dessous).
• la classe "UGalaxie" qui reste à écrire, représentant un ensemble d'étoiles, pas forcément une galaxie (amas globulaire, etc...)

Un moteur "EDOMoteur", dont nous ne sommes pas encore persuadé qu'il doit s'écrire sous forme de classe (c'est probable). Ce moteur prendrait les objets d'une galaxie un à un et les fera évoluer dans le temps. Une méthode de Runge-Kutta simple sera implémentée au début, que l'on fera évoluer vers une méthode avec pas de temps adaptatif.

Conventions

On adoptera une convention de noms pour les classes et variables membres des classes, ainsi que pour les méthodes :

• tous les noms de classe commenceront par un U majuscule suivi d'une majuscule
  • • Exemple : UEtoile, UGalaxie, UEDOMoteur
• tous les noms de variable membres d'une classe commenceront par un "m" minuscule suivi d'une majuscule
  • • Exemple : mX, mMasse
• tous les noms de méthode de classe commenceront par une majuscule suivie de minuscules, sauf pour marque le début d'un mot.
  • • Exemple : Trace() ou TraceUnCarre()
• on verra au fil du temps pour les variables internes aux méthodes

Utilisation de la classe "UEtoile"

Charger [UEtoile.h] et [UEtoile.C] dans ROOT avec :

root[1] .L UEtoile.C

Ensuite, on définit une étoile de masse 12 (unités à définir) avec

root[1] e = new UEtoile(1.2,0.12,1.2345,0,0,0,12)

On peut afficher les valeurs de position masse et vitesse de l'objet e avec

root[1] e->Print()

Et on peut récupérer la valeur d'une variable avec le "Getter" et la modifier avec le "Setter" correspondant :

root[1] e->GetX()
root[1] e->SetX(1.5)
root[1] e->Print()

On peut déclarer et définir un ensemble d'étoiles dans un tableau de pointeurs (si il y en a peu, ça peut aller)

UEtoile* e[10];
for (int i=0;i<10;i++) {
   e[i] = new UEtoile(0,0,0,0,0,0,10);
}

Attention, ce faisant, on a défini toutes les étoiles au même endroit (0,0,0) et avec la même vitesse (0,0,0). Il faut évidemment définir position et vitesses de façon aléatoire (voir la classe TRandom).

On peut ensuite les utiliser avec e[i]. Par exemple :

e[3]->Print();

ou

e[j]->SetX(0.56);

(j est défini par ailleurs)

Affichage graphique

En attente de tests. Premier test : définir quelques étoiles et tenter de les tracer avec les outils ROOT.

Moteur de résolution d'équations différentielles

Méthodes numériques de résolutions d'équations différentielles

   I. Méthode d'Euler

          Cette méthode est la plus simple mais, en même temps, la moins précise. Il peut néanmoins être utile de la détailler car certaines méthodes de résolution numériques ont pour base le calcul d’Euler, en y ajoutant certaines modifications (nous verrons cela plus loin).

        A. Principe.

          Une fonction y quelconque, par exemple la position d’un objet en fonction du temps, admet le développement en série de Taylor suivant :

$y(t+\Delta t)=y(t)+\Delta t{\frac {dy}{dt}}+{\frac {\Delta t^{2}}{2!}}{\frac {d^{2}y}{dt^{2}}}+{\frac {\Delta t^{3}}{3!}}{\frac {d^{3}y}{dt^{3}}}+...$

(I.A_1)

en supposant y(0)=y₀ connu.

A l'ordre 1, on obtient :

$y(t+\Delta t)=y(t)+\Delta t{\frac {dy}{dt}}$

(I.A_2)

Pour Δt constant, on a :

$y(t+\Delta )_{i+1}=y_{i}+\Delta t*f(y_{i},t_{i}){\frac {}{}}$

(I.A_3)

La fonction y est connue à t=0, donc, la fonction y peut être déterminée à tout autre instant ultérieur.

        B. Avantages et inconvénients.

          Cette méthode est la plus simple à mettre en œuvre mais elle conduit très rapidement à des erreurs notamment parce qu’elle effectue l’intégration de y’(xi) à y’(xi+1) à partir de la valeur calculée en xi. Nous ne la retiendrons donc pas.

   II. Méthodes de Runge-Kutta.

          Les méthodes de Runge-Kutta sont basées sur la méthode d’Euler. Elles visent à approximer le résultat du calcul entre xi et xi+1 en estimant des valeurs intermédiaires sur cet intervalle.

      A. Méthode d’ordre 2.

Sur l’intervalle [xi ; xi+1], on calcule la valeur que prend la fonction au milieu de cet intervalle et on calcul la dérivée en ce point. Supposons yi connu. On centre le calcul sur le milieu de l’intervalle. On résoudra donc :

y_(i+1)=y_i+∆t*f(y_(i+1/2)^a,t_(i+1/2)) (II.A_1) On rappelle la formule de Taylor : y(t+∆t)=y(t)+∆t dy/dt+〖∆t〗^2/2! dy/dx+〖∆t〗^3/3! (d^2 y)/〖dx〗^2 +⋯

On estime la valeur y ai+1/2 par : y_(i+1/2)^a=y_i+∆t/2 f(y_i,t_i) (II.A_2) D’autre part, on a évidemment : t_(i+1/2)=t_i+∆t/2 (II.A_3) Alors : y_(i+1)=y_i+∆t*f(y_(i+1/2)^a,t_(i+1/2)) (II.A_4)

      B. Méthode d’ordre 4.

          Le principe est le même qu’à l’ordre deux (‘ordre deux’ parce que deux évaluations successives de f sont effectuées : celle pour calculer yi+1/2 et celle permettant de trouver yi+1). Cette méthode consiste en fait à évaluer plusieurs fois la valeur de f(yi+1/2) entre yi et yi+1 : trois fois pour être précis. On part de la relation suivante :

y_(i+1)=y_i+∆t/6 [f(y_i,t_i )+2f(y_(i+1/2)^a,t_(i+1/2) )+2f(y_(i+1/2)^b,t_(i+1/2) )+f(y_(i+1)^c,t_(i+1) ) ] (II.B_1) Dans laquelle : y_(i+1/2)^a=y_i+∆t/2 f(y_i,t_i ) y_(i+1/2)^b=y_i+∆t/2 f(y_(i+1/2)^a,t_(i+1/2)) y_(i+1/2)^c=y_i+∆t/2 f(y_(i+1/2)^b,t_(i+1/2)) (II.B_2)

      C. Avantages et inconvénients.

          Les méthodes ci-dessus présentent plusieurs avantages : facilité de programmation, stabilité de la solution, modification simple du pas et la connaissance de y(0) suffit pour intégrer l'équation différentielle. Plus l’ordre de la méthode est élevé, meilleure est la précision. La conséquence immédiate est l’augmentation des calculs à effectuer. Les inconvénients de cette méthode se résument donc au temps de calcul lent et à la difficulté de l'estimation de l'erreur locale. A noter néanmoins que la méthode de Runge-Kutta d’ordre 4 semble être assez largement utilisée en simulation de phénomènes physiques (gage d’efficacité ??).

   III. Méthodes d’Adams.   

          Les méthodes suivantes ont également pour base la méthode d’Euler, en particulier elles s’appuient sur le développement de Taylor y(t+∆t)=y(t)+∆t dy/dt+〖∆t〗^2/2! dy/dx+〖∆t〗^3/3! (d^2 y)/〖dx〗^2 +⋯ (III_1) en supposant y(0)=y0 connu. Il existe deux méthodes, ou formules d’Adams : les formules dites « ouvertes » ; les formules dites « fermées ».

      A. Formules d’Adams ouvertes.

          En posant f’(y,t)=dy/dt on peut réécrire immédiatement le développement de Taylor y_(i+1)=y_i+∆t*f_i+〖∆t〗^2/2! f_i^((1))+〖∆t〗^3/3! f_i^((2))+〖∆t〗^4/4! f_i^((3))+⋯ (III.A_1)

         1. Formules d’Adams ouvertes d’ordre 2.

          La résolution s’effectue en tenant compte du DL de Taylor jusqu’à l’ordre 2 (d’où le nom de la méthode). y_(i+1)=y_i+∆t*f_i+〖∆t〗^2/2! f_i^((1))+〖Ο(∆t)〗^3 (III.A.1_1) Or, on obtient la dérivée première d’une fonction par f^' (t)=lim┬(∆t→0)⁡〖(f(t+∆t)-f(t))/Δt〗 (III.A.1_2) On peut donc écrire f_i^((1) )=(f_i-f_(i-1))/Δt (III.A.1_3) où fi = f(yi, ti) et fi-1 = f(yi-1, ti-1). On peut alors calculer yi+1 par y_(i+1)=y_i+∆t*f_i+〖∆t〗^2/2! ((f_i-f_(i-1))/Δt) ⟺y_(i+1)=y_i+2∆t f_i/2+〖∆t〗^2/2 ((f_i-f_(i-1))/Δt) ⟺y_(i+1)=y_i+∆t/2 (〖3f〗_i-f_(i-1) ) (III.A.1_4) Remarquons au passage que le terme y1 ne peut pas être calculé par cette méthode car il suppose la connaissance du terme y-1, ce qui n’est pas le cas. Ce terme doit donc être évalué par une autre méthode.

         2. Formules d’Adams ouvertes d’ordre n.

          Généralisation de la méthode précédente, elle consiste à faire le calcul en prenant l’ordre du DL de Taylor à l’ordre n et à exprimer yi+1 sous forme d’une combinaison linéaire de fi-k, k variant de 0 à n. Pour cela, on utilise la même méthode que ci-dessus : on réinjecte dans le DL l’expression des n-1 dérivées à calculer, évaluées comme la relation (III.A.1_3) le fait. La relation utilisée est la suivante : y_(i+1)=y_i+∆t∑_(k=0)^n▒〖β_nk f_(i-k) 〗 (III.A.2_1) Les βnk sont des coefficients, calculables par une formule qui nous est inconnue, mais déterminables par le calcul explicite. Voici un tableau récapitulatif de ces coefficients pour une formule d’ordre 6.

(Source : http://www.sciences.univ-nantes.fr/physique/perso/aloui/m_numeri/61eqdiff/61eqdiff.htm#1)

      B. Avantages et inconvénients.

          A notre connaissance, il n’y a pas d’avantage particulier déclaré hormis la précision accrue devant la méthode d’Euler. Par contre, comme vu pour la méthode d’ordre 1, il faut évaluer le terme y1 par une méthode autre. De même, pour une méthode d’ordre n, il faut évaluer les yk, k variant de 1 à n-1, par une méthode alternative, à déterminer par ailleurs (Runge-Kutta, ou pourquoi pas une utilisation successive des formules d’Adams pour les évaluer et partir construire ainsi progressivement la formule au rang n que l’on souhaite…). La plus utilisée semble être la méthode d’ordre 4.

   IV. Formules d’Adams fermées.

          Cette méthode est similaire à la précédente, sauf que l’on utilise le DL de Taylor à ‘l’envers’."