Équations différentielles
Les équations différentielles apparaissent naturellement dans de nombreux domaines au-delà des
mathématiques. Elles permettent de modéliser des phénomènes d’évolution en physique, biologie,
économie... Nous expliquons ici comment trouver des solutions approchées de ces équations grâce à
des méthodes numériques de discrétisation.
1. Méthode d’Euler
1.1. Qu’est-ce qu’une équation différentielle ?
Une équation différentielle n’est pas comme une équation classique où le but est de trouver une valeur x.
Dans une équation différentielle l’inconnue est une fonction y(x) et l’équation fait intervenir la fonction
y(x) ainsi que sa dérivée y ′ (x) (et éventuellement sa dérivée seconde y ′′ (x)).
Voici un exemple d’équation différentielle du premier ordre :
y ′ (x) = y(x).
Il s’agit donc de trouver une fonction y(x) telle que sa dérivée est égale à la fonction elle-même. Une
solution est simplement y(x) = e x . Cette solution n’est pas unique, y(x) = ke x est aussi une solution, pour
toute constante k ∈ R.
Le principe fondamental de la mécanique, qui relie l’accélération d’un objet aux forces en présence, conduit
naturellement à des équations différentielles du second ordre car l’accélération est la dérivée seconde de la
position.
Pour certaines équations différentielles il est possible de trouver une solution exacte (on renvoie à un cours
classique de mathématiques), mais ce n’est pas toujours possible. Ici nous allons calculer numériquement
des solutions approchées des équations différentielles.
1.2. Algorithme
L’idée est simple, si on connaît la valeur de y(x) et de y ′ (x) en une valeur x 0 alors on peut calculer une
valeur approchée de y(x 0 + h) par la formule du développement limité à l’ordre 1 :
y(x 0 + h) ≃ y(x 0 ) + h y ′ (x 0 )
On va ensuite calculer, de proche en proche, la valeur de y(x 0 + 2h), y(x 0 + 3h). . . Nous allons considérer
les équations différentielles du premier ordre données par une formule :
y ′ (x) = f (x, y(x))
où f : R2 → R est une fonction donnée. Cette forme d’équation permet de calculer la dérivée y ′ (x 0 ) en
utilisant la valeur de y(x 0 ) et donc par la formule précédente de calculer une valeur approchée de y(x 0 + h).
ÉQUATIONS DIFFÉRENTIELLES 2
Par exemple, pour f (x, y) = x y, l’équation différentielle est y ′ (x) = x y(x) (dont en fait les solutions sont
2
y(x) = ke x /2 (k ∈ R)). Si on considère x 0 = 2 et par exemple on sait que y(2) = 3, alors par l’équation
différentielle on sait que y ′ (2) = x 0 y ′ (x 0 ) = 6. Ainsi y(2+h) ≃ 3+6h. Par exemple y(2.1) ≃ 3+6×0.1 = 3.6.
L’algorithme est donc le suivant :
Algorithme.
Méthode d’Euler.
Entrées :
• une fonction f définissant une équation différentielle y ′ (x) = f (x, y(x)),
• un intervalle [a, b] et un nombre n de subdivisions,
• une valeur initiale y0 pour y(a).
Sortie : une liste yi des valeurs approchées de y(x i ) pour i = 0, . . . , n.
b−a
• h= n
• x0 = a
• y0 = y(a)
• Pour i variant de 0 à n − 1 :
— yi+1 = yi + h f (x i , yi )
— x i+1 = x i + h
1.3. Calculs numériques
Considérons l’équation différentielle
1
y ′ (x) =
2 y(x)
avec la condition initiale y(1) = 1.
p
La solution de ce problème est tout simplement y(x) = x. Nous allons appliquer la méthode d’Euler sur
p
l’intervalle [1, 2] afin d’obtenir une valeur approchée de y(2) = 2.
Voici ce que cela donne pour n = 10 subdivisions :
i xi yi
0 1 1
1 1.1 1.05
2 1.2 1.09
3 1.3 1.14
4 1.4 1.18
5 1.5 1.22
6 1.6 1.26
7 1.7 1.30
8 1.8 1.34
9 1.9 1.38
10 2 1.42
p p
On trouve la valeur approchée de 2 ≃ 1.42. La valeur exacte étant 2 = 1.4142 . . .
ÉQUATIONS DIFFÉRENTIELLES 3
Recommençons avec n = 100 subdivisions :
i xi yi
0 1 1
1 1.01 1.0050
2 1.02 1.0099
···
100 2 1.4148
p
On obtient donc trois décimales exactes pour 2.
1.4. Visualisation
Comparons les graphes de nos solutions approchées avec la solution exacte.
Exemple.
Soit tout d’abord l’équation différentielle avec une condition initiale :
1
y ′ (x) = et y(1) = 1.
2 y(x)
p
La solution exacte est y(x) = x.
Si on ne fixe pas de condition initiale, l’équation différentielle
1
y ′ (x) =
2 y(x)
p
a une infinité de solutions y(x) = x + c, où c est une constante. Dès que l’on fixe une condition initiale
y(x 0 ) = y0 alors la solution devient unique. Pour quelques valeurs initiales, on trace le graphe des
solutions ci-dessous.
ÉQUATIONS DIFFÉRENTIELLES 4
Exemple.
Considérons l’équation différentielle avec une condition initiale :
y ′ (x) = y(x) et y(1) = 1.
On trace le graphe de la solution exacte et des solutions approchées avec n = 10 et n = 20 subdivisions
(pour n = 100 le graphe de la solution approchée serait confondu avec le graphe de la solution exacte).
1.5. Vitesse d’un parachutiste
Déterminons la vitesse d’un parachutiste lors de son saut. Le parachutiste est soumis à deux forces : son
#» #» #»
poids P = m #»g et la force de frottement F = − f ∥ #» v
v ∥2 ∥ #»
v ∥ proportionnelle au carré de la vitesse ( f est le
coefficient de frottement). Le principe fondamental de la mécanique s’écrit :
#» #»
P + F = m #» a
Ce qui conduit à :
dv(t)
mg − f v 2 (t) = m
dt
ÉQUATIONS DIFFÉRENTIELLES 5
Ainsi l’équation différentielle est :
dv(t) f
= g − v 2 (t).
dt m
′ f
Autrement dit avec nos notations précédentes y (t) = g − m y 2 (t) où y(t) serait la vitesse du parachutiste.
Comme condition initiale on choisit v(0) = 0, c’est-à-dire qu’à l’instant initial le parachutiste démarre avec
une vitesse nulle. Il n’est pas évident de déterminer la fonction v solution de cette équation différentielle.
Cependant la méthode d’Euler nous fournit une solution approchée (dessinée ci-dessous pour g = 10,
m = 60 et f = 1).
La solution graphique met en évidence un phénomène vérifié expérimentalement : la vitesse augmente mais
est plafonnée par une vitesse limite.
Si on s’intéresse à la position x du parachutiste (c’est-à-dire sa hauteur) cela revient à trouver une primitive
dx(t)
de v, ce qui est une forme simple de résolution d’une équation différentielle : v(t) = dt . On verra en fin
de chapitre de tels exemples.
2. Méthodes de Runge-Kutta
2.1. Retour sur la méthode d’Euler
Expliquons graphiquement un pas de la méthode d’Euler. On se place en x n , on suppose connue (ou estimée)
la valeur yn de la fonction y en x n . L’équation différentielle nous donne la valeur de la dérivée y ′ (x n ) en
x n . On trace la tangente à la courbe y(x) en x n , c’est-à-dire la droite passant par le point (x n , yn ) et de
pente y ′ (x n ) = f (x n , yn ). On considère que la tangente approche bien la courbe y(x) sur un petit intervalle
[x n , x n + h], donc la valeur de y(x) sur cet intervalle est approximativement donnée par la tangente, en
particulier :
y(x n + h) ≃ yn + h y ′ (x n ).
ÉQUATIONS DIFFÉRENTIELLES 6
yn+1
tangente de pente f (x n )
′
yn
courbe
h
xn x n+1
Méthode d’Euler
La méthode d’Euler itère ce processus, autrement dit, on obtient une formule de récurrence, initialisée par
x 0 et y0 , puis :
x n+1 = x n + h
yn+1 = yn + h f (x n , yn )
2.2. Runge-Kutta d’ordre 2
Pour la méthode de Runge-Kutta d’ordre 2 (RK2), chaque pas se fait en deux étapes préalables, puis un
calcul de moyenne :
1. On calcule d’abord la pente de la tangente au point de départ k1 = f (x n , yn ). On en déduit une
approximation de la valeur de y en x n + h, c’est-à-dire yn + hk1 .
2. On calcule ensuite la pente de la tangente au point d’arrivée : k2 = f (x n + h, yn + hk1 ).
3. On prend pour valeur de y en x n + h la valeur correspondant à la moyenne des ces deux pentes :
yn+1 = yn + 2h (k1 + k2 ).
La méthode de Runge-Kutta d’ordre 2 itère ce processus. La formule combinant ces étapes est :
h
y(x n + h) ≃ yn + f (x n , yn ) + f (x n + h, yn + h f (x n , yn )) .
2
La formule de récurrence est :
k1 = f (x n , yn )
k2 = f (x n + h, yn + hk1 )
x n+1 = x n + h
yn+1 = yn + 2h (k1 + k2 )
ÉQUATIONS DIFFÉRENTIELLES 7
1 tangente au dé-
part de pente k1
2 tangente à l’ar- yn+1
rivée de pente k2
yn
k1 +k2
3 droite de pente 2
h
xn x n+1
Méthode de Runge-Kutta d’ordre 2
Pourquoi cette méthode est-elle meilleure que la méthode d’Euler ? La méthode d’Euler est basée sur la
tangente au point initial et fait partir l’approximation franchement au-dessus ou en-dessous du graphe de
la solution. Par contre avec la méthode de Runge-Kutta d’ordre 2, la direction suivie combine la pente au
départ et à l’arrivée et est beaucoup plus proche de la solution exacte. On fera une comparaison numérique
entre les méthodes dans une section suivante. Il existe une variante de la méthode de Runge-Kutta d’ordre
2 ou la pente utilisée est la pente calculée en x n + 2h .
Si on fait une analogie avec le calcul approché d’intégrale : la méthode d’Euler s’apparente à la méthode des
rectangles, la méthode RK2 avec la méthode des trapèzes, la méthode RK4 avec la méthode de Simpson.
a x
b
α+β
f 2
f (β) f (β) courbe
f (α) f (α) f (α) parabole
α β α β α α+β β
2
Méthode des rectangles Méthode des trapèzes Méthode de Simpson
ÉQUATIONS DIFFÉRENTIELLES 8
2.3. Runge-Kutta d’ordre 4
La méthode de Runge-Kutta d’ordre 4 (RK4) est encore plus précise que la méthode d’ordre 2. C’est l’une
des méthodes les plus populaires. La formule est obtenue à l’aide d’une moyenne pondérée de quatre pentes
de tangentes. Voici sa formule de récurrence x n+1 = x n + h et :
h
yn+1 = yn + k1 + 2k2 + 2k3 + k4
6
avec :
k1 = f (x n , yn )
k2 = f (x n + 2h , yn + 2h k1 )
k3 = f (x n + 2h , yn + 2h k2 )
k4 = f (x n + h, yn + hk3 )
4 tangente à l’ar-
rivée de pente k4
2 tangente au mi-
lieu de pente k2
3 tangente au mi-
lieu de pente k3 yn+1
1 tangente au dé-
yn
part de pente k1
5 droite de pente
k1 +2k2 +2k3 +k4
6
h h
2 2
xn x n+1
Méthode de Runge-Kutta d’ordre 4
2.4. Comparaisons
À chaque étape nos méthodes font des approximations. Plus le pas h est petit, plus l’erreur commise est
petite. Ces erreurs se cumulent lors des itérations. On peut donc comparer les méthodes en fonction de
l’erreur commise à chaque étape et de l’erreur totale commise à la fin de l’itération. Voici un tableau des
ordres de grandeurs des ces erreurs :
Méthode erreur locale erreur globale
Euler O(h2 ) O(h)
RK2 O(h3 ) O(h2 )
RK4 O(h5 ) O(h4 )
Ainsi, si h = 0.1 l’erreur totale commise par la méthode d’Euler est de l’ordre de 0.1, celle de la méthode RK2
est de l’ordre de 0.01 et celle de la méthode RK4 est de l’ordre de 0.0001. On voit donc que les méthodes
RK2 et RK4 sont beaucoup plus précises que la méthode d’Euler. Il faut cependant noter que la méthode
RK2 évalue la fonction f deux fois à chaque étape et est donc deux fois plus lente que la méthode d’Euler,
mais comme l’erreur est quadratique en h, le gain de précision est substantiel.
p
Voici les valeurs numériques pour le calcul de 2 obtenu comme la valeur en 2 de la solution de l’équation
différentielle y ′ = 21y avec y(1) = 1. En gras sont indiquées les décimales exactes.
ÉQUATIONS DIFFÉRENTIELLES 9
pas Euler RK2 RK4
h = 0.1 1.420519799291044 1.4142157109943851 1.414213576569407
h = 0.01 1.4148279673659128 1.4142135644527205 1.414213562374446
Sur le graphe ci-dessous on dessine les solutions exactes et approchées du problème y ′ (x) = y(x) avec
y(0) = 1 pour un pas h = 0.1. Les méthodes de Runge-Kutta sont si efficaces que l’on n’arrive pas à distinguer
le graphe de la solution exacte ( y(x) = e x ) avec les solutions approchées RK2 et RK4.
3. Systèmes d’équations différentielles
3.1. Intégration des champs de vecteurs
#»
Un champ de vecteurs est la donnée pour chaque point (x, y) ∈ R2 du plan d’un vecteur F (x, y).
#»
Ci-dessous le champ de vecteurs F (x, y) = ( y, −x).
ÉQUATIONS DIFFÉRENTIELLES 10
Une courbe intégrale est un arc γ(t) = (x(t), y(t)) tel que pour tout t ∈ R :
#»
(x ′ (t), y ′ (t)) = F (x, y).
#»
Si on note F (x, y) = ( f (x, y), g(x, y)) alors cela revient à trouver une solution du système différentiel
suivant : ′
x (t) = f (x(t), y(t))
y ′ (t) = g(x(t), y(t))
Graphiquement, en partant d’une position initiale (x 0 , y0 ), la courbe intégrale correspond à suivre les flèches
du champ de vecteurs. Ci-dessous : une portion d’une courbe intégrale (à gauche), les courbes intégrales
sont en fait des paramétrisations de cercles (à droite).
La méthode d’Euler s’étend sans difficulté à ce cas.
Méthode d’Euler.
Soit x 0 , y0 ∈ R et h > 0.
x n+1 = x n + h f (x n , yn )
yn+1 = yn + hg(x n , yn )
Voici la méthode de Runge-Kutta d’ordre 2.
Méthode RK2.
Soit x 0 , y0 ∈ R et h > 0.
= f (x n , yn )
k1
ℓ1 = g(x n , yn )
k2 = f (x n + hk1 , yn + hℓ1 )
ℓ2 = g(x n + hk1 , yn + hℓ1 )
= x n + 2h (k1 + k2 )
x
n+1
yn+1 = yn + 2h (ℓ1 + ℓ2 )
#»
Voici un autre exemple de champ de vecteurs : F (x, y) = (x − y, x y).
ÉQUATIONS DIFFÉRENTIELLES 11
3.2. Équations différentielles d’ordre 2
Considérons une équation différentielle d’ordre 2, c’est-à-dire faisant intervenir y(x), y ′ (x) et y ′′ (x) :
y ′′ (x) = f (x, y, y ′ ).
Elle se transforme en un système de deux équations différentielles d’ordre 1. On pose z(x) = y ′ (x) et on a
alors : ′
y (x) = z(x)
z ′ (x) = f (x, y, z)
La solution est généralement unique une fois que l’on a fixé les conditions initiales y(x 0 ) et y ′ (x 0 ). On peut
alors appliquer les méthodes numériques précédentes.
Le principe fondamental de la mécanique conduit naturellement à un système d’équations différentielles
d’ordre 2. En effet la vitesse est la dérivée de la position et l’accélération est la dérivée seconde de la position.
Appliquons ceci au mouvement d’un satellite orbitant dans le plan orbital de la Terre et du Soleil (supposés
tous les deux fixes !). Le satellite est soumis à la force d’attraction de la Terre et à celle du Soleil. Ces forces
s’écrivent :
#» GmM T #» #» GmMS #»
FT = − 2
eT FS = − eS
rT r T2
et le principe fondamental de la mécanique donne :
# » #»
F + F = m #»
T S a
Les données sont :
• m : masse du satellite,
• M T : masse de la Terre,
• MS : masse du Soleil,
• G : constante gravitationnelle,
• r T : distance du satellite à la Terre,
• rS : distance du satellite au Soleil,
#»
• e T : vecteur unitaire issu de la Terre dans la direction du satellite,
#»
• e S : vecteur unitaire issu du Soleil dans la direction du satellite.
ÉQUATIONS DIFFÉRENTIELLES 12
Si (x(t), y(t)) est la position du satellite, alors, le principe fondamental, projeté sur l’axe x, donne :
G M T (x(t) − x S ) G MS (x(t) − x S )
x ′′ (t) = − 23 − 3
(x(t) − x S )2 + ( y(t) − yS )2 (x(t) − x S )2 + ( y(t) − yS )2 2
On obtient une équation similaire pour y ′′ (t). On transforme ces deux équations d’ordre 2 en un système
de 4 équations d’ordre 1.
Voici des exemples d’orbites que l’on peut obtenir.
