1. No category

Chapitre 9 : Méthode d’Euler
MPSI / Informatique commune
Chapitre 9 : Résolution numérique d’une équation
différentielle par la méthode d’Euler
I. Introduction
• De nombreux phénomènes physiques se modélisent naturellement à l’aide d’équations différentielles pour
lesquelles on ne dispose pas de solutions analytiques : un pendule amorti nous amène naturellement à étudier l’équation θ¨ = −k1 sin(θ) − k2 θ˙ , les problèmes de cinétique chimique conduisent naturellement à des
systèmes différentiels non linéaires décrivant les évolutions de différents réactifs au cours du temps, et en
mécanique des fluides, une partie des phénomènes sont décrits par les équations de type Navier-Stokes.
• En mathématiques, ces équations ont leur intérêt propre, et étudier le comportement qualitatif de solutions
est nettement plus aisé si on peut visualiser une approximation raisonnable de celles-ci.
II. Principe de la méthode d’Euler
On considère sur l’intervalle I = [a; b] l’équation différentielle
′
y = F(t, y)
(E) :
y(t0 ) = y0
Exemple 1
y ′ + a(t)y = b(t) ⇐⇒ y ′ = −a(t)y − b(t) = F(t, y)
Exemple 2
2ty
= F(t, y)
t 2 + y2
Lorsqu’il n’est pas possible de résoudre théoriquement (E), c’est-à-dire de trouver une expression explicite de
la solution f de (E) qui vérifie f (t0 ) = y0 , on met en oeuvre un tracé approché de la courbe intégrale sur l’intervalle
I = [a; b].
Pour h suffisamment petit, on peut assimiler la courbe de f à sa tangente avec l’approximation suivante :
(t 2 + y2 )y ′ + 2ty = 0 ⇐⇒ y ′ = −
f (t + h) ≃ f (t) + h f ′ (t)
soit en remplaçant :
f (t + h) ≃ f (t) + hF(t, y(t))
Pour effectuer cette approximation, on va donc dans une première étape découper l’intervalle d’étude [a; b] en
b−a
n sous-intervalles de même longueur h =
: on obtient alors une suite finie de points (t0 ,t1 , . . . ,tn ) qui vérifie :
n
t0 = a,tn = b et pour tout i ∈[[ 0, n − 1 ]], ti+1 = ti + h
1
Chapitre 9 : Méthode d’Euler
MPSI / Informatique commune
Cette suite est appelée une subdivision régulière de [a; b] de pas h.
Remarque
(ti ) est donc une suite géométrique de raison h, on en déduit que pour tout i ∈[[ 0, n ]], ti = t0 + ih
L’objectif est maintenant d’approcher la courbe intégrale sur chaque sous-intervalle [ti ;ti+1 ] par sa tangente en ti .
III. Construction de l’algorithme et codage
1. Algorithme
En chaque point de la subdivision, on approchera donc f (ti ) par yi défini par :
Pour tout i ∈[[ 0, n − 1 ]], yi+1 = yi + h f ′ (ti ) soit
yi+1 = yi + hF(ti , yi ) et y0 = y(t0 ) = a
On approche alors f sur [a; b] par une fonction affine par morceaux et la courbe intégrale de f par une ligne
polygonale.
L’algorithme s’effectue pas à pas puisque la position d’un point au temps ti nécessite la connaissance de la
position au temps ti−1 .
On se donne donc une condition initiale : le point est en y0 à l’instant t0 .
y1 = y0 + hF(t0 , y0 )
y2 = y1 + hF(t1 , y)
..
.
yn = yn−1 + hF(tn−1 , yn−1 )
2. Codage en python
• La fonction euler va prendre en entrée une fonction F, les bornes a et b de l’intervalle d’étude, la condition initiale y0 , et le nombre n de sous-intervalles définies par la subdivision de I. Plus précisément :
avec ces données, la fonction va déterminer les approximations de la solution de y ′ (t) = F(t, y(t)) avec
la condition initiale y(a) = y0 , en rendant un tableau de valeurs approchées de la solution pour chaque
temps a + ih majorés par b. Ce tableau est construit à l’aide d’une liste à laquelle on adjoint les nouveaux termes calculés.
• Codage en python
2
Chapitre 9 : Méthode d’Euler
MPSI / Informatique commune
IV. Erreur de consistance, erreur absolue
1. Définitions
• Au temps ti+1 , on définit l’erreur de consistance par l’écart entre la valeur exacte f (ti+1 ) et la valeur
approchée yi+1 .
Donc l’erreur de consistance ei est définie par
• L’erreur globale est définie par
ei = f (ti+1 ) − yi+1
θn = max |ei |
06i6n−1
pour tout i ∈[[ 0, n − 1 ]].
résultant du calcul des valeurs successives y1 , . . . , yn .
Les fonctions f , f1 , f2 et f3 représentent ici les solutions exactes passant par les points (t0 , y0 ), (t1 , y1 ), (t2 , y2 )
et (t3 , y3 ).
Remarque
On pourrait montrer que l’erreur globale diminue lorsque le pas h diminue
.
2. Compléments : consistance et stabilité d’une méthode numérique
• Définition (méthode consistante)
On dit que la méthode est consistante si pour toute solution exacte f , la somme des erreurs de
n−1
consistance relatives à f , soit ∑ |ei |, tend vers 0 quand h tend vers 0.
i=0
3
Chapitre 9 : Méthode d’Euler
MPSI / Informatique commune
• Une autre notion fondamentale est la notion de stabilité. Dans la pratique, le calcul récurrent des points
yi est en effet entaché d’erreurs d’arrondi εi . Pour que les calculs soient significatifs, il est indispensable
que la propagation de ces erreurs reste contrôlable.
Définition (méthode stable)
On dit que la méthode est stable s’il existe une constante S > 0, appelée constante de stabilité,
telle que pour toutes suites (yi )06i6n et (y˜i )06i6n définies par
∀i ∈[[ 0, n − 1 ]], yi+1 = yi + hF(ti , yi )
y˜i+1 = y˜i + hF(ti , y˜i )
et
on ait
n−1
!
max |y˜i − yi | 6 S |y˜0 − y0 | + ∑ |ei |
06i6n
i=0
Remarque
Autrement dit, une petite erreur initiale ε0 = |y˜0 − y0 | et de petites erreurs d’arrondi εi = |y˜i − yi |
dans le calcul récurrent des y˜i provoquent une erreur finale max |y˜i − yi | contrôlable.
06i6n
V. Choix du pas
• Le choix du pas est une étape obligée lors de la mise en place d’une méthode numérique de résolution. Si
on choisit ce pas trop petit, le temps de calcul sera élevé. Si au contraire h est trop grand, l’erreur de consistance sera trop importante. Il s’agit donc comme toujours en calcul scientifique de faire un compromis, et
ce n’est pas toujours simple.
• Travaillons sur un exemple et déterminons les temps de calculs nécessaires pour résoudre numériquement
y ′ = y sur [0; 1] avec la condition initiale y(0) = 1, pour différents pas h = 1n .
On évalue par ailleurs l’erreur globale θn = max | f (ti ) − yi |. Sachant que la solution de l’équation diffé06i6n−1
rentielle est la fonction exponentielle qui est convexe, les valeurs approchées s’éloignent irrémédiablement
des valeur réelles et on a θn = |e − yn |.
Les temps sont exprimés en secondes et on donne deux chiffres significatifs :
n
103
104
105
106
107
Temps
6, 7.10−4
6, 9.10−3
6, 8.10−2
6, 7.10−1
6, 7
−3
−4
−5
Erreur globale 1, 36.10
1, 36.10
1, 36.10
1, 36.10−6 1, 36.10−7
• Les résultats sont sans surprise : d’une part, le temps est proportionnel au nombre de termes calculés, et
i
d’autre part, on a ici yi = 1 + 1n , donc :
1 n
e
θn = e − 1 +
∼
n n→+∞ 2n
Avec la méthode d’Euler, on ne peut obtenir une précision de l’ordre de 10−6 qu’avec un temps de calcul
assez important.
4
Chapitre 9 : Méthode d’Euler
MPSI / Informatique commune
VI. Méthode d’Euler « vectorialisée »
Un grand nombre de problèmes dans tous les domaines conduisent à des systèmes d’équations différentielles
ou à des équations différentielles d’ordre supérieur. L’idée est d’adapter la méthode précédente afin de résoudre
numériquement ce type de problèmes.
1. Système différentiel
Considérons un système différentiel de la forme
 ′
x (t) = f1 (t, x1 , x2 )


 1′
x2 (t) = f2 (t, x1 , x2 )
x (t ) = y1


 1 0
x2 (t0 ) = y2
On pose Y (t) le vecteur de coordonnées (x1 (t), x2 (t)). Le système est alors équivalent au problème de
Cauchy suivant
′
Y (t) = F(t,Y )
Y (t0 ) = Y0
En posant Y0 = (x1 (t0 ), x2 (t0 )) et F(t,Y (t)) = F(t, x1 (t), x2 (t)) = ( f1 (t, x1 , x2 ), f2 (t, x1 , x2 ))
2. Équations différentielles du second ordre
Considérons une équation différentielle du second ordre de la forme
x ′′ + ax ′ + bx = c avec
x(t0 ) = y0 , x ′ (t0 ) = y′0
L’objectif est de se ramener à un système différentiel pour appliquer la méthode précédente. Pour cela
on pose (x1 (t), x2 (t)) = (x(t), x ′ (t)). L’équation différentielle est alors équivalente au système différentiel
suivant :
 ′
x (t) = x2 (t)


 1′
x2 (t) = −ax2 (t) − bx1 (t) + c
x (t ) = y0


 1 0
x2 (t0 ) = y′0
On pose Y (t) le vecteur de coordonnées (x1 (t), x2 (t)). Le système est alors équivalent au problème de
Cauchy suivant
′
Y (t) = F(t,Y )
Y (t0 ) = Y0
En posant Y0 = (x1 (t0 ), x2 (t0 )) et F(t,Y (t)) = F(t, x1 (t), x2 (t)) = (x2 (t), −ax2 (t) − bx1 (t) + c)
5