Programmation Fonctionnelle, projet

Version & licenses
Creative Commons License

Programmation Fonctionnelle, projet : Problème des n corps.

Guyslain Naves

Information générale

Le projet comptera pour 12 points, décomposés comme suit :

  • 3 points pour pertinence de l'architecture
  • 3 points pour qualité du code
  • 4 points selon les fonctionnalités réalisées
  • 2 points pour la soutenance

Le projet sera réalisé en binôme.

Le planning du projet est :

  • 29 mars : début du projet,
  • 15 avril : rendu de l'architecture prévisionnelle,
  • 25 avril à midi : rendu global,
  • 29 avril : soutenance.

Détails des rendus :

  • architecture prévisionnelle : vous devrez rendre un rapport d'au plus deux pages décrivant l'ensemble des fichiers de code que vous comptez écrire, accompagné de l'ensemble des spécifications sous la forme de fichiers d'interface (d'extension .mli).
  • rendu global : il s'agit du code complet de votre projet, accompagné d'un rapport décrivant :
    • comment utiliser votre programme,
    • quelles sont les fonctionnalités de votre programme,
    • éventuellement les changements apportées à l'architecture.
  • soutenance : présentation de votre solution et démonstration en dix minutes, puis réponse aux questions.

Tous les rendus de fichiers se feront sur Ametice.

Description du projet

Le projet consiste en l'écriture d'un programme pouvant simuler le déplacement de corps célestes selon les lois de la gravité énoncées par Newton. Le programme doit pouvoir afficher la position de l'ensemble des corps à un moment donné, et calculer le mouvement de ces corps sur une petite échelle de temps.

Chaque corps dispose à tout instant de plusieurs caractéristiques :

  • sa masse,
  • sa position dans l'espace,
  • sa vitesse, donnée comme un vecteur.

La principale fonctionnalité du projet est une fonction qui prenant un ensemble de corps à un instant $t$ et une durée $\mathrm{d}t$, calcule la position et la vitesse de chaque corps à l'instant $t + \mathrm{d}t$. Les lois physiques sont rappelées plus bas dans cette page.

Le cœur algorithmique consiste en une structure de donnée, les quad-trees, permettant de calculer les nouvelles positions des $n$ corps en temps $O(n \log n)$, là où la méthode naïve prend un temps $\Theta(n^2)$.

Implémenter la version naïve en 2 dimensions serait le minimum acceptable. Implémenter les quad-trees constituerait un travail honorable. Travailler en 3 dimensions serait fortement appréciés. On peut aussi imaginer ajouter d'autres fonctionnalités, par exemple :

  • rentrer les données précises du système solaire,
  • générer des systèmes stellaires aléatoirement avec des objets en orbite autour d'une étoile,
  • améliorer la sortie graphique en donnant une couleur et un rayon à chaque corps,
  • rendre la sortie graphique interactive, en permettant à la caméra de se déplacer, de ralentir ou accélérer la simulation, …
  • affiner les solutions mathématiques en utilisant une méthode de Runge-Kutta.

Un peu de physique

En l'absence de tout autre corps, un corps céleste se déplace en ligne droite selon sa vitesse. Ainsi, un corps en position $\vec{u}$ à un instant $t$ et ayant un vitesse $\vec{v}$ et une masse $m$, sera à l'instant $t+\mathrm{d}t$ en position $\vec{u} + \mathrm{d}t \cdot \vec{v}$, s'il n'est soumis à aucune autre force. Il gardera aussi sa vitesse $\vec{v}$.

Si le même corps est plutôt soumis à une force constante $\vec{F}$, alors sa vitesse évolue avec le temps. À l'instant $t + \mathrm{d}t$, sa vitesse sera alors $\vec{v} + (\mathrm{d}t/m) \cdot \vec{F}$, et sa position $\vec{u} + \mathrm{d}t \cdot \vec{v} + (\mathrm{d}t^2/(2m)) \cdot \vec{F}$.

Jusque là c'est facile, mais malheureusement les forces qui s'appliquent à un corps sont rarement constantes. Si les forces ne sont plus constantes, il devient très difficile de calculer précisément la position et la vitesse du corps au cours du temps. C'est là que nous allons faire une première approximation : nous considérerons des valeurs $\mathrm{d}t$ telles que $\vec{F}$ ne changent presque pas entre $t$ et $t + \mathrm{d}t$. Ainsi nous pourrons simplement appliquer les formules ci-dessus pour passer de $t$ à $\mathrm{d}t$, puis calculer à nouveau la force d'exerçant sur chaque objet à l'instant $t + \mathrm{d}t$, et recommencer. De telle sorte que pour calculer l'évolution d'un système, nous le ferons simplement évoluer par petits pas de temps. La précision du calcul dépendra donc de la finesse du pas de temps $\mathrm{d}t$ choisi.

En physique Newtonienne, le calcul de la force exercée sur un corps se fait en sommant l'influence de tous les autres corps du système. Pour un corps $A$ en position $\vec{a}$ et de masse $m_A$ et un corps $B$ en position $\vec{b}$ et de masse $m_B$, la force exercée par $B$ sur $A$ est donnée par:

$$ \vec{F_{AB}} = g \frac{m_A \cdot m_B}{\left\|AB\right\|^3} \vec{AB}$$

avec $g = 6.67398 × 10^{-11} \text{m}^3 \text{kg}^{-1} \text{s}^{-2}$ la constante gravitationnelle, et $\left\|AB\right\|$ la distance entre $A$ et $B$.

Pour obtenir la force totale auquelle un corps est soumis, il suffit donc de sommer :

$$ \vec{F_A} = \sum_{B \in \mathcal{S} \setminus \{A\}} \vec{F_{AB}}$$

où $\mathcal{S}$ est l'ensemble des corps du système.

Quad-trees

Ainsi pour calculer l'ensemble des forces, du système, il faut calculer les interactions de chaque paire de corps, donc $\theta(n^2)$calculs élémentaires sont nécessaires. C'est trop pour pouvoir faire tourner une simulation avec plusieurs milliers de corps dans un temps raisonnable.

Nous utilisons alors une deuxième approximation: pour un corps donné $A$, on traite les groupes de corps les plus éloignés comme un seul corps regroupant toute la masse au barycentre du groupe. Par exemple, supposons que $B_1, B_2,\ldots,B_k$ soient des corps loin de $A$ mais suffisamment proche les uns des autres. Notons $m_i$ la masse de $B_i$ et $u_i$ sa position. Alors on calcule la masse totale :

$$M = \sum_{i=1}^k m_i$$

et le barycentre : $\vec{B} = (1/M) \sum_{i=1}^k m_i \vec{u_i}$puis, la force exercée par $B_1,\ldots,B_k$ sur $A$ est approximée par :

$$\vec{F} = g \frac{m_A \cdot M}{\left\|AB\right\|} \vec{AB}$$

Comment trouver les points à grouper ensemble ? Pour cela, nous codons l'ensemble des corps du système par un arbre d'arité 4 (chaque nœud interne a exactement 4 fils) appelé quad-tree. À chaque nœud est associé une région carrée de l'espace (un cube en dimension 3). La région associée à chaque feuille doit contenir au plus un corps. Ainsi, si en ajoutant un corps dans l'arbre, on doit l'ajouter à une feuille dont la région contient déjà un corps, il faudra découper cette région en 4 carrés plus petits et créer les nœuds correspondants dans l'arbre.

Voici l'exemple d'un quad-tree:

un quad-tree

À chaque nœud seront associées les informations suivantes :

  • s'il s'agit d'une feuille qui contient un corps, on lui associe le corps en question,
  • s'il s'agit d'un nœud interne, on lui associe la masse totale des corps contenus dans le sous-arbre, ainsi que le barycentre correspondant.

Un fois le quad-tree correctement construit, pour calculer la force exercée par le système sur un corps $A$, il suffit de parcourir le quadtree depuis la racine :

  • si on parcourt un noeud suffisamment loin du corps $A$ (distance entre le centre de la région et de $A$ supérieure à $\theta$ fois la longueur d'un côté de la région), alors on calcule seulement la force exercée par la masse totale du sous-arbre concentrée au barycentre.
  • si on parcourt un noeud plus proche, on calcule simplement la force exercée par chaque sous-arbre, et on les additionne.
  • si le nœud est une feuille, on calcule juste la force exercée par le corps contenu dans la feuille.

On pourra prendre $\theta = 2$.

Avec cette méthode, le nombre de calcul élémentaire pour obtenir la force exercée par un corps est $O(\log n)$, et donc pour tout le système $O(n \log n)$, si $n$ est le nombre de corps dans le système.