#! /usr/bin/env python3
# -*- coding: utf-8; mode: python -*-
""" 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 :**
"""
from __future__ import print_function # Python 2 compatibility
if __name__ == '__main__':
print(__doc__)
# %% Partie I : Quelques fonctions utilitaires
print("\n\nPartie I : Quelques fonctions utilitaires")
# I.A.1)
if __name__ == '__main__':
print("[1, 2, 3] + [4, 5, 6] = [1, 2, 3, 4, 5, 6]")
assert [1, 2, 3] + [4, 5, 6] == [1, 2, 3, 4, 5, 6]
# I.A.2)
if __name__ == '__main__':
print("2 * [1, 2, 3] = [1, 2, 3, 1, 2, 3]")
assert 2 * [1, 2, 3] == [1, 2, 3, 1, 2, 3]
# I.B)
[docs]def smul(alpha, liste):
r""" 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) <./_modules/Centrale_MP__2015.html#smul>`_, avec une boucle et des ``append``, plus proche de la syntaxe au programme.
"""
n = len(liste)
l = []
for i in range(n):
l.append(alpha * liste[i])
return l
# return [alpha * x for x in liste]
# I.C) Arithmétique des listes
print("\nI.C) Arithmetique des listes")
# I.C.1)
[docs]def vsom(liste1, liste2):
r""" 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) <./_modules/Centrale_MP__2015.html#vsom>`_, avec une boucle et des ``append``, plus proche de la syntaxe au programme.
"""
n1 = len(liste1)
n2 = len(liste2)
assert n1 == n2, "vsom(...) : erreur, les deux listes ont des tailles differentes !"
l = []
for i in range(n1):
l.append(liste1[i] + liste2[i])
return l
# return [x1 + x2 for (x1, x2) in zip(liste1, liste2)]
# I.C.2)
[docs]def vdif(liste1, liste2):
r""" 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) <./_modules/Centrale_MP__2015.html#vdif>`_, avec une boucle et des ``append``, plus proche de la syntaxe au programme.
"""
n1 = len(liste1)
n2 = len(liste2)
assert n1 == n2, "vdif(...) : erreur, les deux listes ont des tailles differentes !"
l = []
for i in range(n1):
l.append(liste1[i] - liste2[i])
return l
# return [x1 - x2 for (x1, x2) in zip(liste1, liste2)]
# %% Partie II : Étude de schémas numériques
print("\n\nPartie II : Etude de schemas numeriques")
# II.A) Mise en forme du problème
print("\nII.A) Mise en forme du probleme")
[docs]def question_ii_a_1():
r""" **Partie II : Étude de schémas numériques.**
Réponse mathématique à la question II.A.1) :
Si :math:`z(t) = y'(t)`, alors l'équation *(II.1)* s'écrit :math:`z'(t) = f(y(t))`,
qu'on peut mettre sous forme vectorielle (en :math:`[y, z]`) :
.. math::
:label: equation_diff_S
\begin{cases}
y'(t) &= \; z(t) \\
z'(t) &= \; f(y(t))
\end{cases}
.. tip::
Cette manipulation est très classique, on transforme une équation différentielle scalaire (en :math:`y`) d'ordre 2 *(II.1)* en une équation différentielle vectorielle (en :math:`[y, z]`) d'ordre 1 *(S)* (:eq:`equation_diff_S`).
Il faut savoir faire ce genre de manipulation sans réfléchir, c'est très classique et ça tombe très souvent.
.. seealso:: `La solution du TP5 sur la méthode d'Euler <../../TPs/solutions/TP5.html#TP5.f>`_.
"""
pass
[docs]def question_ii_a_2():
r""" 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)* (:eq:`equation_diff_S`).
Par définition de :math:`z(t)`, on remarque déjà que
.. math::
\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).
Donc comme demandé pour l'équation *(II.3)*, on a déjà l'égalité de gauche :
.. math::
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)* (:eq:`equation_diff_S`), on a :math:`z' = f(y)`, donc :
.. math::
\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).
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 : :math:`\forall i, \; t_{i+1} - t_i = h`) :
.. math::
:label: approximations_riemann_left
\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}
C'est l'approximation fait par la méthode des rectangles gauches de Riemann.
Voir `cette page (en anglais) pour plus de détails <https://mec-cs101-integrals.readthedocs.io/en/latest/integrals.html#integrals.riemann_left>`_,
ou `cette page Wikipédia <https://fr.wikipedia.org/wiki/Somme_de_Riemann>`_.
.. image:: Riemann_sum_left.gif
:target: https://fr.wikipedia.org/wiki/Somme_de_Riemann
:align: center
:scale: 65%
Notez qu'en (:eq:`approximations_riemann_left`) on a aussi rajouté (à la fin de la ligne) que les valeurs discrétisées :math:`y_i` et :math:`z_i` vérifieront :math:`y_i \simeq y(t_i)` et :math:`z_i \simeq z(t_i) = y'(t_i)`.
.. seealso::
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 <https://fr.wikipedia.org/wiki/M%C3%A9thode_d%27Euler#Principe_du_sch.C3.A9ma_num.C3.A9rique>`_ pour un rappel).
"""
pass
# II.B) Schéma d'Euler explicite
print("\nII.B) Schema d'Euler explicite")
# II.B.1)
[docs]def question_ii_b_1():
r""" Réponse mathématique à la question II.B.1) :
On utilise les deux équations *(II.3)* et les deux approximations données plus haut :eq:`approximations_riemann_left` :
.. math::
:label: schema_euler
\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}
Ces deux équations de récurrence, sont couplées mais vectorielles : on peut écrire :math:`X = [y, z]`, et on aura :math:`X_{i+1} \leftarrow X_i + h \times F(X_i)` avec :math:`F(X) = [X_1, f(X_2)]`.
Dès qu'on connaît :math:`h` le pas de temps, :math:`f` l'opérateur de l'équation différentielle initiale, et les deux valeurs initiales :math:`y_0, z_0`, les équations :eq:`schema_euler` donnent bien deux suites :math:`(y_i)_{i \in J_n}` et :math:`(z_i)_{i \in J_n}`, définies de façon unique.
"""
pass
# II.B.2)
[docs]def euler(n, tmin, tmax, y0, z0, f):
r""" Utilisation du schéma d'Euler :eq:`schema_euler`, pour calculer les deux suites :math:`(y_i)_{i \in J_n}` et :math:`(z_i)_{i \in J_n}` en fonction de leurs valeurs initiales :math:`y_0, z_0`, du nombre de points :math:`n`, et de l'intervalle de temps :math:`[t_{\min}, t_{\max}]`, où :math:`t_{\min} =` ``tmin`` et :math:`t_{\max} =` ``tmax``, pour l'équation différentielle :math:`y''(t) = f(y(t))` *(II.1)*.
- Renvoie deux listes ``y`` et ``z`` qui contiennent :math:`n + 1` points, en utilisant le pas de temps :math:`h = \frac{t_{\max} - t_{\min}}{n - 1} > 0`.
- Complexité en *temps* : :math:`O(n)`,
- Complexité en *mémoire* : :math:`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 :eq:`schema_euler`. Pour un exemple très similaire utilisant un autre schéma numérique, voir :py:func:`verlet`.
.. tip:: 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 :py:func:`euler` a besoin d'une fonction ``f``, soit définie plus haut, soit donnée en argument. En partie III, les deux fonctions :py:func:`force2` et :py:func:`forceN` auront besoin d'une constante :py:const:`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 :math:`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 # doctest: +ELLIPSIS
[3, 3.0, 2.891..., ..., 17.787..., 17.905...]
>>> print(z) # 100 points # doctest: +ELLIPSIS
[0, -3.588..., -7.177..., ..., 3.891..., -17.388...]
>>> # Dernier point, comme vu sur la Figure 1
>>> print(y[-1], z[-1]) # doctest: +ELLIPSIS
17.905... -17.388...
.. seealso:: Si besoin, lisez `cette page Wikipédia sur la méthode d'Euler <https://fr.wikipedia.org/wiki/M%C3%A9thode_d%27Euler>`_ (au moins le début), ça ne peut pas faire de mal...
"""
assert tmax > tmin, "euler(...) : tmax <= tmin n'est pas autorise."
assert n > 0, "euler(...) : Pas assez de points ! n doit etre > 0."
# n = int((tmax - tmin) / h)
# Calcule le pas de temps h
h = (tmax - tmin) / float(n - 1)
# Liste des valeurs y_i,z_i (on ne demandait pas d'utiliser des tableaux numpy ici)
y = [y0]
z = [z0]
for i in range(0, n+1):
# Avec des tableaux numpy on pourrait faire directement X.append(X[i] + h*F(X[i]))
y_i = y[i]
z_i = z[i]
y_i_plus_un = y_i + h * z_i
z_i_plus_un = z_i + h * f(y_i)
y.append(y_i_plus_un)
z.append(z_i_plus_un)
return y, z
from math import pi
#: Exemple de valeur de la pulsation propre du système, :math:`\omega = 2 \pi`.
#: ``pi`` peut être importé depuis le module :py:mod:`math` ou :py:mod:`numpy` (ou défini manuellement ``pi = 3.1415.... = 3.141592653589793`` à la précision que vous voulez).
omega0 = 2 * pi
[docs]def f(y, omega=omega0):
r""" Exemple de fonction :math:`f` définissant l'équation différentielle étudiée : :math:`f(y) = -\omega^2 y`.
"""
return - omega**2 * y
# II.B.3)
[docs]def question_ii_b_3():
r""" On s'intéresse maintenant à l'exemple de :math:`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 :math:`y = \theta \simeq 0`.
- On sait la forme qu'auront les solutions et les portraits de phases (des sinusoïdes pour :math:`y(t)`, resp. des ellipses pour :math:`y(y'))`). Mais notez que *nulle part* le sujet ne demande de résoudre la solution (ie. d'expliciter la forme qu'aura :math:`y(t)`), et qu'on n'en a pas besoin !
"""
pass
# II.B.3.a)
[docs]def question_ii_b_3_a():
r""" Réponse mathématique à la question II.B.3.a) :
- On s'inspire du cours de physique pour définir cette énergie :math:`E(t)`, qu'on suppose tout d'abord dépendante du temps :
.. math::
:label: energie_E
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 :math:`E(t) = E = E(t_0)` est constante, ce qui se fait en deux lignes de calcul en dérivant :math:`E(t)` dans l'équation ci-dessus (:eq:`energie_E`).
- Notez qu'il suffit d'écrire :math:`f(y) = - \omega^2 y` et :math:`g = - f'` (par définition), puis d'intégrer pour obtenir :math:`g(y) = - \frac{1}{2} \omega^2 y^2 + C` (et on se fiche de la constante).
- Cette énergie :math:`E` a exactement la forme connue par le cours de physique : énergie cinétique + énergie potentielle.
.. tip:: 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 !).
"""
pass
# II.B.3.b)
[docs]def question_ii_b_3_b():
r""" Réponse mathématique à la question II.B.3.b) :
- Il suffit d'écrire que :math:`E(t_i) = \frac{1}{2} \left(y'(t_i)\right)^2 + \frac{1}{2} \omega^2 \left(y(t_i)\right)^2`, donc :math:`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 :math:`E_i`).
- Et donc, on peut écrire :math:`E_{i+1}` en utilisant les équations de récurrence :eq:`schema_euler`, et le développer pour obtenir l'expression demandée, en utilisant le fait que :math:`f(y) = -\omega^2 y` :
.. math::
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)
- Un terme se simplifie, et on trouve donc le résultat demandé : :math:`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 :math:`O(h^2)`, ce qui est déjà très rassurant : plus le pas de temps :math:`h` est petit (ie. plus on calcule de point :math:`n`), moins on fera d'erreur d'approximations.
- Le schéma de Verlet, étudié plus bas, sera en :math:`O(h^3)`, et donc sera plus précis (voir :py:func:`question_ii_c_2_a` plus bas).
"""
pass
# II.B.3.c)
[docs]def question_ii_b_3_c():
r""" Réponse mathématique à la question II.B.3.c) :
- Un schéma numérique qui satisfait à la conservation de l'énergie :math:`E` satisferait la conservation des énergies *approchées* :math:`E_i`, et donc :
.. math::
\forall i \in J_n, \;\;\;\; E_{i+1} - E_i = 0.
.. attention::
C'est bien ce terme :math:`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 (:py:func:`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 :math:`h^2 \omega^2 \bigl( 2 E_i \bigr) = 0`", ça n'a pas de sens.
"""
pass
# II.B.3.d)
[docs]def question_ii_b_3_d():
r""" 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 :math:`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 :math:`y` et :math:`z = y'`, elles n'ont pas la même unité et n'ont aucune raison d'avoir le même ordre de grandeur.
"""
pass
# II.B.3.e)
[docs]def question_ii_b_3_e():
r""" 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 :math:`(E_i)_{i \in J_n}` était *une suite strictement croissante*, et donc que la norme du point :math:`(\omega y_i, z_i)`, :math:`R_i = z_i^2 + \omega^2 y_i^2` augmente *strictement au cours des itérations* (on multiplie la coordonnée :math:`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 :math:`(y_0, z_0) = (3, 0)` au centre) :
.. image:: Centrale_MP__2015___figure_1.png
:align: center
:scale: 50%
"""
pass
# II.C) Schéma de Verlet
print("\nII.C) Schema de Verlet")
[docs]def verlet(n, tmin, tmax, y0, z0, f):
r""" Utilisation du schéma de Verlet :eq:`schema_verlet`, pour calculer les deux suites :math:`(y_i)_{i \in J_n}` et :math:`(z_i)_{i \in J_n}` en fonction de leurs valeurs initiales :math:`y_0, z_0`, du nombre de points :math:`n`, et de l'intervalle de temps :math:`[t_{\min}, t_{\max}]`, où :math:`t_{\min} =` ``tmin`` et :math:`t_{\max} =` ``tmax``, pour l'équation différentielle :math:`y''(t) = f(y(t))` *(II.1)*.
.. math::
:label: schema_verlet
\forall i \in J_n, \;\; f_i := f(y_i).
\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}
- À noter que cette fois, la mise à jour de la première coordonnée (ie. le calcul de :math:`y_{i+1}`) doit impérativement se faire avant la mise à jour de la seconde coordonnée : :math:`z_{i+1}` utilise en effet :math:`f_{i+1} = f(y_{i+1})`.
- Renvoie deux listes ``y`` et ``z`` qui contiennent :math:`n + 1` points, en utilisant le pas de temps :math:`h = \frac{t_{\max} - t_{\min}}{n - 1} > 0`.
- Complexité en *temps* : :math:`O(n)` comme pour le schéma d'Euler,
- Complexité en *mémoire* : :math:`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 :eq:`schema_verlet`. Elle est de plus très proche de :py:func:`euler`.
- Mêmes arguments que la fonction :py:func:`euler`, et donc mêmes conseils quant au choix de ces arguments.
- Exemple : comme plus haut, on étudie encore l'exemple de :math:`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 # doctest: +ELLIPSIS
[3, 2.945..., 2.784..., ..., 2.928..., 2.751...]
>>> print(z) # 100 points # doctest: +ELLIPSIS
[0, -3.556..., -6.983..., ..., -4.0816..., -7.478...]
>>> # Dernier point, comme vu sur la Figure 1
>>> print(y[-1], z[-1]) # doctest: +ELLIPSIS
2.751... -7.478...
"""
assert tmax > tmin, "verlet(...) : tmax <= tmin n'est pas autorise."
assert n > 0, "verlet(...) : Pas assez de points ! n doit etre > 0."
# n = int((tmax - tmin) / h)
# Calcule le pas de temps h
h = (tmax - tmin) / float(n - 1)
# Liste des valeurs y_i,z_i (on ne demandait pas d'utiliser des tableaux numpy ici)
y = [y0]
z = [z0]
for i in range(0, n+1):
# Avec des tableaux numpy on pourrait faire directement X.append(X[i] + h*F(X[i]))
y_i = y[i]
f_i = f(y_i) # Pas de vecteur des valeurs de f, juste deux constantes f_i
z_i = z[i]
y_i_plus_un = y_i + h * z_i + (h**2 / 2.0) * f_i
f_i_plus_un = f(y_i_plus_un) # Et f_i_plus_un
z_i_plus_un = z_i + h * (f_i + f_i_plus_un) / 2.0
y.append(y_i_plus_un)
z.append(z_i_plus_un)
return y, z
# II.C.2.a)
[docs]def question_ii_c_2_a():
r""" Réponse mathématique à la question II.C.2.a) :
- Il faut tout d'abord d'écrire que :math:`E(t_i) = \frac{1}{2} \left(y'(t_i)\right)^2 + \frac{1}{2} \omega^2 \left(y(t_i)\right)^2`, donc :math:`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 :math:`E_i`).
- Et donc, on peut écrire :math:`E_{i+1}` en utilisant les équations de récurrence :eq:`schema_verlet`, et le développer pour obtenir l'expression demandée, en utilisant le fait que :math:`f(y) = -\omega^2 y` :
.. math::
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)
- Le terme en :math:`h` se simplifie,
- Puis on redéveloppe :math:`y_{i+1}`, et le terme en :math:`h^2` se simplifie encore.
- Finalement en regroupant les termes selon leur puissance de :math:`h`, on trouve le résultat demandé : :math:`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) :py:func:`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 <http://perso.crans.org/besson/cv.fr.pdf#page=2>`_, une fois *ça me suffit*).
"""
pass
# II.C.2.b)
[docs]def question_ii_c_3_b():
r""" 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, :math:`(E_i)_{i \in J_n}` est presque constante, et donc la norme du vecteur :math:`(y_i, z_i)`, :math:`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 :math:`(y_0, z_0) = (3, 0)`) :
.. image:: Centrale_MP__2015___figure_2.png
:align: center
:scale: 50%
"""
pass
# II.C.2.c)
[docs]def question_ii_c_3_c():
r""" 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** :math:`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** !
.. seealso::
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.
"""
pass
# %% Partie III : Problème à N corps
print("\n\nPartie III : Probleme a N corps")
#: Constante de la gravitation universelle, :math:`G = 6.67 \times 10^{-11} \mathrm{N}\cdot\mathrm{m}^2\cdot\mathrm{kg}^{-2}`.
#:
#: .. attention:: À ne pas confondre avec la valeur :math:`g` (champ d'attraction gravitationnelle de la Terre), qui vaut d'habitude :math:`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) !
#:
G = 6.67e-11
# III.A) Position du problème
print("\nIII.A) Position du probleme")
# III.A.1)
[docs]def question_iii_a_1():
r""" 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 :math:`k (k \neq j)` sur un autre corps :math:`j`, on a directement :
.. math::
:label: sommeforces
\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}
"""
pass
from math import sqrt
[docs]def norme2(u):
r""" Calcule la norme 2 d'un vecteur :math:`u \in \mathbb{R}^3`, selon la formule classique :
.. math::
\|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 :py:func:`math.sqrt` (importée depuis le module :py:mod:`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 = # doctest: +ELLIPSIS
1.73205...
>>> w = [0, 1, 1]
>>> norme2(w) # Racine de 2 = # doctest: +ELLIPSIS
1.41421...
"""
somme_carree = sum(ui ** 2 for ui in u)
return sqrt(somme_carree)
# III.A.2)
[docs]def force2(m1, p1, m2, p2):
r""" Calcule le vecteur force exercée par un corps :math:`k` sur un autre corps :math:`j`, d'après la formule :
.. math::
\overrightarrow{F}_{k / j} = G \frac{m_j m_k}{r_{jk}^3} \overrightarrow{P_j P_k}
- Le premier corps :math:`j` est le corps 1 (cible), décrit par sa masse ``m1`` et son vecteur position ``p1``.
- Le second corps :math:`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 <https://fr.wikipedia.org/wiki/Syst%C3%A8me_international_d%27unit%C3%A9s>`_).
- On pensera à utiliser la fonction préliminaire :py:func:`vdif` pour calculer le vecteur :math:`\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ù :math:`\mathcal{O}` désigne l'origine, ``[0, 0, 0]``, et ``p1,p2`` les positions des corps 1 et 2, ie. :math:`k` et :math:`j`).
- On pensera aussi à utiliser la fonction préliminaire :py:func:`smul` pour multiplier ce vecteur par le scalaire :math:`\alpha = G \frac{m_1 m_2}{r_{jk}^3}`.
- Enfin, on a définit plus haut la fonction :py:func:`norme2` pour calculer cette distance :math:`r_{jk} = \|\overrightarrow{P_j P_k}\|`.
- Complexité en *temps* : :math:`O(1)`, puisqu'on doit effectue juste quelques calculs sur des vecteurs de de :math:`\mathbb{R}^3` (ie. une liste de taille 3), pas sur les :math:`N` positions ou vitesses.
- Complexité en *mémoire* : :math:`O(1)`, parce qu'on n'utilise pour mémoire supplémentaire que des variables scalaires, et un vecteur de :math:`\mathbb{R}^3`.
"""
p_1_2 = vdif(p2, p1) # Vecteur P_j P_k
r_1_2 = norme2(p_1_2) # Distance r_j_k
alpha = (G * m1 * m2 / (r_1_2 ** 3)) # Scalaire
return smul(alpha, p_1_2)
# III.A.3)
[docs]def forceN(j, m, pos):
r""" Calcule le vecteur force :math:`\overrightarrow{F}_j`, somme des vecteurs forces exercées par les autres corps :math:`k (k \neq j)` sur le corps :math:`j`, d'après la formule :eq:`sommeforces`.
- On pensera à utiliser la fonction préliminaire :py:func:`vsom` pour calculer la somme des vecteurs forces, chacun donnés par :py:func:`force2`.
- Il suffit de faire une somme sur ``k``, et un test ``if k != j: ...`` (le corps :math:`j` n'applique aucune attirance sur lui-même -- il n'y pas d'amour propre en physique !).
- Complexité en *temps* : :math:`O(N)`, puisqu'on doit calculer une somme de :math:`N - 1` forces.
- Complexité en *mémoire* : :math:`O(1)`, parce qu'on n'utilise pour mémoire supplémentaire que des variables scalaires, et un vecteur de :math:`\mathbb{R}^3` ie. une liste de taille 3.
"""
N = len(m)
mj = m[j]
pj = pos[j]
force_totale = [0] * 3
for k in range(N):
if k != j:
mk = m[k]
pk = pos[k]
force_k_sur_j = force2(mj, pj, mk, pk)
force_totale = vsom(force_totale, force_k_sur_j)
# Très inefficace, on ferait mieux d'utiliser des tableaux numpy !
return force_totale
# III.B) Approche numérique
print("\nIII.B) Approche numerique")
# III.B.1)
[docs]def question_iii_b_1():
r""" 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 :math:`N` corps sera exécuté durant ``n`` étapes, pour produire des listes de ``n+1`` positions et vitesses approchées (comme plus haut, par :py:func:`euler` et :py:func:`verlet`).
- Soit :math:`i \in \{1, \dots, N \}` l'indice d'un des :math:`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 :math:`i`. Comme expliqué dans l'énoncé, ``position[i][j]`` sera un vecteur de :math:`\mathbb{R}^3`, sous la forme déjà utilisée plus haut : ``position[i][j]`` :math:`= [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]`` :math:`= [v_{xij}, v_{yij}, v_{zij}] \simeq [\dot{x}_{i}(t_j), \dot{y}_{i}(t_j), \dot{z}_{i}(t_j)]`.
"""
pass
[docs]def arrondi_liste(pos):
r""" Fonction pour arrondir les :math:`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.
"""
return [[round(s, 3) for s in p] for p in pos]
#: Liste artificielle de masses de ``N = 3`` corps (valeurs arbitraires).
#: Notez que, bien évidemment, les masses doivent toutes être positives !
masse = m = [10, 20, 30]
#: Constante pour écrire les distantes en `unités astronomiques <https://fr.wikipedia.org/wiki/Unit%C3%A9_astronomique>`_.
#: Il faut juste multiplier les distances par une constante (:math:`1 \;\mathrm{ua} = 1.5 \times 10^{8} \;\mathrm{km}`).
ua = 1.5e8
#: Positions artificielles de nos trois points (choix arbitraire pour l'exemple numérique montré plus bas).
pos = [
[0, ua*10, 0],
[-ua*10, 0, 0],
[0, 0, -ua*10]
]
#: Vitesses artificielles de nos trois points (choix arbitraire).
vit = [
[5e5, 0, 0],
[0, 5e5, 0],
[0, 0, 5e5]
]
#: Pas de temps (valeur arbitraire).
h = 100
# III.B.2)
[docs]def pos_suiv(m, pos, vit, h):
r""" Utilise le schéma de Verlet :eq:`schema_verlet` pour mettre à jour la liste des positions des :math:`N` corps (ie. calcule les positions au temps suivant :math:`t_{i+1}`), dans le but de résoudre numériquement le système de :math:`N` équations différentielles :
.. math::
:label: equation_dynamique_complete
\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}
- Cette équation différentielle :eq:`equation_dynamique_complete`, d'ordre 2, non-linéaire, couplée à :math:`N` corps, vient de votre cours de physique (seconde loi de Newton). Consultez `cette page Wikipédia si besoin <https://fr.wikipedia.org/wiki/Principe_fondamental_de_la_dynamique>`_.
- En discrétisant selon les temps :math:`t_i` comme avant, et en appliquant le schéma de Verlet pour la mise à jour des positions, cela donne :
.. math::
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 (:math:`P_j(t_i), P_j(t_{i+1})`) par les valeurs discrètes approchées (:math:`P_{j,i}, P_{j,i+1}`).
- Cela donne les équations suivantes, qu'il suffit alors d'écrire en Python (grâce à :py:func:`forceN` appliqué avec chaque indice :math:`j` pour calculer l'opérateur :math:`\overrightarrow{F}_j(P_{j,i})`) :
.. math::
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 :py:func:`smul` et :py:func:`vsom` pour faciliter l'arithmétique sur les vecteurs de :math:`\mathbb{R}^3`.
- Chaque mise à jour des vitesses coûte :math:`O(3 \times N) = O(N)` en mémoire, et :math:`O(N^2)` en temps (pour chaque point, on doit calculer une somme de :math:`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)) # doctest: +ELLIPSIS
[[0.333, 10.0, -0.0], [-10.0, 0.333, -0.0], [-0.0, 0.0, -9.667]]
"""
# print("pos_suiv(m, pos, vit, h)")
# print("m =", m, "pos =", pos, "vit =", vit, "h =", h)
N = len(m)
pos_i_plus_un = [[] for j in range(N)]
for j in range(N):
delta_1_pos_j = smul(h, vit[j])
F_j_t_i = forceN(j, m, pos)
delta_2_pos_j = smul(h**2 / (2.0 * m[j]), F_j_t_i)
pos_i_plus_un[j] = vsom(pos[j], vsom(delta_1_pos_j, delta_2_pos_j))
return pos_i_plus_un
# III.B.3)
[docs]def etat_suiv(m, pos, vit, h):
r""" Utilise le schéma de Verlet :eq:`schema_verlet` pour mettre à jour la liste des vitesses des :math:`N` corps, après avoir utiliser :py:func:`pos_suiv` pour calculer les positions au temps suivant :math:`t_{i+1}`.
- Cela donne les équations suivantes, qu'il suffit alors d'écrire en Python (grâce à :py:func:`forceN`) :
.. math::
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 :math:`P_{j,i+1}`) doit impérativement se faire avant la mise à jour des vitesses : :math:`P'_{j,i+1}` utilise en effet :math:`\overrightarrow{F}_j(P_{j,i+1}) = \overrightarrow{F}_{j,i+1}`.
- Chaque mise à jour des vitesses coûte :math:`O(3 \times N) = O(N)` en mémoire, et :math:`O(N^2)` en temps (pour chaque point, on doit calculer une somme de :math:`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)) # doctest: +ELLIPSIS
[[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 :py:func:`pos_suiv` doit renvoyer ``pos_i_plus_un`` et :py:func:`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.
"""
# print("etat_suiv(m, pos, vit, h)")
# print("m =", m, "pos =", pos, "vit =", vit, "h =", h)
N = len(m)
# On calcule les positions suivantes
pos_i_plus_un = pos_suiv(m, pos, vit, h)
# Puis les vitesses suivantes, l'une après l'autre.
vit_i_plus_un = [[] for j in range(N)]
for j in range(N):
F_j_t_i = forceN(j, m, pos)
F_j_t_i_plus_un = forceN(j, m, pos_i_plus_un)
delta_vit_j = vsom(F_j_t_i, F_j_t_i_plus_un)
vit_i_plus_un[j] = vsom(vit[j], smul(h/2.0, delta_vit_j))
return pos_i_plus_un, vit_i_plus_un
# III.B.4.a)
[docs]def question_iii_b_4_a():
r""" Réponse mathématique à la question III.B.4.a) :
.. image:: Centrale_MP__2015___figure_3.png
:align: center
:scale: 80%
- On voit *vaguement* une relation linéaire entre :math:`\log (\tau_N)` et :math:`\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*.
"""
pass
# III.B.4.b)
[docs]def question_iii_b_4_b():
r""" Réponse mathématique à la question III.B.4.b) :
- D'après la question III.B.4.a) précédente (:py:func:`question_iii_b_4_a`), on a :
.. math::
\log(\tau_N) \simeq \text{constante} + (2 \times \log(N))
- Ce qui se traduit donc par une relation exponentielle entre :math:`\tau_N` et :math:`N`, où on peut ignorer l'ordonnée à l'origine (ie. le terme :math:`\text{constante}`) :
.. math::
:label: complexite_etat_suiv_exper
\tau_N = O(N^2)
"""
pass
# III.B.5)
[docs]def question_iii_b_5():
r""" 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 :math:`j`, soit deux sommes de :math:`N-1` vecteurs de :math:`\mathbb{R}^3`, ce qui est en temps :math:`O(N)` pour chaque point, et donc :math:`O(N^2)` en tout. On trouve donc :
.. math::
:label: complexite_etat_suiv
\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 :eq:`complexite_etat_suiv_exper`.
.. tip:: À croire que les sujets sont conçus pour que tout fonctionne bien ! *C'est dingue, non ?*
"""
pass
# %% Partie IV : Exploitation d'une base de données
print("Partie IV : Exploitation d'une base de donnees")
# IV.A)
[docs]def question_iv_a():
r""" **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.
.. image:: Centrale_MP__2015___figure_4.png
:align: center
:scale: 120%
.. 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``) :
.. code-block:: SQL
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 <http://perso.crans.org/besson/contact>`_ 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.
"""
pass
# IV.B.1)
[docs]def question_iv_b_1():
r""" 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) :
.. code-block:: SQL
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``) :
.. code-block:: SQL
SELECT COUNT(DISTINCT id_corps)
FROM etat
WHERE datem <= tmin()
"""
pass
# IV.B.2)
[docs]def question_iv_b_2():
r""" 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``) :
.. code-block:: SQL
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``.
"""
pass
# IV.B.3)
[docs]def question_iv_b_3():
r""" 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 :
.. code-block:: SQL
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) :
.. code-block:: SQL
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 :math:`-\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 <http://pygments.org/languages/>`_. Notez que `toutes les solutions <http://perso.crans.org/besson/infoMP/>`_ sont en fait générées grâce à `Sphinx <http://www.sphinx-doc.org/>`_, lui aussi un outil codé en Python et pour Python.
"""
pass
#: Nombre de corps.
N = 3
#: 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.
masse = m = [10, 20, 30]
#: Date en nombre entier de secondes du temps :math:`t_{\min}` ou ``tmin``.
t0 = 1000
#: Positions artificielles de nos trois points.
p0 = [
[0, 10, 0],
[-10, 0, 0],
[0, 0, -10]
]
#: Vitesses artificielles de nos trois points.
v0 = [
[5e2, 0, 0],
[0, 5e2, 0],
[0, 0, 5e2]
]
# IV.C)
[docs]def simulation_verlet(deltat, n):
r""" Fonction Python pour la question IV.C).
- On utilise les deux fonctions :py:func:`pos_suiv` et :py:func:`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 <https://fr.wikipedia.org/wiki/Unit%C3%A9_astronomique>`_ et les vitesses en :math:`\mathrm{km}/\mathrm{h}`. Mais ce n'est pas si dur : il faut juste multiplier les distances par une constante (:math:`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 :math:`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 ! # doctest: +ELLIPSIS
[[3.367, 10.0, -0.0], [-10.0, 3.367, -0.0], [-0.0, 0.0, -6.633]]
.. warning:: 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 !
"""
ua = 1.5e8
tmin = t0
tmax = t0 + n * deltat
h = (tmax - tmin) / (n - 1)
p0_km = [smul(ua, p) for p in p0]
v0_kmbys = [smul(1e3, v) for v in v0]
positions = [p0_km]
vitesses = [v0_kmbys]
for i in range(1, n+1):
pos = positions[i-1]
vit = vitesses[i-1]
pos_i_plus_un, vit_i_plus_un = etat_suiv(masse, pos, vit, h)
positions.append(pos_i_plus_un)
vitesses.append(vit_i_plus_un)
# Fin de la simulation, on reconvertit en u.a. :
return [[smul(1.0/ua, p) for p in pos] for pos in positions]
if __name__ == '__main__':
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)
print(positions[0]) # Premières positions (= p0)
print("positions[-1] :") # Dernières positions !
print(arrondi_liste(positions[-1])) # Dernières positions !
# %% Fin de la correction.
if __name__ == '__main__':
from doctest import testmod
print("\nTest automatique de toutes les doctests écrites dans la documentation (docstring) de chaque fonction :")
testmod(verbose=True)
testmod()
print("\nPlus de détails sur ces doctests sont dans la documentation de Python:\nhttps://docs.python.org/3/library/doctest.html (en anglais)")
# Fin de Centrale_MP__2015.py