#! /usr/bin/env python3
# -*- coding: utf-8; mode: python -*-
""" TP5 au Lycée Lakanal (sujet rédigé par Arnaud Basson), sur le pendule pesant.
- *Date :* mercredi 13-01 et jeudi 14-01 (*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).
"""
from __future__ import print_function, division # Python 2 compatibility
import numpy as np
import matplotlib.pyplot as plt
# Exercice I : Le pendule pesant
# %% Exercice I.1
if __name__ == '__main__':
print("\nTP5 au Lycee Lakanal\nExercice I : Le pendule pesant")
print("I.1. On se ramene à un systeme du premier ordre.")
print("I.1.a. Mettre cette equation sous forme d'un systeme differentiel d'ordre 1 :")
#: Variable globale :math:`l = 30` cm (longueur de la tige).
l = 0.30
#: Variable globale :math:`g = 9.80665 \; m/s^2` (accélération de la pesanteur à la surface de la Terre).
g = 9.80665
[docs]def f(Y, l=0.30, g=9.80665):
r""" L'opérateur :math:`f: Y \mapsto f(Y)` :eq:`operateur_f`, écrit comme une fonction Python.
.. math::
:label: operateur_f
\frac{dY}{dt} &= \Big[\frac{d\theta}{dt}, \frac{d^2\theta}{dt^2}\Big] \\
&= [\dot\theta, - (g/l) \sin(\theta)] \\
f(Y) &= [Y_1, - (g/l) \sin(Y_0)]
- *Arguments :*
- ``Y`` est tableau ``numpy`` monodimensionnel (à deux éléments) et renvoie un tableau numpy de même taille. :math:`Y = [\theta, \dot\theta]` donc :math:`Y_0 = \theta, Y_1 = \dot\theta`,
- (optionnel) ``g`` est l'accélération de la pesanteur à la surface de la Terre (:math:`g = 9.80665 \; m/s^2`),
- (optionnel) ``l`` est la longueur de la tige (:math:`l = 30` cm dans les applications).
- Exemple :
>>> Y = np.array([0, 10]); print(Y)
[ 0 10]
>>> l = 0.30; g = 9.80665
>>> print(f(Y, l, g))
[10 0]
>>> Y = np.array([10, 5])
>>> print(f(Y, l, g))
[ 5 17]
>>> l = 0.60 # Une corde plus longue => plus faible vitesse
>>> print(f(Y, l, g))
[5 8]
"""
fY = np.zeros_like(Y)
fY[0] = Y[1]
fY[1] = - (g / l) * np.sin(Y[0])
return fY
# %% Exercice I.2
if __name__ == '__main__':
print("I.2. Resolution par la methode d'Euler")
print("I.2.a. Voir votre cours pour un rappel de la definition du schema d'Euler...")
#: Choix de la méthode, pour changer le titre des graphiques.
methode = u"par methode d'Euler"
#: Angle initial :math:`\theta_0 = \theta(0)` (en radian).
theta0 = 0
#: Vitesse angulaire initiale :math:`\dot{\theta_0} = \dot\theta(0)` (en radian par seconde).
thetaPoint0 = 6
#: Temps maximum :math:`t_{\max}` (en seconde).
tmax = 5
#: Nombre de points calculés :math:`n`. Doit être de l'ordre du milier pour un graphique correct.
nbPoints = 2000
#: Pas de temps, :math:`h = t_{\max} / n > 0` (en seconde). Doit être *"relativement petit"* par rapport à :math:`t_{\max}`.
h = tmax / nbPoints
[docs]def euler(h, tmax, theta0, thetaPoint0):
r"""
Résoud numériquement l'équation différentielle d'ordre 2 du pendule simple :eq:`pendule` par `la méthode d'Euler <https://fr.wikipedia.org/wiki/M%C3%A9thode_d%27Euler>`_ pour :math:`t \in [0, t_{\max}]`, avec un pas de temps :math:`h > 0`.
.. math:: \frac{d^2 \theta}{d t^2}(t) + \frac{g}{l} \sin(\theta(t)) = 0.
:label: pendule
Pour utiliser un schéma d'Euler, il nous faut écrire cette équation (ordinaire d'ordre 2, mais *non linéaire !*) sous la forme :math:`Y'(t) = f(t, Y)`, où :math:`Y = [\theta, \dot\theta]` :
.. math:: f(Y) = [Y_1, - (g/l) \sin(Y_0)]
Et dès lors, le schéma de mise à jour d'Euler s'écrit :
.. math:: Y_{i+1} = Y_i + (t_{i+1} - t_i) \times f(t_i, Y_i)
Mais dans notre cas, les temps discrets :math:`t_i` sont *uniformément* espacés,
donc :math:`\forall i \in \{0,\dots,n-1\},\; t_{i+1} - t_{i} = h` (le *pas de temps*), donc le schéma se simplifie en :eq:`schemaeuler` :
.. math:: Y_{i+1} = Y_i + h \times f(t_i, Y_i)
:label: schemaeuler
- L'opérateur :math:`f: (t, Y) \mapsto f(t, Y)` est vectoriel (de taille :math:`2 \times 2`), et il est donné par la fonction :func:`f`.
- Les conditions initiales sont données par ``theta0`` et ``thetaPoint0``.
- Renvoit un tableau (``list``) contenant les valeurs de :math:`u_c` calculées aux instants :math:`t_i = i \times h` (pour :math:`0 \leq i \leq t_{\max} / h)`.
- *Arguments* :
- ``h`` est un pas de temps (:math:`h > 0`),
- ``tmax`` est la durée totale de la simulation numérique (:math:`t_{\max} > 0`),
- ``theta0``:math:`= \theta_0` est la valeur de :math:`\theta(t=0)`,
- et ``thetaPoint0``:math:`= \dot{\theta_0}` est la valeur de :math:`\dot\theta(t=0) = \left(\frac{d\theta}{dt}\right)(t=0)`.
- *Hypothèse* : On résoudra l'équation pour :math:`t \geq 0`.
- Exemple (question I.2.c) :
>>> theta0 = 0
>>> thetaPoint0 = 6
>>> tmax = 5
>>> nbPoints = 2000 # On choisit 2000 points
>>> h = tmax / nbPoints
>>> thetas = euler(h, tmax, theta0, thetaPoint0)
>>> np.shape(thetas)
(2000, 2)
>>> thetas[0] # Vérifions les conditions initiales
array([ 0., 6.])
>>> thetas[-1] # Dernière valeur ? # doctest: +ELLIPSIS
array([ 0.988..., 4.510...])
>>> maxTheta = max(thetas[:, 0]) # Plus grand angle atteint
>>> maxTheta # doctest: +ELLIPSIS
1.277...
>>> maxThetaPoint = max(thetas[:, 1]) # Plus grande vitesse atteinte
>>> maxThetaPoint # doctest: +ELLIPSIS
7.013...
- Courbe :math:`(t, \theta(t))`, on constate déjà que ce n'est pas exactement une solution harmonique (il y a une légère amplification due aux erreurs numériques introduites par l'approximation de Taylor qui justifie le schéma de mise à jour de la méthode d'Euler) :
.. image:: TP5__euler_ordre2_position.png
:align: center
:scale: 80%
- Portrait de phase (:math:`(\theta(t), \dot\theta(t))`), qui n'est pas un cercle et donc on constate que la méthode d'Euler n'est pas exacte pour cette simulation :
.. image:: TP5__euler_ordre2_vitesse.png
:align: center
:scale: 80%
- Complexité en *temps* : :math:`O(n)`, avec :math:`n = \lceil t_{\max} / h \rceil`.
- Complexité en *mémoire* (question I.2.c) : des constantes, plus un tableau numpy de flottants sur 64 bits (``numpy.float64``, soit 4 octets), de taille exactement :math:`(n, 2)`, avec :math:`n = \lceil t_{\max} / h \rceil` (la taille d'un tableau numpy peut être obtenue avec la fonction :func:`numpy.shape`).
- Exemple :
>>> h = 0.001; tmax = 60
>>> n = int(np.ceil(tmax / h))
>>> memoire = 4 * 2 * n # Pour le numpy array de taille (n, 2)
>>> print("La procedure euler consommera de l'ordre de %g octets..." % memoire)
La procedure euler consommera de l'ordre de 480000 octets...
>>> print("Soit environ %g Ko (ca reste tres raisonnable !)." % (memoire / 1024))
Soit environ 468.75 Ko (ca reste tres raisonnable !).
"""
i_max = int(np.ceil(tmax / h))
# On créé un tableau qu'on va remplir petit à petit
thetas = np.zeros((i_max, 2))
# Conditions initiales
thetas[0, 0] = theta0
thetas[0, 1] = thetaPoint0
# Iteration
for i in range(0, i_max - 1):
# On applique la relation d'Euler vectorielle
thetas[i + 1] = thetas[i] + h * f(thetas[i])
return thetas
if __name__ == '__main__':
print("I.2.b. Pulsation propre")
[docs]def pulsation_propre(l=0.30, g=9.80665):
r""" :math:`\omega_0 = \sqrt{g / l}` est la **pulsation propre du système harmonique** (voir `cette page pour plus de détails <https://fr.wikipedia.org/wiki/Oscillateur_harmonique>`_).
- Exemple :
>>> g = 9.80665; l = 0.30 # 30 cm
>>> w0 = pulsation_propre(l, g)
>>> w0 # doctest: +ELLIPSIS
5.717...
"""
return np.sqrt(g / l)
if __name__ == '__main__':
print("I.2.c Trace de courbe et de portrait de phase")
[docs]def tracerCourbe(listeTemps, theta, methode="par méthode d'Euler"):
r""" Trace la courbe :math:`\theta(t)`, valeur de l'angle :math:`\theta` en fonction du temps :math:`t`.
- *Arguments*:
- ``listeTemps`` est un *array-like* (liste ou numpy array) contenant les valeurs :math:`t_i`,
- ``theta`` contient les valeurs (approchées ou non) de :math:`\theta_i = \theta(t_i)`.
.. image:: TP5__euler_ordre2_position.png
:align: center
:scale: 120%
"""
X = listeTemps
Y = theta
plt.figure()
plt.plot(X, Y, 'r-')
plt.grid()
plt.xlabel(r"Valeur discrète du temps $t$")
plt.ylabel(r"Angle $\theta$ (en radian)")
plt.title("Résolution du pendule pesant simple {}.".format(methode))
plt.show()
[docs]def portraitPhase(thetas, methode="par méthode d'Euler"):
r""" Trace le portrait de phase, soit la dérivée :math:`\dot\theta` en fonction de la valeur de l'angle :math:`\theta`.
- *Argument*:
- ``theta`` est un *array-like* (liste ou numpy array) de taille :math:`(n, 2)`, qui contient les valeurs (approchées ou non) de :math:`Y = [\theta, \dot\theta]` aux instants :math:`t_i`.
.. image:: TP5__euler_ordre2_vitesse.png
:align: center
:scale: 120%
"""
plt.figure()
X = thetas[:, 0]
Y = thetas[:, 1]
plt.plot(X, Y, 'b-')
plt.grid()
plt.xlabel(r"Angle $\theta$ (en radian)")
plt.ylabel(r"Vitesse angulaire $\dot\theta$ (en radian par seconde)")
plt.title("Portrait de phase du pendule pesant simple {}.".format(methode))
plt.show()
# %% Exercice I.2
# if __name__ == '__main__':
if False and __name__ == '__main__':
print("Resolution approche du pendule pesant simple par methode d'Euler, exemple I.2.c")
methode = "par méthode d'Euler"
theta0 = 0
thetaPoint0 = 6
tmax = 5
nbPoints = 2000
h = tmax / nbPoints
# On calcule la solution
thetas = euler(h, tmax, theta0, thetaPoint0)
# Affichages
listeTemps = np.linspace(0, tmax, nbPoints)
print("Affichage de la courbe theta(t)...")
tracerCourbe(listeTemps, thetas[:, 0], methode=methode)
print("Affichage du portrait de phase (theta, thetaPoint)... \n ==> Ce n'est pas un vrai cercle !")
portraitPhase(thetas, methode=methode)
print("Remarque : la methode d'Euler n'est pas exacte, elle apporte une erreur a la solution calculee numeriquement !")
print("Termine pour l'exemple I.2.c...")
# %% Exercice I.3
if __name__ == '__main__':
print("I.3. Des methodes de resolution plus precises.")
print("I.3.a. Methode du point milieu.")
[docs]def pointMilieu(h, tmax, theta0, thetaPoint0):
r"""
Résoud numériquement l'équation différentielle d'ordre 2 du pendule simple :eq:`pendule2` par `la méthode du point milieu <http://uel.unisciel.fr/mathematiques/eq_diff/eq_diff_ch04/co/apprendre_ch4_03.html>`_ pour :math:`t \in [0, t_{\max}]` avec un pas de temps :math:`h`.
.. math:: \frac{d^2 \theta}{d t^2}(t) + \frac{g}{l} \sin(\theta(t)) = 0.
:label: pendule2
Pour utiliser cet autre schéma, il nous faut aussi écrire cette équation différentielle (ordinaire d'ordre 2, mais *non linéaire !*) sous la forme :math:`Y'(t) = f(t, Y)`, où :math:`Y = [\theta, \dot\theta]` :
.. math:: f(Y) = [Y_1, - (g/l) \sin(Y_0)]
C'est le même opérateur :math:`f` (calculé par la fonction :func:`f` ci-dessus).
Et dès lors, le *schéma de mise à jour* de *la méthode du point milieu* s'écrit :
.. math::
\begin{cases}
\Delta_i &= \frac{(t_{i+1} - t_i)}{2} f(t_i, Y_i) \\
Y_{i+1} &= Y_i + (t_{i+1} - t_i) \times f(t_i, Y_i + \Delta_i) \\
\end{cases}
Mais dans notre cas, les temps discrets :math:`t_i` sont uniformément espacés,
donc :math:`\forall i \in \{0,\dots,n-1\},\; t_{i+1} - t_{i} = h` (le *pas de temps*), donc le schéma se simplifie en :eq:`schemamilieu` :
.. math::
:label: schemamilieu
\begin{cases}
\Delta_i &= \frac{h}{2} f(t_i, Y_i) \\
Y_{i+1} &= Y_i + h \times f(t_i, Y_i + \Delta_i) \\
\end{cases}
- L'opérateur :math:`f: (t, Y) \mapsto f(t, Y)` est vectoriel (de taille :math:`2 \times 2`), et c'est le même que la méthode d'Euler (donné par la fonction :func:`f`).
- Les conditions initiales sont données par ``theta0`` et ``thetaPoint0``.
- Renvoit un tableau (``list``) contenant les valeurs de :math:`u_c` calculées aux instants :math:`t_i = i \times h` (pour :math:`0 \leq i \leq t_{\max} / h)`.
- *Arguments* :
- ``h`` est un pas de temps (:math:`h > 0`),
- ``tmax`` est la durée totale de la simulation numérique (:math:`t_{\max} > 0`),
- ``theta0``:math:`= \theta_0` est la valeur de :math:`\theta(t=0)`,
- et ``thetaPoint0``:math:`= \dot{\theta_0}` est la valeur de :math:`\dot\theta(t=0) = \left(\frac{d\theta}{dt}\right)(t=0)`.
- *Hypothèse* : On résoudra l'équation pour :math:`t \geq 0`.
- *Résultat* : Elle est d'ordre 2, donc plus précise que la méthode d'Euler qui est seulement d'ordre 1.
- Exemple (question I.3.a) :
>>> theta0 = 0
>>> thetaPoint0 = 6
>>> tmax = 5
>>> nbPoints = 2000 # On choisit 2000 points
>>> h = tmax / nbPoints
>>> thetas = pointMilieu(h, tmax, theta0, thetaPoint0)
>>> thetas[0] # Vérifions les conditions initiales
array([ 0., 6.])
>>> thetas[-1] # Dernière valeur ? # doctest: +ELLIPSIS
array([ 1.058..., 1.631...])
>>> maxTheta = max(thetas[:, 0]) # Plus grand angle atteint
>>> maxTheta # doctest: +ELLIPSIS
1.104...
>>> maxThetaPoint = max(thetas[:, 1]) # Plus grande vitesse atteinte
>>> maxThetaPoint # doctest: +ELLIPSIS
6.0...
- Courbe :math:`(t, \theta(t))`, cette fois ça ressemble beaucoup plus à une solution harmonique idéale (ni amortie ni amplifiée) :
.. image:: TP5__ptmilieu_position.png
:align: center
:scale: 80%
- Portrait de phase (:math:`(\theta(t), \dot\theta(t))`), sur lequel on constate avec plaisir que la méthode du point milieu est, elle, parfaitement exacte pour la résolution de cette équation différentielle :eq:`pendule2` :
.. image:: TP5__ptmilieu_vitesse.png
:align: center
:scale: 80%
- Complexités : :math:`O(n)` en temps, :math:`O(n)` en mémoire, avec :math:`n = \lceil t_{\max} / h \rceil`.
.. seealso:: D'autres références sur la méthode du point milieu :
- `Ces slides <http://wwwens.aero.jussieu.fr/lefrere/master/mni/mncs/cours/equa-diff.pdf#page=13>`_ (slide 13),
- `Ce sujet de TP <http://www.math.u-bordeaux1.fr/~cdossal/Enseignements/L3-AnaNum/TPEDO1.pdf>`_ (§ 3),
- `Ces autres slides <http://maverick.inria.fr/~Nicolas.Holzschuch/cours/class5.pdf>`_,
- (et rien sur Wikipédia, mais... n'hésitez pas à participer à Wikipédia et créer ou modifier une page sur la méthode du point milieu).
"""
i_max = int(np.ceil(tmax / h))
# On créé un tableau qu'on va remplir petit à petit
thetas = np.zeros((i_max, 2))
# Conditions initiales
thetas[0, 0] = theta0
thetas[0, 1] = thetaPoint0
# Iteration
for i in range(0, i_max - 1):
# La méthode du point milieu est définie par les relations
delta_i = (h / 2) * f(thetas[i])
thetas[i + 1] = thetas[i] + h * f(thetas[i] + delta_i)
return thetas
# if __name__ == '__main__':
if False and __name__ == '__main__':
print("Resolution exacte du pendule pesant simple par methode du point milieu, exemple I.3.a")
methode = "par méthode du point milieu"
theta0 = 0
thetaPoint0 = 6
tmax = 5
nbPoints = 2000
h = tmax / nbPoints
# On calcule la solution
thetas = pointMilieu(h, tmax, theta0, thetaPoint0)
# Affichages
listeTemps = np.linspace(0, tmax, nbPoints)
print("Affichage de la courbe theta(t)...")
tracerCourbe(listeTemps, thetas[:, 0], methode=methode)
print("Affichage du portrait de phase (theta, thetaPoint)... \n ==> C'est un vrai cercle !")
portraitPhase(thetas, methode=methode)
print("Remarque : la methode du point milieu est ici exacte, elle n'apporte aucune erreur a la solution calculee numeriquement !")
print("Termine pour l'exemple I.3.a...")
if __name__ == '__main__':
print("I.3.b. La bibliothèque scipy (scientific Python)")
from scipy.integrate import odeint
[docs]def G(Y, t):
r""" Opérateur :math:`g: (Y, t) \mapsto g(Y, t)`, ici simplement égal à :math:`g(Y, t) = f(Y)` (donné en :func:`f`).
En utilisant ``scipy.integrate.odeint`` (`voir sa doc <https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.odeint.html>`_, issue du `module scipy.integrate <https://docs.scipy.org/doc/scipy/reference/integrate.html>`_), on obtient la courbe et le portrait de phase suivant :
.. image:: TP5__scipy_integrate_odeint_position.png
:align: center
:scale: 80%
.. image:: TP5__scipy_integrate_odeint_vitesse.png
:align: center
:scale: 80%
"""
return f(Y)
# if __name__ == '__main__':
if False and __name__ == '__main__':
print("Resolution exacte du pendule pesant simple par scipy.integrate.odeint, exemple I.3.b")
theta0 = 0
thetaPoint0 = 6
thetas0 = (theta0, thetaPoint0)
tmax = 5
# On choisit 1000 points
nbPoints = 1000
listeTemps = np.linspace(0, tmax, nbPoints)
# On calcule la solution
thetas = odeint(G, thetas0, listeTemps)
# Affichages
methode = "par scipy.integrate.odeint (%i points)" % nbPoints
print("Affichage de la courbe theta(t)...")
tracerCourbe(listeTemps, thetas[:, 0], methode=methode)
print("Affichage du portrait de phase (theta, thetaPoint)... \n ==> C'est un vrai cercle !")
portraitPhase(thetas, methode=methode)
print("Remarque : la methode scipy.integrate.odeint est ici exacte, elle n'apporte aucune erreur a la solution calculee numeriquement !")
print("Termine pour l'exemple I.3.b...")
# %% Exercice I.4
if __name__ == '__main__':
print("I.4. Portrait de phase.")
[docs]def portraitsPhase(kmin=1, kmax=15, tmax=2):
r""" Utilise :func:`scipy.integrate.odeint` pour afficher un portrait de phase du pendule simple comportant plusieurs trajectoires.
- On prendra pour cela :math:`t_{\max} = 2 s`, :math:`\theta_0 = 0`, et :math:`\dot{\theta_0} = k \; \text{rad} / s`, pour :math:`k \in [1,.., 15]`.
.. note:: Remarque sur les trois régimes (sous-critique, sur-critique, critique).
- Lorsque :math:`\theta_0 = 0` est assez petit (:math:`\dot{\theta_0} < \dot{\theta}_{\text{critique}}`), le pendule effectue des oscillations autour de la position d'équilibre stable :math:`\theta_{\infty} = 0` (**régime sous-critique**),
- Lorsque :math:`\theta_0 = 0` est assez grand (:math:`\dot{\theta_0} > \dot{\theta}_{\text{critique}}`), il a un mouvement de révolution autour de son point d'attache (**régime sur-critique**),
- Dans le **régime critique** (qui se produit lorsque :math:`\theta_0 = 0 = 2\omega_0` ), le pendule converge infiniment lentement vers la position d'équilibre instable (:math:`\theta = \pi`). Ce dernier cas est théorique, *il ne s'observe jamais en pratique*.
- Dans nos exemples, :math:`2\omega_0 \simeq 11.4`, donc si on choisit des valeurs entières entre 8 et 15 (inclus), 4 seront en régime sous-critique (8,9,10,11), soit 4 cercles ou ellipses, et 4 seront en régime sur-critique (12,13,14,15), soit 4 trajectoires sinusoïdales dans le portrait de phase.
>>> g = 9.80665; l = 0.30 # 30 cm
>>> w0 = pulsation_propre(l, g)
>>> k_critique = 2 * w0
>>> k_critique # doctest: +ELLIPSIS
11.434...
- En plus pour ce graphique, j'ai ajouté le régime critique pour bien le visualiser (en vert sombre). On constate effectivement que le pendule converge vers l'équilibre (instable) où :math:`\theta = \pi` :
.. image:: TP5__15_portraits_de_phase.png
:align: center
:scale: 90%
"""
g = 9.80665
l = 0.30 # 30 cm
w0 = pulsation_propre(l, g)
theta0 = 0
nbPoints = 1000 # On choisit 1000 points
listeTemps = np.linspace(0, tmax, nbPoints)
plt.figure()
plt.hold(True)
# Amélioration du graphique
plt.grid()
plt.xlabel(r"Angle $\theta$ (en radian)")
plt.xticks(np.arange(-np.pi, 8 * np.pi, np.pi), [(r'$%i\pi$' % i) for i in range(-1, 9)])
plt.ylabel(r"Vitesse angulaire $\dot\theta$ (en radian par seconde)")
methode = r"par odeint (%i points, $\dot{\theta_0} = %i..%i$)" % (nbPoints, kmin, kmax)
plt.title("Portraits de phase du pendule pesant simple, {}.".format(methode))
# 15 portraits de phases, cas sous-critique (cercle) puis sur-critique (sinus)
for k in np.arange(kmin, kmax + 1):
thetas = odeint(G, (theta0, k), listeTemps)
regime = "sous" if k < 2 * w0 else "sur"
plt.plot(thetas[:, 0], thetas[:, 1], label=(r"$\dot{\theta_0} = %i$ (%s-critique)" % (k, regime)))
# Cas critique
k = 2 * w0
thetas = odeint(G, (theta0, k), listeTemps)
plt.plot(thetas[:, 0], thetas[:, 1], label=(r"$\dot{\theta_0} = %g$ (critique)" % k), linewidth=3)
# Tous les portraits de phases sont affichés, on termine le graphique
plt.legend(loc='lower right')
plt.hold(False)
plt.show()
# if __name__ == '__main__':
if False and __name__ == '__main__':
print("Affichage de 15 portraits de phases du pendule pesant simple (par scipy.integrate.odeint), exemple I.4")
portraitsPhase(kmin=8, kmax=15, tmax=2)
# %% Exercice I.5
if __name__ == '__main__':
print("I.5. Determination de la periode des oscillations.")
[docs]def periode1(h, thetas):
r"""
Calcule une valeur approchée de la période :math:`T` du système dynamique supposé en régime sous-critique.
À cette fin, on observera que :math:`\theta(T /2) = 0` et :math:`\theta(t) > 0` pour :math:`t \in ]0, T/2[`.
Il suffira alors de déterminer de façon *approchée* le plus petit :math:`t > 0` tel que :math:`\theta(t) = 0` : ::
angles = thetas[1:, 0] # On ignore theta0 = 0
indice_min_theta = np.argmin(np.abs(angles))
T = 2 * imin * h
- Exemple :
>>> g = 9.80665; l = 0.30 # 30 cm
>>> w0 = pulsation_propre(l, g)
>>> k_critique = 2 * w0
>>> theta0 = 0; thetaPoint0 = 6
>>> assert thetaPoint0 < k_critique, "Erreur le regime n'est pas sous-critique !"
>>> # C'est bon, on est en regime sous-critique.
>>> tmax = 5
>>> nbPoints = 1000 # On choisit 2000 points
>>> h = tmax / nbPoints
>>> thetas = pointMilieu(h, tmax, theta0, thetaPoint0)
>>> T1 = periode1(h, thetas)
>>> T1 # doctest: +ELLIPSIS
1.1799...
>>> T_ideal = (2 * np.pi) / w0
>>> T_ideal # doctest: +ELLIPSIS
1.0989...
"""
angles = thetas[1:, 0]
imin = np.argmin(np.abs(angles))
return 2 * imin * h
from scipy.integrate import quad
[docs]def periode2(theta0, thetaPoint0, l=0.30, g=9.80665):
r"""
Calcule une valeur approchée de la période :math:`T` du système dynamique supposé en régime sous-critique.
On calcule :math:`\omega_0`, puis :math:`\cos(\theta_{\max})` par :eq:`cos_theta_max` qu'on inverse pour avoir :math:`\theta_{\max}` l'amplitude maximale des oscillations :
.. math:: \cos(\theta_{max}) = \cos(\theta_0) - \frac{\dot{\theta_0}^2}{2 \omega_0^2}.
:label: cos_theta_max
Puis calcule :math:`T` via cette formule :eq:`integraleT` utilisant une intégrale généralisée :
.. math::
:label: integraleT
T = \sqrt{\frac{8 l}{g}} \left(\int_{0}^{\theta_{\max}} \frac{1}{\sqrt{\cos(\theta) - \cos(\theta_{\max})}} \mathrm{d}\theta\right).
.. note::
On utilisera la fonction :func:`scipy.integrate.quad` (`voir sa documentation ? <https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.quad.html>`_), du `même module <https://docs.scipy.org/doc/scipy/reference/integrate.html>`_ :mod:`scipy.integrate`, dont `voici quelques explications <https://docs.scipy.org/doc/scipy/reference/tutorial/integrate.html#integration-scipy-integrate>`_).
Cette fonction :func:`quad` est une façon très générique de calculer une intégrale numériquement.
Pour des fonctions réelles ou multi-variées, continues ou non, :func:`quad` fonctionne très bien (pour des intégrales définies, mais pas uniquement des domaines bornés).
Voir `cet article Wikipédia sur les méthodes de quadrature de Gauss <https://fr.wikipedia.org/wiki/M%C3%A9thodes_de_quadrature_de_Gauss>`_ pour plus de détails si cela vous intéresse,
ou `cette page en anglais <http://mec-cs101-integrals.readthedocs.io/en/latest/integrals.html#integrals.gaussian_quad>`_, issue de la correction d'un projet que j'avais donné à mes `élèves indiens <http://www.mahindraecolecentrale.edu.in/>`_ l'an dernier (avril 2015).
On peut même intégrer, par ces méthodes, sur des domaines plus compliqués que des simples hyper-rectangles, voir `ces exemples là <http://mec-cs101-integrals.readthedocs.io/en/latest/integrals.html#integrals.nd_quad>`_ (`l'idée est toute bête mais très astucieuse ! <http://scikit-monaco.readthedocs.io/en/latest/tutorials/getting_started.html#complex-integration-volumes>`_).
- Exemple :
>>> theta0 = 0; thetaPoint0 = 6
>>> T2 = periode2(theta0, thetaPoint0)
>>> T2 # doctest: +ELLIPSIS
1.1891...
>>> g = 9.80665; l = 0.30 # 30 cm
>>> w0 = pulsation_propre(l, g)
>>> T_ideal = (2 * np.pi) / w0
>>> T_ideal # doctest: +ELLIPSIS
1.0989...
"""
w0 = pulsation_propre(l, g)
cos_thetaMax = np.cos(theta0) - (thetaPoint0**2) / (2 * w0**2)
thetaMax = np.arccos(cos_thetaMax)
def fonction_f(theta):
r""" Fonction :math:`f(\theta)` qu'on veut intégrer."""
return 1.0 / np.sqrt(np.cos(theta) - np.cos(thetaMax))
val = quad(fonction_f, 0, thetaMax)[0]
T = np.sqrt((8 * l) / g) * val
return T
# %% Exercice I.6
if __name__ == '__main__':
print("I.6. Réalisation d'une animation.")
from matplotlib.animation import FuncAnimation
[docs]def animationq6(regime='critique', show=False):
r"""
Calcule (et affiche si ``show=True``) une animation du mouvement d'un pendule pesant.
- *Arguments :*
- ``regime`` peut valoir ``"sous-critique"``, ``"critique"``, ``"sur-critique"`` ou ``"nul"`` (:math:`\dot{\theta_0} = 0` donc aucun mouvement).
- ``show=False`` par défaut, mais l'animation est affichée avant d'être sauvegardée si ``True``.
- Exemple de **régime sous-critique** (:math:`\dot{\theta_0} < \dot{\theta}_{\text{critique}}`), avec des oscillations perpétuelles, et un angle toujours borné (:math:`0 \leq \theta \leq \theta_{\max} < \pi/2`) et en dessous de l'horizontale :
.. image:: TP5__animationq6__sous-critique.gif
:align: center
:scale: 85%
- Exemple de **régime presque critique**, malgré le choix de :math:`\dot{\theta_0} = \dot{\theta}_{\text{critique}} = 2 \omega_0`, le pendule n'atteint pas l'équilibre instable :math:`\theta = \pi` (à la première oscillation on y croit presque, mais non) :
.. image:: TP5__animationq6__critique.gif
:align: center
:scale: 85%
- Exemple de **régime sur-critique** (:math:`\dot{\theta_0} > \dot{\theta}_{\text{critique}}`), avec des oscillations perpétuelles, toujours dans le même sens, :math:`0 \leq \theta \to +\infty` qui diverge (au-delà de :math:`2\pi`, ça compte les tours) :
.. image:: TP5__animationq6__sur-critique.gif
:align: center
:scale: 85%
- Exemple de régime *"nul"*, avec :math:`\theta_0 = 0` et :math:`\dot{\theta_0} = 0` , il n'y aucun mouvement (et c'est donc l'animation la plus excitante du jour...) :
.. image:: TP5__animationq6__nul.gif
:align: center
:scale: 85%
"""
# Constants de la simulation
theta0 = 0
if regime == 'critique':
thetaPoint0 = 11.434829834034844
elif regime == 'sous-critique':
thetaPoint0 = 4
elif regime == 'sur-critique':
thetaPoint0 = 13
else:
print("Regime inconnu...")
thetaPoint0 = 0
print("Animation du mouvement du pendule en regime %s, thetaPoint0 = %g." % (regime, thetaPoint0))
tmax = 5
h = .05
# Résolution approchée
solution = pointMilieu(h, tmax, theta0, thetaPoint0)
# Affichage du pendule mouvant
fig = plt.figure()
plt.grid()
tige = plt.plot([0, 0], [0, 0], lw=4)
boule = plt.plot([0], [0], 'ko', ms=20)
plt.axis('equal')
plt.xlim([-1.5 * l, 1.5 * l])
plt.ylim([-1.5 * l, 1.5 * l])
plt.title(r"Mouvement d'un pendule pesant idéal en régime %s Ni amorti ni excité, vitesse initiale $\dot{\theta_0} = %g$." % (regime + ".\n", thetaPoint0))
# la variable l désigne la longueur de la tige
def bouger(i):
""" Mise à jour de la tige. """
if i >= len(solution):
return
# le tableau solution stocke les valeurs
theta = solution[i, 0]
# approchées de theta et thetaPoint
tige[0].set_data([0, l * np.sin(theta)], [0, -l * np.cos(theta)])
boule[0].set_data([l * np.sin(theta)], [-l * np.cos(theta)])
# Lancer l'animation !
anim = FuncAnimation(fig, bouger, interval=(1000 * h), repeat=True)
if show:
plt.show()
nomFichier = 'TP5__animationq6__%s' % regime
anim.save(nomFichier + '.gif', writer='imagemagick') # , fps=60
# anim.save(nomFichier+'.mp4', writer='ffmpeg') # , fps=60
print("L'animation est finie, elle a ete sauvegardee dans %s.gif !" % nomFichier)
# if __name__ == '__main__':
if False and __name__ == '__main__':
print("Affichage de l'animation, exemple I.5")
animationq6(regime='nul')
animationq6(regime='sous-critique')
animationq6(regime='critique')
animationq6(regime='sur-critique')
# %% Exercice I.7
if __name__ == '__main__':
print("I.7. Un pendule excité : le botafumeiro.")
[docs]def botafumeiro():
r"""
Désormais, on considère le **pendule excité**.
TODO !
Et une vidéo montrant le Botafumeiro en action (via YouTube) :
.. youtube:: https://www.youtube.com/watch?v=S_s2Rf0Z0eE
"""
raise ValueError("TODO")
# %% Fin de la correction.
if __name__ == '__main__':
from doctest import testmod
print("Test automatique de toutes les doctests ecrites dans la documentation (docstring) de chaque fonction :")
testmod(verbose=True)
# testmod()
print("\nPlus de details sur ces doctests peut etre trouve dans la documentation de Python:\nhttps://docs.python.org/3/library/doctest.html (en anglais)")
# Fin de TP5.py