Projet: Evolution d'une galaxie

[Code source courant (23/11/09), version 0.3]

[Code source ancien (23/11/09), version 0.2]

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)

Conditions initiales

Le choix des conditions initiales, c'est à dire des positions et vecteurs vitesse initiales des étoiles, est un problème en soi. On commence par un cas simple en répartissant uniformément les étoiles dans une sphère.

Répartition uniforme dans une sphère

Voici une petite macro qui répartit aléatoirement et uniformément des étoiles dans une sphère :

{
   double x,y,z;
   double r,theta,phi;

   double PI = 4*atan(1);

   h = new TH3D("h","hhh",30,-1.2,1.2,30,-1.2,1.2,30,0,1.2);

   for (int i=0;i<1000000;i++) {
      r = pow(gRandom->Uniform(),1./3.);
      theta = acos(1-gRandom->Uniform()*2);
      phi = 2*PI*gRandom->Uniform();

      x = r*sin(theta)*cos(phi);
      y = r*sin(theta)*sin(phi);
      z = r*cos(theta);

      h->Fill(x,y,z);
   }

   h->Draw("");
   h->SetFillColor(kRed);
}

En partant de cette macro, essayez de comprendre pourquoi on obtient bien une distribution uniforme et démontrez-le.

Indication : on reviendra au cours de PHYS504 de l'année dernière. On part de trois variables aléatoires uniformes entre 0 et 1. Appelons-les $(A,B,C)$. On veut une distribution uniforme dans des coordonnées sphériques $(R,\Theta ,\Phi )$. Les éléments de volume donc les probabilités correspondantes doivent être les mêmes :

$dA.dB.dC=\left({\frac {3}{4\pi }}\right)r^{2}\sin \theta dr.d\theta .d\phi$

Le facteur ${\frac {3}{4\pi }}$ est un facteur de normalisation correspondant au volume d'une sphère de diamètre 1. Rappel : la somme des probabilités intégrées doit être égale à 1...

Au travail !

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

          La section suivante vise à rapporter de manière détaillée différentes méthodes de résolution numériques d'équation différentielles ordinaires (ODE) existantes. Cette liste n'est bien sûr pas exhaustive.

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\left(y_{i},t_{i}\right)$

(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}+\Delta t*f\left(y_{i+{\frac {1}{2}}}^{a},t_{i+{\frac {1}{2}}}\right)$

(II.A_1)

On rappelle la formule de Taylor :

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

On estime la valeur y ai+1/2 par :

$y_{i+1/2}^{a}=y_{i}+{\frac {\Delta t}{2}}f(y_{i},t_{i})$

(II.A_2)

D’autre part, on a évidemment :

$t_{i+{\frac {1}{2}}}=t_{i+{\frac {\Delta t}{2}}}$

(II.A_3)

Alors :

$y_{i+1}=y_{i}+\Delta t*f\left(y_{i+{\frac {1}{2}}}^{a},t_{i+{\frac {1}{2}}}\right)$

(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 y_i+1/2 et celle permettant de trouver y_i+1). Cette méthode consiste en fait à évaluer plusieurs fois la valeur de f(y_i+1/2) entre y_i et y_i+1 : trois fois pour être précis. On part de la relation suivante :

$y_{i+1}=y_{i}+{\frac {\Delta t}{6}}\left[f(y_{i},t_{i})+2f\left(y_{i+{\frac {1}{2}}}^{a},t_{i+{\frac {1}{2}}}\right)+2f\left(y_{i+{\frac {1}{2}}}^{b},t_{i+{\frac {1}{2}}}\right)+f\left(y_{i+1}^{c},t_{i+1}\right)\right]$

(II.B_1)

Dans laquelle :

${\color {red}y_{i+{\frac {1}{2}}}^{a}}=y_{i}+{\frac {\Delta t}{2}}f(y_{i},t_{i})$

${\color {green}y_{i+{\frac {1}{2}}}^{b}}=y_{i}+{\frac {\Delta t}{2}}f\left({\color {red}y_{i+{\frac {1}{2}}}^{a}},t_{i+{\frac {1}{2}}}\right)$

$y_{i+{\frac {1}{2}}}^{c}=y_{i}+{\frac {\Delta t}{2}}f\left({\color {green}y_{i+{\frac {1}{2}}}^{b}},t_{i+{\frac {1}{2}}}\right)$

(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+\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}}}+...$

(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)={\frac {dy}{dt}}$, on peut réécrire immédiatement le développement de Taylor de la façon suivante :

$y(t+\Delta t)=y(t)+\Delta tf_{i}+{\frac {\Delta t^{2}}{2!}}f_{i}^{1}+{\frac {\Delta t^{3}}{3!}}f_{i}^{2}+{\frac {\Delta t^{4}}{4!}}f_{i}^{3}+...$

(III.A_1)

où f_i⁽ⁿ⁾ désigne la dérivée n-ième de f(y_i,t).

            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(t+\Delta t)=y(t)+\Delta tf_{i}+{\frac {\Delta t^{2}}{2!}}f_{i}^{1}+{\frac {\Delta t^{3}}{3!}}f_{i}^{2}+O(t)^{3}$

(III.A.1_1)

Or, on obtient la dérivée première d’une fonction par :

$f^{'}(t)=\lim _{\Delta t\rightarrow 0}{\frac {f(t+\Delta t)-f(t)}{\Delta t}}\,$

(III.A.1_2)

On peut donc écrire :

$f_{i}^{(1)}={\frac {f_{i}-f_{i-1}}{\Delta t}}$

(III.A.1_3)

où fi = f(y_i, t_i) et f_i-1 = f(y_i-1, t_i-1). On peut alors calculer y_i+1 par :

$y_{i+1}=y_{i}+\Delta t*f_{i}+{\frac {\Delta t^{2}}{2!}}{\frac {f_{i}-f_{i-1}}{\Delta t}}$

$\Leftrightarrow y_{i+1}=y_{i}+\Delta t{\frac {f_{i}}{2}}+{\frac {\Delta t^{2}}{2}}{\frac {f_{i}-f_{i-1}}{\Delta t}}$

$\Leftrightarrow y_{i+1}=y_{i}+{\frac {\Delta t}{2}}(3f_{i}-f_{i-1})$

(III.A.1_4)

Remarquons au passage que le terme y₁ 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 y_i+1 sous forme d’une combinaison linéaire de f_i-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 le fait la relation (III.A.1_3). La relation utilisée est la suivante :

$y_{i+1}=y_{i}+\Delta t\sum _{k=0}^{n}\beta _{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 des formules allant jusqu'à l’ordre 6.

$k\rightarrow$ $n\downarrow$	0	1	2	3	4	5	Ordre de la méthode
0	1						1
1	${\frac {3}{2}}$	$-{\frac {1}{2}}$					2
2	${\frac {23}{12}}$	$-{\frac {16}{12}}$	${\frac {5}{12}}$				3
3	${\frac {55}{24}}$	$-{\frac {59}{24}}$	${\frac {37}{24}}$	$-{\frac {9}{24}}$			4
4	${\frac {1901}{720}}$	$-{\frac {2774}{720}}$	${\frac {2616}{720}}$	$-{\frac {1274}{720}}$	${\frac {251}{720}}$		5
5	${\frac {4277}{1440}}$	$-{\frac {7923}{1440}}$	${\frac {9982}{1440}}$	$-{\frac {7298}{140}}$	${\frac {2877}{1440}}$	$-{\frac {475}{1440}}$	6

        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 y₁ par une méthode autre. De même, pour une méthode d’ordre n, il faut évaluer les y_k, 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.

        A. Principe.

          Cette méthode est similaire à la précédente, sauf que l’on utilise le DL de Taylor à ‘l’envers’. En effet, on peut toujours écrire $y(t)=y(t+\Delta t-\Delta t)$(curieux non??). Dès lors, le DL de Taylor s'écrit :

$y(t)=y(t+\Delta t)-\Delta tf+{\frac {\Delta t^{2}}{2!}}f^{1}-{\frac {\Delta t^{3}}{3!}}f^{2}+{\frac {\Delta t^{4}}{4!}}f^{3}+...+(-1)^{n}{\frac {\Delta t}{n!}}f^{n-1}$

(VI.A_1)

D'où l'écriture simplifiée suivante de y_i+1 en fonction de y_i:

$y_{i+1}=y_{i}-\Delta tf_{i+1}+{\frac {\Delta t^{2}}{2!}}f_{i+1}^{1}-{\frac {\Delta t^{3}}{3!}}f_{i+1}^{2}+{\frac {\Delta t^{4}}{4!}}f_{i+1}^{3}+...+(-1)^{n}{\frac {\Delta t}{n!}}f_{i+1}^{n-1}$

(VI.A_2)

Il est donc évident que, pour connaître y_i+1, il faut connaître y_i+1! On doit donc estimer une première valeur de y_i+1 que l'on notera y⁽⁰⁾_i+1.

$y_{i+1}^{(1)}=y_{i}-\Delta tf_{i+1}\left(y_{i+1}^{(0)},t\right)+{\frac {\Delta t^{2}}{2!}}f_{i+1}^{1}\left(y_{i+1}^{(0)},t\right)-{\frac {\Delta t^{3}}{3!}}f_{i+1}^{2}\left(y_{i+1}^{(0)},t\right)+...+(-1)^{n}{\frac {\Delta t}{n!}}f_{i+1}^{n-1}\left(y_{i+1}^{(0)},t\right)$

(VI.A_3)

Par la formule (III.A.1_3) et la méthode qui lui fait suite, on peut réécrire le développement précédent :

$y_{i+1}=y_{i}+\Delta t\sum _{k=0}^{n}\alpha _{nk}f_{i+1-k}$

(VI.A_3)

Cette méthode est dite itérative car elle s'effectue par approximations successives de la valeur réelle visée. En effet, elle suppose également l'utilisation d'un contrôle de la précision pour interrompre l'itération. Cela nécessite donc la définition d'une précision ε telle que quand $\epsilon \geq {\begin{vmatrix}{y_{i+1}^{(r+1)}-y_{i+1}^{(r)}}\end{vmatrix}}$.
Les valeur des α_nk sont données dans le tableau suivant pour de développement allant jusqu'à l'ordre 6.

$k\rightarrow$ $n\downarrow$	0	1	2	3	4	5	Ordre de la méthode
0	1						1
1	${\frac {1}{2}}$	${\frac {1}{2}}$					2
2	${\frac {5}{12}}$	${\frac {8}{12}}$	$-{\frac {1}{12}}$				3
3	${\frac {9}{24}}$	${\frac {19}{24}}$	$-{\frac {5}{24}}$	${\frac {1}{24}}$			4
4	${\frac {251}{720}}$	${\frac {646}{720}}$	$-{\frac {264}{720}}$	${\frac {106}{720}}$	$-{\frac {19}{720}}$		5
5	${\frac {1475}{1440}}$	${\frac {1427}{1440}}$	$-{\frac {798}{1440}}$	${\frac {482}{140}}$	$-{\frac {173}{1440}}$	${\frac {27}{1440}}$	6

        B. Avantages et inconvénients.

          Cette méthode itérative a l'avantage d'être plus précise et plus stable que la méthode ouverte, à ordre égaux. Néanmoins, à cause du calcul de convergence qu'elle comporte, elle demande plus de temps que la méthode ouverte. De même que dans la méthode précédente, il faut pouvoir évaluer les yk, k variant de 1 à n-1. D'autre part, à chaque nouvelle évaluation, il faut estimer une première valeur de y_i+1, ce que cette méthode ne permet pas : il faut en utiliser une autre. La section suivante est très adaptée à ce problème.

   V. Méthode de Prédicteur-Correcteur.

        A. Principe.

          La méthode de prédicteur-correcteur est basée sur la combinaison des méthodes d'Adams ouverte et fermée. Le prédicteur du terme y⁽⁰⁾_i+1 utilisé dans la méthode fermée, à l'ordre n, est évalué à l'aide d'une méthode ouverte, à l'ordre n également.
Explicitons nos propos. Le terme y⁽⁰⁾_i+1 prend la forme suivante :

${\color {red}y_{i+1}^{(0)}}=y_{i}+\Delta t\sum _{k=0}^{n}\beta _{nk}f_{i-k}$

(V.A_1)

C'est le prédicteur.

Ce terme est injecté dans la méthode fermée et l'itération est amorcée. Du coup, on obtient :

$y_{i+1}^{(1)}=y_{i}+\Delta t\left(\alpha _{n0}f({\color {red}y_{i+1}^{(0)}},t_{i+1})+\sum _{k=1}^{n}\alpha _{nk}f_{i+1-k}\right)$

(V.A_2)

C'est le correcteur.
L'itération de la méthode fermée est alors lancée jusqu'à ce que la précision requise est atteinte.

        B. Avantages et inconvénients.

          Tout comme dans toutes les méthodes d'Adams, il faut évaluer les yk, k variant de 1 à n-1 par une méthode différente pour initialiser le calcul.
Néanmoins, par cette approche on peut accélérer la convergence de la méthode fermée sans perdre en précision.

   VI. Méthode de Tir.

         

Article à venir.

   VII. Méthode de Newmark.

         

Article à venir.

   VIII. Méthode matricielle.

         

Article à venir.

   IX. Sources.

          Voilà les sources qui nous ont servies à rédiger cette section.

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