Simulation d'un système à N corps — Wikipédia

En physique et en astronomie, une simulation à N corps est une simulation d'un système dynamique de particules, généralement sous l'influence de forces physiques, telles que la gravité (voir problème à N corps pour d'autres applications). Les simulations à N corps sont des outils largement utilisés en astrophysique^[1], depuis l'étude de la dynamique de systèmes à quelques corps comme le système Terre - Lune - Soleil, jusqu'à la compréhension de l'évolution de la structure à grande échelle de l'univers observable^[2]. En cosmologie, les simulations à N corps sont utilisées pour étudier les processus de formation de structures non linéaires tels que les filaments de galaxies et les halos de galaxies sous l'influence de la matière noire^[3]. Des simulations directes à N corps sont utilisées pour étudier l'évolution dynamique des amas d'étoiles.

Les «  particules » traitées par la simulation peuvent correspondre ou non à des objets physiques. Par exemple, dans une simulation à N corps d'un amas d'étoiles chaque particule peut avoir une signification physique et correspondre à une étoile. Mais dans une simulation d'un nuage interstellaire, il n'est pas possible d'avoir une particule pour chaque atome ou molécule de gaz, car cela nécessiterait de l'ordre de 10²³ particules pour chaque mole de matière (voir constante d'Avogadro) ; donc une seule « particule » représentera une quantité beaucoup plus grande de gaz (souvent mise en œuvre à l'aide de l'hydrodynamique des particules lissées ). Cette quantité n'a pas besoin d'avoir une signification physique, mais doit être choisie comme un compromis entre la précision et les ressources informatiques disponibles.

Dans les simulations gravitationnelles directes à N corps, les équations du mouvement d'un système de N particules sous l'influence de leurs forces gravitationnelles mutuelles sont intégrées numériquement sans aucune approximation simplificatrice. Ces calculs sont utilisés dans des situations où les interactions entre des objets individuels, tels que des étoiles ou des planètes, sont importantes pour l'évolution du système.

Les premières simulations gravitationnelles directes d'un système à N corps ont été réalisées par Erik Holmberg à l'Observatoire de Lund en 1941, déterminant les forces entre les étoiles lors de la rencontre des galaxies via l'équivalence mathématique entre la propagation de la lumière et l'interaction gravitationnelle : la méthode consistait à mettre des ampoules à la place de chaque étoile et à mesurer les flux lumineux directionnels vers les autres étoiles par une cellule photoélectrique ; la complexité pour l'intégration des équations du mouvement est proportionnelle à N ^[4]. Les premières simulations purement numériques ont été réalisées par Sebastian von Hoerner à l'Institut de calcul astronomique, à Heidelberg en Allemagne. Sverre Aarseth (en) de l'Université de Cambridge (Royaume-Uni) a consacré toute sa vie scientifique au développement d'une série de codes hautement efficaces pour le problème à N corps dans des applications astrophysiques ; ces codes utilisent des pas de temps adaptatifs (hiérarchiques), un schéma de voisinage d'Ahmad-Cohen et la régularisation des collisions entre étoiles (voir infra). La régularisation est une astuce mathématique pour supprimer la singularité dans la loi newtonienne de la gravitation lorsque deux particules deviennent arbitrairement proches^[5]. Les codes de Sverre Aarseth sont utilisés pour étudier la dynamique des amas d'étoiles, des systèmes planétaires et des noyaux galactiques^[6].

De nombreuses simulations sont de tailles suffisamment importantes pour que les effets de la relativité générale ne soient pas négligeables dans l'établissement d'une cosmologie basée sur la métrique Friedmann-Lemaitre-Robertson-Walker. On en tient compte dans la simulation par une mesure des distances évolutive (ou facteur d'échelle) dans un système de coordonnées co-mobiles ; ceci se traduit par un ralentissement des particules dans les coordonnées co-mobiles (ainsi que le décalage vers le rouge de leur énergie physique). Cependant, d'autres contributions de la relativité générale (comme la vitesse finie de la gravité) peuvent être ignorées, car les échelles de temps dynamiques pour la simulation sont longues par rapport au temps de traversée de la lumière, et la courbure spatio-temporelle induite par les particules et leurs vitesses est petite. Les conditions aux limites de ces simulations cosmologiques sont généralement périodiques (ou toroïdales), de sorte qu'un bord du volume de simulation correspond au bord opposé.

Les simulations à N corps sont simples dans leur principe, car elles consistent simplement à intégrer les 6 N équations différentielles ordinaires définissant les mouvements des particules dans la gravité newtonienne. En réalité, le nombre N de particules physiques impliquées est généralement très grand (les simulations typiques en incluent plusieurs millions, la simulation du Millennium (en) en inclut dix milliards^[7]) et le nombre d'interactions particule-particule à calculer est proportionnel à ${\textstyle N^{2}}$ ; l'intégration directe des équations différentielles peut être d'un coût de calcul prohibitif. Par conséquent, un certain nombre de raffinements sont couramment utilisés^[1].

L'intégration numérique est généralement effectuée par des petits pas de temps, à l'aide d'une méthode comme celle du saute-mouton (en anglais « leapfrog »). Cependant toute intégration numérique conduit à des erreurs. Diminuer les pas d'intégration diminue les erreurs, mais augmente le temps de calcul. Le « leapfrog » est une méthode d'intégration du 2e ordre en pas de temps, d'autres méthodes comme celle de Runge – Kutta peuvent avoir une précision du 4e ordre ou plus élevée.

L'un des raffinements les plus simples consiste à poser que chaque particule se meut avec son propre pas de temps, de sorte que les particules avec des temps dynamiques très différents n'ont pas toutes à évoluer au rythme de celles avec le temps le plus court.

Il existe deux schémas d'approximation de base pour réduire le temps de calcul pour de telles simulations. Ils peuvent ramener la complexité de calcul à l'ordre ${\textstyle N\log N}$, ou même mieux, au détriment de la précision :

Dans les méthodes arborescentes^[1], telles qu'une simulation de Barnes-Hut, une arborescence (octree) est généralement utilisée pour diviser le volume en cellules cubiques et seules les interactions entre les particules des cellules voisines doivent être traitées individuellement ; les particules dans les cellules distantes peuvent être traitées collectivement comme une seule grande particule centrée sur le centre de masse de la cellule distante (ou comme une expansion multipolaire d'ordre inférieur). Cela peut réduire considérablement le nombre d'interactions de paires de particules qui doivent être calculées.

Pour éviter que la simulation ne soit submergée par le calcul des interactions particule-particule, les cellules doivent être subdivisées en cellules plus petites dans les parties plus denses de la simulation, celles qui contiennent de nombreuses particules par cellule. Pour les simulations où les particules ne sont pas uniformément réparties, les méthodes de de Callahan et Kosaraju donnent un temps de calcul d'ordre ${\textstyle N\log N}$ par itération, pour une dimension fixe.

Une autre possibilité est la méthode des particules sur mailles (en)^[2], dans laquelle l'espace est discrétisé sur un maillage et pour le calcul du potentiel gravitationnel, les particules sont supposées être réparties entre les sommets les plus proches du maillage. Le calcul de l'énergie potentielle ${\displaystyle \Phi }$ est facile, car l'équation de Poisson :

${\displaystyle \nabla ^{2}\Phi =4\pi G{\rho },\,}$

est triviale à résoudre en utilisant la transformée de Fourier rapide : G est la constante de gravitation et ${\displaystyle {\rho }}$ la densité (nombre de particules aux points du maillage). Par la transformation de Fourier, l'équation de Poisson a la forme simple dans le domaine fréquentiel :

${\displaystyle {\hat {\Phi }}=-4\pi G{\frac {\hat {\rho }}{k^{2}}},\,}$

où ${\displaystyle {\vec {k}}}$ est le nombre d'onde, ${\displaystyle {\hat {\Phi }}\,}$et ${\displaystyle {\hat {\rho }}\,}$désignent les transformées de Fourier du potentiel et de la densité. Comme ${\displaystyle {\vec {g}}=-{\vec {\nabla }}\Phi }$, le champ gravitationnel est obtenu en multipliant ${\displaystyle {\hat {\Phi }}\,}$ par ${\displaystyle -i{\vec {k}}}$ et en calculant la transformée de Fourier inverse. Étant donné que cette méthode est limitée par la taille du maillage, en pratique, un maillage plus petit ou une autre technique (comme la combinaison avec un arbre ou un algorithme simple particule-particule) sont utilisés pour calculer les forces à petite échelle. Parfois, un maillage adaptatif est utilisé, dans lequel les mailles sont beaucoup plus petites dans les régions les plus denses de la simulation (voir figure).

Plusieurs algorithmes différents de perturbation gravitationnelle sont utilisés pour obtenir des estimations assez précises de la trajectoire des objets dans le système solaire.

Il est souvent nécessaire de placer un satellite sur une orbite stationnaire. La trajectoire d'un satellite en orbite proche de la Terre peut être modélisée avec précision à partir de l'orbite elliptique à 2 corps autour du centre de la Terre, et en ajoutant de petites corrections dues à la géométrie de la Terre, à l'attraction gravitationnelle du Soleil et de la Lune, à la traînée atmosphérique, etc. Il est possible de trouver une orbite stationnaire sans calculer la trajectoire réelle du satellite.

La trajectoire d'une petite planète, d'une comète ou d'un vaisseau spatial à longue portée peut souvent être modélisée avec précision à partir de l'orbite elliptique à 2 corps autour du soleil, et en ajoutant de petites corrections dues à l'attraction gravitationnelle des plus grandes planètes sur leurs orbites connues.

Certaines caractéristiques des trajectoires d'un système de particules sur des temps très longs peuvent être calculées directement. La trajectoire réelle d'une particule particulière n'a pas besoin d'être calculée comme étape intermédiaire. Les caractéristiques que l'on peut estimer incluent la stabilité de Liapounov, le durée de Liapounov, diverses grandeurs mesurables de la théorie ergodique, etc.

Bien qu'il y ait des millions ou des milliards de particules réelles dans les simulations classiques, elles sont représentées par des particules virtuelles avec une très grande masse, typiquement 10 ⁹ masses solaires. Cela peut introduire des problèmes d'interactions à courte portée entre les particules telles que la formation de systèmes binaires à deux particules. Comme les particules sont censées représenter un grand nombre de particules de matière noire ou de groupes d'étoiles, ces binaires ne sont pas physiques. Pour éviter cela, une loi de force newtonienne régularisée est utilisée^[5], qui ne diverge pas comme l'inverse du rayon au carré aux courtes distances. La plupart des simulations implémentent cela assez naturellement en exécutant les simulations sur des cellules de taille finie. Il est important de mettre en œuvre la procédure de discrétisation de manière que les particules exercent toujours une force nulle sur elles-mêmes.

La régularisation est une astuce numérique utilisée dans les techniques à N corps pour éviter les divergences numériques lorsqu'une particule s'approche trop près d'une autre (et donc que l'intensité de la force tende vers l'infini). Ceci est obtenu en régularisant le potentiel gravitationnel de chaque particule, sous la forme

${\displaystyle \Phi =-{\frac {1}{\sqrt {r^{2}+\epsilon ^{2}}}},}$

(plutôt que 1/r) où ε est un paramètre ajustable. La valeur de ce paramètre doit être suffisamment petite pour que les simulations restent réalistes.

De nombreuses simulations simulent uniquement la matière noire froide, et n'incluent donc que la force gravitationnelle. L'incorporation de baryons, de leptons et de photons dans les simulations augmente considérablement leur complexité^[8] et des simplifications souvent radicales de la physique sous-jacente doivent être faites. Cependant, il s'agit d'un domaine extrêmement important et de nombreuses simulations modernes tentent maintenant de comprendre les processus qui se produisent lors de la formation des galaxies.

Reif et Tate ^[9] ont étudié la classe de complexité des simulations à N-corps ; par exemple dans un problème de N corps satisfaisant une loi de potentiel électrostatique fixe, déterminer si un corps atteint une zone (sphérique) de destination dans une limite de temps donnée, est un problème qui appartient à la classe de complexité PSPACE.

D'un autre côté, si la question est de savoir si le corps atteint finalement la sphère de destination, le problème est PSPACE-difficile (voir Problème_NP-complet).

Cet article est issu d'une traduction de la page en anglais N-body simulation, sans le code en langage C donné en annexe, merci aux auteurs.

• Problème à N corps
• Simulation informatique
• Millenium simulation
• Structure à grande échelle du cosmos
• Formation et évolution des galaxies
• Simulation Barnes–Hutte

• (en) Mark Vogelsberger, Federico Marinacci, Paul Torrey et Ewald Puchwein, "Cosmological simulations of galaxy formation" Nature Reviews Physics, vol.  2, n^o  1 janvier 2020, p.  42–66 (ISSN 2522-5820, DOI 10.1038/s42254-019-0127-2
• (en) Sverre J. Aarseth, Gravitational N-body Simulations: Tools and Algorithms, Cambridge University Press, 2003 (ISBN 978-0-521-12153-8)
• Bertschinger, « Simulations of structure formation in the universe », Annual Review of Astronomy and Astrophysics, vol. 36, n^o 1,‎ 1998, p. 599–654 (DOI 10.1146/annurev.astro.36.1.599, Bibcode 1998ARA&A..36..599B)
• James Binney et Scott Tremaine, Galactic Dynamics, Princeton University Press, 1987 (ISBN 978-0-691-08445-9)
• Paul B. Callahan et Sambasiva Rao Kosaraju, STOC '92: Proc. ACM Symp. Theory of Computing, ACM, 1992, « A decomposition of multidimensional point sets with applications to k-nearest-neighbors and n-body potential fields (preliminary version) ».

• Portail de la physique