DM n°7 : Résoudre des équations différentielles en Python, tracé de courbes
A rendre en trinômes pour le lundi 14 mars 2022
Vous me rendrez sur cahier de prépa un unique script Python contenant votre réponse pour les trois
exercices. Je ne vous demande rien d’autre que le code, mais chaque ligne devra être commentée.
Je vous conseille de chercher ces exercices vous-même (pendant les vacances par exemple). Aidez-vous
entre vous, et si vous récupérer le script de quelqu’un d’autre, essayez de le comprendre et de l’utiliser avec
des paramètres différents, pour vous l’approprier.
N’hésitez pas à chercher des renseignements dans l’aide de Python et sur internet. Allez consulter l’annexe
du TP sur les battements pour une liste de fonctions pour tracer les courbes (ou sinon l’aide de pyplot).
Python propose un outil pour résoudre numériquement les équations différentielles. La fonction odeint de la
bibliothèque scipy.integrate permet de résoudre des équations différentielles du premier ordre qui s’écrivent
sous la forme y’=f(y,t). La syntaxe complète de la fonction est décrite dans l’aide Python.
L’appel de cette fonction se fait de la façon suivante :
from scipy.integrate import odeint
y=odeint(f,y0,t)
• f : fonction qui définit l’équation différentielle telle que y’=f(y,t) (attention à l’ordre des paramètres !
y doit être avant t). Il faut se ramener à une équation du premier ordre, qui peut être un système de
k équations. Dans ce cas, y doit être un tableau et f doit retourner un tableau de même dimension.
• y0 : la condition initiale (à t=0). C’est une liste de k éléments
• t : liste de taille N comportant les temps sur lesquels on veut résoudre l’équation (exemple : de 0 à 3s)
• y sera une matrice de k lignes et N colonnes comportant la solution de l’équation différentielle évaluée
aux temps de la liste t.
Exemples :
𝑑𝑦
• Si 𝑑𝑡 = 𝐴𝑦 est l’équation à résoudre, on prend
def f(y,t):
return A*x
𝑑2 𝑦 𝑑𝑦 𝑑𝑌
• Si 𝑎 𝑑𝑡 2 + 𝑏 𝑑𝑡 + 𝑐 = 0 est l’équation à résoudre, on va résoudre le système 𝑑𝑡
= 𝑓(𝑌) avec 𝑌 =
𝑦
(𝑑𝑦)
𝑑𝑡
def f(Y,t):
dYdt=[Y[1],-c/a-b/a*X[0]]
return dYdt
1. Résolution de l’équation différentielle y’=y (dont la solution analytique est simple).
a. Ecrire dans un script Python une fonction nommée f qui à deux variables d’entrée, t et y,
associe une variable de sortie ydot soit ydot=f(y,t). Pour répondre au problème ici : ydot=y.
b. Définir un vecteur temps de 101 valeurs comprises entre 0 et 1 par pas de 0.1. On utilisera
np.linspace
c. Calculer le vecteur ycalcul résultat de l’équation différentielle aux instants définis par le
vecteur temps en considérant comme conditions initiales y(0)=1 et en utilisant la fonction
odeint().
d. Afficher les courbes ycalcul et exponentielle (np.exp) sur un même graphique en fonction du
vecteur temps.
2. Le modèle « proies-prédateurs »
Les équations différentielles de Lotka-Volterra décrivent le modèle de relation de type prédateur-proie entre
deux espèces. Volterra introduisit ce modèle en 1926 pour expliquer l’évolution périodique de deux espèces
de poisson de la mer Adriatique, des équations de la même forme furent introduite en 1922 par Lotka en
cinétique chimique.
Ces équations s’écrivent :
𝑑𝑥(𝑡)
= 𝑥(𝑡)(𝑎 − 𝑏𝑦(𝑡))
𝑑𝑡
𝑑𝑦(𝑡)
= −𝑦(𝑡)(𝑑 − 𝑐𝑥(𝑡))
𝑑𝑡
Avec : x(t) l’effectif des proies
y(t) l’effectif des prédateurs
t le temps
a le taux de reproduction des proies (constant et indépendant du nombre de prédateurs)
b le taux de mortalité des proies dû aux prédateurs rencontrés
c le taux de reproduction des prédateurs
d le taux de mortalité des prédateurs (constant et indépendant du nombre de proies)
Prendre pour débuter a=5/an, b=1/an, c=1/an et d=3/an.
• Définir un vecteur temps qui va de 0 à 10an par pas de 0.1an
• Définir dans un script python la fonction f avec Y=f(X) telle que
Y=[X[0]*(a-b*X[1]),-X[1]*(d-c*X[0])]
• En utilisant cette fonction, résoudre le système d’équations différentielles du modèle proie prédateur
avec la condition initiale y0=[x(0),y(0)]=[2,2].
• Afficher dans la même fenêtre : le graphique représentant les deux populations proies/prédateurs en
fonctions du temps et le graphique représentant la population de prédateurs en fonction de la
population de proies.
NB : Quand on a une matrice M de n lignes et m colonnes, M[ :,i] représente la colonne i et M[j, :]
représente la ligne j.
Pour la représentation graphique, on pourra utiliser plt.subplot(nbre de lignes, nombre de colonnes,
n° de plot). Ici par exemple 2 lignes et 1 colonnes :
plt.subplot(2,1,1)
Plt.plot(…)
plt.subplot(2,1,2)
Plt.plot(…)
On pourra jouer avec les conditions initiales et les paramètres
3. Résolution d’une équation différentielle du 2nd ordre
On va maintenant chercher à utiliser odeint pour résoudre une équation différentielle du 2nd ordre, comme
nous l’avons vu dans le cours pour trouver la solution de l’équation différentielle non linéaire du pendule
simple.
On considère un plateau P de masse m monté sur un ressort de raideur
k=150 N/m et de longueur à vide 𝑙0 = 20𝑐𝑚. L’axe vertical Oz est dirigé
vers le haut, et l’origine z=0 se trouve à la position d’équilibre du
système. On note z(t) la position du plateau à l’instant t.
Lors de son mouvement, le plateau est soumis à une force de frottement
𝑓⃗ = −𝑎𝑣⃗
Le plateau est compressé à t=0 à 𝑧 = −𝑧0 = −5𝑐𝑚
1. Montrer que l’équation différentielle décrivant le système s’écrit :
𝑧̈ + 2𝜉𝜔0 + 𝜔02 𝑧 = 0
𝑘 𝑎
Avec 𝜔0 = √ et 𝜉 =
𝑚 2√𝑘𝑚
2. Dans un script Python, définir :
• Les différents paramètres du problème avec leurs valeurs numériques (on prendre 𝑎 = 1 𝑆. 𝐼. pour
commencer, on pourra l’ajuster par la suite).
2𝜋
• Un vecteur temps de 500 valeurs entre t0=0 et t0=10T avec 𝑇 = 𝜔
0
• La fonction f permettant de décrire l’équation différentielle (on s’inspirera de l’exemple vu en cours
et du script Python correspondant que l’on trouvera sur cahier de prépa dans le dossier documents
de cours/mécanique)
3. Toujours en s’inspirant du script vu en cours, utiliser odeint pour résoudre l’équation différentielle en
prenant comme conditions initiales 𝑧(0) = 𝑧0 et 𝑧̇ (0) = 0 .
4. Tracer la position et la vitesse en fonction du temps sur deux graphes différents pendant le régime
transitoire (on utilisera de nouveau subplot). On adaptera le vecteur temps pour observer l’ensemble du
régime transitoire.
5. Pour quelle valeur du coefficient de frottement a est-ce que l’on passe du régime pseudopériodique au
régime apériodique ? Est-ce que cela correspond à ce qu’on attend théoriquement ?