Code source de TP5

#! /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