TP5 : Résolution numérique des équations différentielles

On s’intéresse à la méthode d’Euler pour la résolution approchée du pendule pesant harmonique ou amorti, ainsi qu’à d’autres méthodes.

_images/TP5__euler_ordre2_vitesse.png _images/TP5__ptmilieu_vitesse.png

Plusieurs portraits de phases, montrant les trois régimes possibles (sous-critique, critique, sur-critique) :

_images/TP5__15_portraits_de_phase.png

Et une jolie animation d’un pendule pesant idéal en régime sous-critique :

_images/TP5__animationq6__sous-critique.gif

Documentation

TP5 au Lycée Lakanal (sujet rédigé par Arnaud Basson), sur le pendule pesant.

TP5.l = 0.3

Variable globale \(l = 30\) cm (longueur de la tige).

TP5.g = 9.80665

Variable globale \(g = 9.80665 \; m/s^2\) (accélération de la pesanteur à la surface de la Terre).

TP5.f(Y, l=0.3, g=9.80665)[source]

L’opérateur \(f: Y \mapsto f(Y)\) (1), écrit comme une fonction Python.

(1)\[\begin{split}\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)]\end{split}\]
  • Arguments :
    • Y est tableau numpy monodimensionnel (à deux éléments) et renvoie un tableau numpy de même taille. \(Y = [\theta, \dot\theta]\) donc \(Y_0 = \theta, Y_1 = \dot\theta\),
    • (optionnel) g est l’accélération de la pesanteur à la surface de la Terre (\(g = 9.80665 \; m/s^2\)),
    • (optionnel) l est la longueur de la tige (\(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]
TP5.euler(h, tmax, theta0, thetaPoint0)[source]

Résoud numériquement l’équation différentielle d’ordre 2 du pendule simple (2) par la méthode d’Euler pour \(t \in [0, t_{\max}]\), avec un pas de temps \(h > 0\).

(2)\[\frac{d^2 \theta}{d t^2}(t) + \frac{g}{l} \sin(\theta(t)) = 0.\]

Pour utiliser un schéma d’Euler, il nous faut écrire cette équation (ordinaire d’ordre 2, mais non linéaire !) sous la forme \(Y'(t) = f(t, Y)\), où \(Y = [\theta, \dot\theta]\) :

\[f(Y) = [Y_1, - (g/l) \sin(Y_0)]\]

Et dès lors, le schéma de mise à jour d’Euler s’écrit :

\[Y_{i+1} = Y_i + (t_{i+1} - t_i) \times f(t_i, Y_i)\]

Mais dans notre cas, les temps discrets \(t_i\) sont uniformément espacés, donc \(\forall i \in \{0,\dots,n-1\},\; t_{i+1} - t_{i} = h\) (le pas de temps), donc le schéma se simplifie en (3) :

(3)\[Y_{i+1} = Y_i + h \times f(t_i, Y_i)\]
  • L’opérateur \(f: (t, Y) \mapsto f(t, Y)\) est vectoriel (de taille \(2 \times 2\)), et il est donné par la fonction f().
  • Les conditions initiales sont données par theta0 et thetaPoint0.
  • Renvoit un tableau (list) contenant les valeurs de \(u_c\) calculées aux instants \(t_i = i \times h\) (pour \(0 \leq i \leq t_{\max} / h)\).
  • Arguments :
    • h est un pas de temps (\(h > 0\)),
    • tmax est la durée totale de la simulation numérique (\(t_{\max} > 0\)),
    • theta0\(= \theta_0\) est la valeur de \(\theta(t=0)\),
    • et thetaPoint0\(= \dot{\theta_0}\) est la valeur de \(\dot\theta(t=0) = \left(\frac{d\theta}{dt}\right)(t=0)\).
  • Hypothèse : On résoudra l’équation pour \(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 ?  
array([ 0.988...,  4.510...])
>>> maxTheta = max(thetas[:, 0])  # Plus grand angle atteint
>>> maxTheta  
1.277...
>>> maxThetaPoint = max(thetas[:, 1])  # Plus grande vitesse atteinte
>>> maxThetaPoint  
7.013...
  • Courbe \((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) :
_images/TP5__euler_ordre2_position.png
  • Portrait de phase (\((\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 :
_images/TP5__euler_ordre2_vitesse.png
  • Complexité en temps : \(O(n)\), avec \(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 \((n, 2)\), avec \(n = \lceil t_{\max} / h \rceil\) (la taille d’un tableau numpy peut être obtenue avec la fonction 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 !).
TP5.pulsation_propre(l=0.3, g=9.80665)[source]

\(\omega_0 = \sqrt{g / l}\) est la pulsation propre du système harmonique (voir cette page pour plus de détails).

  • Exemple :
>>> g = 9.80665; l = 0.30  # 30 cm
>>> w0 = pulsation_propre(l, g)
>>> w0  
5.717...
TP5.tracerCourbe(listeTemps, theta, methode="par méthode d'Euler")[source]

Trace la courbe \(\theta(t)\), valeur de l’angle \(\theta\) en fonction du temps \(t\).

  • Arguments:
    • listeTemps est un array-like (liste ou numpy array) contenant les valeurs \(t_i\),
    • theta contient les valeurs (approchées ou non) de \(\theta_i = \theta(t_i)\).
_images/TP5__euler_ordre2_position.png
TP5.portraitPhase(thetas, methode="par méthode d'Euler")[source]

Trace le portrait de phase, soit la dérivée \(\dot\theta\) en fonction de la valeur de l’angle \(\theta\).

  • Argument:
    • theta est un array-like (liste ou numpy array) de taille \((n, 2)\), qui contient les valeurs (approchées ou non) de \(Y = [\theta, \dot\theta]\) aux instants \(t_i\).
_images/TP5__euler_ordre2_vitesse.png
TP5.pointMilieu(h, tmax, theta0, thetaPoint0)[source]

Résoud numériquement l’équation différentielle d’ordre 2 du pendule simple (4) par la méthode du point milieu pour \(t \in [0, t_{\max}]\) avec un pas de temps \(h\).

(4)\[\frac{d^2 \theta}{d t^2}(t) + \frac{g}{l} \sin(\theta(t)) = 0.\]

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 \(Y'(t) = f(t, Y)\), où \(Y = [\theta, \dot\theta]\) :

\[f(Y) = [Y_1, - (g/l) \sin(Y_0)]\]

C’est le même opérateur \(f\) (calculé par la fonction f() ci-dessus). Et dès lors, le schéma de mise à jour de la méthode du point milieu s’écrit :

\[\begin{split}\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}\end{split}\]

Mais dans notre cas, les temps discrets \(t_i\) sont uniformément espacés, donc \(\forall i \in \{0,\dots,n-1\},\; t_{i+1} - t_{i} = h\) (le pas de temps), donc le schéma se simplifie en (5) :

(5)\[\begin{split}\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}\end{split}\]
  • L’opérateur \(f: (t, Y) \mapsto f(t, Y)\) est vectoriel (de taille \(2 \times 2\)), et c’est le même que la méthode d’Euler (donné par la fonction f()).
  • Les conditions initiales sont données par theta0 et thetaPoint0.
  • Renvoit un tableau (list) contenant les valeurs de \(u_c\) calculées aux instants \(t_i = i \times h\) (pour \(0 \leq i \leq t_{\max} / h)\).
  • Arguments :
    • h est un pas de temps (\(h > 0\)),
    • tmax est la durée totale de la simulation numérique (\(t_{\max} > 0\)),
    • theta0\(= \theta_0\) est la valeur de \(\theta(t=0)\),
    • et thetaPoint0\(= \dot{\theta_0}\) est la valeur de \(\dot\theta(t=0) = \left(\frac{d\theta}{dt}\right)(t=0)\).
  • Hypothèse : On résoudra l’équation pour \(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 ?  
array([ 1.058...,  1.631...])
>>> maxTheta = max(thetas[:, 0])  # Plus grand angle atteint
>>> maxTheta  
1.104...
>>> maxThetaPoint = max(thetas[:, 1])  # Plus grande vitesse atteinte
>>> maxThetaPoint  
6.0...
  • Courbe \((t, \theta(t))\), cette fois ça ressemble beaucoup plus à une solution harmonique idéale (ni amortie ni amplifiée) :
_images/TP5__ptmilieu_position.png
  • Portrait de phase (\((\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 (4) :
_images/TP5__ptmilieu_vitesse.png
  • Complexités : \(O(n)\) en temps, \(O(n)\) en mémoire, avec \(n = \lceil t_{\max} / h \rceil\).

Voir aussi

D’autres références sur la méthode du point milieu :

  • Ces slides (slide 13),
  • Ce sujet de TP (§ 3),
  • Ces autres slides,
  • (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).
TP5.h = 0.0025

Pas de temps, \(h = t_{\max} / n > 0\) (en seconde). Doit être “relativement petit” par rapport à \(t_{\max}\).

TP5.G(Y, t)[source]

Opérateur \(g: (Y, t) \mapsto g(Y, t)\), ici simplement égal à \(g(Y, t) = f(Y)\) (donné en f()).

En utilisant scipy.integrate.odeint (voir sa doc, issue du module scipy.integrate), on obtient la courbe et le portrait de phase suivant :

_images/TP5__scipy_integrate_odeint_position.png _images/TP5__scipy_integrate_odeint_vitesse.png
TP5.theta0 = 0

Angle initial \(\theta_0 = \theta(0)\) (en radian).

TP5.thetaPoint0 = 6

Vitesse angulaire initiale \(\dot{\theta_0} = \dot\theta(0)\) (en radian par seconde).

TP5.tmax = 5

Temps maximum \(t_{\max}\) (en seconde).

TP5.nbPoints = 2000

Nombre de points calculés \(n\). Doit être de l’ordre du milier pour un graphique correct.

TP5.methode = "par methode d'Euler"

Choix de la méthode, pour changer le titre des graphiques.

TP5.portraitsPhase(kmin=1, kmax=15, tmax=2)[source]

Utilise scipy.integrate.odeint() pour afficher un portrait de phase du pendule simple comportant plusieurs trajectoires.

  • On prendra pour cela \(t_{\max} = 2 s\), \(\theta_0 = 0\), et \(\dot{\theta_0} = k \; \text{rad} / s\), pour \(k \in [1,.., 15]\).

Note

Remarque sur les trois régimes (sous-critique, sur-critique, critique).

  • Lorsque \(\theta_0 = 0\) est assez petit (\(\dot{\theta_0} < \dot{\theta}_{\text{critique}}\)), le pendule effectue des oscillations autour de la position d’équilibre stable \(\theta_{\infty} = 0\) (régime sous-critique),
  • Lorsque \(\theta_0 = 0\) est assez grand (\(\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 \(\theta_0 = 0 = 2\omega_0\) ), le pendule converge infiniment lentement vers la position d’équilibre instable (\(\theta = \pi\)). Ce dernier cas est théorique, il ne s’observe jamais en pratique.
  • Dans nos exemples, \(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  
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ù \(\theta = \pi\) :
_images/TP5__15_portraits_de_phase.png
TP5.periode1(h, thetas)[source]

Calcule une valeur approchée de la période \(T\) du système dynamique supposé en régime sous-critique.

À cette fin, on observera que \(\theta(T /2) = 0\) et \(\theta(t) > 0\) pour \(t \in ]0, T/2[\). Il suffira alors de déterminer de façon approchée le plus petit \(t > 0\) tel que \(\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  
1.1799...
>>> T_ideal = (2 * np.pi) / w0
>>> T_ideal  
1.0989...
TP5.periode2(theta0, thetaPoint0, l=0.3, g=9.80665)[source]

Calcule une valeur approchée de la période \(T\) du système dynamique supposé en régime sous-critique.

On calcule \(\omega_0\), puis \(\cos(\theta_{\max})\) par (6) qu’on inverse pour avoir \(\theta_{\max}\) l’amplitude maximale des oscillations :

(6)\[\cos(\theta_{max}) = \cos(\theta_0) - \frac{\dot{\theta_0}^2}{2 \omega_0^2}.\]

Puis calcule \(T\) via cette formule (7) utilisant une intégrale généralisée :

(7)\[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 scipy.integrate.quad() (voir sa documentation ?), du même module scipy.integrate, dont voici quelques explications). Cette fonction 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, 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 pour plus de détails si cela vous intéresse, ou cette page en anglais, issue de la correction d’un projet que j’avais donné à mes élèves indiens 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à (l’idée est toute bête mais très astucieuse !).

  • Exemple :
>>> theta0 = 0; thetaPoint0 = 6
>>> T2 = periode2(theta0, thetaPoint0)
>>> T2  
1.1891...
>>> g = 9.80665; l = 0.30  # 30 cm
>>> w0 = pulsation_propre(l, g)
>>> T_ideal = (2 * np.pi) / w0
>>> T_ideal  
1.0989...
TP5.animationq6(regime='critique', show=False)[source]

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" (\(\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 (\(\dot{\theta_0} < \dot{\theta}_{\text{critique}}\)), avec des oscillations perpétuelles, et un angle toujours borné (\(0 \leq \theta \leq \theta_{\max} < \pi/2\)) et en dessous de l’horizontale :

_images/TP5__animationq6__sous-critique.gif
  • Exemple de régime presque critique, malgré le choix de \(\dot{\theta_0} = \dot{\theta}_{\text{critique}} = 2 \omega_0\), le pendule n’atteint pas l’équilibre instable \(\theta = \pi\) (à la première oscillation on y croit presque, mais non) :
_images/TP5__animationq6__critique.gif
  • Exemple de régime sur-critique (\(\dot{\theta_0} > \dot{\theta}_{\text{critique}}\)), avec des oscillations perpétuelles, toujours dans le même sens, \(0 \leq \theta \to +\infty\) qui diverge (au-delà de \(2\pi\), ça compte les tours) :
_images/TP5__animationq6__sur-critique.gif
  • Exemple de régime “nul”, avec \(\theta_0 = 0\) et \(\dot{\theta_0} = 0\) , il n’y aucun mouvement (et c’est donc l’animation la plus excitante du jour...) :
_images/TP5__animationq6__nul.gif
TP5.botafumeiro()[source]

Désormais, on considère le pendule excité.

TODO !

Et une vidéo montrant le Botafumeiro en action (via YouTube) :

Sortie du script

$ python TP5.py

TP5 au Lycee Lakanal
Exercice I : Le pendule pesant
I.1. On se ramene à un systeme du premier ordre.
I.1.a. Mettre cette equation sous forme d'un systeme differentiel d'ordre 1 :
I.2. Resolution par la methode d'Euler
I.2.a. Voir votre cours pour un rappel de la definition du schema d'Euler...
I.2.b. Pulsation propre
I.2.c Trace de courbe et de portrait de phase
I.3. Des methodes de resolution plus precises.
I.3.a. Methode du point milieu.
I.3.b. La bibliothèque scipy (scientific Python)
I.4. Portrait de phase.
I.5. Determination de la periode des oscillations.
I.6. Réalisation d'une animation.
I.7. Un pendule excité : le botafumeiro.
Test automatique de toutes les doctests ecrites dans la documentation (docstring) de chaque fonction :
Trying:
    theta0 = 0
Expecting nothing
ok
Trying:
    thetaPoint0 = 6
Expecting nothing
ok
Trying:
    tmax = 5
Expecting nothing
ok
Trying:
    nbPoints = 2000  # On choisit 2000 points
Expecting nothing
ok
Trying:
    h = tmax / nbPoints
Expecting nothing
ok
Trying:
    thetas = euler(h, tmax, theta0, thetaPoint0)
Expecting nothing
ok
Trying:
    np.shape(thetas)
Expecting:
    (2000, 2)
ok
Trying:
    thetas[0]  # Vérifions les conditions initiales
Expecting:
    array([ 0.,  6.])
ok
Trying:
    thetas[-1]  # Dernière valeur ?  # doctest: +ELLIPSIS
Expecting:
    array([ 0.988...,  4.510...])
ok
Trying:
    maxTheta = max(thetas[:, 0])  # Plus grand angle atteint
Expecting nothing
ok
Trying:
    maxTheta  # doctest: +ELLIPSIS
Expecting:
    1.277...
ok
Trying:
    maxThetaPoint = max(thetas[:, 1])  # Plus grande vitesse atteinte
Expecting nothing
ok
Trying:
    maxThetaPoint  # doctest: +ELLIPSIS
Expecting:
    7.013...
ok
Trying:
    h = 0.001; tmax = 60
Expecting nothing
ok
Trying:
    n = int(np.ceil(tmax / h))
Expecting nothing
ok
Trying:
    memoire = 4 * 2 * n  # Pour le numpy array de taille (n, 2)
Expecting nothing
ok
Trying:
    print("La procedure euler consommera de l'ordre de %g octets..." % memoire)
Expecting:
    La procedure euler consommera de l'ordre de 480000 octets...
ok
Trying:
    print("Soit environ %g Ko (ca reste tres raisonnable !)." % (memoire / 1024))
Expecting:
    Soit environ 468.75 Ko (ca reste tres raisonnable !).
ok
Trying:
    Y = np.array([0, 10]); print(Y)
Expecting:
    [ 0 10]
ok
Trying:
    l = 0.30; g = 9.80665
Expecting nothing
ok
Trying:
    print(f(Y, l, g))
Expecting:
    [10  0]
ok
Trying:
    Y = np.array([10, 5])
Expecting nothing
ok
Trying:
    print(f(Y, l, g))
Expecting:
    [ 5 17]
ok
Trying:
    l = 0.60  # Une corde plus longue => plus faible vitesse
Expecting nothing
ok
Trying:
    print(f(Y, l, g))
Expecting:
    [5 8]
ok
Trying:
    g = 9.80665; l = 0.30  # 30 cm
Expecting nothing
ok
Trying:
    w0 = pulsation_propre(l, g)
Expecting nothing
ok
Trying:
    k_critique = 2 * w0
Expecting nothing
ok
Trying:
    theta0 = 0; thetaPoint0 = 6
Expecting nothing
ok
Trying:
    assert thetaPoint0 < k_critique, "Erreur le regime n'est pas sous-critique !"
Expecting nothing
ok
Trying:
    tmax = 5
Expecting nothing
ok
Trying:
    nbPoints = 1000  # On choisit 2000 points
Expecting nothing
ok
Trying:
    h = tmax / nbPoints
Expecting nothing
ok
Trying:
    thetas = pointMilieu(h, tmax, theta0, thetaPoint0)
Expecting nothing
ok
Trying:
    T1 = periode1(h, thetas)
Expecting nothing
ok
Trying:
    T1  # doctest: +ELLIPSIS
Expecting:
    1.1799...
ok
Trying:
    T_ideal = (2 * np.pi) / w0
Expecting nothing
ok
Trying:
    T_ideal  # doctest: +ELLIPSIS
Expecting:
    1.0989...
ok
Trying:
    theta0 = 0; thetaPoint0 = 6
Expecting nothing
ok
Trying:
    T2 = periode2(theta0, thetaPoint0)
Expecting nothing
ok
Trying:
    T2  # doctest: +ELLIPSIS
Expecting:
    1.1891...
ok
Trying:
    g = 9.80665; l = 0.30  # 30 cm
Expecting nothing
ok
Trying:
    w0 = pulsation_propre(l, g)
Expecting nothing
ok
Trying:
    T_ideal = (2 * np.pi) / w0
Expecting nothing
ok
Trying:
    T_ideal  # doctest: +ELLIPSIS
Expecting:
    1.0989...
ok
Trying:
    theta0 = 0
Expecting nothing
ok
Trying:
    thetaPoint0 = 6
Expecting nothing
ok
Trying:
    tmax = 5
Expecting nothing
ok
Trying:
    nbPoints = 2000  # On choisit 2000 points
Expecting nothing
ok
Trying:
    h = tmax / nbPoints
Expecting nothing
ok
Trying:
    thetas = pointMilieu(h, tmax, theta0, thetaPoint0)
Expecting nothing
ok
Trying:
    thetas[0]  # Vérifions les conditions initiales
Expecting:
    array([ 0.,  6.])
ok
Trying:
    thetas[-1]  # Dernière valeur ?  # doctest: +ELLIPSIS
Expecting:
    array([ 1.058...,  1.631...])
ok
Trying:
    maxTheta = max(thetas[:, 0])  # Plus grand angle atteint
Expecting nothing
ok
Trying:
    maxTheta  # doctest: +ELLIPSIS
Expecting:
    1.104...
ok
Trying:
    maxThetaPoint = max(thetas[:, 1])  # Plus grande vitesse atteinte
Expecting nothing
ok
Trying:
    maxThetaPoint  # doctest: +ELLIPSIS
Expecting:
    6.0...
ok
Trying:
    g = 9.80665; l = 0.30  # 30 cm
Expecting nothing
ok
Trying:
    w0 = pulsation_propre(l, g)
Expecting nothing
ok
Trying:
    k_critique = 2 * w0
Expecting nothing
ok
Trying:
    k_critique  # doctest: +ELLIPSIS
Expecting:
    11.434...
ok
Trying:
    g = 9.80665; l = 0.30  # 30 cm
Expecting nothing
ok
Trying:
    w0 = pulsation_propre(l, g)
Expecting nothing
ok
Trying:
    w0  # doctest: +ELLIPSIS
Expecting:
    5.717...
ok
6 items had no tests:
    __main__
    __main__.G
    __main__.animationq6
    __main__.botafumeiro
    __main__.portraitPhase
    __main__.tracerCourbe
7 items passed all tests:
  18 tests in __main__.euler
   7 tests in __main__.f
  13 tests in __main__.periode1
   7 tests in __main__.periode2
  12 tests in __main__.pointMilieu
   4 tests in __main__.portraitsPhase
   3 tests in __main__.pulsation_propre
64 tests in 13 items.
64 passed and 0 failed.
Test passed.

Plus de details sur ces doctests peut etre trouve dans la documentation de Python:
https://docs.python.org/3/library/doctest.html (en anglais)

Le fichier Python se trouve ici : TP5.py.