Résolution numérique

Équations Différentielles Ordinaires (EDO)

A. Fassi Fihri, M. El Ossmani, A. Taakili

Ecole Nationale Supérieure d’Arts et Métiers de Meknès

Université Moulay Ismail

November 21, 2023

Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 1 / 64
Sommaire

1 Problème de Cauchy et Équations de premier ordre

2 Les systèmes du 1er ordre

3 Les équations différentielles d’ordre >1

4 Approximation numérique des ODE d’ordre 1

Analyse de la méthode d’Euler
Schéma de Heun
Schémas de Runge-Kutta

5 Schémas numériques à pas liés

Principe général
Méthodes d’ADAMS-BASHFORTH
Méthodes d’ADAMS-MOULTON

Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 2 / 64
Introduction
On s’intéresse ici à des problèmes physiques où l’inconnue est une fonction :
Chute libre : Trouver la vitesse v (t) satisfaisant

mv 0 (t) = mg − kv (t)2 ∀t

Pendule : Trouver l’angle θ(t) satisfaisant

ml θ00 (t) = −mg sin θ(t) − lcf θ0 (t) ∀t

Masse-ressort : Trouver le déplacement x (t) satisfaisant

x 00 (t) + cx 0 (t) + ω 2 x (t) = 0 ∀t

Croissance de population : Trouver le nombre d’individus N(t)

N 0 (t) = (n1 − m1 N(t))N(t) ∀t

Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 3 / 64
Introduction : modèle de Lokta-Volterra

Le modèle de Lokta-Volterra s’écrit :

(
y10 (t) = y1 (t)(α − βy2 (t)),
(1)
y20 (t) = y2 (t)(δy1 (t) − γ),

y1 (t)) est l’effectif des proies, y2 (t) l’effectif des prédateurs

t le temps, y10 (t) et y20 (t) représentent la variation des populations au cours de temps.
α taux de reproduction des proies en l’absence des prédateurs;
β taux de mortalité des proies due aux prédateur;
γ taux de mortalité des prédateurs;
δ taux de reproduction des prédateurs en fonction des proies rencontrées et mangées;

Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 4 / 64
Equation différentielle
Une équation différentielle est une relation faisant intervenir les dérivées d’une fonction inconnue y (t)
à déterminer sur un intervalle [t0 , T ] :

F (t, y (t), y 0 (t), y 00 (t), . . . , y (n) (t)) = 0, ∀t ∈ [t0 , T ].

Définition
L’ordre d’une équation différentielle est déterminé par l’ordre maximal de la dérivée de la fonction
inconnue (variable dépendante).

Exemple
Ordre 1 : y 0 (t) = ty (t) alors F (t, y , y 0 ) = y 0 (t) − ty (t).
Ordre 2 : y 00 (x ) = y 0 (x ) − y (x ) alors F (x , y , y 0 , y 00 ) = y 00 (x ) − y 0 (x ) + y (x ).

Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 5 / 64
Problème de Cauchy et Équations de premier ordre

Définition
Soit f une fonction définie de [t0 , T ] × R dans R. Le problème de Cauchy consiste à trouver une
fonction y : [t0 , T ] 7−→ R solution de
(
y 0 (t) = f (t, y (t)), t ∈ [t0 , T ],
(2)
y (t0 ) = y0 .

La condition y (t0 ) = y0 est une condition initiale ou la condition de Cauchy.

Théorème de Cauchy
Si on suppose que la fonction f est continue par rapport aux deux variables t et y et que f est
uniformément Lipschitzienne par rapport à y , c’est à dire que

∃L > 0, ∀t ∈ [t0 , T ], ∀y1 , y2 ∈ R, |f (t, y1 ) − f (t, y2 | ≤ L|y1 − y2 |,

alors le problème de Cauchy (2) admet une unique solution y ∈ C 1 ([0, T ]).
Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 6 / 64
Remarque
Le problème de Cauchy est un problème d’évolution, c’est à dire à partir de la condition initiale, on
peut calculer la solution à l’instant t comme suit :
Z t
y (t) = y (t0 ) + f (s, y (s)) ds, (3)
t0

et donc la solution à l’instant t dépend uniquement de la solution aux instants t0 ≤ s ≤ t.

La solution (3) ne donne pas y de façon explicite sauf dans des cas simples.
Système du 1er ordre
Il s’écrit sous la forme suivante :
(
Y 0 (t) = F (t, Y (t)), t ∈ [t0 , T ],
(4)
Y (t0 ) = Y0 ,

avec Y = (y1 , y2 , . . . , yN ), F : [t0 , T ] × RN 7−→ RN définie par

F (t, Y ) = (f1 (t, Y ), f2 (t, Y ), . . . , fN (t, Y )).

Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 7 / 64
Réduction à l’ordre 1

EDO linéaire : a0 (t)u + a1 (t)u 0 + · · · + an (t)u (n) = b(t) , dite homogène si b(t) = 0
Exemple : u 00 − u 0 + u = t + 1 EDO linéaire d’ordre 2
On pose z1 (t) = u(t), z2 (t) = u 0 (t) et on écrit

z10 = z2
⇔ z 0 (t) = F (t, z(t))
z20 = z2 − z1 + t + 1

   
z1 z2 (t)
avec z = et F (t, z(t)) = ,
z2 z2 (t) − z1 (t) + t + 1

Forme matricielle  0     
z1 (t) 0 1 z1 (t) 0
= +
z2 (t) −1 1 z2 (t) t +1

Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 8 / 64
Réduction à l’ordre 1
EDO non linéaire : Exemple : équation du pendule
( g
θ00 (t) = − sin θ(t),
L (5)
θ(0) = θ0 , θ0 (0) = θ1 .

Cette équation différentielle non linéaire d’ordre 2 se ramène à un système d’ordre 1 en posant y1 (t) = θ(t) et
y2 (t) = θ0 (t). En effet :

g g
y10 (t) = θ0 (t) et y20 (t) = θ00 (t) = − sin θ(t) = − sin y1 (t).
L L
   
y1 (t) y2 (t)
On note Y (t) = et F (t, Y (t)) = ,
y2 (t) − gL sin y1 (t)

alors l’équation (5) est équivalente à 

 Y 0 (t) = 
F (t,Y (t),
0 θ (6)
 Y (0) = θ1 .

Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 9 / 64
Les équations différentielles d’ordre > 1 suite

Dans le cas général d’une équation différentielle de la forme :

u (p) (t) = h(t, u(t), u 0 (t), u 00 (t), . . . , u (p−1) (t))


(7)
u (j) (0) = uj0 , j = 0, . . . , p − 1.

On se ramène à un système d’ordre 1 en posant : y1 = u, y2 = u 0 , . . . , yp = u (p−1) , soit

0

y1 = u 0 ⇒ y1 0 = y2 00,

y2 = u ⇒ y2 = u = y3 ,


. (8)

 ..


yp = u (p−1) ⇒ yp0 = u (p) = h(t, u(t), u 0 (t), u 00 (t), . . . , u (p−1) (t)).


Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 10 / 64
Exemples

1 Le problème du pendule satisfait les hypothèses du Théorème de Cauchy.

En effet : la fonction F de (6) vérifie pour tous Y et Z de R2 :
|F1 (Y ) − F1 (Z )| = |y2 − z2 | ≤ kY − Z k,
|F2 (Y ) − F2 (Z )| = ω 2 |sin(y1 ) − sin(z1 )| ≤ ω 2 kY − Z k,
puisque
y 1 − z1 y 1 + z1
|sin(y1 ) − sin(z1 )| = 2|sin( )cos( )| ≤ |y1 − z1 | ≤ kY − Z k.
2 2
Donc il admet une solution unique. Cependant, on ne peut pas expliciter sa solution à l’aide de fonctions usuelles. Cette solution peut
s’exprimer à l’aide d’une fonction spéciale appelée fonction “sinus amplitude” de Jacobi, mais c’est beaucoup plus compliqué que de rechercher une
approximation numérique de la solution en utilisant une des méthodes que nous allons exposer dans ce cours.

y 0 (t) = y (t)1/2 , t ∈ [0, 1],


2 Le problème ne satisfait pas les hypothèses du Théorème de
y (0) = y0 ,
Cauchy. Il admet au moins deux solutions distinctes y1 (t) = 0 et y2 (t) = (t/2)2 .

Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 11 / 64
Objectif et Vocabulaire

Objectif
Le but est d’étudier des méthodes numériques consistantes et stables permettant le calcul de
bonnes approximations de la solution exacte de l’EDO (2).
Au niveau numérique, on ne peut pas évaluer la solution y (t) en tout point t.
On va chercher une approximation de la solution pour certaines valeurs de la variable indépendante
t que l’on notera ti .
La méthode la plus célèbre est celle de Euler.

Vocabulaire
La distance entre ti+1 et ti est notée hi = ti+1 − ti , on l’appelle le pas du temps.
On utilisera généralement un pas de temps constant noté h.
On note y (ti ) la solution analytique évaluée en t = ti .
On note yi la solution approchée de y (ti ) obtenue par une méthode numérique.

Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 12 / 64
Discrétisation
Supposons que l’on désire calculer la solution y (t) d’une équation différentielle pour t > t0 . Une discrétisation de l’équation définissant y (t) est définie par
T −t0
la donnée d’un ensemble discret de valeur pour la variable indépendante t : en général, on fixe le pas h = N
avec N est le nombre de points de la
discrétisation. On pose
ti = t0 + ih ∀i = 0, 1, . . . , N.

Produisant la suite
t0 < t1 < t2 < · · · < ti < · · · < tN = T = t0 + Nh,

et d’une méthode numérique (schéma) fournissant une approximation de la solution exacte au points ti :

yi ≈ y (ti ).

Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 13 / 64
Approximation numérique des ODE d’ordre 1
Méthode de Euler. On se donne une subdivision de I = [t0 , T ] en N intervalles de pas

T − t0
h = tn+1 − tn = .
N
La méthode d’Euler consiste à approcher y 0 (tn ) par la formule de Taylor comme suit :

y (tn+1 ) = y (tn ) + hy 0 (tn ) + O(h2 ).

Soit
y (tn+1 ) − y (tn )
y 0 (tn ) = + O(h) = f (tn , y (tn )).
h
Soit yn une approximation de y (tn ) (yn ≈ y (tn )), le schéma d’Euler s’écrit

yn+1 = yn + hf (tn , yn ), n = 0, . . . , N − 1,
y (0) = y0 .

Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 14 / 64
Interprétation de la méthode d’Euler

Définitions
On appelle schéma à un pas si yn+1 est une fonction de tn et yn uniquement (ex. Euler).
On appelle schéma à deux pas si yn+1 est une fonction de tn , yn , tn−1 , yn−1 uniquement.
On appelle schéma à k pas si yn+1 est une fonction de tn , yn , tn−1 , yn−1 , . . . , tn−(k−1) , yn−(k−1)
On appelle schéma à k pas est explicite si on peut exprimer yn+1 sous la forme

yn+1 = E (tn , tn−1 , . . . , tn−(k−1) , yn , yn−1 , . . . , yn−(k−1) ).

Exemple. Le schéma d’Euler est explicite.

Un schéma à k pas est implicite si yn+1 est solution d’une équation non linéaire de la forme

yn+1 = I(tn+1 , tn , tn−1 , . . . , tn−(k−1) , yn+1 , yn , yn−1 , . . . , yn−(k−1) ).

Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 15 / 64
Schémas classiques

Une façon d’obtenir une multitude de schémas numériques est d’intégrer l’EDO sur [tn , tn+1 ] :
Z tn+1
y (tn+1 ) − y (tn ) = f (t, y (t))dt,
tn

et ensuite d’approcher l’intégrale. Par exemple

Intégration par la méthode des rectangles à gauche :

Z tn+1
f (t, y (t))dt ≈ hf (tn , y (tn )),
tn

ce qui donne le schéma d’Euler explicite : yn+1 = yn + hf (tn , yn )

Intégration par la méthode des rectangles à droite :

Z tn+1
f (t, y (t))dt ≈ hf (tn+1 , y (tn+1 )),
tn

ce qui donne le schéma d’Euler implicite: yn+1 = yn + hf (tn+1 , yn+1 ).

yn+1 est solution d’une équation non linéaire.
Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 16 / 64
Schémas classiques

Une façon d’obtenir une multitude de schémas numériques est d’intégrer l’EDO sur [tn , tn+1 ] :
Z tn+1
y (tn+1 ) − y (tn ) = f (t, y (t))dt,
tn

et ensuite d’approcher l’intégrale. Par exemple

Intégration par la méthode des rectangles à gauche :

Z tn+1
tn tn+1
f (t, y (t))dt ≈ hf (tn , y (tn )),
tn

ce qui donne le schéma d’Euler explicite : yn+1 = yn + hf (tn , yn )

Intégration par la méthode des rectangles à droite :

Z tn+1
f (t, y (t))dt ≈ hf (tn+1 , y (tn+1 )),
tn

ce qui donne le schéma d’Euler implicite: yn+1 = yn + hf (tn+1 , yn+1 ).

yn+1 est solution d’une équation non linéaire.
Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 16 / 64
Schémas classiques

Une façon d’obtenir une multitude de schémas numériques est d’intégrer l’EDO sur [tn , tn+1 ] :
Z tn+1
y (tn+1 ) − y (tn ) = f (t, y (t))dt,
tn

et ensuite d’approcher l’intégrale. Par exemple

Intégration par la méthode des rectangles à gauche :

Z tn+1
tn tn+1
f (t, y (t))dt ≈ hf (tn , y (tn )),
tn

ce qui donne le schéma d’Euler explicite : yn+1 = yn + hf (tn , yn )

Intégration par la méthode des rectangles à droite :

Z tn+1
f (t, y (t))dt ≈ hf (tn+1 , y (tn+1 )),
tn

ce qui donne le schéma d’Euler implicite: yn+1 = yn + hf (tn+1 , yn+1 ).

yn+1 est solution d’une équation non linéaire.
Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 16 / 64
Schémas classiques

Une façon d’obtenir une multitude de schémas numériques est d’intégrer l’EDO sur [tn , tn+1 ] :
Z tn+1
y (tn+1 ) − y (tn ) = f (t, y (t))dt,
tn

et ensuite d’approcher l’intégrale. Par exemple

Intégration par la méthode des rectangles à gauche :

Z tn+1
tn tn+1
f (t, y (t))dt ≈ hf (tn , y (tn )),
tn

ce qui donne le schéma d’Euler explicite : yn+1 = yn + hf (tn , yn )

Intégration par la méthode des rectangles à droite :

Z tn+1
f (t, y (t))dt ≈ hf (tn+1 , y (tn+1 )),
tn

ce qui donne le schéma d’Euler implicite: yn+1 = yn + hf (tn+1 , yn+1 ).

yn+1 est solution d’une équation non linéaire.
Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 16 / 64
Schémas classiques

Une façon d’obtenir une multitude de schémas numériques est d’intégrer l’EDO sur [tn , tn+1 ] :
Z tn+1
y (tn+1 ) − y (tn ) = f (t, y (t))dt,
tn

et ensuite d’approcher l’intégrale. Par exemple

Intégration par la méthode des rectangles à gauche :

Z tn+1
tn tn+1
f (t, y (t))dt ≈ hf (tn , y (tn )),
tn

ce qui donne le schéma d’Euler explicite : yn+1 = yn + hf (tn , yn )

Intégration par la méthode des rectangles à droite :

Z tn+1
f (t, y (t))dt ≈ hf (tn+1 , y (tn+1 )),
tn

ce qui donne le schéma d’Euler implicite: yn+1 = yn + hf (tn+1 , yn+1 ).

yn+1 est solution d’une équation non linéaire.
Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 16 / 64
Schémas classiques

Une façon d’obtenir une multitude de schémas numériques est d’intégrer l’EDO sur [tn , tn+1 ] :
Z tn+1
y (tn+1 ) − y (tn ) = f (t, y (t))dt,
tn

et ensuite d’approcher l’intégrale. Par exemple

Intégration par la méthode des rectangles à gauche :

Z tn+1
tn tn+1
f (t, y (t))dt ≈ hf (tn , y (tn )),
tn

ce qui donne le schéma d’Euler explicite : yn+1 = yn + hf (tn , yn )

Intégration par la méthode des rectangles à droite :

Z tn+1
f (t, y (t))dt ≈ hf (tn+1 , y (tn+1 )),
tn tn tn+1
ce qui donne le schéma d’Euler implicite: yn+1 = yn + hf (tn+1 , yn+1 ).
yn+1 est solution d’une équation non linéaire.
Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 16 / 64
Schémas classiques

Une façon d’obtenir une multitude de schémas numériques est d’intégrer l’EDO sur [tn , tn+1 ] :
Z tn+1
y (tn+1 ) − y (tn ) = f (t, y (t))dt,
tn

et ensuite d’approcher l’intégrale. Par exemple

Intégration par la méthode des rectangles à gauche :

Z tn+1
tn tn+1
f (t, y (t))dt ≈ hf (tn , y (tn )),
tn

ce qui donne le schéma d’Euler explicite : yn+1 = yn + hf (tn , yn )

Intégration par la méthode des rectangles à droite :

Z tn+1
f (t, y (t))dt ≈ hf (tn+1 , y (tn+1 )),
tn tn tn+1
ce qui donne le schéma d’Euler implicite: yn+1 = yn + hf (tn+1 , yn+1 ).
yn+1 est solution d’une équation non linéaire.
Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 16 / 64
Schémas classiques ...

Intégration par la méthode des trapèzes

Z tn+1
h
f (t, y (t))dt ≈ (f (tn , y (tn )) + f (tn+1 , y (tn+1 ))),
tn 2

ce qui donne le schéma de Cranck-Nicholson semi-implicite :

tn tn+1
h
yn+1 = yn + (f (tn , yn ) + f (tn+1 , yn+1 )).
2
Intégration par la méthode du point milieu
Z tn+1
h h
f (t, y (t))dt ≈ hf (tn + , y (tn + ))),
tn 2 2
h
Pour approcher la solution au point tn + 2
, on utilise le schéma d’Euler explicite :

y (tn + h
2
,
) ≈ y (tn ) + h
2
f (tn , y (tn ))

ce qui donne le schéma d’Euler modifié explicite :

h h
yn+1 = yn + hf (tn + , yn + f (tn , yn )).
2 2
Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 17 / 64
Schémas classiques ...

Intégration par la méthode des trapèzes

Z tn+1
h
f (t, y (t))dt ≈ (f (tn , y (tn )) + f (tn+1 , y (tn+1 ))),
tn 2

ce qui donne le schéma de Cranck-Nicholson semi-implicite :

tn tn+1
h
yn+1 = yn + (f (tn , yn ) + f (tn+1 , yn+1 )).
2
Intégration par la méthode du point milieu
Z tn+1
h h
f (t, y (t))dt ≈ hf (tn + , y (tn + ))),
tn 2 2
h
Pour approcher la solution au point tn + 2
, on utilise le schéma d’Euler explicite :

y (tn + h
2
,
) ≈ y (tn ) + h
2
f (tn , y (tn ))

ce qui donne le schéma d’Euler modifié explicite :

h h
yn+1 = yn + hf (tn + , yn + f (tn , yn )).
2 2
Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 17 / 64
Schémas classiques ...

Intégration par la méthode des trapèzes

Z tn+1
h
f (t, y (t))dt ≈ (f (tn , y (tn )) + f (tn+1 , y (tn+1 ))),
tn 2

ce qui donne le schéma de Cranck-Nicholson semi-implicite :

tn tn+1
h
yn+1 = yn + (f (tn , yn ) + f (tn+1 , yn+1 )).
2
Intégration par la méthode du point milieu
Z tn+1
h h
f (t, y (t))dt ≈ hf (tn + , y (tn + ))),
tn 2 2
h
Pour approcher la solution au point tn + 2
, on utilise le schéma d’Euler explicite :

y (tn + h
2
,
) ≈ y (tn ) + h
2
f (tn , y (tn ))

ce qui donne le schéma d’Euler modifié explicite :

h h
yn+1 = yn + hf (tn + , yn + f (tn , yn )).
2 2
Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 17 / 64
Schémas classiques ...

Intégration par la méthode des trapèzes

Z tn+1
h
f (t, y (t))dt ≈ (f (tn , y (tn )) + f (tn+1 , y (tn+1 ))),
tn 2

ce qui donne le schéma de Cranck-Nicholson semi-implicite :

tn tn+1
h
yn+1 = yn + (f (tn , yn ) + f (tn+1 , yn+1 )).
2
Intégration par la méthode du point milieu
Z tn+1
h h
f (t, y (t))dt ≈ hf (tn + , y (tn + ))),
tn 2 2
h
Pour approcher la solution au point tn + 2
, on utilise le schéma d’Euler explicite :

y (tn + h
2
,
) ≈ y (tn ) + h
2
f (tn , y (tn ))
tn tn+1
ce qui donne le schéma d’Euler modifié explicite :

h h
yn+1 = yn + hf (tn + , yn + f (tn , yn )).
2 2
Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 17 / 64
Schémas classiques ...

Intégration par la méthode des trapèzes

Z tn+1
h
f (t, y (t))dt ≈ (f (tn , y (tn )) + f (tn+1 , y (tn+1 ))),
tn 2

ce qui donne le schéma de Cranck-Nicholson semi-implicite :

tn tn+1
h
yn+1 = yn + (f (tn , yn ) + f (tn+1 , yn+1 )).
2
Intégration par la méthode du point milieu
Z tn+1
h h
f (t, y (t))dt ≈ hf (tn + , y (tn + ))),
tn 2 2
h
Pour approcher la solution au point tn + 2
, on utilise le schéma d’Euler explicite :

y (tn + h
2
,
) ≈ y (tn ) + h
2
f (tn , y (tn ))
tn tn+1
ce qui donne le schéma d’Euler modifié explicite :

h h
yn+1 = yn + hf (tn + , yn + f (tn , yn )).
2 2
Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 17 / 64
Forme générale d’un schéma à un pas

La forme générale d’un schéma à un pas

yn+1 = yn + hΦ(tn , yn ; h)
où Φ ne dépend que de f et est appelée fonction incrément
Le choix de la fonction incrément Φ détermine le schéma
Exemples :
1 Euler explicite : yn+1 = yn + hΦ(tn , yn ; h) avec

Φ = f (tn , yn ).

2 Euler modifié : yn+1 = yn + hΦ(tn , yn ; h) avec

f (tn + h2 , yn + h2 f (tn , yn )
Φ= .
2

Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 18 / 64
Une mise en oeuvre réussie

Le problème
y 0 (t) = y (t),

t ∈ [0, 1],
(9)
y (0) = 1,
admet comme unique solution y (t) = e t .
La méthode d’Euler explicite consiste à calculer, pour n ≥ 0,

yn+1 = yn + hyn
= (1 + h)yn .

Partons de y0 = 1, on obtient yn = (1 + h)n . Avec h = N1 , on retrouve pour yN une valeur approchée de

1
y (1) = e, puisque lim yN = lim (1 + )N = e.
N→∞ N→∞ N
La méthode implicite fournit quant à elle yN = (1−11 )N qui converge également vers e lorsque N tend
N
vers l’infini.

Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 19 / 64
 Figure: Les méthodes d’Euler pour le problème (9)

Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 20 / 64
En fait, ces comportements s’expliquent par un développement limité :
 1
yN,exp = e N ln(1+ N ) = e(1 − 1
+ O( N12 )),
1
2N (10)
yN,imp = e −N ln(1− N ) = e(1 + 1
2N
+ O( N12 )).

Cette estimation est confortée par les valeurs données dans le tableau suivant, où εN désigne l’erreur absolue au point
t = 1 pour chacune des méthodes.

Euler Explicite Euler Implicite

εN εN
N εN h
εN h
10 0.1245 1.2454 0.1497 1.4969
25 0.0524 1.3111 0.0564 1.4110
50 0.0267 1.3347 0.0277 1.3845
100 0.0135 1.3468 0.0137 1.3717
250 0.0054 1.3542 0.0055 1.3641
500 0.0027 1.3567 0.0027 1.3616
1000 0.0014 1.3579 0.0014 1.3604
εN e
Les valeurs du rapport h
encadrent de mieux en mieux 2
≈ 1.3591.

Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 21 / 64
Une mise en oeuvre défaillante (exemple surprenant)

Considérons le problème suivant : 

y 0 (t) = 3y (t) − 3t, t ∈ [0, T ],
y (0) = 13 ,
(11)

dont la solution est y (t) = 13 + t. Appliquons la méthode d’Euler pour en chercher la solution approchée à l’instant
T , pour différentes valeurs de T et différents choix du pas h. Les résultats obtenus peuvent surprendre. Le tableau
suivant donne les erreurs dans le calcul de yN :

h T =5 T = 10 T = 20
1 1.8651 10−14 1.9403 10−11 2.0345 10−05
0.5 4.9560 10−13 4.7278 10−09 4.2999 10−01
0.1 5.2402 10−14 2.6318 10−08 6.5252 10+03
0.05 1.5525 10−11 1.8232 10−05 2.5142 10+07

On constate que l’erreur semble croître avec le pas, contrairement à ce que l’on espérait. Par ailleurs, pour T = 20,
cette erreur atteint des valeurs astronomiques !!!

Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 22 / 64
Analyse de l’exemple

Calculons la solution générale de l’équation différentielle : y 0 (t) = 3y (t) − 3t.

La solution de l’équation homogène est y (t) = Ke 3t . La méthode de variation de la constante fournit pour notre
équation

1 −3t 1
K 0 (t) = −3te −3t ⇒ K (t) = (t + )e + k ⇒ y (t) = t + + ke 3t .
3 3

La constante k peut être évaluée à l’aide de la valeur de y en 0, avec k = 0 si y (0) = 13 . Mais si y (0) = 13 + ε. On
obtient k = ε implique que la solution réelle du problème devient yε (t) = 31 + t + εe 3t . Les valeurs de e 3T sont
respectivement :
pour T = 5, e 3T ' 3.269 106 ,
pour T = 10, e 3T ' 1.068 1013 ,
pour T = 20, e 3T ' 1.142 1026 .
Avec une unité d’arrondi de l’ordre de 1.1102 10−16 , on retrouve l’ordre de grandeur des erreurs
constatées ! Le problème (11) est typiquement un problème instable; une petite perturbation de la donnée initiale
entraîne une grande perturbation de la solution.

Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 23 / 64
Sommaire

1 Problème de Cauchy et Équations de premier ordre

2 Les systèmes du 1er ordre

3 Les équations différentielles d’ordre >1

4 Approximation numérique des ODE d’ordre 1

Analyse de la méthode d’Euler
Schéma de Heun
Schémas de Runge-Kutta

5 Schémas numériques à pas liés

Principe général
Méthodes d’ADAMS-BASHFORTH
Méthodes d’ADAMS-MOULTON

Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 24 / 64
Analyse des méthodes à un pas

Le but est que l’écart entre la solution exacte et la solution approchée

en = u(tn ) − un , n = 0, · · · N

soit petit pour h assez petit. (convergence)

La solution exacte u(tn+1 ) est généralement n’est pas connue
∗
On introduit un+1 solution de
∗
un+1 = u(tn ) + hφ(tn , u(tn ), h)
∗ ∗
L’erreur en+1 = (u(tn+1 ) − un+1 ) + (un+1 − un+1 )
∗
L’écart εn+1 = u(tn+1 ) − un+1 est évalué par développement limité (consistance)
∗
L’écart un+1 − un+1 représente la propagation de l’erreur accumulé antre tn et tn+1 (stabilité)

Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 25 / 64
Analyse des méthodes à un pas
Consistance et ordre

L’erreur de consistance du schéma (∗) à l’instant tn+1 est définie par

∗
u(tn+1 ) − un+1 u(tn+1 ) − u(tn )
τn+1 (h) = = − φ(tn , u(tn ), h)
h h
ou encore
u(tn+1 ) = u(tn ) + hφ(tn , u(tn ), h) + hτn+1 (h)
L’erreur de consistance globale
τ (h) = max |τn+1 (h)|.
0≤n≤N−1

Le schéma (∗) est alors consistant si on a lim |τ (h)| = 0.

h→0
Il est d’ordre p si |τ (h)| ≤ C hp (schéma p-consistant).

Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 26 / 64
Analyse des méthodes à un pas
Consistance et ordre

Définition de la consistance indépendamment des points de la discrétisation

Le schéma (∗) est dit consistant si ∀u ∈ C 1 , ∀t ∈ [a, b], l’erreur de consistance

u(t + h) − u(t)
R(t, u, h) = − φ(t, u(t), h)
h

tends vers 0 lorsque h tends vers 0.

ou encore
u(t + h) = u(t) + hφ(t, u(t), h) + hR(t, u, h)
Le schéma (∗) est consistant d’ordre p > 0 si kR(t, u, h)k ≤ C (h)p
On remarque que τn+1 (h) = R(tn , u, h)

Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 27 / 64
Analyse des méthodes à un pas
Consistance et ordre

Pour déterminer la consistance et l’ordre d’un schéma, nous utilisons la formule de Taylor.

Théorème (Taylor)
Supposons u de classe C p ([a, b]), alors ∀t ∈ [a, b] et ∀h ∈ R tels que t + h ∈ [a, b] on a
p
u(t + h) = u(t) + hu 0 (t) + · · · + hp! u (p) (t) + O(hp+1 )

On peut écrire aussi

(−h)p
u(t) = u(t + h) − hu 0 (t + h) + · · · + p!
u (p) (t + h) + O(hp+1 )
Exemples :
1 u(t + h) = u(t) + hu 0 (t) + O(h2 ) (pour u ∈ C 1 )
2 u(t) = u(t + h) − hu 0 (t + h) + O(h2 ) (pour u ∈ C 1 )
3 u(t + h) = u(t) + O(h1 ) (pour u ∈ C 0 )

Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 28 / 64
Consistance des schémas d’Euler explicite et implicite
Le schéma d’Euler explicite est consistant d’ordre 1
u(t+h)−u(t)
L’erreur de consistance : R(t, u, h) = h
− f (t, u(t))

u(t + h) = u(t) + hu 0 (t) + O(h2 )

= u(t) + hf (t, u(t)) + O(h2 )
u(t+h)−u(t)
R(t, u, h) = h − f (t, u(t)) = O(h)

Le schéma d’Euler implicite est consistant d’ordre 1

u(t+h)−u(t)
L’erreur de troncature : R(t, u, h) = h
− f (t + h, u(t + h))

u(t) = u(t + h) − hu 0 (t + h) + O(h2 )

= u(t + h) − hf (t + h, u(t + h)) + O(h2 )
u(t+h)−u(t)
R(t, u, h) = h − f (t + h, u(t + h)) = O(h)

Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 29 / 64
Consistance du schéma semi-implicite

Le schéma de Cranck Nicholson est consistant d’ordre 2

u(t+h)−u(t) f (t+h,u(t+h))+f (t,u(t))
L’erreur de troncature : R(t, u, h) = h − 2
h2 00
u(t + h) = u(t) + hu 0 (t) + 2 u (t) + O(h3 ) (1)
h2 00
u(t) = u(t + h) − hu 0 (t + h) + 2 u (t + h) + O(h3 ) (2)
Alors (1)- (2) ⇒
O(h3 )
z }| {
h2
u(t + h) = u(t) + h2 (u 0 (t) + u 0 (t + h)) + (u 00 (t) − u 00 (t + h)) +O(h3 )
4
h
= u(t) + 2 (f (t, u(t)) + f (t + h, u(t + h))) + O(h3 )
Au final R(t, u, h) = O(h2 ).

Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 30 / 64
Consistance du schéma d’Euler modifié

Le schéma d’Euler modifié est consistant d’ordre 2

u(t+h)−u(t)
 
L’erreur de troncature R(t, u, h) = h
− f t + h2 , u(t) + h2 f (u(t), t)
h2 00
u(t + h) = u(t + h2 ) + h2 u 0 (t + h2 ) + 8 u (t + h2 ) + O(h3 )
h2 00
u(t) = u(t + h2 ) − h2 u 0 (t + h2 ) + 8 u (t + h2 ) + O(h3 )
donc u(t + h) = u(t) + hu 0 (t + h2 ) + O(h3 ) = u(t) + hf (t + h2 , u(t + h2 )) + O(h3 )
Et u(t + h2 ) = u(t) + h2 u 0 (t) + O(h2 ) = u(t) + h2 f (t, u(t)) + O(h2 )
donc
h i
u(t + h) = u(t) + hf t + h2 , u(t) + h2 f (u(t), t) + O(h3 )
h i
u(t+h)−u(t)
D’où R(t, u, h) = h − f t + h2 , u(t) + h2 f (u(t), t) = O(h2 )

Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 31 / 64
Consistance et ordre (astuce)
La formule de Taylor donne

u(t + h) − u(t) h h2 (3) hp−1 (p)

= u 0 (t) + u 00 (t) + u (t) + . . . + u (t) + . . .
h 2! 3! p!

∂Φ h2 ∂ 2 Φ
Φ(t, u, h) = Φ(t, u(t), 0) + h (t, u, 0) + (t, u, 0)
∂h 2! ∂ 2 h
h3 ∂ 3 Φ hp−1 ∂ p−1 Φ
+ (t, u, 0) + . . . + (t, u, 0) + . . .
3! ∂ 3 h (p − 1)! ∂ p−1 h

alors

 
h2 ∂2Φ
h1 ∂Φ
i 1 (3)
R(t, u, h) = [u 0 (t) − Φ(t, u, 0)] + h u 00 (t) − (t, u, 0) + u (t) − 2 (t, u, 0)
2 ∂h 2! 3 ∂ h
p−1
 
h 1 (p) ∂ p−1 Φ
+ ... + u (t) − p−1 (t, u, 0) .
(p − 1)! p ∂ h

Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 32 / 64
Consistance et ordre (astuce)
Le schéma est au moins d’ordre 1 (consistant) si

Φ(t, u, 0) = u 0 (t) = f (t, u(t))

Le schéma est au moins d’ordre 2 si de plus

∂Φ 1 1 d 1
(t, u, 0) = u 00 (t) = f (t, u(t)) = [∂t f (t, u(t)) + u 0 (t)∂u f (t, u(t))]
∂h 2 2 dt 2
1
= [∂t f (t, u(t)) + f (t, u(t))∂u f (t, u(t))] .
2
Le schéma est au moins d’ordre 3 si de plus

∂2Φ 1 1 d2
(t, u, 0) = u (3) (t) = f (t, u(t)).
∂2h 3 3 dt 2
Le schéma est au moins d’ordre p si de plus

∂ p−1 Φ 1 1 d p−1
(t, u, 0) = u (p) (t) = f (t, u(t)).
∂ p−1 h p p dt p−1

Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 33 / 64
Consistance et ordre (astuce)
Montrer que le schéma semi-implicite est consistant d’ordre 2

f (t+h,u(t+h))+f (t,u(t))
On a Φ(t, u, h) = 2 , donc Φ(t, u, 0) = f (t, u(t))
Et ∂Φ
∂h (t, u, h) = 1
2 [∂t f (t + h, u(t + h)) + u 0 (t + h)∂u f (t + h, u(t + h))] .
Alors ∂Φ
∂h (t, u, 0) = 1
2 [∂t f (t, u(t)) + u 0 (t)∂u f (t, u(t))] .
Le schéma est alors consistant d’ordre 2.

Montrer que le schéma d’Euler modifié est consistant d’ordre 2

On a Φ(t, u, h) = f [t + h2 , u(t) + h2 f (t, u(t))], donc Φ(t, u, 0) = f (t, u(t))

∂Φ
Et ∂h (t, u, h) = 12 ∂t f [t + h2 , u(t) + h2 f (t, u(t))] + 12 f (t, u(t))∂u f [t + h2 , u(t) + h2 f (t, u(t))].
∂Φ 1
Alors ∂h (t, u, 0) = 2 [∂t f (t, u(t)) + f (t, u(t))∂u f (t, u(t))] .
Le schéma est alors consistant d’ordre 2.
Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 34 / 64
Stabilité

Soient (un )n∈{0,··· ,N} solution du schéma (∗) et (vn )n∈{0,··· ,N} solution du schéma perturbé

vn+1 = vn + hφ(tn , vn , h) + εn

Définition (Stabilité)
Le schéma (∗) est dit stable si

N−1
!
X
max |un − vn | ≤ C |u0 − v0 | + |εn |
n=0,··· ,N
n=0

Ceci signifie que le schéma est peu sensible aux erreurs (de méthode, de données, de troncature)

Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 35 / 64
Stabilité

La stabilité d’un schéma est évaluée comme suit :

1 Introduire une perturbation à l’instant n : ũn = un + εn
2 Calculer l’évolution à l’instant n + 1 : εn+1
ũn+1 = ũn + hφ(tn , ũn , h)
un+1 + εn+1 = un + εn + hφ(tn , un + εn , h)
On utilise la formule de Taylor
φ(tn , un + εn , h) = φ(tn , un , h) + εn ∂φ(t∂u
n ,un ,h)
+ O(ε2n )
On aura  
∂φ(tn , un , h)
εn+1 = εn 1 + h + hO(εn )
∂u

Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 36 / 64
Stabilité

Soit l’équation différentielle du

dt = λu, u(0) = u0
 Le schéma d’Euler explicite :un+1 = (1 + λh)un
un+1 + εn+1 = (1 + λh)(un + εn )
εn+1 = (1 + λh)εn
1
 Le schéma d’Euler implicite :un+1 = 1−λh un
1
εn+1 = 1−λh εn
1+ λh
 Le schéma de Crank Nicholson :un+1 = 2
1− λh
un
2
1+ λh
εn+1 = 2
1− λh
εn
2
 La stabilité de ces schémas dépend de λ et h.

Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 37 / 64
Convergence

Définition (Convergence)
Le schéma (∗) est dit convergent si pour toute solution exacte u et toute solution approchée un , on a

lim max |en | = 0.

h→0 0≤n≤N

Proposition
Un schéma numérique à un pas qui est stable et consistant est convergent.

Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 38 / 64
Convergence du schéma d’Euler

Proposition
Supposons u ∈ C 2 ([a, b]). Alors la méthode d’Euler explicite est convergente, et on a la majoration de l’erreur
suivante
M e (b−a)L − 1
|en | ≤ h
2 L
où M = maxt∈[a;b] |u 00 (t)| et L est la constante de Lipschitz.

D’après Lagrange, il existe ξn ∈ [tn ; tn+1 ]

h2 0
u(tn+1 ) = u(tn ) + hf (tn , u(tn )) + f (ξn , u(ξn ))
2

D’autre part un+1 = un + hf (tn , un )

h2
donc u(tn+1 ) − un+1 = u(tn ) − un + h[f (tn , u(tn )) − f (tn , un )] + 2
f 0 (ξn , u(ξn ))
h2 h2
Et f Lipschitzienne, alors |en+1 | ≤ |en | + hL|en | + 2
M = |en |(1 + hL) + 2
M
nδ
−1
On utilise le résultat suivant pour conclure : si dn est une suite vérifiant dn+1 ≤ dn (1 + δ) + µ, avec δ, µ > 0 alors dn ≤ e nδ d0 + µ e δ
avec ici d0 = e0 = 0, δ = hL
2
h
et µ = 2
M

Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 39 / 64
Consistance, Stabilité, Convergence

Discrétisation Schéma
EDO Consistance numérique

Résolution
Résolution? Stabilité

Solution Représentativité Solution

théorique Convergence approchée

Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 40 / 64
Consistance, Stabilité, Convergence

Discrétisation Schéma
EDO Consistance numérique

Résolution
Résolution? Stabilité

Solution Représentativité Solution

théorique Convergence approchée

Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 40 / 64
Consistance, Stabilité, Convergence

Discrétisation Schéma
EDO Consistance numérique

Résolution
Résolution? Stabilité

Solution Représentativité Solution

théorique Convergence approchée

Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 40 / 64
Sommaire

1 Problème de Cauchy et Équations de premier ordre

2 Les systèmes du 1er ordre

3 Les équations différentielles d’ordre >1

4 Approximation numérique des ODE d’ordre 1

Analyse de la méthode d’Euler
Schéma de Heun
Schémas de Runge-Kutta

5 Schémas numériques à pas liés

Principe général
Méthodes d’ADAMS-BASHFORTH
Méthodes d’ADAMS-MOULTON

Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 41 / 64
Schéma de Heun - schéma prédicteur/correcteur

On écrit la forme intégrale sur [tn , tn+1 ]

Z tn+1
u(tn+1 ) = u(tn ) + f (t, u(t))dt
tn

Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 42 / 64
Schéma de Heun - schéma prédicteur/correcteur
R tn+1
L’intégrale tn
f (t, u(tn ) est approchée par la méthode des trapèzes, on aura

h
un+1 = un + (f (tn , un ) + f (tn+1 , un+1 ))
2

Cette forme est implicite en un+1 qui sera approchée par pn+1 grace au schéma d’Euler (schéma prédicteur)

pn+1 = un + hf (un , tn ) ≈ un+1

Le schéma de Heun (schéma correcteur)

h
un+1 = un + [f (tn , un ) + f (tn+1 , un + hf (un , tn ))]
2

La fonction incrément dans ce cas est

1
φ(t, u; h) = [f (tn , un ) + f (tn+1 , un + hf (un , tn ))]
2

Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 43 / 64
Sommaire

1 Problème de Cauchy et Équations de premier ordre

2 Les systèmes du 1er ordre

3 Les équations différentielles d’ordre >1

4 Approximation numérique des ODE d’ordre 1

Analyse de la méthode d’Euler
Schéma de Heun
Schémas de Runge-Kutta

5 Schémas numériques à pas liés

Principe général
Méthodes d’ADAMS-BASHFORTH
Méthodes d’ADAMS-MOULTON

Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 44 / 64
Schémas de Runge-Kutta (Principe général)

On introduit q points intermédiaires dans l’intervalle [tn , tn+1 ] :

tn,i = tn + αi h, αi ∈ [0, 1], i = 1, · · · , q
On écrit la forme Zintégrale sur [tn , tn,i ]
tn,i Z αi
u(tn,i ) = u(tn ) + f (t, u(t))dt = u(tn ) + h f [tn + s h, u(tn + s h)] ds
tn 0
i−1 i−1
X X
≈ u(tn ) + h aij f [tn,j , u(tn,j ] aij = αi
j=1 j=1

La forme intégrale sur [tn , tn+1 ]

Z 1
u(tn+1 ) = u(tn ) + h f [tn + s h, u(tn + s h)] ds
0
q q
X X
≈ u(tn ) + h ωi f [tn,i , u(tn,i )] wi = 1
i=1 " i=1 #
q i−1
X X
≈ u(tn ) + h ωi f tn,i , u(tn ) + h aij f [tn,j , u(tn,j )]
i=1 j=1

Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 45 / 64
Schémas de Runge-Kutta (Principe général)

En résumé, les méthodes de Runge-Kutta de rang q s’écrit

q
X
un+1 = un + h ωi ki
i=1
Pq
où i=1 ωi = 1 et les fonctions kj sont définies par

 k1 (tn , un ) =
 f (tn , un ) !
i−1
P
 ki (tn , un )
 = f tn,i , u(tn ) + h aij f [tn,j , u(tn,j )] , ∀i ≥ 2
j=1

Les poids ωi , et les coefficients αi et aij sont à choisir de tel sort à assurer la consistance.

Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 46 / 64
Schémas de Runge-Kutta (tableau de Butcher)

On résume la méthode de Runge Kutta par le tableau de Butcher

α2 a21
α3 a31 a32
.. .. .. ..
. . . .
αq−1 aq−1,1 aq−1,2 · · · aq−1,q−2
αq aq,1 aq,2 ··· aq,q−2 aq,q−1
ω1 ω2 ω3 ··· ωq−1 ωq
La méthode de Runge-Kutta est consistante si
i−1
X
ai,j = αi , pour i = 2, · · · , q
j=1

Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 47 / 64
Schéma de Runge-Kutta d’ordre 2 (RK2)

Pour un schéma RK2, on écrit

un+1 = un + h(ω1 k1 + ω2 k2 )

où 
k1 = f (tn , un )
k2 = f (tn + α2 h, un + a2,1 hf (tn , un ))
1
On a ω1 + ω2 = 1 et ω2 α2 = ω2 a2,1 = 2
(Formule de Taylor)
En effet, on a :
un+1 = un + h [ω1 f (tn , un ) + ω2 f (tn + α2 h, un + a2,1 hf (tn , un ))]
= un + h [ω1 f (tn , un ) + ω2 f (tn , un ) + ω2 α2 h∂t f (tn , un )
+ω2 a2,1 hf (tn , un )∂u f (tn , un )]
Et d’après le développement de Taylor :
2
u(tn + h) = u(tn ) + hf (tn , u(tn )) + h2 df dt
(tn , un ) + O(h3 )
2
h2
= u(tn ) + hf (tn , u(tn )) + h2 ∂t f (tn , un ) + 2
f (tn , un )∂u f (tn , un ) + O(h3 )
Donc par identification, on trouve le résultat.

Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 48 / 64
Schéma de Runge-Kutta d’ordre 2 (RK2)
En posant λ = ω2 , pour le schéma RK2, on obtient :

h h
 
φ(t, u; h) = (1 − λ)f (t, u) + λf t+ ,u + f (t, u)
2λ 2λ
1 1
Le tableau de Butcher 2λ 2λ
1−λ λ
Cas particuliers :
h hf (t,u)


1 Schéma d’Euler modifié : λ = 1, φ(t, u; h) = f t+ 2
,u + 2

(
k1 = f (tn , un ),
1 1
k2 = f (tn + h2 , un + h
2
k1 )) 2 2
0 1
un+1 = un + hk2

1 f (t,u)+f (t+h,u+hf (t,u))

2 Schéma de Heun : λ = 2
, φ(t, u; h) = 2


k1 = f (tn , un ),
 k2 = f (tn+1 , un + hk1 )) 1 1
h 1 1
 un+1 = un + (k1 + k2 ) 2 2
2

Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 49 / 64
Schéma de Runge-Kutta d’ordre 3 (RK3)

Un schéma de Runge-Kutta d’ordre 3 s’écrit

1
un+1 = un + h (k1 + 4k2 + k3 )
6

où 
 k1 = f (tn , un )
f tn + h2 , un + h2 k1


k =
 2
k3 = f (tn + h, un − hk1 + 2hk2 )
1 1
2 2
Le tableau de Butcher de RK3 1 −1 2
1/6 2/3 1/6

Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 50 / 64
Schéma de Runge-Kutta d’ordre 4 (RK4)

Le schéma classique de Runge-Kutta d’ordre 4 s’écrit

1
un+1 = un + h (k1 + 2k2 + 2k3 + k4 )
6

où 

 k1 = f (tn , un )
f tn + h2 , un + h2 k1



k2 =

k = f tn + h2 , un + h2 k2
 3


k4 = f (tn + h, un + hk3 )
1/2 1/2
1/2 0 1/2
Le tableau de Butcher de RK4
1 0 0 1
1/6 2/6 2/6 1/6

Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 51 / 64
Sommaire

1 Problème de Cauchy et Équations de premier ordre

2 Les systèmes du 1er ordre

3 Les équations différentielles d’ordre >1

4 Approximation numérique des ODE d’ordre 1

Analyse de la méthode d’Euler
Schéma de Heun
Schémas de Runge-Kutta

5 Schémas numériques à pas liés

Principe général
Méthodes d’ADAMS-BASHFORTH
Méthodes d’ADAMS-MOULTON

Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 52 / 64
Schémas numériques à pas liés
Principe général

Dans un schéma à un pas (ou pas séparé), le terme un+1 dépend explicitement uniquement du
terme un .
Pour les méthodes multi-pas (ou à pas liés), le terme un+1 dépend de plusieurs valeurs précédentes
de la solution (tk , uk ), k ≤ n.

Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 53 / 64
Schémas numériques à pas liés
Principe général

Dans un schéma à un pas (ou pas séparé), le terme un+1 dépend explicitement uniquement du
terme un .
Pour les méthodes multi-pas (ou à pas liés), le terme un+1 dépend de plusieurs valeurs précédentes
de la solution (tk , uk ), k ≤ n.

Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 53 / 64
Schémas numériques à pas liés
Principe général

Définition (Méthode à p + 1 pas)

Une méthode à p + 1 pas (p ≥ 0) est telle que ∀n ≥ p, le terme un+1 dépend directement de un−p , mais
d’aucun termes uk tel que k < n − p.

Exemples :
 Méthode explicite à 2 pas (Méthode du point milieu)

un+1 = un−1 + 2hfn , n≥1

u0 donné et u1 à déterminer.
 Méthode implicite à 2 pas (Méthode de Simpson)

h
un+1 = un−1 + [fn−1 + 4fn + fn+1 ], n≥1
3
u0 donné et u1 à déterminer.
Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 54 / 64
Schémas numériques à pas liés
Principe général

Les méthodes multi-pas linéaires à p + 1 pas (p ≥ 0) sont définies par

p
X p
X
un+1 = aj un−j + h bj fn−j + hb−1 fn+1 , n = p, p + 1, · · ·
j=0 j=0

les coefficients aj , bj caractérisent le schéma (on suppose que ap 6= 0 ou bp 6= 0)

u0 est donné, Mais les uk , k = 1, · · · , p doivent être initialisées (par exemple par une méthode à un
pas)
Si b−1 = 0 le schéma est explicite. Sinon, le schéma est implicite, il est nécessaire, dans ce cas, de
résoudre un problème de point fixe ou d’utiliser une technique de prédiction-correction

Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 55 / 64
Schémas numériques à pas liés
Principe général

Les méthodes multi-pas linéaires à p + 1 pas (p ≥ 0) sont définies par

p
X p
X
un+1 = aj un−j + h bj fn−j + hb−1 fn+1 , n = p, p + 1, · · ·
j=0 j=0

les coefficients aj , bj caractérisent le schéma (on suppose que ap 6= 0 ou bp 6= 0)

u0 est donné, Mais les uk , k = 1, · · · , p doivent être initialisées (par exemple par une méthode à un
pas)
Si b−1 = 0 le schéma est explicite. Sinon, le schéma est implicite, il est nécessaire, dans ce cas, de
résoudre un problème de point fixe ou d’utiliser une technique de prédiction-correction

Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 55 / 64
Schémas numériques à pas liés
Principe général

Les méthodes multi-pas linéaires à p + 1 pas (p ≥ 0) sont définies par

p
X p
X
un+1 = aj un−j + h bj fn−j + hb−1 fn+1 , n = p, p + 1, · · ·
j=0 j=0

les coefficients aj , bj caractérisent le schéma (on suppose que ap 6= 0 ou bp 6= 0)

u0 est donné, Mais les uk , k = 1, · · · , p doivent être initialisées (par exemple par une méthode à un
pas)
Si b−1 = 0 le schéma est explicite. Sinon, le schéma est implicite, il est nécessaire, dans ce cas, de
résoudre un problème de point fixe ou d’utiliser une technique de prédiction-correction

Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 55 / 64
Sommaire

1 Problème de Cauchy et Équations de premier ordre

2 Les systèmes du 1er ordre

3 Les équations différentielles d’ordre >1

4 Approximation numérique des ODE d’ordre 1

Analyse de la méthode d’Euler
Schéma de Heun
Schémas de Runge-Kutta

5 Schémas numériques à pas liés

Principe général
Méthodes d’ADAMS-BASHFORTH
Méthodes d’ADAMS-MOULTON

Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 56 / 64
Schémas numériques à pas liés
Méthodes d’ADAMS-BASHFORTH

La forme générale des méthodes d’ADAMS-BASHFORTH à p + 1 pas est

p
X
un+1 = un + h βi fn−i
i=0

Les βi sont choisis de sorte que l’ordre de la méthode soit maximal

Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 57 / 64
Schémas numériques à pas liés
Méthodes d’ADAMS-BASHFORTH

La forme générale des méthodes d’ADAMS-BASHFORTH à p + 1 pas est

p
X
un+1 = un + h βi fn−i
i=0

Les βi sont choisis de sorte que l’ordre de la méthode soit maximal

Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 57 / 64
Schémas numériques à pas liés
Méthodes d’ADAMS-BASHFORTH

La forme intégrale entre tn et tn+1 Z tn+1

u(tn+1 ) = u(tn ) + f (t, u(t))dt

tn

L’idée est d’approcher f (t, u(t)) dans l’intégrale par un polynôme d’interpolation (Ln,p ) aux points (tn−i )0≤i≤p

Z tn+1

u(tn+1 ) ≈ u(tn ) + Ln,p (t)dt

tn

L’intégrale est calculée analytiquement

Pour p = 0 (Ln,0 (t) = fn ), on retrouve le schéma à un pas d’Euler explicite (ordre 1)

Z tn+1

un+1 = un + Ln,0 (t)dt = un + hfn

tn

Pour p = 1, 2 et 3, on retrouve les schémas d’ADAMS- BASHFORTH à 2, 3 et 4 pas respectivement.

Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 58 / 64
Schémas numériques à pas liés
Méthodes d’ADAMS-BASHFORTH

Méthode d’ADAMS-BASHFORTH à 2 pas (ordre 2)

h
un+1 = un + (3fn − fn−1 )
2

Méthode d’ADAMS-BASHFORTH à 3 pas (ordre 3)

h
un+1 = un + (23fn − 16fn−1 + 5fn−2 )
12

Méthode d’ADAMS-BASHFORTH à 4 pas (ordre 4)

h
un+1 = un + (55fn − 59fn−1 + 37fn−2 − 9fn−3 )
24

Ces schémas sont explicites et leur ordre correspond au nombre de pas

Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 59 / 64
Sommaire

1 Problème de Cauchy et Équations de premier ordre

2 Les systèmes du 1er ordre

3 Les équations différentielles d’ordre >1

4 Approximation numérique des ODE d’ordre 1

Analyse de la méthode d’Euler
Schéma de Heun
Schémas de Runge-Kutta

5 Schémas numériques à pas liés

Principe général
Méthodes d’ADAMS-BASHFORTH
Méthodes d’ADAMS-MOULTON

Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 60 / 64
Schémas numériques à pas liés
Méthodes d’ADAMS-MOULTON

La forme générale des méthodes d’ADAMS-MOULTON à p + 1 pas

p
X
un+1 = un + h βi fn+1−i , avec β0 6= 0
i=0

La fonction f (t, u(t)) est approchée, dans ce cas, par un polynôme d’interpolation (L∗n,p ) aux points
(tn+1−i )0≤i≤p
Z tn+1
u(tn+1 ) ≈ u(tn ) + L∗n,p (t)dt
tn

Pour p = 0 (L∗n,0 (t) = fn+1 ), on retrouve le schéma à un pas d’Euler implicite (ordre 1)

Z tn+1
un+1 = un + L∗n,0 (t)dt = un + hfn+1
tn

Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 61 / 64
Schémas numériques à pas liés
Méthodes d’ADAMS-MOULTON

Pour p = 1, L∗n,1 (t) est un polynôme d’interpolation de f aux points (tn , fn ) et (tn+1 , fn+1 ), on retrouve le schéma
d’ADAMS-MOULTON à un pas (ordre 2)

h
un+1 = un + (fn+1 + fn ) ,
2

Le schéma d’ADAMS-MOULTON à 2 pas (ordre 3)

h
un+1 = un + (5fn+1 + 8fn − fn−1 ) ,
12

Le schéma d’ADAMS-MOULTON à 3 pas (ordre 4)

h
un+1 = un + (9fn+1 + 19fn − 5fn−1 + fn−2 ) ,
24

Ces schémas sont implicites et leur ordre correspond au nombre de pas plus un.

Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 62 / 64
Références

André Fortin, "Analyse numérique pour ingénieurs", Montréal : Presses internationales

Polytechnique, 2011.
Gabriel Nagy, "Ordinary differential equations" Mathematics Department, Michigan State
University, East Lansing, MI, 48824.
Alfio Maria Quarteroni,Fausto Saleri,Paola Gervasio, "Calcul Scientifique: Cours, exercices corrigés
et illustrations en Matlab et Octave" Springer-Verlag Italia 2008

Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 63 / 64
 Fin Chapitre EDO

Fassi Fihri, El Ossmani, Taakili (ENSAM) Résolution numérique November 21, 2023 64 / 64