C HAPITRE 4

R ÉSOLUTION NUMÉRIQUE DES

ÉQUATIONS DIFFÉRENTIELLES
ORDINAIRES

Matière : Analyse Numérique

Sommaire
4.1 Équations différentielles ordinaires (EDO). . . . . . . . . . . . . . . . 38
4.2 Méthodes (ou schémas) numériques à un pas. . . . . . . . . . . . . . . 39
4.2.1 Maillage du domaine. . . . . . . . . . . . . . . . . . . . . . . . . 39
4.2.2 Méthode d’Euler. . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
4.2.3 Méthode de Crank-Nicolson. . . . . . . . . . . . . . . . . . . . . 42
4.2.4 Méthode de Runge-Kutta. . . . . . . . . . . . . . . . . . . . . . . 42
4.3 Étude générale des méthodes numériques à un pas. . . . . . . . . . . 43
4.3.1 Convergence d’une méthode numérique. . . . . . . . . . . . . . 43
4.3.2 Stabilité absolue d’une méthode numérique. . . . . . . . . . . . 46

Les équations différentielles décrivent l’évolution de nombreux phénomènes dans

des domaines variés. Elles constituent l’outil le plus fréquemment utilisé dans la mo-
délisation des problèmes des sciences physiques et ceux de l’ingénieur. Une équation
différentielle est une équation impliquant une ou plusieurs dérivées d’une fonction in-
connue. Si toutes les dérivées sont prises par rapport à une seule variable, on parle
d’équation différentielle ordinaire (EDO). Une équation mettant en jeu des dérivées
partielles est appelée équation aux dérivées partielles (EDP).

37
38 Chapitre 4. Résolution numérique des équations différentielles ordinaires

On dit qu’une équation différentielle (ordinaire ou aux dérivées partielles) est d’ordre
k si elle implique des dérivées d’ordre au plus k.

Dans ce chapitre, on considère des équations différentielles ordinaires (EDO) d’ordre

un. On présente quelques méthodes (ou schémas) numériques pour la résolution des
EDO munies des conditions initiales.

4.1 Équations différentielles ordinaires (EDO).

Une équation différentielle ordinaire admet généralement une infinité de solutions.
Pour en sélectionner une, on doit imposer une condition supplémentaire qui corres-
pond à la valeur prise par la solution en un point de l’intervalle d’intégration.
Par conséquent, toute équation différentielle ordinaire (EDO) munie d’une condition
initiale construit un problème appelé problème de Cauchy.

Définition 4.1 (Problème de Cauchy)

Soit F : I × R → R une fonction donnée et y 0 la dérivée de y par rapport à t.
On appelle problème de Cauchy, le problème :

 Trouver y : I ⊂
R → R tel que
y 0 (t) = F t, y(t) , ∀t ∈ I (4.1)
y(t0 ) = y0


avec t0 un point de I et y0 une valeur appelée donnée initiale. L’intervalle I s’appelle intervalle
de vie de la solution.

L’étude mathématique de l’existence et de l’unicité de solution d’un problème de Cau-

chy peut être une affaire délicate. Dans ce chapitre, nous nous contentons de rappeler
un résultat d’existence et d’unicité globale, au sens où on peut intégrer le problème de
Cauchy jusqu’à t = ∞ :

Théorème 4.2 (Théorème de Cauchy-Lipschitz)

On suppose que la fonction F : I × R → R est :
1. continue par rapport à ses deux variables,
2. lipschitzienne par rapport à sa deuxième variable :
c-à-d. qu’il existe une constante positive L (appelée constante de Lipschitz) telle que
|F (t, y1 ) − F (t, y2 )| ≤ L |y1 − y2 |, ∀t ∈ I, ∀y1 , y2 ∈ R.
Alors la solution y = y(t) du problème de Cauchy (4.1) existe, est unique et appartient à C 1 (I).
4.2. Méthodes (ou schémas) numériques à un pas. 39

En pratique, on ne peut expliciter les solutions que pour des équations différentielles
ordinaires très particulières. Dans certains cas, on ne peut exprimer la solution que
sous forme implicite. Dans d’autres cas, on ne parvient même pas à représenter la
solution sous forme implicite. Pour ces raisons, on cherche des méthodes numériques
capables d’approcher la solution de toutes les équations différentielles ordinaires qui
admettent une solution.

4.2 Méthodes (ou schémas) numériques à un pas.

On considère le problème de Cauchy et on suppose que l’on ait montré l’existence
d’une solution y. On s’intéresse à la détermination d’une approximation de la solution y
du problème de Cauchy.

4.2.1 Maillage du domaine.

Avant d’approcher numériquement la solution de l’EDO étudiée, une approximation
du domaine est requise en premier lieu. Le domaine sera représenté par un maillage,
qui est un ensemble de points isolés.
En dimension 1, on appelle maillage : toute subdivision du domaine en segments.
Pour un intervalle I = [t0 , T ], avec t0 , T < +∞, la seule subdivision est une décompo-
sition en segments successifs :
On subdivise l’intervalle I en N sous-intervalles :
N
[ −1
I= [tn , tn+1 ]
n=0

avec une collection de points tn (N + 1 points) tel que

t0 < t1 < . . . < tn < tn+1 < . . . < tN = T

 Chaque point du maillage définit un noeud.

 La distance entre deux noeuds successifs définit un pas de discrétisation (∆t ou h) :

∆t = tn+1 − tn

Ce pas de discrétisation pourrait être constant dans le cas d’un maillage équidistant,
on parlera dans ce cas d’un maillage régulier, ou d’un maillage irrégulier dans le cas
où le pas est variable :
40 Chapitre 4. Résolution numérique des équations différentielles ordinaires

Définition 4.3 (Maillage 1D)

• Le maillage est dit uniforme (ou régulier), si les noeuds tn sont équidistants, c-à-d :
T − t0
∆t = et tn = t0 + n ∆t pour 0 ≤ n ≤ N
N
• Le maillage est dit non-uniforme (ou irrégulier), si le pas de discrétisation ∆t est variable.
Dans ce cas ∆t correspond à la taille (ou longueur) maximale des sous-intervalles du maillage.

Ainsi, pour chaque nœud tn , 0 ≤ n ≤ N , on cherche la valeur inconnue yn qui approche

y(tn ) : yn ≈ y(tn )

L’ensemble des valeurs {y0 , y1 , . . . , yN } représente la solution numérique.

Les schémas numériques permettent de calculer yn+1 à partir de yn et il est donc pos-
sible de calculer successivement y1 , y2 , . . ., en partant de y0 .

Si on intègre l’EDO : y 0 (t) = F t, y(t) entre tn et tn+1 , on obtient

Z tn+1


y(tn+1 ) − y(tn ) = F t, y(t) dt
tn

Soit yn une approximation de y(tn ) et yn+1 une approximation de y(tn+1 ). On peut

construire différents schémas numériques selon la formule de quadrature utilisée pour
approcher le membre de droite.

Formules de quadratures / Formules d’intégration numérique :

• Formule de quadrature du rectangle à gauche :

Z tn+1



F t, y(t) dt ≈ ∆t F tn , y(tn )
tn

• Formule de quadrature du rectangle à droite :

Z tn+1



F t, y(t) dt ≈ ∆t F tn+1 , y(tn+1 )
tn

• Formule de quadrature du point milieu :

Z tn+1

∆t ∆t

F t, y(t) dt ≈ ∆t F tn + , y(tn + )
tn 2 2
• Formule du trapèze :
Z tn+1



F tn , y(tn ) + F tn+1 , y(tn+1 )
F t, y(t) dt ≈ ∆t
tn 2
4.2. Méthodes (ou schémas) numériques à un pas. 41

4.2.2 Méthode d’Euler.

• Si on utilise la formule de quadrature du rectangle à gauche, on obtient alors le
schéma d’Euler explicite (ou schéma d’Euler progressif) :

y0 = y(t0 )
yn+1 = yn + ∆t F (tn , yn ), n = 0, . . . , N − 1

• Si on utilise la formule de quadrature du rectangle à droite, on obtient alors le schéma

d’Euler implicite (ou schéma d’Euler rétrograde) :

y0 = y(t0 )
yn+1 = yn + ∆t F (tn+1 , yn+1 ), n = 0, . . . , N − 1

Ces deux méthodes d’Euler sont dites à un pas : pour calculer la solution numérique
yn+1 au noeud tn+1 , on a besoin des informations disponibles au noeud précédent tn .
Plus précisément, pour la méthode d’Euler progressive, yn+1 ne dépend que de la va-
leur yn calculée précédemment, tandis que pour la méthode d’Euler rétrograde, yn+1
dépend aussi "de lui-même" à travers la valeur de F (tn+1 , yn+1 ). C’est pour cette rai-
son que la méthode d’Euler progressive est dite explicite tandis que pour la méthode
d’Euler rétrograde est dite implicite.

• Si on utilise la formule de quadrature du point milieu, on obtient un nouveau schéma :


y0 = y(t0 )
∆t
yn+1 = yn + ∆t F (tn + 2
, yn+ 1 ), n = 0, . . . , N − 1
2

∆t
où yn+ 1 est une approximation de y(tn + 2
).
2

On peut utiliser le schéma d’Euler progressif pour approcher yn+ 1 dans le terme F (tn +
2
∆t
2
, yn+ 1 ) par :
2
∆t
yen+ 1 = yn + F (tn , yn ).
2 2
On construit alors un nouveau schéma appelé schéma d’Euler modifié (ou schéma du
point milieu explicite) :

y = y(t0 )

 0

∆t
yen+ 1 = yn + F (tn , yn )
 2 2
yn+1 = yn + ∆t F (tn + ∆t , yen+ 1 ), n = 0, . . . , N − 1

2 2

Il s’agit d’un schéma explicite car il permet d’expliciter yn+1 en fonction de yn .

42 Chapitre 4. Résolution numérique des équations différentielles ordinaires

4.2.3 Méthode de Crank-Nicolson.

• Si on utilise la formule du trapèze, on obtient alors le schéma de Crank-Nicolson (ou
schéma du trapèze) :
(
y0 = y(t0 )
∆t ∆t
yn+1 = yn + F (tn , yn ) + F (tn+1 , yn+1 ), n = 0, . . . , N − 1
2 2
Il s’agit à nouveau d’un schéma implicite car il ne permet pas d’expliciter directement
yn+1 en fonction de yn lorsque la fonction F n’est pas triviale. En fait, ce schéma fait la
moyenne des schémas d’Euler progressif et rétrograde.

• Pour éviter le calcule implicite de yn+1 dans le schéma de Crank-Nicolson, on peut

utiliser le schéma d’Euler explicite pour approcher yn+1 dans le terme F (tn+1 , yn+1 )
par :
yen+1 = yn + ∆t F (tn , yn ).
On construit ainsi un nouveau schéma appelé schéma de Heun (ou schéma du trapèze
explicite) :

 y0 = y(t0 )

yen+1 = yn + ∆t F (tn , yn )
 yn+1 = yn + ∆t F (tn , yn ) + F (tn+1 , yen+1 ),

n = 0, . . . , N − 1
2
Il s’agit alors d’un schéma explicite car il permet d’expliciter yn+1 en fonction de yn .

4.2.4 Méthode de Runge-Kutta.

Les méthodes de Runge-Kutta (RK) permettent d’obtenir une plus grande précision
que les méthodes précédentes.
La méthode de Runge-Kutta classique est donnée par :
(
y0 = y(t0 )
∆t ∆t ∆t ∆t
yn+1 = yn + y1n + y2n + y3n + y4n , n = 0, . . . , N − 1
6 3 3 6
où
 y1n = F (tn , yn )
∆t ∆t
 y2n = F (tn + 2
, yn + 2
y1n )
∆t ∆t
 y3n = F (tn + 2
, yn + 2
y2n )
 y4n = F (tn+1 , yn + ∆t y3n )
Cette méthode de Runge-Kutta d’ordre 4 est une méthode explicite très populaire. Elle
calcule la valeur de la fonction F en quatre points intermédiaires.
4.3. Étude générale des méthodes numériques à un pas. 43

4.3 Étude générale des méthodes numériques à un pas.

En générale, il faut qu’un schéma numérique soit convergent pour qu’il donne des bons
résultats sur n’importe quelle équation différentielle.

Deux questions naturelles se posent :

1. Que se passe-t-il lorsqu’on fait tendre le pas de discrétisation h vers 0 ?
2. Que se passe-t-il lorsqu’on fixe le pas h > 0, mais on fait tendre T vers l’infini ?
(où T est défini dans l’intervalle I = [t0 , T ])
Dans les deux cas, le nombre de nœuds (N ) tend vers l’infini, mais dans le 1er cas
on s’intéresse à l’erreur en chaque point, dans le 2ème cas il s’agit du comportement
asymptotique de la solution et de son approximation.

4.3.1 Convergence d’une méthode numérique.

Dans les méthodes numériques à un pas, le calcul de yn+1 se fait à partir de tn , yn et h
(où h ≡ ∆t est le pas de discrétisation).
Écrivons la forme générale d’un schéma numérique à un pas :

y0 = y(t0 )
(4.2)
yn+1 = yn + h Φ(tn , yn , h), n = 0, . . . , N − 1

où Φ est une fonction continue sur I × R × [0, hmax ], avec hmax > 0 donné. (Le cas de
T − t0
maillage uniforme, hmax = h ≡ ∆t = , pour I = [t0 , T ].)
N
Ainsi, yn+1 est calculé à partir de yn par l’intermédiaire de la fonction Φ. Cette méthode
peut être explicite ou implicite.
On va étudier les conditions qu’il faut imposer à Φ ainsi que le lien entre Φ et F , pour
que la méthode numérique soit "bonne".

Dans ce paragraphe, on fait une théorie générale des méthodes numériques à un pas
en vue de l’étude de l’erreur de discrétisation :

en = y(tn ) − yn

En fait, il nécessaire d’analyser dans quelle mesure la valeur calculer yn approche suf-
fisamment la valeur exacte y(tn ) et donc, évaluer cette erreur de discrétisation en . Ceci
nous amène alors à introduire les notions de consistance, de stabilité et de convergence.
44 Chapitre 4. Résolution numérique des équations différentielles ordinaires

Définition 4.4 (Consistance)

La méthode numérique (4.2) est dite consistante avec l’équation différentielle (4.1), si
y(tn+1 ) − y(tn )

lim max − Φ tn , y(tn ), h = 0
h→0 0≤n≤N h
pour toute solution y(t) de l’équation différentielle (4.1).
y(tn+1 ) − y(tn )

La quantité en+1 = − Φ tn , y(tn ), h représente dans un certain sens
h
l’erreur que l’on fait au nième pas, en remplaçant l’équation différentielle (4.1) par le
schéma numérique (4.2).

Définition 4.5 (Stabilité)

Soient (yn )0≤n≤N et (zn )0≤n≤N les solutions des systèmes suivants :

y0 fixé
yn+1 = yn + h Φ(tn , yn , h)
et 
z0 = y0 + ξ0  
zn+1 = zn + h Φ(tn , yn , h) + ξn+1
où ξn est quelconque.
La méthode numérique est dite stable s’il existe une constante C indépendante de h, telle que

max |yn − zn | ≤ C max |ξn |

0≤n≤N 0≤n≤N

Cette notion de stabilité implique qu’une petite perturbation sur les données n’entraîne
qu’une petite perturbation sur la solution approchée et ceci indépendamment de h, ce
qui est absolument nécessaire pour le traitement numérique du problème. Un schéma
"instable" ne présente aucun intérêt pratique. Notons que cette notion de stabilité est
intrinsèque au méthode de résolution numérique.

Définition 4.6 (Convergence)

La méthode numérique (4.2) est dite convergente si

lim max |y(tn ) − yn | = 0

h→0 0≤n≤N

Cette notion de convergence indique que l’erreur de discrétisation en = y(tn ) − yn

diminue lorsque h diminue et à la limite, yn doit tendre vers y(tn ) lorsque h tend vers
zéro.

En fait, les trois notions précédentes : consistance, stabilité et convergence, ne sont pas
indépendantes.
4.3. Étude générale des méthodes numériques à un pas. 45

Théorème 4.7
Si la méthode numérique à un pas donnée par (4.2) est consistante et stable, alors elle est
convergente.

Théorème 4.8 (Condition nécessaire et suffisante de consistance)

La méthode numérique (4.2) est consistante si et seulement si

∀t ∈ I, ∀y ∈ R, Φ(t, y, 0) = F (t, y)

Théorème 4.9 (Condition suffisante de stabilité)

Si la fonction Φ est lipschitzienne par rapport à sa deuxième variable :

∃M > 0, ∀t ∈ I, ∀y, y ∈ R, ∀h ∈ [0, hmax ], |Φ(t, y, h) − Φ(t, y, h)| ≤ M |y − y|,

alors la méthode numérique (4.2) est stable

Pouvoir évaluer l’erreur de discrétisation en fonction de h, ceci nous permettra d’obte-

nir l’ordre de convergence de la méthode numérique.

Définition 4.10 (Ordre de convergence)

La méthode numérique (4.2) est dite d’ordre p, si pour toute solution y de l’équation différen-
tielle (4.1), on a
y(tn+1 ) − y(tn )
− Φ tn , y(tn ), h ≤ K hp


max
0≤n≤N h
où K est une constante indépendante de h.
 y(tn+1 ) − y(tn ) 
− Φ tn , y(tn ), h = O(hp ) .


ou encore si max
0≤n≤N h
Théorème 4.11
Si la méthode (4.2) est stable et d’ordre p, alors il existe une constante K > 0 telle que

max yn − y(tn ) ≤ K hp
0≤n≤N

Théorème 4.12 (Condition nécessaire et suffisante pour que la méthode soit d’ordre p)
∂Φ ∂ p−1 Φ
On suppose que F est de classe C p sur I × R et les fonctions Φ, , . . . , p−1 existent et sont
∂h ∂h
continues sur I × R × [0, hmax ]. Alors la méthode (4.2) est d’ordre p, si et seulement si,
pour tout (t, y) ∈ I × R :


 Φ(t, y, 0) = F (t, y)

 ∂Φ 1
(t, y, 0) = F (1) (t, y)



∂h 2

..

 .
 p−1
 ∂ Φ (t, y, 0) = 1 F (p−1) (t, y)



∂hp−1

p
46 Chapitre 4. Résolution numérique des équations différentielles ordinaires

où les fonctions F (k) sont définies par la relation de récurrence :

 (0)
 F =F
∂F (k) ∂F (k)
 F (k) = +F
∂t ∂y

4.3.2 Stabilité absolue d’une méthode numérique.

Dans la section précédente, on a considéré la résolution du problème de Cauchy (4.1)
sur des intervalles bornés. Dans ce cadre, le nombre N de sous-intervalles ne tend vers
l’infini que quand h tend vers zéros. Il existe cependant de nombreuses situations dans
lesquelles le problème de Cauchy doit être intégré sur des intervalles en temps très
grands ou même infini. Dans ce cas, même pour h fixé, N tend vers l’infini. On s’in-
téresse donc à des méthodes capables d’approcher la solution pour des intervalles en
temps arbitrairement grands, même pour des pas h "assez grands".

On considère le problème de Cauchy dans le cas particulier suivant :

 0
y (t) = −β y(t), ∀t > 0
(4.3)
y(0) = y0
avec β > 0 un nombre réel positif donné et y0 6= 0 une donnée initiale. La solution de
ce problème de Cauchy (4.3) est y(t) = y0 e−β t . Puisque β est positif, ce problème est
numériquement bien posé : y(t) décroît exponentiellement lorsque t tend vers l’infini :

lim y(t) = 0
t→+∞

Définition 4.13 (A-Stabilité)

Soit h > 0 un pas de discrétisation donné, tn = n h pour n ∈ N et notons yn ≈ y(tn ) une
approximation de la solution y au point tn . Si, sous d’éventuelles conditions sur h, on a

lim yn = 0
n→+∞

alors on dit que le schéma numérique, appliqué au problème de Cauchy (4.3), est A-stable (ou
absolument stable).

La stabilité d’un schéma numérique décrit la façon dont les erreurs s’accumulent sur
un intervalle borné I = [t0 , T ]. La A-stabilité décrit la façon dont les erreurs traduisent
le comportement de la solution lorsque t tend vers l’infini. En effet, la A-stabilité ex-
prime le fait que, pour h fixé, la solution numérique reproduit le comportement de la
solution exacte quand t tend vers l’infini.
4.3. Étude générale des méthodes numériques à un pas. 47

On peut tirer des conclusions analogues quand β est un complexe ou une fonction
positive de t.

À première vue, il semble que le schéma d’Euler explicite soit préférable aux schémas
d’Euler implicite et de Crank-Nicolson puisque ces derniers ne sont pas explicites.
Cependant, on va voir sur un exemple que le schéma explicite peut engendrer des dif-
ficultés que les schémas implicites n’engendrent pas.

• Le schéma d’Euler explicite est A-stable sous certaine condition sur h, appelée condi-
tion de stabilité ou condition CFL (Courant-Friedrichs-Lewy) : elle limite le pas de dis-
crétisation h d’avance en t. Sinon, on dit que le schéma d’Euler explicite est "instable".
• Contrairement au schéma d’Euler explicite, le schéma d’Euler implicite et le schéma
de Crank-Nicolson sont toujours A-stables, sans condition sur h.

On en conclut que la méthode d’Euler explicite est conditionnellement absolument stable,

tandis que les méthodes d’Euler implicite et de Crank-Nicolson sont inconditionnelle-
ment absolument stables.

De nombreux autres schémas implicites sont inconditionnellement absolument stables.

Cette propriété rend les méthodes implicites attractives, bien qu’elles soient plus coû-
teuses que les méthodes explicites.