n-body problem - numerique-exos

pour faire cette activité sur votre ordi localement, commencez par télécharger le zip

dans ce TP on vous invite à écrire un simulateur de la trajectoire de n corps qui interagissent entre eux au travers de leurs masses, pour produire des sorties de ce genre

on suppose:

on se place dans un monde en 2 dimensions
on fixe au départ le nombre de corps N
chacun est décrit par une masse constante
et au début du monde chacun possède une position et une vitesse

imports¶

on pourra utiliser le mode ipympl de matplotlib

import numpy as np
import matplotlib.pyplot as plt

%matplotlib ipympl

initialisation aléatoire¶

en fixant arbitrairement des limites dans l’espace des positions, des vitesses et des masses, la fonction init_problem() tire au hasard une configuration de départ pour la simulation

# les bornes pour le tirage au sort initial
mass_max = 3.
    
x_min, x_max = -10., 10.
y_min, y_max = -10., 10.

speed_max = 1.

# votre code

def init_problem(N):
    """
    retourne un tuple masses, positions, speeds
    de formes resp.   (N,)    (2, N)     (2, N)
    tiré au sort dans les bornes définies ci-dessus
    """
    return None, None, None

# pour tester

# normalement vous devez pouvoir faire ceci

masses, positions, speeds = init_problem(10)

# et ceci devrait afficher OK
try:
    masses.shape == (10,) and positions.shape == speeds.shape == (2, 10)
    print("OK")
except:
    print("KO")

KO

initialisation reproductible¶

par commodité on vous donne la fonction suivante qui crée 3 objets:

le premier - pensez au soleil - de masse 3, an centre de la figure, de vitesse nulle
et deux objets de masse 1, disposés symétriquement autour du soleil
- position initiale (5, 1) et vitesse initiale (-1, 0)
- symétrique en (-5, -1) et vitesse initiale (1, 0)

# for your convenience

def init3():
    # first element is sun-like: heavy, at the center, and no speed
    masses = np.array([3, 1, 1], dtype=float)
    positions = np.array([
        [0, 5, -5], 
        [0, 1, -1]], dtype=float)
    speeds = np.array([
        [0, -1, 1], 
        [0, 0, 0]], dtype=float)
    return masses, positions, speeds

les forces¶

à présent, on va écrire une fonction qui va calculer les influences de toutes les particules entre elles, suivant la loi de Newton

\vec{F}_i = \sum_{\substack{j=1 \\ j \neq i}}^N G \, m_i m_j \, \frac{\vec{r}_j - \vec{r}_i}{\lvert \vec{r}_j - \vec{r}_i \rvert^3}

(1)

pour cela on se propose d’écrire la fonction suivante

# votre code

def forces(masses, positions, G=1.0):
    """
    returns an array of shape (2,  N)
    that contains the force felt by each mass from all the others
    G is the Newton constant and by default is set to 1 
    (we have abstract units anyway)
    """
    pass

# pour tester, voici les valeurs attendues avec la config prédéfinie

masses, positions, speeds = init3()

f = forces(masses, positions)

## uncomment, should be true
# np.all(np.isclose(f, np.array([
#     [ 0.        , -0.12257258,  0.12257258],
#     [ 0.        , -0.02451452,  0.02451452]])))

le simulateur¶

à présent il nous reste à utiliser cette brique de base pour “faire avancer” le modèle depuis son état initial et sur un nombre fixe d’itérations

cela pourrait se passer dans une fonction qui ressemblerait à ceci

# votre code

def simulate(masses, positions, speeds, dt=0.1, nb_steps=100):
    """
    should return the positions across time
    so an array of shape (nb_steps, 2, N)
    optional dt is the time step
    """
    pass

# pour tester

SMALL_STEPS = 4

s = simulate(masses, positions, speeds, nb_steps=SMALL_STEPS)

try:
    if s.shape == (SMALL_STEPS, 2, 3):
        print("shape OK")
except Exception as exc:
    print(f"OOPS {type(exc)} {exc}")

OOPS <class 'AttributeError'> 'NoneType' object has no attribute 'shape'

# pour tester
# uncomment - should be true

# # isolate first step
# positions1 = s[1]

# np.all(np.isclose(positions1, np.array([
#     [ 0.        ,  4.89877427, -4.89877427],
#     [ 0.        ,  0.99975485, -0.99975485]
# ])))

dessiner¶

ne reste plus qu’à dessiner; quelques indices potentiels:

1. chaque corps a une couleur; l’appelant peut vous passer un jeu de couleurs, sinon en tirer un au hasard
2.a pour l’épaisseur de chaque point, on peut imaginer utiliser la masse de l’objet
2.b ou peut-être aussi, à tester, la vitesse de l’objet (plus c’est lent et plus on l’affiche en gros ?)

# votre code

def draw(simulation, masses, colors=None, scale=10.):
    """
    takes as input the result of simulate() above,
    and draws the nb_steps positions of each of the N bodies
    ideally it should return a matplotlib Axes object

    one can provide a collection of N colors to use for each body
    if not provided this is randomized

    also the optional scale parameter is used as a constant
    multiplier to obtain the final size of each dot on the figure
    """
    pass

un jeu de couleurs¶

# for convenience

colors3 = np.array([
    [32, 32, 32],
    (228, 90, 146),
    (111, 0, 255),
]) / 255

on assemble le tout¶

pour commencer et tester, on se met dans l’état initial reproductible

# décommentez ceci pour tester votre code

# masses, positions, speeds = init3()
# draw(simulate(masses, positions, speeds), masses, colors3)

et avec ces données vous devriez obtenir plus ou moins une sortie de ce genre
mais voyez aussi la discussion ci-dessous sur les diverses stratégies possibles

après vous avez le droit de vous enhardir avec des scénarii plus compliqués par exemple avec ce code

m5, p5, s5 = init_problem(5)
sim5 = simulate(m5, p5, s5, nb_steps=1000)
draw(sim5, m5, scale=3);
plt.savefig("random5.png")

j’ai pu obtenir ceci

partie optionnelle¶

option 1: la 3D¶

modifiez votre code pour passer à une simulation en 3D

et appliquez-le aux conditions initiales suivantes

# des conditions initiales pour la 3D
def init3_3d():
    # masses: 3 bodies with the same mass
    masses = np.array([1, 1, 1], dtype=float)
    positions = np.array([
        [1, 0, 0],        
        [0, 1, 0],        
        [0, 0, 1]], dtype=float)
    speeds = np.array([
        [0, 0, 1],        
        [1, 0, 0],        
        [0, 1, 0]], dtype=float)
    return masses, positions, speeds

option 2: un rendu plus interactif¶

le rendu sous forme de multiples scatter plots donne une idée du résultat mais c’est très améliorable
voyez un peu si vous arrivez à produire un outil un peu plus convivial pour explorer les résultats de manière interactive; avec genre

une animation qui affiche les points au fur et à mesure du temps
qu’on peut controler un peu comme une vidéo avec pause / backward / forward
l’option de laisser la trace du passé
et si vous avez un code 3d, la possibilité de changer le point de vue de la caméra sur le monde
etc etc...

voici une possibilité avec matplotlib; mais cela dit ne vous sentez pas obligé de rester dans Jupyter Lab ou matplotlib, il y a plein de technos rigolotes qui savent se décliner sur le web, vous avez l’embarras du choix...

# prune-remove-input

# credit: Damien Corral
# with good old matplotlib FuncAnimation

from matplotlib.animation import FuncAnimation
from IPython.display import HTML

def animate(simulation, masses, colors=None, scale=5., interval=50):
    nb_steps, _, N = simulation.shape
    colors = (colors if colors is not None
              else np.random.uniform(0.3, 1., size=(N, 3)))

    fig, ax = plt.subplots()
    ax.set_title(f"we have {N} bodies over {nb_steps} steps")

    ax.set_xlim(simulation[:, 0].min() - 1, simulation[:, 0].max() + 1)
    ax.set_ylim(simulation[:, 1].min() - 1, simulation[:, 1].max() + 1)

    scat = ax.scatter(np.zeros(N), np.zeros(N), c=colors, s=(masses*scale)**2)

    def init():
        scat.set_offsets(np.zeros((nb_steps, N)))
        return scat

    def update(step):
        x, y = simulation[step]
        scat.set_offsets(np.c_[x, y])
        return scat

    animation = FuncAnimation(
        fig, update, frames=nb_steps,
        init_func=init, blit=True, interval=interval
    )
    plt.close()
    return animation

def animate_from_file(filename):
    simulation = np.loadtxt(filename).reshape((100, 2, 3))
    animation = animate(simulation, masses, colors=colors3)
    return HTML(animation.to_jshtml())

animate_from_file("data/init3-simu-1.txt")

plusieurs stratégies¶

Pour les curieux, vous avez sans doute observé qu’il y a plusieurs façons possibles d’écrire la fonction simulate()

dans ce qui suit, on note l’accélération $a$ , la vitesse $s$ et la position $p$

Bref, vous voyez qu’il y a énormément de liberté sur la façon de s’y prendre
Ce qui peut expliquer pourquoi vous n’obtenez pas la même chose que les illustrations avec pourtant les mêmes données initiales

D’autant que, c’est bien connu, ce problème des n-corps est l’exemple le plus célèbre de problème instable, et donc la moindre divergence entre deux méthodes de calcul peut entrainer de très sérieux écarts sur les résultats obtenus

Voici d’ailleurs les résultats obtenus avec ces deux approches alternatives, et vous pouvez constater qu’effectivement les résultats sont tous très différents !

imports¶

initialisation aléatoire¶

initialisation reproductible¶

les forces¶

le simulateur¶

dessiner¶

un jeu de couleurs¶

on assemble le tout¶

partie optionnelle¶

option 1: la 3D¶

option 2: un rendu plus interactif¶

plusieurs stratégies¶

approche 1bis¶

approche 2¶