n-body: elaborate on several possible strategies

Author: parmentelat

notebooks/numpy-tps/n-body/.teacher/README-n-body-corrige-nb.md

@@ -1,7 +1,4 @@
 ---
-authors:
-- name: Aubin Geoffre
-- name: Thierry Parmentelat
 jupytext:
   formats: md:myst
   text_representation:
@@ -25,7 +22,7 @@ pour faire cette activité sur votre ordi localement, {download}`commencez par t
 
 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
 
-```{image} media/init3.png
+```{image} media/init3-1.png
 :align: center
 :width: 600px
 ```
@@ -372,13 +369,14 @@ masses, positions, speeds = init3()
 simulation3 = simulate(masses, positions, speeds, nb_steps=steps)
 draw(simulation3, masses, colors3);
 
-plt.savefig("init3.png")
-np.savetxt("simulation3.txt", simulation3.reshape(-1))
+plt.savefig("init3-1.png")
+np.savetxt("init3-simu-1.txt", simulation3.reshape(-1))
 print(simulation3.shape)
 ```
 
-et avec ces données vous devriez obtenir plus ou moins une sortie de ce genre
-```{image} media/init3.png
+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](label-n-body-strategies)
+```{image} media/init3-1.png
 ```
 
 +++
@@ -466,13 +464,174 @@ def animate(simulation, masses, colors=None, scale=5., interval=50):
         scat.set_offsets(np.c_[x, y])
         return scat
 
-    ani = FuncAnimation(fig, update, frames=nb_steps,
-                        init_func=init, blit=True, interval=interval)
+    animation = FuncAnimation(
+        fig, update, frames=nb_steps,
+        init_func=init, blit=True, interval=interval
+    )
     plt.close()
-    return ani
+    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())
 
-simulation3 = np.loadtxt("data/simulation3.txt").reshape((100, 2, 3))
-animation = animate(simulation3, masses, colors=colors3)
-HTML(animation.to_jshtml())
+animate_from_file("data/init3-simu-1.txt")
 ```
+
+(label-n-body-strategies)=
+## 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$
+
+````{admonition} Approche 1
+:class: note
+dans l'implémentation qui a servi à calculer l'illustration ci-dessus, on a écrit principalement ceci:
+- on calcule l'accélération
+- ce qui permet d'extrapoler les vitesses  
+  $s = s + a.dt$
+- et ensuite d'extrapoler les positions  
+  $p = p + s.dt$
+```{admonition} 1bis
+:class: tip
+une variante consiste à intervertir les deux, en arguant du fait que la vitesse à l'instant $t$ agit sur la position à l'instant $t$
+- on extrapole d'abord les positions
+- puis seulement on calcule l'accélération
+- enfin on extrapole les vitesses
+```
+````
+````{admonition} Approche 2
+:class: tip
+dans cette approche plus fine, on utiliserait deux versions de l'accélération (l'instant présent et l'instant suivant), et un dévelopmment du second ordre, ce qui conduirait à
+- calculer la position comme $p = p + s.dt + \frac{a}{2}.dt^2$
+- calculer les accélérations $a_+$ sur la base de cette nouvelle position
+- estimer l'accélération sur l'intervalle comme la demie-somme entre les deux accélérations
+  $a_m = (a+a_+)/2$
+- mettre à jour les vitesses
+  $s = s + a_m.dt$
+- ranger $a_+$ dans $a$ pour le prochain instant
+````
+
++++
+
+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 !
+
++++
+
+### approche 1bis
+
+```{code-cell} ipython3
+# prune-cell
+
+def simulate_1bis(masses, positions, speeds, dt=0.1, nb_steps=100):
+    # initialize result
+    N = masses.shape[0]
+    result = np.empty((nb_steps, 2, N))
+    # store initial position
+    result[0] = positions
+
+    # because we have stored the initial state
+    # we need to start at 1 here, we dont have space for more
+    for step_counter in range(1, nb_steps):
+        # apply on positions
+        positions += dt*speeds
+
+        # compute forces and accelerations
+        accelerations = forces(masses, positions) / masses
+
+        # apply on speeds
+        speeds += dt*accelerations
+
+        # store current position
+        result[step_counter] = positions
+    return result
+```
+
+```{code-cell} ipython3
+# prune-cell
+
+steps = 100
+
+masses, positions, speeds = init3()
+simulation3 = simulate_1bis(masses, positions, speeds, nb_steps=steps)
+draw(simulation3, masses, colors3);
+
+plt.savefig("init3-1bis.png")
+np.savetxt("init3-simu-1bis.txt", simulation3.reshape(-1))
+print(simulation3.shape)
+```
+
+````{grid} 2 2 2 2 
+```{image} media/init3-1bis.png
+```
+```{code-cell} python
+:tags: [remove-input]
+animate_from_file("data/init3-simu-1bis.txt")
+```
+````
+
++++
+
+### approche 2
+
+```{code-cell} ipython3
+# prune-cell
+
+def simulate_2(masses, positions, speeds, dt=0.1, nb_steps=100):
+    # initialize result
+    N = masses.shape[0]
+    result = np.empty((nb_steps, 2, N))
+    # store initial position
+    result[0] = positions
+
+    # initialize accelerations
+    accelerations = forces(masses, positions) / masses
+    
+    # because we have stored the initial state
+    # we need to start at 1 here, we dont have space for more
+    for step_counter in range(1, nb_steps):
+        # compute positions
+        positions = positions + speeds * dt + accelerations * dt**2 * 0.5
+        # compute forces and accelerations for next instant
+        next_accelerations = forces(masses, positions) / masses
+        # use half-mean as an approximation of acceleration over the period
+        speeds += 0.5 * (accelerations + next_accelerations)/2 * dt
+        # save for next iteration
+        accelerations = next_accelerations
+
+        # store current position
+        result[step_counter] = positions
+    return result
+```
+
+```{code-cell} ipython3
+# prune-cell
+
+steps = 100
+
+masses, positions, speeds = init3()
+simulation3 = simulate_2(masses, positions, speeds, nb_steps=steps)
+draw(simulation3, masses, colors3);
+
+plt.savefig("init3-2.png")
+np.savetxt("init3-simu-2.txt", simulation3.reshape(-1))
+print(simulation3.shape)
+```
+
+````{grid} 2 2 2 2 
+```{image} media/init3-2.png
+```
+```{code-cell} python
+:tags: [remove-input]
+animate_from_file("data/init3-simu-2.txt")
+```
+````

notebooks/numpy-tps/n-body/ARTEFACTS

@@ -0,0 +1,3 @@
+README*
+data/
+media/

notebooks/numpy-tps/n-body/ARTEFACTS-n-body.zip

@@ -1,7 +1,4 @@
 ---
-authors:
-- name: Aubin Geoffre
-- name: Thierry Parmentelat
 jupytext:
   formats: md:myst
   text_representation:
@@ -25,7 +22,7 @@ pour faire cette activité sur votre ordi localement, {download}`commencez par t
 
 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
 
-```{image} media/init3.png
+```{image} media/init3-1.png
 :align: center
 :width: 600px
 ```
@@ -258,8 +255,9 @@ pour commencer et tester, on se met dans l'état initial reproductible
 # draw(simulate(masses, positions, speeds), masses, colors3)
 ```
 
-et avec ces données vous devriez obtenir plus ou moins une sortie de ce genre
-```{image} media/init3.png
+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](label-n-body-strategies)
+```{image} media/init3-1.png
 ```
 
 +++
@@ -347,13 +345,93 @@ def animate(simulation, masses, colors=None, scale=5., interval=50):
         scat.set_offsets(np.c_[x, y])
         return scat
 
-    ani = FuncAnimation(fig, update, frames=nb_steps,
-                        init_func=init, blit=True, interval=interval)
+    animation = FuncAnimation(
+        fig, update, frames=nb_steps,
+        init_func=init, blit=True, interval=interval
+    )
     plt.close()
-    return ani
+    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())
 
-simulation3 = np.loadtxt("data/simulation3.txt").reshape((100, 2, 3))
-animation = animate(simulation3, masses, colors=colors3)
-HTML(animation.to_jshtml())
+animate_from_file("data/init3-simu-1.txt")
 ```
+
+(label-n-body-strategies)=
+## 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$
+
+````{admonition} Approche 1
+:class: note
+dans l'implémentation qui a servi à calculer l'illustration ci-dessus, on a écrit principalement ceci:
+- on calcule l'accélération
+- ce qui permet d'extrapoler les vitesses  
+  $s = s + a.dt$
+- et ensuite d'extrapoler les positions  
+  $p = p + s.dt$
+```{admonition} 1bis
+:class: tip
+une variante consiste à intervertir les deux, en arguant du fait que la vitesse à l'instant $t$ agit sur la position à l'instant $t$
+- on extrapole d'abord les positions
+- puis seulement on calcule l'accélération
+- enfin on extrapole les vitesses
+```
+````
+````{admonition} Approche 2
+:class: tip
+dans cette approche plus fine, on utiliserait deux versions de l'accélération (l'instant présent et l'instant suivant), et un dévelopmment du second ordre, ce qui conduirait à
+- calculer la position comme $p = p + s.dt + \frac{a}{2}.dt^2$
+- calculer les accélérations $a_+$ sur la base de cette nouvelle position
+- estimer l'accélération sur l'intervalle comme la demie-somme entre les deux accélérations
+  $a_m = (a+a_+)/2$
+- mettre à jour les vitesses
+  $s = s + a_m.dt$
+- ranger $a_+$ dans $a$ pour le prochain instant
+````
+
++++
+
+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 !
+
++++
+
+### approche 1bis
+
++++
+
+````{grid} 2 2 2 2 
+```{image} media/init3-1bis.png
+```
+```{code-cell} python
+:tags: [remove-input]
+animate_from_file("data/init3-simu-1bis.txt")
+```
+````
+
++++
+
+### approche 2
+
++++
+
+````{grid} 2 2 2 2 
+```{image} media/init3-2.png
+```
+```{code-cell} python
+:tags: [remove-input]
+animate_from_file("data/init3-simu-2.txt")
+```
+````