Une approche algorithmique naïve du problème à N corps en interactions gravitationnelles implique une complexité de l'ordre de NxN. Pour 1000 corps, la mise à jour des vecteurs vitesse sera donc obtenue au prix inacceptable d'un million de calculs à chaque frame. Cela est à l'évidence un problème rédhibitoire, qui empêche en pratique de simuler un nombre appréciable de particules.
Cependant, en tolérant une approximation raisonnable du problème, on peut faire diminuer ce nombre de calculs très sensiblement : en 2D, c'est la méthode des Quad-trees (ou Octrees en 3D) encore appelée Simulation de Barnes-Hut.
Si on s'intéresse à un corps et que l'on cherche la somme des forces gravitationnelles que les autres corps exercent sur lui, plutôt que de calculer chaque couple de corps, on peut considérer sans trop faire d'erreur qu'un groupe de point, s'il est suffisamment éloigné, pourra être assimilé à un point unique situé au barycentre de ces corps et de masse égale à la somme des masses de ses constituants.
Toute la difficulté est de identifier ces zones possibles à "résumer". Pour y parvenir, l'espace va être récursivement découpé en sous-parties, organisées dans un arbre (graphe) :
L'arbre est recaclulé à chaque image, ce qui n'est pas si coûteux en performances.
Enfin, lors du calcul de la somme des forces appliquées à un corps, on commence par la racine de l'arbre et on calcule la distance entre ce corps et le centre de gravité du quad. Si cette distance est inférieure à un seuil alors le quad est "résumé", sinon, on explore les quads "enfant" de ce quad.
Ça ne forme pas une galaxie très convaincante mais c'est un début ! Comme la simulation avait tendance à déporter le groupe de particules hors du cadre, j'ai ajouté un attracteur au centre. Ce n'est peut-être pas si arbitraire que cela étant donné que l'on estime que la plupart des grandes galaxies possèdent un trou noir supermassif en leur centre. J'apprécie toutefois une esquisse de bras spiraux dans les premières secondes de la simulation !
À propos des paramètres, le réglage de seuil est celui qui autorise ou non le "résumé" d'un quadtree à son centre de gravité. 0.5 est la valeur préconisée. Plus ce chiffre est élevé, plus des groupes d'étoiles proches seront résumés (moins exact mais plus rapide). Le paramètre d'occupation permet de limiter la profondeur de l'arbre (et accélérer au prix d'erreurs).
Reste à faire ça en WebGL désormais !