Chapitre IV.

Résolution numérique des équations

différentielles ordinaires
Introduction: La plupart des phénomènes physiques sont modélisés par des équations
différentielles. Le modèle mathématique est complètement définit lorsque les conditions
initiales où aux limites sont spécifiées.
EXP: Oscillations amorties d'un ressort.

L'équation différentielle du mouvement:

𝑑𝑑 2 𝑥𝑥 𝑑𝑑𝑑𝑑
+ 2𝛼𝛼 + 𝜔𝜔02 𝑥𝑥 = 0
𝑑𝑑𝑡𝑡 2 𝑑𝑑𝑑𝑑
avec les conditions initiales:
𝑑𝑑𝑑𝑑
𝑡𝑡 = 0 → 𝑥𝑥 = 𝑥𝑥0 et 𝑡𝑡 = 0 → =0
𝑑𝑑𝑑𝑑

Ces équations différentielles peuvent être transformer par changement de variables, en un

système d'équations différentielles ordinaires d'ordre 1.
Méthode d'Euler:
Le problème est de résoudre l'équation différentielle:
𝑑𝑑𝑑𝑑
= 𝑓𝑓(𝑥𝑥, 𝑦𝑦)
� 𝑑𝑑𝑑𝑑
𝑥𝑥 = 𝑥𝑥0 → 𝑦𝑦 = 𝑦𝑦0
Sur un intervalle donné [x 0 , x n ].
Les étapes de calcul:
1. Discrétiser où diviser l'intervalle du domaine de calcul [x 0 , x n ] en "n" sous intervalles.
? ? ? ? ?
y0 y1 y2 y3 y n-1 yn
x0 x1 x2 x3 x n-1 xn
Le nombre de nœuds (points) égale au nombre d'intervalles plus un (n+1).
La précision des calculs augmente avec le nombre de divisions n.
2. Résoudre l'équation différentielle revient à calculer les valeurs "y i " en chaque point "x i ". Le
Développement limité en série de Taylor est utilisé dans la méthode d'Euler:

yi-1 yi-1 yi-1

x x x
xi-1 xi-1 xi-1

𝑥𝑥 𝑛𝑛 −𝑥𝑥 0
Le pas de calcul est constant: h=x 1 - x0 =x 2 -x 1 =.........=.x n -x n-1 =
𝑛𝑛
ℎ 𝑑𝑑𝑑𝑑 ℎ2 𝑑𝑑 2 𝑦𝑦 ℎ𝑘𝑘 𝑑𝑑 𝑘𝑘 𝑦𝑦
𝑦𝑦𝑖𝑖 = 𝑦𝑦𝑖𝑖−1 + � � + � 2 � + ⋯ + � 𝑘𝑘 � + 𝜃𝜃(ℎ𝑘𝑘+1 )
�����
1! 𝑑𝑑𝑑𝑑 𝑖𝑖−1 2! 𝑑𝑑𝑥𝑥 𝑖𝑖−1 𝑘𝑘! 𝑑𝑑𝑥𝑥 𝑖𝑖−1
𝐸𝐸𝐸𝐸𝐸𝐸𝐸𝐸𝐸𝐸𝐸𝐸 𝑑𝑑 ′ 𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜 𝑘𝑘+1

Dans la méthode d'Euler, la série de Taylor est tronquée à l'ordre 1 (l'erreur est d'ordre 2):
𝑑𝑑𝑑𝑑
𝑦𝑦𝑖𝑖 = 𝑦𝑦𝑖𝑖−1 + ℎ � �
𝑑𝑑𝑑𝑑 𝑖𝑖−1
𝑑𝑑𝑑𝑑
Or l'équation différentielle à résoudre est: = 𝑓𝑓(𝑥𝑥, 𝑦𝑦)
𝑑𝑑𝑑𝑑

Donc: 𝒚𝒚𝒊𝒊 = 𝒚𝒚𝒊𝒊−𝟏𝟏 + 𝒉𝒉𝒉𝒉(𝒙𝒙𝒊𝒊−𝟏𝟏 , 𝒚𝒚𝒊𝒊−𝟏𝟏 ) avec: i=1,n

Donc à partir de l'information disponible au début du domaine de calcul (condition initiale
où au limite), la solution est avancée par la formule d'Euler, pour calculer toutes les valeurs
"y i " en chaque nœud "x i " du domaine de calcul.
Exemple (1): (calcul manuel)
Résoudre l'équation différentielle suivante par la méthode d'Euler:
𝑑𝑑𝑑𝑑
= −2𝑥𝑥𝑥𝑥
�𝑑𝑑𝑑𝑑
𝑥𝑥 = 0 → 𝑦𝑦 = 1
Le domaine de calcul est [0 , 0.05] et le pas de calcul est ∆x=0.1
Solution:
?
y 0 =1 y1 y2? y3? y4? y5?

x 0 =0 x 1 =0.1 x 2 =0.2 x 3 =0.3 x 4 =0.4 x 5 =0.5

𝑦𝑦𝑖𝑖 = 𝑦𝑦𝑖𝑖−1 + ∆𝑥𝑥(−2𝑥𝑥𝑖𝑖−1 𝑦𝑦𝑖𝑖−1 ) avec: i=1,5
i=1: x 1 =x 0 +0.1000=0.1000 y 1 =y 0 +∆x* (-2*x 0 *y 0 )=1+0.1000*(-2*0.0000*1.000) y 1 =1.0000
i=2: x 2 =x 1 +0.1000=0.2000 y 2 =y 1 +∆x* (-2*x 1 *y 1 )=1+0.1000*(-2*0.1000*1.000) y 2 =0.9800
i=3: x 3 =x 2 +0.1000=0.3000 y 3 =y 2 +∆x* (-2*x 2 *y 2 )=0.9800+0.1000*(-2*0.2000*0.9800) y 3 =0.9408
i=4: x 4 =x 3 +0.1000=0.4000 y 4 =y 3 +∆x* (-2*x 3 *y 3 )=0.9408+0.1000*(-2*0.3000*0.9408) y 4 =0.8844
i=5: x 5 =x 4 +0.1000=0.5000 y 5 =y 4 +∆x* (-2*x 4 *y 4 )=0.8844+0.1000*(-2*0.4000*0.8844) y 5 =0.8140
Comparaison des résultats obtenus par la méthode d'Euler avec la solution exacte:
𝑑𝑑𝑑𝑑 𝑑𝑑𝑑𝑑
= −2𝑥𝑥𝑥𝑥 ⟹ = −2𝑥𝑥𝑥𝑥𝑥𝑥
𝑑𝑑𝑑𝑑 𝑦𝑦
2
en intégrant cette équation on obtient la solution exacte: 𝑦𝑦 = 𝐶𝐶 𝑒𝑒 −𝑥𝑥
2
avec la condition au limite imposée, on aboutit à la solution suivante: 𝑦𝑦 = 𝑒𝑒 −𝑥𝑥

x y Euler Y exacte
0.0000 1.0000 1.0000
0.1000 1.0000 0.9900
0.2000 0.9800 0.9608
0.3000 0.9408 0.9139
0.4000 0.8844 0.8521
0.5000 0.8140 0.7788

Pour mieux comparer les résultats, on trace les deux courbes sur le même graphe (Origin):

On peut améliorer les résultats obtenus par la méthode d'Euler en diminuant le pas de Calcul
∆x. Dans ce cas la courbe de Euler devient plus proche à la courbe des résultats exactes.
Exemple (2): Décharge d'un condensateur -Voir cours Phys 2-

𝑑𝑑𝑑𝑑 𝑑𝑑𝑑𝑑
𝑈𝑈 = 𝑅𝑅𝑅𝑅 = − = −𝑅𝑅𝑅𝑅
𝑑𝑑𝑑𝑑 𝑑𝑑𝑑𝑑
On cherche à déterminer la variation de la différence de potentiel U en fonction du temps t
𝑑𝑑𝑑𝑑 1
=− 𝑈𝑈
𝑑𝑑𝑑𝑑 𝑅𝑅𝑅𝑅
avec la condition initiale: t =0 U=5 V

Les données: RC=1 , T max =1.5 et ∆t=0.05

Indications:
• Commencer par calculer le nombre de nœuds (points) du domaine de calcul.
• Utiliser Excel où écrire le programme Fortran pour calculer U dans tous les nœuds:
1
𝑈𝑈𝑖𝑖 = 𝑈𝑈𝑖𝑖−1 + Δ𝑡𝑡 �− 𝑈𝑈𝑖𝑖−1 �
𝑅𝑅𝑅𝑅
• On trouve que le condensateur se décharge dans la résistance avec la variation suivante:
Exemple (3): Chute d'un corps avec résistance du milieu - Voir cours Phys 1-

∑ �F⃗ = m γ
�⃗ y'
��⃗
�P⃗ + �F⃗ = m dV �F⃗
dt

En projetant sur l’axe yy’ :

dV
mg − KV 2 = m ...........(I) �⃗
P
dt

Condition initiale: à t=0, V=0 y

K
Données du problème: g=9.81 m/s2 ; = 1 m−1 , ∆t=0.05 ; domaine de calcul [0 , 1.5]
m

Solution: (Voir: Exemples du cours traités par Excel)

On écrit l'équation (I) sous la forme suivante:
dV K
= g − V2
� dt m

t=0→V=0

𝑡𝑡 𝑚𝑚𝑚𝑚𝑚𝑚 −𝑡𝑡 0 𝑡𝑡 𝑚𝑚𝑚𝑚𝑚𝑚 −𝑡𝑡 0 1.5−0

Calcul du nombre de points: Δ𝑡𝑡 = 𝑛𝑛 −1
⇒ 𝑛𝑛 =
Δ𝑡𝑡
+1, 𝑛𝑛 =
0.05
+ 1 = 31

V0 V1 Vi-1 Vi

t0 t1 ti-1 ti

𝐾𝐾 2
𝑡𝑡𝑖𝑖 = 𝑡𝑡𝑖𝑖−1 + ∆𝑡𝑡 et 𝑉𝑉𝑖𝑖 = 𝑉𝑉𝑖𝑖−1 + ∆𝑇𝑇 �𝑔𝑔 −
𝑚𝑚
𝑉𝑉𝑖𝑖−1 � avec: i=1,n

La courbe représentative de la variation de la vitesse en fonction du temps (Origin):

3,5
v=f(t)
3,0

2,5
vitesse(m/s)

2,0

1,5

1,0

0,5

0,0
0,0 0,2 0,4 0,6 0,8 1,0 1,2 1,4 1,6
temps(s)
 Résolution d'une équation différentielle d'ordre 2 par la méthode d'Euler

Il suffit de transformer l'équation différentielle d'ordre 2, en un système de deux équations

d'ordre 1 . La transformation est effectuée par changement de variable:

𝑑𝑑 2 𝑦𝑦 𝑑𝑑𝑑𝑑
𝛼𝛼 2 + 𝛽𝛽 + 𝛾𝛾 = 0
𝑑𝑑𝑥𝑥 𝑑𝑑𝑑𝑑

𝑑𝑑𝑑𝑑
On pose: = 𝑍𝑍
𝑑𝑑𝑑𝑑

𝑑𝑑𝑑𝑑 𝑑𝑑𝑑𝑑 𝛽𝛽 𝛾𝛾
donc: 𝛼𝛼 + 𝛽𝛽𝛽𝛽 + 𝛾𝛾 = 0 d'ou = − 𝑍𝑍 −
𝑑𝑑𝑑𝑑 𝑑𝑑𝑑𝑑 𝛼𝛼 𝛼𝛼

On obtient un système d'équations résolubles par la méthode d'Euler:

𝑑𝑑𝑑𝑑 𝛽𝛽 𝛾𝛾
= − 𝑍𝑍 −
�𝑑𝑑𝑑𝑑 𝛼𝛼 𝛼𝛼
𝑑𝑑𝑑𝑑
= 𝑍𝑍
𝑑𝑑𝑑𝑑

Exemple: Pendule amortie - Voir cours Phys3-

La variation de θ en fonction de t, est donnée par l'équation différentielle suivante:

𝑑𝑑 2 𝜃𝜃 𝐶𝐶𝑓𝑓 𝑑𝑑𝑑𝑑 𝑔𝑔
+ + sin(𝜃𝜃 ) = 0
𝑑𝑑𝑡𝑡 2 𝑚𝑚 𝑑𝑑𝑑𝑑 𝑙𝑙
θ
�𝑡𝑡 = 0 → 𝜃𝜃 = 0.1
𝑑𝑑𝑑𝑑
𝑡𝑡 = 0 → = 0
𝑑𝑑𝑑𝑑

Les données du problème: g=9.81 m/s2 , m=0.1 kg , l=50 cm , C f = (0 ; 0.01 ; 0.05 ; 0.1) ,
∆t=0.1 s et l’intervalle de résolution [0 , 20] s

Solution: (Voir: Exemples du cours traités par Excel)

Les résultats obtenus sont présentés sur les graphes suivants:
 Cf=0.0
0,10 0,10 Cf=0.01

0,05 0,05

angle(rd)
angle(rd)

0,00 0,00

-0,05 -0,05

-0,10 -0,10

0 5 10 15 20 0 5 10 15 20
temps(s) temps(s)

0,12
0,10 Cf=0.05 0,10
C=0.1
f
0,08
0,06
0,05
0,04
angle(rd)

angle(rd)

0,02
0,00 0,00
-0,02
-0,05 -0,04
-0,06
-0,08
-0,10
0 5 10 15 20 0 5 10 15 20
temps(s) temps(s)