DS 3 : Centrale MP 2015

3ème TD/DS, algorithmique sur des listes, puis schémas numériques (Euler, Verlet), appliqué à un problème de mécanique à N corps, et enfin un dernier exercice avec une base de données (et des requêtes SQL).

Correction du TD/DS du mercredi 27 et jeudi 28 janvier (2016).

Voici ci-dessous la documentation de la correction complète du sujet d’écrit d’informatique (Centrale, 2015), pour MP, PC, PSI et TSI.

Ce sujet d’écrit a été posé en devoir écrit non noté (TD/DS) pour le cours d’informatique pour tous en prépa MP au Lycée Lakanal.

Veuillez consulter le code source pour tous les détails svp. Ci dessous se trouve simplement la documentation de ce programme.

• Date : Samedi 16 janvier 2016,
• Auteur : Lilian Besson, pour le cours d’informatique pour tous en prépa MP (http://perso.crans.org/besson/infoMP/),
• Licence : MIT Licence (http://lbesson.mit-license.org).

---

Correction :

Centrale_MP__2015.smul(alpha, liste)

Prend un nombre alpha et une liste de nombres, multiple chaque élément de la liste par le nombre et renvoie une nouvelle liste (il ne faut pas modifier la liste argument !).

• Exemple :

>>> print(smul(2, [1, 2, 3]))
[2, 4, 6]

• Version courte :

  return [alpha * x for x in liste]
  

• Version plus longue (smul), avec une boucle et des append, plus proche de la syntaxe au programme.

Centrale_MP__2015.vsom(liste1, liste2)

Prend deux listes de nombres de même longueur liste1 et liste2 et qui renvoie une nouvelle liste (il ne faut pas modifier aucune des listes données en argument !), constituée de la somme terme à terme de ces deux listes.

• Exemple :

>>> print(vsom([1, 2, 3], [4, 5, 6]))
[5, 7, 9]

• Version courte :

  return [x1 + x2 for (x1, x2) in zip(liste1, liste2)]
  

• Version plus longue (vsom), avec une boucle et des append, plus proche de la syntaxe au programme.

Centrale_MP__2015.vdif(liste1, liste2)

Prend deux listes de nombres de même longueur liste1 et liste2 et qui renvoie une nouvelle liste (il ne faut pas modifier aucune des listes données en argument !), constituée de la différence terme à terme de ces deux listes.

• Exemple :

>>> print(vdif([1, 2, 3], [4, 5, 6]))
[-3, -3, -3]

• Version courte :

  return [x1 - x2 for (x1, x2) in zip(liste1, liste2)]
  

• Version plus longue (vdif), avec une boucle et des append, plus proche de la syntaxe au programme.

Centrale_MP__2015.question_ii_a_1()

Partie II : Étude de schémas numériques.

Réponse mathématique à la question II.A.1) :

Si \(z(t) = y'(t)\), alors l’équation (II.1) s’écrit \(z'(t) = f(y(t))\), qu’on peut mettre sous forme vectorielle (en \([y, z]\)) :

(1)\[\begin{split}\begin{cases} y'(t) &= \; z(t) \\ z'(t) &= \; f(y(t)) \end{cases}\end{split}\]

Astuce

Cette manipulation est très classique, on transforme une équation différentielle scalaire (en \(y\)) d’ordre 2 (II.1) en une équation différentielle vectorielle (en \([y, z]\)) d’ordre 1 (S) ((1)).

Il faut savoir faire ce genre de manipulation sans réfléchir, c’est très classique et ça tombe très souvent.

Voir aussi

La solution du TP5 sur la méthode d’Euler.

Centrale_MP__2015.question_ii_a_2()

Réponse mathématique à la question II.A.2) :

Pour les deux calculs avec les intégrales, il suffit de calculer l’intégrale en utilisant une primitive, les deux étant données par le système différentiel (S) ((1)).

Par définition de \(z(t)\), on remarque déjà que

\[\begin{split}\int_{t_i}^{t_{i+1}} z(t) \mathrm{d}t &= \int_{t_i}^{t_{i+1}} y'(t) \mathrm{d}t \\ &= \bigl[ y(t) \bigr]_{t_i}^{t_{i+1}} = y(t_{i+1}) - y(t_i).\end{split}\]

Donc comme demandé pour l’équation (II.3), on a déjà l’égalité de gauche :

\[y(t_{i+1}) = y(t_i) + \int_{t_i}^{t_{i+1}} z(t) \mathrm{d}t.\]

Ensuite, pour le terme de gauche, on procède exactement pareil, puisque d’après (S) ((1)), on a \(z' = f(y)\), donc :

\[\begin{split}\int_{t_i}^{t_{i+1}} f(y(t)) \mathrm{d}t &= \int_{t_i}^{t_{i+1}} z'(t) \mathrm{d}t \\ &= \bigl[ z(t) \bigr]_{t_i}^{t_{i+1}} = z(t_{i+1}) - z(t_i).\end{split}\]

L’énoncé nous demande ensuite d’approcher chaque terme sous le signe intégrale est remplacé par sa valeur prise en la borne inférieure (et on utilise le fait que les pas de temps sont constants : \(\forall i, \; t_{i+1} - t_i = h\)) :

(2)\[\begin{split}\begin{cases} \int_{t_i}^{t_{i+1}} f(y(t)) \mathrm{d}t &\simeq (t_{i+1} - t_i) \times f(y(t_i)) \simeq h \times f(y_i) = h \times f_i, \\ \int_{t_i}^{t_{i+1}} z(t) \mathrm{d}t &\simeq (t_{i+1} - t_i) \times z(t_i) \simeq h \times z_i. \end{cases}\end{split}\]

C’est l’approximation fait par la méthode des rectangles gauches de Riemann. Voir cette page (en anglais) pour plus de détails, ou cette page Wikipédia.

Notez qu’en ((2)) on a aussi rajouté (à la fin de la ligne) que les valeurs discrétisées \(y_i\) et \(z_i\) vérifieront \(y_i \simeq y(t_i)\) et \(z_i \simeq z(t_i) = y'(t_i)\).

Voir aussi

Pour retrouver le schéma d’Euler gauche (dit explicite), on utilise l’approximation de Riemann gauche. Et inversement, si chaque terme sous le signe intégrale est remplacé par sa valeur prise en la borne supérieure, on retrouve le schéma d’Euler implicite (relisez votre cours de l’an dernier, ou ce paragraphe sur Wikipédia pour un rappel).

Centrale_MP__2015.question_ii_b_1()

Réponse mathématique à la question II.B.1) :

On utilise les deux équations (II.3) et les deux approximations données plus haut (2) :

(3)\[\begin{split}\begin{cases} y_{i+1} \; \leftarrow& y_i + h \times z_i \\ z_{i+1} \; \leftarrow& z_i + h \times f(y_i) \end{cases}\end{split}\]

Ces deux équations de récurrence, sont couplées mais vectorielles : on peut écrire \(X = [y, z]\), et on aura \(X_{i+1} \leftarrow X_i + h \times F(X_i)\) avec \(F(X) = [X_1, f(X_2)]\). Dès qu’on connaît \(h\) le pas de temps, \(f\) l’opérateur de l’équation différentielle initiale, et les deux valeurs initiales \(y_0, z_0\), les équations (3) donnent bien deux suites \((y_i)_{i \in J_n}\) et \((z_i)_{i \in J_n}\), définies de façon unique.

Centrale_MP__2015.euler(n, tmin, tmax, y0, z0, f)

Utilisation du schéma d’Euler (3), pour calculer les deux suites \((y_i)_{i \in J_n}\) et \((z_i)_{i \in J_n}\) en fonction de leurs valeurs initiales \(y_0, z_0\), du nombre de points \(n\), et de l’intervalle de temps \([t_{\min}, t_{\max}]\), où \(t_{\min} =\) tmin et \(t_{\max} =\) tmax, pour l’équation différentielle \(y''(t) = f(y(t))\) (II.1).

• Renvoie deux listes y et z qui contiennent \(n + 1\) points, en utilisant le pas de temps \(h = \frac{t_{\max} - t_{\min}}{n - 1} > 0\).

• Complexité en temps : \(O(n)\),

• Complexité en mémoire : \(O(n)\) aussi.

• Symétriquement, on aurait pu prendre le pas de temps h, et calculer le nombre de point n grâce à la relation :

  n = int((tmax - tmin) / h)
  

Note

À part le choix des paramètres, cette fonction est très simple à écrire, il s’agit juste d’écrire en Python les équations mathématiques (3). Pour un exemple très similaire utilisant un autre schéma numérique, voir verlet().

Astuce

Devoir choisir soi-même les paramètres peut être déstabilisant, voici quelques conseils :

• Don’t panic!
• Astuce : écrire le corps de la fonction, et ensuite mettre en argument tout ce dont elle a besoin ! Ici, on a au moins besoin de y0, z0, et ensuite soit n et h, soit n, tmin et tmax. Et f.
• Si vous n’êtes pas sûr si une variable doit être un argument ou supposée définie avant (ie. plus haut dans le programme Python), faites un choix, et spécifiez le ! Par exemple, cette fonction euler() a besoin d’une fonction f, soit définie plus haut, soit donnée en argument. En partie III, les deux fonctions force2() et forceN() auront besoin d’une constante G par exemple, qu’on peut soit supposée définie avant, soit recalculée/redéfinie dans le corps de la fonction, il faut juste être clair.
• Et rappelez-vous que l’ordre des arguments n’a pas vraiment d’importance.

• Exemple : voir aussi plus bas, où on continue d’étudier l’exemple de \(f(y) = -\omega^2 y\).

>>> y0 = 3; z0 = 0
>>> tmin = 0.0; tmax = 3.0
>>> n = 100
>>> h = (tmax - tmin) / float(n - 1)
>>> print("n =", n, "et h =", h)
n = 100 et h = 0.030303030303030304
>>> from math import *
>>> omega = 2 * pi
>>> f = lambda y: - omega**2 * y
>>> y, z = euler(n, tmin, tmax, y0, z0, f)
>>> # Verification, conditions initiales bien respectees
>>> print(y0 == y[0] and z0 == z[0])
True
>>> print(y)  # 100 points  
[3, 3.0, 2.891..., ..., 17.787..., 17.905...]
>>> print(z)  # 100 points  
[0, -3.588..., -7.177..., ..., 3.891..., -17.388...]
>>> # Dernier point, comme vu sur la Figure 1
>>> print(y[-1], z[-1])  
17.905... -17.388...

Voir aussi

Si besoin, lisez cette page Wikipédia sur la méthode d’Euler (au moins le début), ça ne peut pas faire de mal...

Centrale_MP__2015.omega0 = 6.283185307179586

Exemple de valeur de la pulsation propre du système, \(\omega = 2 \pi\). pi peut être importé depuis le module math ou numpy (ou défini manuellement pi = 3.1415.... = 3.141592653589793 à la précision que vous voulez).

Centrale_MP__2015.f(y, omega=6.283185307179586)

Exemple de fonction \(f\) définissant l’équation différentielle étudiée : \(f(y) = -\omega^2 y\).

Centrale_MP__2015.question_ii_b_3()

On s’intéresse maintenant à l’exemple de \(f(y) = -\omega^2 y\).

• On reconnaît une équation différentielle classique, souvent rencontrée en physique, celle du pendule pesant idéal, sans frottement ni excitation, linéarisée pour de petites oscillations \(y = \theta \simeq 0\).
• On sait la forme qu’auront les solutions et les portraits de phases (des sinusoïdes pour \(y(t)\), resp. des ellipses pour \(y(y'))\)). Mais notez que nulle part le sujet ne demande de résoudre la solution (ie. d’expliciter la forme qu’aura \(y(t)\)), et qu’on n’en a pas besoin !

Centrale_MP__2015.question_ii_b_3_a()

Réponse mathématique à la question II.B.3.a) :

• On s’inspire du cours de physique pour définir cette énergie \(E(t)\), qu’on suppose tout d’abord dépendante du temps :

(4)\[E(t) := \frac{1}{2} \left(y'(t)\right)^2 + \frac{1}{2} \omega^2 \left(y(t)\right)^2\]

• Ensuite, le but de la question est simplement de montrer que \(E(t) = E = E(t_0)\) est constante, ce qui se fait en deux lignes de calcul en dérivant \(E(t)\) dans l’équation ci-dessus ((4)).
• Notez qu’il suffit d’écrire \(f(y) = - \omega^2 y\) et \(g = - f'\) (par définition), puis d’intégrer pour obtenir \(g(y) = - \frac{1}{2} \omega^2 y^2 + C\) (et on se fiche de la constante).
• Cette énergie \(E\) a exactement la forme connue par le cours de physique : énergie cinétique + énergie potentielle.

Astuce

Tout le sujet manipule des quantités physiques, donc vérifier les unités permet de vérifier qu’on écrit pas n’importe quoi (faire à l’œil une analyse dimensionnelle rapide peut aider !).

Centrale_MP__2015.question_ii_b_3_b()

Réponse mathématique à la question II.B.3.b) :

• Il suffit d’écrire que \(E(t_i) = \frac{1}{2} \left(y'(t_i)\right)^2 + \frac{1}{2} \omega^2 \left(y(t_i)\right)^2\), donc \(E(t_i) \simeq \frac{1}{2} \left(z_i\right)^2 + \frac{1}{2} \omega^2 \left(y_i\right)^2 =: E_i\) (c’est ainsi qu’on définit \(E_i\)).
• Et donc, on peut écrire \(E_{i+1}\) en utilisant les équations de récurrence (3), et le développer pour obtenir l’expression demandée, en utilisant le fait que \(f(y) = -\omega^2 y\) :

\[\begin{split}2 E_{i+1} &= \left(z_{i+1}\right)^2 + \omega^2 \left(y_{i+1}\right)^2 \\ &= \left(z_i - \omega^2 h y_i\right)^2 + \omega^2 \left(y_i + h z_i\right)^2 \\ &= \bigl(\left(z_i\right)^2 - 2 \omega^2 h z_i y_i + \left(h \omega^2 y_i\right)^2\bigr) + \omega^2 \bigl(\left(y_i\right)^2 + 2 h y_i z_i + \left(h z_i\right)^2\bigr) \\ &= 2 E_i + h^2 \omega^2 \bigl( 2 E_i \bigr)\end{split}\]

• Un terme se simplifie, et on trouve donc le résultat demandé : \(E_{i+1} - E_i = h^2 \omega^2 E_i > 0\) (l’énergie est strictement croissante au cours des itérations).

Note

• Notez que cette erreur d’approximation dans les énergies est en \(O(h^2)\), ce qui est déjà très rassurant : plus le pas de temps \(h\) est petit (ie. plus on calcule de point \(n\)), moins on fera d’erreur d’approximations.
• Le schéma de Verlet, étudié plus bas, sera en \(O(h^3)\), et donc sera plus précis (voir question_ii_c_2_a() plus bas).

Centrale_MP__2015.question_ii_b_3_c()

Réponse mathématique à la question II.B.3.c) :

• Un schéma numérique qui satisfait à la conservation de l’énergie \(E\) satisferait la conservation des énergies approchées \(E_i\), et donc :

\[\forall i \in J_n, \;\;\;\; E_{i+1} - E_i = 0.\]

Attention

C’est bien ce terme \(E_{i+1} - E_i\) qui vaut zéro, pas le terme de droite calculé avant en II.B.3.b). Le calcul précédent (question_ii_b_3_b()) n’est vrai que pour le schéma d’Euler, qui n’est précisément pas un schéma conservant les énergies. Il ne faut surtout pas écrire “ça implique que \(h^2 \omega^2 \bigl( 2 E_i \bigr) = 0\)”, ça n’a pas de sens.

Centrale_MP__2015.question_ii_b_3_d()

Réponse mathématique à la question II.B.3.d) :

• Le portrait de phase d’un schéma numérique qui satisferait à la conservation de l’énergie \(E\) serait une ellipse parfaite.

Attention

Attention, aucune raison que le portrait de phase ne soit un cercle, à moins de normaliser les deux quantités \(y\) et \(z = y'\), elles n’ont pas la même unité et n’ont aucune raison d’avoir le même ordre de grandeur.

Centrale_MP__2015.question_ii_b_3_e()

Réponse mathématique à la question II.B.3.e) :

• Le graphe de portrait de phase obtenu avec le schéma d’Euler confirme le calcul effectué plus haut, on avait montré que \((E_i)_{i \in J_n}\) était une suite strictement croissante, et donc que la norme du point \((\omega y_i, z_i)\), \(R_i = z_i^2 + \omega^2 y_i^2\) augmente strictement au cours des itérations (on multiplie la coordonnée \(y\) pour être homogène).
• Et... c’est exactement ce qu’on observe dans ce portrait de phase (à gauche), où il faut comprendre que la spirale diverge au cours des itérations (le point initial est \((y_0, z_0) = (3, 0)\) au centre) :

Centrale_MP__2015.verlet(n, tmin, tmax, y0, z0, f)

Utilisation du schéma de Verlet (5), pour calculer les deux suites \((y_i)_{i \in J_n}\) et \((z_i)_{i \in J_n}\) en fonction de leurs valeurs initiales \(y_0, z_0\), du nombre de points \(n\), et de l’intervalle de temps \([t_{\min}, t_{\max}]\), où \(t_{\min} =\) tmin et \(t_{\max} =\) tmax, pour l’équation différentielle \(y''(t) = f(y(t))\) (II.1).

(5)\[ \begin{align}\begin{aligned}\forall i \in J_n, \;\; f_i := f(y_i).\\\begin{split}\begin{cases} y_{i+1} \; \leftarrow& y_i + h \times z_i + \frac{h^2}{2} \times f_i \\ z_{i+1} \; \leftarrow& z_i + h \times \bigl( \frac{f_i + f_{i+1}}{2} \bigr) \end{cases}\end{split}\end{aligned}\end{align} \]

• À noter que cette fois, la mise à jour de la première coordonnée (ie. le calcul de \(y_{i+1}\)) doit impérativement se faire avant la mise à jour de la seconde coordonnée : \(z_{i+1}\) utilise en effet \(f_{i+1} = f(y_{i+1})\).

• Renvoie deux listes y et z qui contiennent \(n + 1\) points, en utilisant le pas de temps \(h = \frac{t_{\max} - t_{\min}}{n - 1} > 0\).

• Complexité en temps : \(O(n)\) comme pour le schéma d’Euler,

• Complexité en mémoire : \(O(n)\) aussi.

• Symétriquement, on aurait pu prendre le pas de temps h, et calculer le nombre de point n grâce à la relation :

  n = int((tmax - tmin) / h)
  

Note

À part le choix des paramètres, cette fonction est très simple à écrire, il s’agit juste d’écrire en Python les équations mathématiques (5). Elle est de plus très proche de euler().

• Mêmes arguments que la fonction euler(), et donc mêmes conseils quant au choix de ces arguments.
• Exemple : comme plus haut, on étudie encore l’exemple de \(f(y) = -\omega^2 y\).

>>> y0 = 3; z0 = 0
>>> tmin = 0.0; tmax = 3.0
>>> n = 100
>>> h = (tmax - tmin) / float(n - 1)
>>> print("n =", n, "et h =", h)
n = 100 et h = 0.030303030303030304
>>> from math import *
>>> omega = 2 * pi
>>> f = lambda y: - omega**2 * y
>>> y, z = verlet(n, tmin, tmax, y0, z0, f)
>>> # Verification, conditions initiales bien respectees
>>> print(y0 == y[0] and z0 == z[0])
True
>>> print(y)  # 100 points  
[3, 2.945..., 2.784..., ..., 2.928..., 2.751...]
>>> print(z)  # 100 points  
[0, -3.556..., -6.983..., ..., -4.0816..., -7.478...]
>>> # Dernier point, comme vu sur la Figure 1
>>> print(y[-1], z[-1])  
2.751... -7.478...

Centrale_MP__2015.question_ii_c_2_a()

Réponse mathématique à la question II.C.2.a) :

• Il faut tout d’abord d’écrire que \(E(t_i) = \frac{1}{2} \left(y'(t_i)\right)^2 + \frac{1}{2} \omega^2 \left(y(t_i)\right)^2\), donc \(E(t_i) \simeq \frac{1}{2} \left(z_i\right)^2 + \frac{1}{2} \omega^2 \left(y_i\right)^2 =: E_i\) (là encore, c’est ainsi qu’on définit \(E_i\)).
• Et donc, on peut écrire \(E_{i+1}\) en utilisant les équations de récurrence (5), et le développer pour obtenir l’expression demandée, en utilisant le fait que \(f(y) = -\omega^2 y\) :

\[\begin{split}2 E_{i+1} &= \left(z_{i+1}\right)^2 + \omega^2 \left(y_{i+1}\right)^2 \\ &= \left(z_i - \omega^2 \frac{h}{2} (y_i + y_{i+1})\right)^2 + \omega^2 \left(y_i + h z_i - \omega^2 \frac{h^2}{2} y_i \right)^2 \\ &= \left(z_i^2 + \omega^4 \frac{h^2}{4} (y_i + y_{i+1})^2 - 2 z_i \omega^2 \frac{h}{2} (y_i + y_{i+1})\right) \dots \\ & \dots + \omega^2 \left(y_i^2 + h^2 z_i^2 + \omega^4 \frac{h^4}{4} y_i^2 + 2 y_i h z_i - 2 y_i \omega^2 \frac{h^2}{2} y_i - 2 h z_i \omega^2 \frac{h^2}{2} y_i \right) \\ &= 2 E_i + \dots\\ &\text{A finir de developper vous-meme !} \\ &= 2 E_i + h^3 \bigl( \text{constante} + h \times \dots \bigr)\\ &= 2 E_i + O(h^3)\end{split}\]

• Le terme en \(h\) se simplifie,
• Puis on redéveloppe \(y_{i+1}\), et le terme en \(h^2\) se simplifie encore.
• Finalement en regroupant les termes selon leur puissance de \(h\), on trouve le résultat demandé : \(E_{i+1} - E_i = O(h^3)\) (le schéma est presque conservatif en énergie).

Note

Ces calculs sont très similaires à ceux effectués plus haut, pour la question II.B.3.b) question_ii_b_3_b(), ils sont juste plus longs (et pas forcément palpitants, j’espère que vous ne m’en voudrez pas de ne pas détailler entièrement ici – j’ai passé les concours Centrale quand j’étais à votre place, une fois ça me suffit).

Centrale_MP__2015.question_ii_c_3_b()

Réponse mathématique à la question II.C.2.b) :

• Le graphe de portrait de phase obtenu avec le schéma de Verlet confirme le calcul effectué plus haut, \((E_i)_{i \in J_n}\) est presque constante, et donc la norme du vecteur \((y_i, z_i)\), \(R_i = y_i^2 + \omega^2 z_i^2\) augmente très peu.
• C’est exactement ce qu’on observe dans ce portrait de phase (à droite), où il faut comprendre que les points décrivent presque un cercle au cours des itérations (le point initial est \((y_0, z_0) = (3, 0)\)) :

Centrale_MP__2015.question_ii_c_3_c()

Réponse mathématique à la question II.C.2.c) :

• On peut conclure des remarques et des observations précédentes que le schéma de Verlet est plus précis et introduit moins d’erreurs d’approximations que le schéma d’Euler pour cette équation particulière \(y'' = - \omega^2 y\).

Attention

On ne peut rien déduire de l’efficacité d’aucun de ces deux schémas dans un cadre plus général. On les a juste comparé sur une Équation Différentielle particulière !

Voir aussi

Malgré tout, vous serez ravis d’apprendre que dans un cadre plus général, on peut montrer que le schéma de Verlet est plus précis, et s’applique avec succès à un large spectre d’ED, linéaires ou non. Mais toutes ces considérations sont bien-sûr hors programme.

Centrale_MP__2015.G = 6.67e-11

Constante de la gravitation universelle, \(G = 6.67 \times 10^{-11} \mathrm{N}\cdot\mathrm{m}^2\cdot\mathrm{kg}^{-2}\).

Attention

À ne pas confondre avec la valeur \(g\) (champ d’attraction gravitationnelle de la Terre), qui vaut d’habitude \(g \simeq 9.81 \mathrm{m}\cdot\mathrm{s}^{-2}\).

Note

Vous devez connaître la valeur de cette constante (pour les écrits comme pour les oraux) !

Centrale_MP__2015.question_iii_a_1()

Réponse mathématique à la question III.A.1) :

• Par le principe de sommation des forces et l’expression de chaque force exercée par un corps \(k (k \neq j)\) sur un autre corps \(j\), on a directement :

(6)\[\begin{split}\overrightarrow{F}_j &= \sum_{k, k \neq j} \overrightarrow{F}_{k / j} \\ &= \sum_{k, k \neq j} G \frac{m_j m_k}{r_{jk}^3} \overrightarrow{P_j P_k}\end{split}\]

Centrale_MP__2015.norme2(u)

Calcule la norme 2 d’un vecteur \(u \in \mathbb{R}^3\), selon la formule classique :

\[\|u\|^2 = u_1^2 + u_2^2 + u_3^2 \implies \|u\| = \sqrt{u_1^2 + u_2^2 + u_3^2}\]

• Ici, dans un sujet écrit, on peut utiliser au choix la fonction math.sqrt() (importée depuis le module math), ou bien mettre à la puissance 0.5 :

>>> from math import sqrt  # Premiere facon
>>> sqrt(4)
2.0
>>> import math  # Deuxieme facon
>>> math.sqrt(4)
2.0
>>> 4 ** 0.5
2.0

• Exemples :

>>> u = [1, 0, 0]
>>> norme2(u)
1.0
>>> v = [1, 1, 1]
>>> norme2(v)  # Racine de 3 =  
1.73205...
>>> w = [0, 1, 1]
>>> norme2(w)  # Racine de 2 =  
1.41421...

Centrale_MP__2015.force2(m1, p1, m2, p2)

Calcule le vecteur force exercée par un corps \(k\) sur un autre corps \(j\), d’après la formule :

\[\overrightarrow{F}_{k / j} = G \frac{m_j m_k}{r_{jk}^3} \overrightarrow{P_j P_k}\]

• Le premier corps \(j\) est le corps 1 (cible), décrit par sa masse m1 et son vecteur position p1.
• Le second corps \(k\) est le corps 2 (acteur), décrit par sa masse m2 et son vecteur position p2.
• Renvoie un vecteur à trois coordonnées, représentant les composantes de la force (dans la base de référence), exprimées en Newtons (unité du Système International).
• On pensera à utiliser la fonction préliminaire vdif() pour calculer le vecteur \(\overrightarrow{P_j P_k} = \overrightarrow{\mathcal{O} P_k} - \overrightarrow{\mathcal{O} P_j}\), qui sera un vecteur p_1_2 = vdif(p2, p1) (où \(\mathcal{O}\) désigne l’origine, [0, 0, 0], et p1,p2 les positions des corps 1 et 2, ie. \(k\) et \(j\)).
• On pensera aussi à utiliser la fonction préliminaire smul() pour multiplier ce vecteur par le scalaire \(\alpha = G \frac{m_1 m_2}{r_{jk}^3}\).
• Enfin, on a définit plus haut la fonction norme2() pour calculer cette distance \(r_{jk} = \|\overrightarrow{P_j P_k}\|\).
• Complexité en temps : \(O(1)\), puisqu’on doit effectue juste quelques calculs sur des vecteurs de de \(\mathbb{R}^3\) (ie. une liste de taille 3), pas sur les \(N\) positions ou vitesses.
• Complexité en mémoire : \(O(1)\), parce qu’on n’utilise pour mémoire supplémentaire que des variables scalaires, et un vecteur de \(\mathbb{R}^3\).

Centrale_MP__2015.forceN(j, m, pos)

Calcule le vecteur force \(\overrightarrow{F}_j\), somme des vecteurs forces exercées par les autres corps \(k (k \neq j)\) sur le corps \(j\), d’après la formule (6).

• On pensera à utiliser la fonction préliminaire vsom() pour calculer la somme des vecteurs forces, chacun donnés par force2().
• Il suffit de faire une somme sur k, et un test if k != j: ... (le corps \(j\) n’applique aucune attirance sur lui-même – il n’y pas d’amour propre en physique !).
• Complexité en temps : \(O(N)\), puisqu’on doit calculer une somme de \(N - 1\) forces.
• Complexité en mémoire : \(O(1)\), parce qu’on n’utilise pour mémoire supplémentaire que des variables scalaires, et un vecteur de \(\mathbb{R}^3\) ie. une liste de taille 3.

Centrale_MP__2015.question_iii_b_1()

Réponse mathématique à la question III.B.1) :

• On suppose que l’algorithme itératif de résolution approchée de l’équation différentielle qui régit ces \(N\) corps sera exécuté durant n étapes, pour produire des listes de n+1 positions et vitesses approchées (comme plus haut, par euler() et verlet()).
• Soit \(i \in \{1, \dots, N \}\) l’indice d’un des \(N\) corps.
• Alors après la résolution approchée, la liste position[i] contiendra n+1 listes de taille 3, contenant les positions successives du corps \(i\). Comme expliqué dans l’énoncé, position[i][j] sera un vecteur de \(\mathbb{R}^3\), sous la forme déjà utilisée plus haut : position[i][j] \(= [x_{ij}, y_{ij}, z_{ij}] \simeq [x_{i}(t_j), y_{i}(t_j), z_{i}(t_j)]\).
• Et on a la même chose pour les vitesses : vitesse[i][j] \(= [v_{xij}, v_{yij}, v_{zij}] \simeq [\dot{x}_{i}(t_j), \dot{y}_{i}(t_j), \dot{z}_{i}(t_j)]\).

Centrale_MP__2015.arrondi_liste(pos)

Fonction pour arrondir les \(N \times 3\) valeurs d’une liste de positions / vitesses.

Note

Pas du tout demandée par l’énoncée, c’est juste pour améliorer la lisibilité des exemples traités ensuite.

Centrale_MP__2015.ua = 150000000.0

Constante pour écrire les distantes en unités astronomiques. Il faut juste multiplier les distances par une constante (\(1 \;\mathrm{ua} = 1.5 \times 10^{8} \;\mathrm{km}\)).

Centrale_MP__2015.pos = [[0, 1500000000.0, 0], [-1500000000.0, 0, 0], [0, 0, -1500000000.0]]

Positions artificielles de nos trois points (choix arbitraire pour l’exemple numérique montré plus bas).

Centrale_MP__2015.vit = [[500000.0, 0, 0], [0, 500000.0, 0], [0, 0, 500000.0]]

Vitesses artificielles de nos trois points (choix arbitraire).

Centrale_MP__2015.h = 100

Pas de temps (valeur arbitraire).

Centrale_MP__2015.pos_suiv(m, pos, vit, h)

Utilise le schéma de Verlet (5) pour mettre à jour la liste des positions des \(N\) corps (ie. calcule les positions au temps suivant \(t_{i+1}\)), dans le but de résoudre numériquement le système de \(N\) équations différentielles :

(7)\[\begin{split}\begin{cases} \forall j \in \{1, \dots, N \} : \\ P_j'' = \frac{1}{m_j} \overrightarrow{F}_j = \frac{1}{m_j} \sum_{k, k \neq j} \overrightarrow{F}_{k / j} \end{cases}\end{split}\]

• Cette équation différentielle (7), d’ordre 2, non-linéaire, couplée à \(N\) corps, vient de votre cours de physique (seconde loi de Newton). Consultez cette page Wikipédia si besoin.
• En discrétisant selon les temps \(t_i\) comme avant, et en appliquant le schéma de Verlet pour la mise à jour des positions, cela donne :

\[P_j(t_{i+1}) = P_j(t_i) + h P_j'(t_i) + \frac{h^2}{2} \overrightarrow{F}_j(P_j(t_i))\]

• Et dès lors, on fait comme précédemment dans les schémas d’Euler et de Verlet, on remplace les valeurs réelles aux temps discrétisés (\(P_j(t_i), P_j(t_{i+1})\)) par les valeurs discrètes approchées (\(P_{j,i}, P_{j,i+1}\)).
• Cela donne les équations suivantes, qu’il suffit alors d’écrire en Python (grâce à forceN() appliqué avec chaque indice \(j\) pour calculer l’opérateur \(\overrightarrow{F}_j(P_{j,i})\)) :

\[P_{j,i+1} \; \leftarrow \; P_{j,i} + h P'_{j,i} + \frac{h^2}{2} \overrightarrow{F}_j(P_{j,i})\]

Attention

On pensera à utiliser les fonctions préliminaires smul() et vsom() pour faciliter l’arithmétique sur les vecteurs de \(\mathbb{R}^3\).

• Chaque mise à jour des vitesses coûte \(O(3 \times N) = O(N)\) en mémoire, et \(O(N^2)\) en temps (pour chaque point, on doit calculer une somme de \(N-1\) forces).
• Exemple (avec les valeurs artificielles définies plus haut) :

>>> pos_i_plus_un = pos_suiv(m, pos, vit, h)
>>> pos_i_plus_un = [smul(1.0/ua, p) for p in pos_i_plus_un]
>>> print(arrondi_liste(pos_i_plus_un))  
[[0.333, 10.0, -0.0], [-10.0, 0.333, -0.0], [-0.0, 0.0, -9.667]]

Centrale_MP__2015.etat_suiv(m, pos, vit, h)

Utilise le schéma de Verlet (5) pour mettre à jour la liste des vitesses des \(N\) corps, après avoir utiliser pos_suiv() pour calculer les positions au temps suivant \(t_{i+1}\).

• Cela donne les équations suivantes, qu’il suffit alors d’écrire en Python (grâce à forceN()) :

\[P'_{j,i+1} \; \leftarrow P'_{j,i} + h \times \Bigl( \frac{\overrightarrow{F}_j(P_{j,i}) + \overrightarrow{F}_j(P_{j,i+1})}{2} \Bigr)\]

• Cette fois aussi, la mise à jour des positions (ie. le calcul de \(P_{j,i+1}\)) doit impérativement se faire avant la mise à jour des vitesses : \(P'_{j,i+1}\) utilise en effet \(\overrightarrow{F}_j(P_{j,i+1}) = \overrightarrow{F}_{j,i+1}\).
• Chaque mise à jour des vitesses coûte \(O(3 \times N) = O(N)\) en mémoire, et \(O(N^2)\) en temps (pour chaque point, on doit calculer une somme de \(N-1\) forces).
• Exemple (avec les valeurs artificielles définies plus haut) :

>>> pos_i_plus_un, vit_i_plus_un = etat_suiv(m, pos, vit, h)
>>> print(arrondi_liste(vit_i_plus_un))  
[[500000.0, -0.0, -0.0], [0.0, 500000.0, -0.0], [-0.0, 0.0, 500000.0]]

Note

Le sujet n’était pas parfaitement clair quant à savoir si pos_suiv() doit renvoyer pos_i_plus_un et etat_suiv() doit renvoyer (pos_i_plus_un, vit_i_plus_un) ou juste vit_i_plus_un. J’ai choisi de renvoyer à la fois les nouvelles positions et les nouvelles vitesses.

Centrale_MP__2015.question_iii_b_4_a()

Réponse mathématique à la question III.B.4.a) :

• On voit vaguement une relation linéaire entre \(\log (\tau_N)\) et \(\log(N)\)...
• Le coefficient directeur de la droite semble être environ deux. Vous avez plus l’habitude que moi, et sur papier on peut prendre sa règle et mesurer.
• ... On trouve presque deux, trust me.

Centrale_MP__2015.question_iii_b_4_b()

Réponse mathématique à la question III.B.4.b) :

• D’après la question III.B.4.a) précédente (question_iii_b_4_a()), on a :

\[\log(\tau_N) \simeq \text{constante} + (2 \times \log(N))\]

• Ce qui se traduit donc par une relation exponentielle entre \(\tau_N\) et \(N\), où on peut ignorer l’ordonnée à l’origine (ie. le terme \(\text{constante}\)) :

(8)\[\tau_N = O(N^2)\]

Centrale_MP__2015.question_iii_b_5()

Réponse mathématique à la question III.B.5), a) et b) :

• a) Comme signalé plus haut, chaque mise à jour des vitesses demande une mise à jour des positions, et le calculs de deux vecteurs de forces pour chaque point \(j\), soit deux sommes de \(N-1\) vecteurs de \(\mathbb{R}^3\), ce qui est en temps \(O(N)\) pour chaque point, et donc \(O(N^2)\) en tout. On trouve donc :

(9)\[\tau_N = O(N^2)\]

• b) ... Sans surprise, c’est exactement le même résultat que celui trouvé expérimentalement à la question d’avant (8).

Astuce

À croire que les sujets sont conçus pour que tout fonctionne bien ! C’est dingue, non ?

Centrale_MP__2015.question_iv_a()

Partie IV : Exploitation d’une base de données.

Attention

Vous avez moins l’habitude de ce genre de question, je vous conseille tout particulièrement d’y passer du temps !

• Ci-dessous se trouve le schéma de la base de données utilisée dans cette dernière partie IV.

Note

Notez que le sujet écrivait les deux noms de tables CORPS et ETAT en majuscules, puis date_mesure en minuscule, j’ai choisi de tous les écrire en minuscules, pour plus de cohérence et de lisibilité.

Réponse SQL à la question IV.A) :

• Écrire une requête SQL qui renvoie la liste des masses (ie. la clé masse) de tous les corps étudiés (dans la table corps) :

SELECT masse
FROM corps

Attention

Je ne suis pas très expérimenté en SQL, la correction donnée ici aux questions de la partie IV demandant d’écrire des requêtes SQL peut être fausse. J’ai fait de mon mieux pour apprendre et proposer une solution à ces questions, mais je ne garanti par leur exactitude. Veuillez redoubler d’attention, et, si besoin, n’hésitez pas à me contacter pour me signaler une erreur !

Note

J’ai utilisé la solution rédigée par Arnaud Basson afin de corriger les quelques erreurs, normalement les réponses suivantes doivent être correctes.

Centrale_MP__2015.question_iv_b_1()

Réponse SQL à la question IV.B.1) :

1. (Fourni par l’énoncé) Le nombre de corps présents dans la base est obtenu grâce à la requête (on utilise la fonction d’agrégation COUNT pour compter) :

SELECT COUNT(*)
FROM corps

2. Écrire une requête SQL qui renvoie le nombre de corps qui ont au moins un état connu antérieur à tmin() (on utilise la clé datem, la fonction d’agrégation COUNT et le filtre DISTINCT) :

SELECT COUNT(DISTINCT id_corps)
FROM etat
WHERE datem <= tmin()

Centrale_MP__2015.question_iv_b_2()

Réponse SQL à la question IV.B.2) :

• Écrire une requête SQL qui renvoie, pour chaque corps (dans CORPS), son identifiant (id_corps) et la date (datem) de son dernier état antérieur à tmin() (on utilise la fonction d’agrégation MAX pour choisir la date la plus récente parmi celles antérieure à tmin(), et on regroupe par la clé primaire id_corps) :

SELECT id_corps, MAX(datem)
FROM etat
WHERE datem <= tmin()
GROUP BY id_corps

• Le résultat de cette requête est stocké dans une nouvelle table appellée date_mesure.

Centrale_MP__2015.question_iv_b_3()

Réponse SQL à la question IV.B.3) :

• On peut utiliser à la fois les fonctions masse_min() et arete().
• Il faut récupérer les positions et vitesses des corps à la date date_der donc on a besoin de la jointure des tables etat et date_mesure (cette jointure porte sur les attributs id_corps et date, car on ne veut qu’un seul état par corps, celui à la date date_der); on fait aussi la jointure avec la table corps pour récupérer les masses.

1. Écrire une requête SQL qui renvoie la masse et l’état initial (sous la forme masse, x, y, z, vx, vy, vz) de chaque corps retenu pour participer à la simulation :

SELECT masse, x, y, z, vx, vy, vz
FROM corps
    JOIN date_mesure
        ON corps.id_corps = date_mesure.id_corps
    JOIN etat
        ON corps.id_corps = etat.id_corps
        AND etat.datem = date_mesure.date_der
WHERE masse >= masse_min()
    AND (x BETWEEN -arete()/2 AND arete()/2)
    AND (y BETWEEN -arete()/2 AND arete()/2)
    AND (z BETWEEN -arete()/2 AND arete()/2)
ORDER BY x*x + y*y + z*z

2. Pour classer les corps dans l’ordre croissant par rapport à leur distance à l’origine du référentiel, il suffit d’ajouter cette ligne à la fin de la requête (inutile de prendre une racine carrée) :

ORDER BY x*x + y*y + z*z

Attention

Le sujet parle d’un cube d’arête donnée par arete(), donc le test à faire est \(-\mathrm{arete}()/2 \leq x,y,z \leq \mathrm{arete}()/2\), en considérant la moitié de la longueur arete().

Note

Les réponses SQL sont colorées joliment, grâce à un outil codé en Python, pygments. Notez que toutes les solutions sont en fait générées grâce à Sphinx, lui aussi un outil codé en Python et pour Python.

Centrale_MP__2015.N = 3

Nombre de corps.

Centrale_MP__2015.masse = [10, 20, 30]

Liste des masses des corps sélectionnés par la requête précédente. On utilise des données artificielles, avec N = 3 points.

Centrale_MP__2015.m = [10, 20, 30]

Liste des masses des corps sélectionnés par la requête précédente. On utilise des données artificielles, avec N = 3 points.

Centrale_MP__2015.t0 = 1000

Date en nombre entier de secondes du temps \(t_{\min}\) ou tmin.

Centrale_MP__2015.p0 = [[0, 10, 0], [-10, 0, 0], [0, 0, -10]]

Positions artificielles de nos trois points.

Centrale_MP__2015.v0 = [[500.0, 0, 0], [0, 500.0, 0], [0, 0, 500.0]]

Vitesses artificielles de nos trois points.

Centrale_MP__2015.simulation_verlet(deltat, n)

Fonction Python pour la question IV.C).

• On utilise les deux fonctions pos_suiv() et vit_suiv() écrites à la partie III.
• On suppose avoir des variables globales pour t0, masse, p0, et v0 (définies plus haut).

Attention

Il faut prendre en compte le fait que les positions sont en unités astronomiques et les vitesses en \(\mathrm{km}/\mathrm{h}\). Mais ce n’est pas si dur : il faut juste multiplier les distances par une constante (\(1 \;\mathrm{ua} = 1.5 \times 10^{8} \;\mathrm{km}\)) au début, les garder en kilomètres durant le calcul, et les reconvertir en ua à la fin.

Attention

De même, les vitesses sont en kilomètres par seconde, mais doivent être en mètres par seconde pour utiliser le schéma de Verlet, on doit les multiplier par \(10^{3}\).

• Exemple. Pour commencer, on affiche les valeurs définies plus haut :

>>> print(t0)
1000
>>> print(masse)
[10, 20, 30]
>>> print(p0)
[[0, 10, 0], [-10, 0, 0], [0, 0, -10]]
>>> print(v0)
[[500.0, 0, 0], [0, 500.0, 0], [0, 0, 500.0]]

• Exemple de simulations par le schéma de Verlet :

>>> n = 100  # Nombre de points à calculer !
>>> deltat = 10  # Nombre de secondes entre chaque simulation (inutilisé)
>>> positions = simulation_verlet(deltat, n)
>>> print(positions[0])  # Premières positions (= p0)
[[0.0, 10.0, 0.0], [-10.0, 0.0, 0.0], [0.0, 0.0, -10.0]]
>>> print(arrondi_liste(positions[-1]))  # Dernières positions !  
[[3.367, 10.0, -0.0], [-10.0, 3.367, -0.0], [-0.0, 0.0, -6.633]]

Avertissement

Il était difficile de montrer des exemples pour les requêtes SQL, mais j’ai essayé de montrer des exemples pour toutes les fonctions Python demandées par l’énoncé.

Note

Fin du sujet ! Bon courage pour les écrits !

---

Conseils

• Soignez la présentation de vos copies,
• Essayez d’être très attentif à la syntaxe de Python (il y a eu trop d’erreurs d’indentation et de : manquants),
• Vous devez être plus attentif aux consignes de l’énoncé (certains élèves oublient de donner la complexité, ou répondent seulement à la moitié de la question),
• Comme dans chaque concours/DS, vous devez essayer de “grapiller” des points là où vous peuvez.

Note

Ce sujet durait 3h et on n’a eu que 2h en TD/DS, donc c’est normal de ne pas avoir pu tout faire !

---

Sortie du script

$ python3 Centrale_MP__2015.py
 Correction du TD/DS du mercredi 27 et jeudi 28 janvier (2016).

Voici ci-dessous la documentation de la correction complète du sujet d'écrit d'informatique (Centrale, 2015), pour MP, PC, PSI et TSI.

Ce sujet d'écrit a été posé en devoir écrit non noté (TD/DS) pour le cours d'informatique pour tous en prépa MP au Lycée Lakanal.


Veuillez `consulter le code source pour tous les détails <_modules/Centrale_MP__2015.html>`_ svp.
Ci dessous se trouve simplement la *documentation* de ce programme.


- *Date :* Samedi 16 janvier 2016,
- *Auteur :* Lilian Besson, pour le cours d'informatique pour tous en prépa MP (http://perso.crans.org/besson/infoMP/),
- *Licence :* MIT Licence (http://lbesson.mit-license.org).

-----

**Correction :**



Partie I : Quelques fonctions utilitaires
[1, 2, 3] + [4, 5, 6] = [1, 2, 3, 4, 5, 6]
2 * [1, 2, 3] = [1, 2, 3, 1, 2, 3]

I.C) Arithmetique des listes


Partie II : Etude de schemas numeriques

II.A) Mise en forme du probleme

II.B) Schema d'Euler explicite

II.C) Schema de Verlet


Partie III : Probleme a N corps

III.A) Position du probleme

III.B) Approche numerique
Partie IV : Exploitation d'une base de donnees
positions[0] :
[[0.0, 10.0, 0.0], [-10.0, 0.0, 0.0], [0.0, 0.0, -10.0]]
positions[-1] :
[[3.367, 10.0, -0.0], [-10.0, 3.367, -0.0], [-0.0, 0.0, -6.633]]

Test automatique de toutes les doctests écrites dans la documentation (docstring) de chaque fonction :
Trying:
    pos_i_plus_un, vit_i_plus_un = etat_suiv(m, pos, vit, h)
Expecting nothing
ok
Trying:
    print(arrondi_liste(vit_i_plus_un))  # doctest: +ELLIPSIS
Expecting:
    [[500000.0, -0.0, -0.0], [0.0, 500000.0, -0.0], [-0.0, 0.0, 500000.0]]
ok
Trying:
    y0 = 3; z0 = 0
Expecting nothing
ok
Trying:
    tmin = 0.0; tmax = 3.0
Expecting nothing
ok
Trying:
    n = 100
Expecting nothing
ok
Trying:
    h = (tmax - tmin) / float(n - 1)
Expecting nothing
ok
Trying:
    print("n =", n, "et h =", h)
Expecting:
    n = 100 et h = 0.030303030303030304
ok
Trying:
    from math import *
Expecting nothing
ok
Trying:
    omega = 2 * pi
Expecting nothing
ok
Trying:
    f = lambda y: - omega**2 * y
Expecting nothing
ok
Trying:
    y, z = euler(n, tmin, tmax, y0, z0, f)
Expecting nothing
ok
Trying:
    print(y0 == y[0] and z0 == z[0])
Expecting:
    True
ok
Trying:
    print(y)  # 100 points  # doctest: +ELLIPSIS
Expecting:
    [3, 3.0, 2.891..., ..., 17.787..., 17.905...]
ok
Trying:
    print(z)  # 100 points  # doctest: +ELLIPSIS
Expecting:
    [0, -3.588..., -7.177..., ..., 3.891..., -17.388...]
ok
Trying:
    print(y[-1], z[-1])  # doctest: +ELLIPSIS
Expecting:
    17.905... -17.388...
ok
Trying:
    from math import sqrt  # Premiere facon
Expecting nothing
ok
Trying:
    sqrt(4)
Expecting:
    2.0
ok
Trying:
    import math  # Deuxieme facon
Expecting nothing
ok
Trying:
    math.sqrt(4)
Expecting:
    2.0
ok
Trying:
    4 ** 0.5
Expecting:
    2.0
ok
Trying:
    u = [1, 0, 0]
Expecting nothing
ok
Trying:
    norme2(u)
Expecting:
    1.0
ok
Trying:
    v = [1, 1, 1]
Expecting nothing
ok
Trying:
    norme2(v)  # Racine de 3 =  # doctest: +ELLIPSIS
Expecting:
    1.73205...
ok
Trying:
    w = [0, 1, 1]
Expecting nothing
ok
Trying:
    norme2(w)  # Racine de 2 =  # doctest: +ELLIPSIS
Expecting:
    1.41421...
ok
Trying:
    pos_i_plus_un = pos_suiv(m, pos, vit, h)
Expecting nothing
ok
Trying:
    pos_i_plus_un = [smul(1.0/ua, p) for p in pos_i_plus_un]
Expecting nothing
ok
Trying:
    print(arrondi_liste(pos_i_plus_un))  # doctest: +ELLIPSIS
Expecting:
    [[0.333, 10.0, -0.0], [-10.0, 0.333, -0.0], [-0.0, 0.0, -9.667]]
ok
Trying:
    print(t0)
Expecting:
    1000
ok
Trying:
    print(masse)
Expecting:
    [10, 20, 30]
ok
Trying:
    print(p0)
Expecting:
    [[0, 10, 0], [-10, 0, 0], [0, 0, -10]]
ok
Trying:
    print(v0)
Expecting:
    [[500.0, 0, 0], [0, 500.0, 0], [0, 0, 500.0]]
ok
Trying:
    n = 100  # Nombre de points à calculer !
Expecting nothing
ok
Trying:
    deltat = 10  # Nombre de secondes entre chaque simulation (inutilisé)
Expecting nothing
ok
Trying:
    positions = simulation_verlet(deltat, n)
Expecting nothing
ok
Trying:
    print(positions[0])  # Premières positions (= p0)
Expecting:
    [[0.0, 10.0, 0.0], [-10.0, 0.0, 0.0], [0.0, 0.0, -10.0]]
ok
Trying:
    print(arrondi_liste(positions[-1]))  # Dernières positions !  # doctest: +ELLIPSIS
Expecting:
    [[3.367, 10.0, -0.0], [-10.0, 3.367, -0.0], [-0.0, 0.0, -6.633]]
ok
Trying:
    print(smul(2, [1, 2, 3]))
Expecting:
    [2, 4, 6]
ok
Trying:
    print(vdif([1, 2, 3], [4, 5, 6]))
Expecting:
    [-3, -3, -3]
ok
Trying:
    y0 = 3; z0 = 0
Expecting nothing
ok
Trying:
    tmin = 0.0; tmax = 3.0
Expecting nothing
ok
Trying:
    n = 100
Expecting nothing
ok
Trying:
    h = (tmax - tmin) / float(n - 1)
Expecting nothing
ok
Trying:
    print("n =", n, "et h =", h)
Expecting:
    n = 100 et h = 0.030303030303030304
ok
Trying:
    from math import *
Expecting nothing
ok
Trying:
    omega = 2 * pi
Expecting nothing
ok
Trying:
    f = lambda y: - omega**2 * y
Expecting nothing
ok
Trying:
    y, z = verlet(n, tmin, tmax, y0, z0, f)
Expecting nothing
ok
Trying:
    print(y0 == y[0] and z0 == z[0])
Expecting:
    True
ok
Trying:
    print(y)  # 100 points  # doctest: +ELLIPSIS
Expecting:
    [3, 2.945..., 2.784..., ..., 2.928..., 2.751...]
ok
Trying:
    print(z)  # 100 points  # doctest: +ELLIPSIS
Expecting:
    [0, -3.556..., -6.983..., ..., -4.0816..., -7.478...]
ok
Trying:
    print(y[-1], z[-1])  # doctest: +ELLIPSIS
Expecting:
    2.751... -7.478...
ok
Trying:
    print(vsom([1, 2, 3], [4, 5, 6]))
Expecting:
    [5, 7, 9]
ok
26 items had no tests:
    __main__
    __main__.arrondi_liste
    __main__.f
    __main__.force2
    __main__.forceN
    __main__.question_ii_a_1
    __main__.question_ii_a_2
    __main__.question_ii_b_1
    __main__.question_ii_b_3
    __main__.question_ii_b_3_a
    __main__.question_ii_b_3_b
    __main__.question_ii_b_3_c
    __main__.question_ii_b_3_d
    __main__.question_ii_b_3_e
    __main__.question_ii_c_2_a
    __main__.question_ii_c_3_b
    __main__.question_ii_c_3_c
    __main__.question_iii_a_1
    __main__.question_iii_b_1
    __main__.question_iii_b_4_a
    __main__.question_iii_b_4_b
    __main__.question_iii_b_5
    __main__.question_iv_a
    __main__.question_iv_b_1
    __main__.question_iv_b_2
    __main__.question_iv_b_3
9 items passed all tests:
   2 tests in __main__.etat_suiv
  13 tests in __main__.euler
  11 tests in __main__.norme2
   3 tests in __main__.pos_suiv
   9 tests in __main__.simulation_verlet
   1 tests in __main__.smul
   1 tests in __main__.vdif
  13 tests in __main__.verlet
   1 tests in __main__.vsom
54 tests in 35 items.
54 passed and 0 failed.
Test passed.

Plus de détails sur ces doctests sont dans la documentation de Python:
https://docs.python.org/3/library/doctest.html (en anglais)

Le fichier Python se trouve ici : Centrale_MP__2015.py.