Méthodes numériques
Utilisation de odeint our les équations
différentielles du 2nd ordre
Objectif : résoudre une équation différentielle d’ordre 2 à l’aide de la fonction odeint
de Python.
Précédemment, nous avons vu la méthode d’Euler explicite pour résoudre une équation
différentielle d’ordre 1. Cette méthode n’est cependant pas toujours fiable. Nous allons
ici utiliser la fonction odeint, de la librairie [Link]. Elle permet de résoudre
des équations différentielles de manière fiable et rapide.
De plus, nous nous intéressons à une équation d’ordre 2, comme celle du pendule simple :
θ̈ + ω02 sin θ = 0, soit aussi θ̈ = −ω02 sin θ. (1)
a/ Reformulation de l’équation en une équation d’ordre 1
La fonction odeint a une syntaxe particulière, qui oblige
à transformer l’équation d’ordre
θ
2 sur θ, en une équation d’ordre 1 sur le couple y = .
θ̇
On parle parfois de vecteur à deux composantes pour désigner y :
– sa première composante est y0 = θ,
– sa seconde composante est y1 = θ̇.
1 Dériver le couple y pour trouver l’équation différentielle qu’il suit.
dy d θ
=
dt dt θ̇
θ̇
=
θ̈
θ̇
=
−ω02 sin θ
y1
=
−ω02 sin y0
Méthodes numériques | odeint 1/7 Raoul Follereau | PTSI
On arrive donc à une équation du premier ordre sur y :
dy
y1
= F (y) avec F (y) = . (2)
dt −ω02 sin y0
Bilan : on
estpassé de l’équation d’ordre 2 sur θ à l’équation ci-dessus, d’ordre 1 sur le
θ
couple y = .
θ̇
Autres exemples pour bien comprendre
(à faire sur feuille à part).
2 On considère l’équation du pendule, mais avec des frottements :
θ̈ = −ω02 sin θ − αθ̇.
Avec la même démarche que ci-dessus, transformer cette équation d’ordre 2 sur θ en une
dy
θ
équation d’ordre 1 sur le couple y = (donc en une équation du type = F (y)).
θ̇ dt
dy d θ
=
dt dt θ̇
θ̇
=
θ̈
θ̇
=
−ω 2 sin θ − αθ̇
0
y1
=
−ω02 sin y0 − α y1
dy
y1
Donc on a = F (y) avec F (y) = .
dt −ω02 sin y0 − α y1
3On considère l’équation suivie par la tension u aux bornes du condensateur dans un
circuit RLC en régime sinusoïdal forcé :
ω0
ü + u̇ + ω02 u = α cos ωt.
Q
Avec la même démarche que ci-dessus, transformer cette équation d’ordre 2 sur u en une
dy
u
équation d’ordre 1 sur le couple y = (donc en une équation du type = F (y,t)).
u̇ dt
Méthodes numériques | odeint 2/7 Raoul Follereau | PTSI
dy d u
=
dt dt u̇
u̇
=
ü
u̇
=
− ωQ0 u̇ − ω02 u + α cos ωt
y1
=
− ωQ0 y1 − ω02 y0 + α cos ωt
dy
y1
Donc on a = F (y,t) avec F (y,t) = (cette fois il y a
dt − ωQ0 y1 − ω02 y0 + α cos ωt.
dépendance explicite en t)
Méthodes numériques | odeint 3/7 Raoul Follereau | PTSI
b/ Utilisation de odeint
On repart sur l’exemple du pendule (équations 1 et 2). on utilise la fonction odeint de
la librairie [Link] dont voici la syntaxe :
import scipy . integrate as sp
y_sol = sp . odeint (F , y_ini , t )
Paramètres de odeint :
I F est une fonction du type F(y,t), qui renvoie la valeur de la dérivée de y à l’instant t.
Pour résoudre une équation d’ordre 1, y est un scalaire.
y
Mais si l’équation est d’ordre 2, alors y est de type vecteur : y = 0 .
y1
Attention, même si F ne dépend pas de t, il faut tout de même que t apparaisse
comme second argument : F(y,t).
I y_ini : valeur initiale de y. C’est un couple qui a autant de composantes que y.
I t est un tableau qui contient les instants auxquels la solution sera calculée.
Par exemple, pour le problème du pendule, on peut définir F, y_ini et t de la façon
suivante (avant l’appel à odeint) :
# Fonction F telle que d / dt ( theta , theta_point ) = F ( theta , theta_point )
def F (y , t ):
return ( y [1] , - w0 **2 * np . sin ( y [0]))
t = np . linspace (0 ,5* T0 ,1000) # tableau des temps , de 0 à 5 T0 , 1000 points .
theta0 = 50 * np . pi /180. # angle initial ( rad ) , ici correspond à 50 degr é s .
theta_prime0 = 0 # vitesse angulaire initiale ( rad / s )
y_ini = ( theta0 , theta_prime0 ) # valeur initiale de y = ( theta , theta_point )
C’est seulement une fois F, y_ini et t définis comme ci-dessus qu’on peut exécuter
odeint :
y_sol = sp . odeint (F , y_ini , t )
Objets retournés par odeint :
I Après exécution de y_sol = [Link](F, y_ini, t), y_sol est une matrice dont
la première colonne contient les valeurs de y0 (t) (donc dans notre exemple, de θ(t)),
et la seconde colonne contient les valeurs de y1 (t) (donc dans notre exemple, de θ̇(t)).
On peut y accéder avec la syntaxe suivante :
theta = y_sol [: ,0] # 1 è re colonne de y_sol
theta_prime = y_sol [: ,1] # 2 è me colonne de y_sol
Il suffit ensuite de tracer theta en fonction de t avec les méthodes usuelles :
[Link](t,theta, label="solution de l’équation non linéaire avec odeint"),
etc.
Méthodes numériques | odeint 4/7 Raoul Follereau | PTSI
c/ Exemple d’utilisation dans le cas du pendule
1 - Reprendre tous les éléments précédents pour écrire un algorithme qui permet de
résoudre l’équation non linéaire du pendule (équation 1, que l’on a transformée en
l’équation 2). On tracera la solution.
Des éléments de code sont déjà présents dans le fichier Capytale numéro dac9-1273229.
2 - Écrire une fonction solution_petits_angles(t,theta0) qui prend en argument un
instant t et un angle initial theta0, et qui retourne la valeur theta0 * [Link](w0*t)
(où w0 est la pulsation, déjà définie par ailleurs).
3 - Ajouter, sur le même graphe que celui de la question 1, la solution pour les petits
angles. Faire pour cela une nouvelle ligne avec [Link].
4 - Explorer avec différents angles initiaux. Tester par exemple 5°, 45° et 140°.
À partir de quel angle voit-on clairement un désaccord entre la solution sans ap-
proximation (équation avec sin θ) et la solution avec petits angles (approximation
sin θ ' θ) ?
5 - À l’aide de votre simulation numérique, mesurer la période T des oscillations pour
des angles initiaux de 10°, 30°, 60° et 80°. Calculer le rapport T /T0 dans chaque cas
(où T0 = 2π/ω0 est la période pour les petits angles) et consigner les résultats dans
un tableau sur votre feuille.
6 - On souhaite vérifier ceci expérimentalement. À l’aide du pendule à votre disposi-
tion.
Faire d’abord la mesure de T0 en prenant un angle petit (10° par exemple).
Faire ensuite des mesures à 30°, 60° et 80°, calculer le rapport T /T0 , et comparer avec
les résultats numériques : y a-t-il accord ?
Vous pouvez essayer pour des angles plus grands, mais vous verrez que le capteur
pose problème et qu’il faut ruser.
7 - Si le temps le permet, caractériser expérimentalement le pendulep: la période aux
petits angles est-elle en accord avec la formule théorique T = 2π l/g ? Comment
mesurer le facteur de qualité ?
θ0 = 5o
1.5
θ0 = 45o
θ0 = 140o
1.0
0.5 Exemple de solutions numériques obtenues avec
θ/ θmax
odeint, sans approximation, pour trois condi-
0.0
tions initiales différentes. Cf question 1 ci-
dessus.
−0.5
−1.0
0 1 2 3 4
t
Méthodes numériques | odeint 5/7 Raoul Follereau | PTSI
Exemples de résultats :
θ0 T /T0 (simulation) T /T0 (exp) T (exp)
10° 1,004 1,003 1,366 s
30° 1,024 1,018 1,386 s
60° 1,080 1,079 1,470 s
90° 1,180 1,184 1,613 s
Pour l’expérience, on a mesuré T0 = 1,362 s pour des angles petits.
Code complet :
# ------ Import des bibliotheques
import matplotlib . pyplot as plt # pour les graphiques
import numpy as np # pour les tableaux
import scipy . integrate as sp # pour les outils numeriques comme odeint
# ------ Parametres physiques -- on ne changera pas ces valeurs
g = 9.81 # pesanteur , m / s **2
l = 1.00 # longeur du pendule , m
w0 = ( g / l )**0.5 # calcul de la pulsation , rad / s
T0 = 2.* np . pi / w0 # calcul de la periode , s , ici on obtient T0 =1 s
print ( " Periode du pendule ( si petits angles ) : T0 = " + str ( T0 ) + " s " )
# ------
# ------ Solution lorsqu ' on fait l ' approximation des petits angles ( oscillateur harmonique
# t : instant t ou tableau d ' instants
# theta0 : amplitude maximale
# Il est sous - entendu que la vitesse initiale est nulle .
def s o l u t i o n _ p e t i t s _ a n g l e s (t , theta0 ):
return theta0 * np . cos ( w0 * t )
# ------ Solution de l 'é quation g é n é rale , sans approximation , à l ' aide d ' odeint
# D é finition de la fonction f telle que d / dt ( theta , theta_point ) = F ( theta , theta_point )
def F (y , t ):
return ( y [1] , - w0 **2* np . sin ( y [0]))
t = np . linspace (0 ,5* T0 ,1000) # tableau des temps , va de 0 à 5 xla p é riode , avec 1000 points
theta0 = 50 * np . pi /180. # angle initial ( rad )
theta_prime0 = 0 # vitesse angulaire initiale ( rad / s )
y_ini = ( theta0 , theta_prime0 ) # valeur initiale de y = ( theta , theta_point )
y_sol = sp . odeint (F , y_ini , t )
theta = y_sol [: ,0]
theta_prime = y_sol [: ,1]
# ------
# ------ Trac é s
plt . figure ( " figure " )
plt . plot (t , s o l u t i o n _ p e t i t s _ a n g l e s (t , theta0 ) , label = " solution cas des petits angles " )
plt . plot (t , theta , label = " solution de l 'é quation non lin é aire avec odeint " )
plt . xlabel ( " t ( s ) " )
plt . ylabel ( " theta ( degr é s ) " )
plt . legend ()
plt . grid ()
plt . show ()
# ------
Méthodes numériques | odeint 6/7 Raoul Follereau | PTSI
# ------ Pour trouver automatiquement la p é riode de la solution
L =[]
for k in range ( len ( t ) -1):
if theta [ k ]* theta [ k +1] <0:
L . append ( t [ k ])
print ( L [2] - L [0])
# ------
Méthodes numériques | odeint 7/7 Raoul Follereau | PTSI