packages = ["matplotlib", "numpy", "scipy"]

lundi 02 octobre 2023 Commencé à 16:07 Terminé à .........

EXERCICE 1 : Circuit RC série soumis à un échelon de tension

On a \( E=5 \text{ V} \), \( R=1000 \text{ Ω} \) et \(C=100 \text{ nF} \).
Initialement le condensateur est déchargé, on a donc \( u_C(0^-)=0 \). A l'instant \( t=0 \), on ferme l'interrupteur.

1Quelle est l'équation différentielle vérifiée par la tension \( u_C \) à partir de cet instant ?

2Que vaut la tension \(u_C(0)\) aux bornes du condensateur juste après avoir fermé l'interrupteur ? Pourquoi ?

3Au bout de combien de temps environ le régime permanent sera atteint (en fonction des paramètres du problème) ?

Nous allons résoudre numériquement cette équation, en utilisant la méthode d'Euler.
On commence par importer les bibliothèques nécessaires (penser à exécuter le code de chaque cellule en cliquant sur le bouton en bas à droite !) :

import numpy as np import matplotlib.pyplot as plt

4• Définir les constantes nécessaires avec leur valeur dans le Système International d'unités : E, R et C.
• Définir également la durée totale tf sur laquelle on veut résoudre, en fonction de R et C (voir question 2).
• Définir un nombre N de divisions de cette intervalle de temps pour appliquer la méthode d'Euler (on pourra par exemple choisir \(N=10\) pour commencer).
• Définir alors le pas de calul dt en fonction de tf et N (ne pas donner de valeur numérique, écrire le calcul en fonction de tf et N, pour pouvoir ajuster les valeurs de tf ou de N facilement dans la suite de l'activité).

E = R = C = tf = N = dt =

5• Définir une liste numpy t contenant les \(N\) instants \(t_i\) régulièrement espacés entre 0 et tf.^[Note]
• Définir une liste numpy Uc contenant pour l'instant \(N\) zéros. On la remplira par la suite : elle contiendra les valeurs de la tension \(u_C(t_i)\) aux différents instants \(t_i\) de calcul.^[Note]
• Définir une deuxième liste numpy dUc_dt contenant elle aussi pour l'instant \(N\) zéros. On la remplira par la suite : elle contiendra les valeurs de la dérivée \( \frac{\text{d}u_C}{\text{d}t}(t_i) \) aux différents instants \(t_i\) de calcul^[Note].

t = Uc = dUc_dt =

6C'est parti, on arrive au coeur de l'algorithme d'Euler !
• Que vaut la dérivée \( \frac{\text{d}u_C}{\text{d}t} \) en fonction de \( u_C \) (voir question 1) ?
\frac{\text{d}u_c}{\text{d}t}=

• Initialisation : fixer dans Uc[0] la valeur de la condition initiale connue \( u_C(0) \).
• Dans une boucle^[Note], on calcule les valeurs suivantes de \( u_C \) de proche en proche :
   - calculer dUc_dt[i] en fonction de Uc[i] conformément à l'équation différentielle à résoudre,
   - puis calculer Uc[i+1] en fonction de Uc[i] suivant la méthode d'Euler.

... ... ... ...

7Il ne reste plus qu'à tracer le graphe des Uc en ordonnées et des t en abscisses, pour visualiser la solution approchée ainsi obtenue !

fig, ax = plt.subplots() plt.plot(t, Uc, 'o-', markersize=3) # avec des ronds (pas trop gros: taille 3) reliés par des segments plt.title("Tension aux bornes du condensateur") plt.xlabel("Temps en s") plt.ylabel("Tension en V") fig # affichage de la figure

8Puisque l'équation différentielle sur laquelle on travaille est en fait linéaire, il existe une solution exacte. Il serait donc intéressant de comparer la solution approchée obtenue avec la solution exacte !
• Quelle est la solution mathématique exacte de cette équation différentielle ?

u_c(t)=

• Dans la cellule de la question 5, ajouter une ligne pour calculer UcExacte en fonction de t, E, R et C.^[Note]
Puis, dans le même graphe, superposer cette solution exacte en ajoutant la ligne plt.plot(t,UcExacte).
• Exécuter à nouveau le programme en choisissant des nombres d'itérations N différents. Commenter.

---

EXERCICE 2 : Chute d'une goutte de pluie avec frottements quadratiques

Intéressons-nous maintenant un exemple où la résolution numérique a un plus grand intérêt, c'est-à-dire à un problème physique décrit par une équation différentielle qui n'est pas linéaire.
On étudie une goutte de pluie de masse \(m\) qui chute verticalement dans le champ de pesanteur. L'air exerce sur cette goutte des frottements fluides proportionnels au carré de la vitesse : \( \overrightarrow{f}=-\alpha v^2 \overrightarrow{u_x} \).

1Quelle est l'équation différentielle vérifiée par la vitesse \(v(t)\) de la goutte ?

2Il s'agit d'une équation non linéaire, sa résolution analytique est difficile... Notre résolution numérique approchée a donc un intérêt.
Néanmoins, sans résoudre, on peut déterminer la vitesse limite atteinte par la goutte en régime permanent. Quelle est son expression ?

v_{lim}=

A vous de jouer : vous allez écrire un programme Python utilisant la méthode d'Euler pour résoudre cette équation différentielle et connaître la vitesse de la goutte pendant le régime transitoire.
Avant tout, exécuter cette ligne pour effacer les variables définies à l'exercice précédent :

del E, R, C, tf, N, dt, t, Uc, dUc_dt, i

Comme dans le problème ouvert du TD du chapitre C3, on prendra une grosse goutte de pluie supposée sphérique de rayon \(R=4\text{ mm}\).
Sa masse est donc de \( m = \rho\times\frac{4}{3}\pi R^3 = 2,7\cdot 10^{-4}\text{ kg} \).
Le coefficient de frottement vaut \( \alpha = \frac{1}{2}\rho\pi R^2 C_x = \frac{1}{2}\times 1000 \times 3,14 \times \left(0,004\right)^2 \times 0,5 = 1,3 \cdot 10^{-2}\text{ N.s}^2\text{.m}^{-2} \)

3Ecrire le programme Python.
On appellera g la pesanteur terrestre, alpha le coefficient de frottement, v le tableau numpy des vitesses, et a celui des accélérations.
Pour l'intervalle total de temps tf sur lequel résoudre, vous pourrez estimer au doigt mouillé (avec votre sens physique) la durée au bout de laquelle une goutte de pluie initialement immobile atteint sa vitesse de chute limite. Puis, après avoir tracé le graphe, vous pourrez réajuster cette valeur.

... ... ... ... ...

4Afficher le graphe de la vitesse approchée \(v(t)\) obtenue :

fig, ax = plt.subplots() plt.plot(t,v, 'o-', markersize=3) # avec des ronds (pas trop gros: taille 3) reliés par des segments plt.title("Chute avec frottements quadratiques") plt.xlabel("Temps en s") plt.ylabel("Vitesse en m/s") fig # affichage de la figure

5Commenter. (La vitesse limite est elle correcte ? Qu'avez-vous constaté pour certains paramétrages de la discrétisation ?....)

Remarque : La solution exacte, difficile à trouver par le calcul (!) est : \[ v(t) = \sqrt{\frac{mg}{\alpha}} \tanh \left(\frac{gt}{\sqrt{\frac{mg}{\alpha}}}\right) \]

---

EXERCICE 3 : La fonction solve_ivp

En fait, il existe une fonction solve_ivp déjà implémentée dans la bibliothèque scipy qui permet de résoudre numériquement des équations différentielles (mais toujours d'ordre 1 seulement !), avec un algorithme plus performant que celui d'Euler.

1 Pour s'entraîner à l'utiliser sur un exemple simple, reprendre le cas du circuit RC de l'exercice 1, et résoudre numériquement l'équation différentielle avec solve_ivp pour déterminer la tension \(u_C(t)\) aux bornes du condensateur entre \(t=0\) et \(t=5RC\).

# On efface d'abord toutes les variables définies dans l'exercice précédent : del tf, N, dt, t, i, v, a, g, alpha # On importe solve_ivp depuis la bibliothèque scipy.integrate : from scipy.integrate import solve_ivp # On redéfinit les paramètres physiques du problème : E=5 # volts R=1000 # ohms C=100e-9 # farads # A vous de jouer ...... ...... ...... ...... ...... ...... ...... # On trace le graphe (A COMPLETER) fig, ax = plt.subplots() plt.plot(............ , ............ , 'o-' , markersize=3) plt.title("Tension aux bornes du condensateur") plt.xlabel("Temps en s") plt.ylabel("Tension en V") fig

---

EXERCICE 4 : Pendule simple aux grands angles

Cette fois-ci, nous allons nous intéresser à un système physique régit par une équation différentielle qui n'admet pas de solution analytique : le pendule simple.
Les conditions initiales seront les suivantes : \[ \begin{Bmatrix} \theta(0)=\theta_0 \\ \dot{\theta}(0)=0 \end{Bmatrix} \]L'angle initial \(\theta_0\) ne sera pas nécessairement petit devant 1 radian.
La longueur du fil sera de \(L=20\text{ cm}\). La pesanteur vaut \(g=9,8\text{ m.s}^{-2}\). On pourra poser \(\omega_0=\sqrt{\frac{g}{L}}\).

1Quelle est l'équation différentielle vérifiée par l'angle \(\theta(t)\) ?

2En vous aidant de l'énoncé papier, utiliser solve_ivp pour résoudre numériquement cette équation différentielle entre \(t=0\) et \(t=6\times\frac{2\pi}{\omega_0}\).

import numpy as np import matplotlib.pyplot as plt from scipy.integrate import solve_ivp # Fixer les valeurs des paramètres : ...... ...... ...... ...... ...... ...... # Définir le système d'E.D. à résoudre dans une fonction "systeme(t,Y)" ...... ...... ...... ...... ...... ...... ...... # On trace le graphe fig, ax = plt.subplots() plt.plot(............ , ............) # liste des absisses, liste des ordonnées (à compléter) plt.title("........") plt.xlabel("........") plt.ylabel("........") fig

On pourra tester plusieurs valeurs d'angles initial \(\theta_0\), éventuellement très proches de \(\pi\) (\(\theta_0 = 0,5\pi\), \(0,8\pi\), \(0,9\pi\), \(0,99\pi\)...).

3Que pensez-vous de l'allure des solutions pour des "grands" angles de départ ?

4Dans l'approximation des petits angles, on rappelle qu'il y a isochronisme des oscillations. Qu'est-ce que cela signifie ?

Que vaut la période des oscillations, dans ce cas ?

T=

5Est-ce encore le cas quand on n'est plus aux "petits angles" ? Pour le savoir, nous allons superposer les courbes de \(\theta(t)\) sur le même graphe, pour 10 valeurs différentes d'angle initial \(\theta_0\) échelonnées entre \(0\) et \(\pi\).
Quelques conseils pour vous aider :

• Faire un copier-coller de la cellule précédente dans la cellule vide ci-dessous.
• Remplacer la valeur unique la condition initiale theta0 par un tableau numpy contenant 10 valeurs de conditions initiales échelonnées entre \(0\) et \(\pi\).^[Note]
• Insérer l'instruction solution = solve_ivp() dans une boucle, en l'adaptant.
• Déplacer l'instruction plt.plot() à l'intérieur de cette boucle, pour qu'une nouvelle courbe s'affiche à chaque valeur de \(\theta_0\).

Attention : l'instruction fig,ax=plt.subplots() (qui crée la nouvelle zone de graphique) doit être située avant les plt.plot(). La déplacer au début du programme.
Pour mieux voir, on pourra tracer sur une durée entre \(t=0\) et \(t_f=3\frac{2\pi}{\omega_0}\), au lieu de \(t_f=6\frac{2\pi}{\omega_0}\) précédemment.

...... ...... ......

Commenter : y a-t-il isochronisme des oscillations ?

FIN DE LA SEANCE