Chapitre 5: Méthodes à pas pour la résolution

numérique des EDO

Table des matières

1 Méthodes à un pas : 1
1.1 Généralités : . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Consistance, stabilité, convergence et ordre d’une méthode . . . . . . . . . . . . . . . . . . 2

2 Les méthodes d’Euler 5

2.1 Méthode d’Euler explicite . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.2 Méthode d’Euler implicite . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

3 Les méthodes de Runge-Kutta 7

3.1 Stabilité des méthodes de Runge-Kutta explicites . . . . . . . . . . . . . . . . . . . . . . . 8
3.2 Exemples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

4 Mise en œuvre et difficultés 10

4.1 Problème bien posé mathématiquement et numériquement . . . . . . . . . . . . . . . . . . 10
4.2 Conditionnement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

L’objectif de ce chapitre est de décrire des méthodes permettant de résoudre numériquement un problème
de Cauchy de condition initiale y (t0 ) = y0 pour une équation différentielle

y 0 = f (t, y), (1)

où f : [t0 , t0 + T ] × R → R est une fonction suffisamment régulière. Le cas traité ici est celui des équations
unidimensionnelles afin de simplifier les notations ; le cas des systèmes dans Rm est tout à fait identique, à
condition de considérer y comme une variable vectorielle et f comme une fonction vectorielle dans les algorithmes
qui vont être décrits.
Considérons une subdivision t0 < t1 < . . . < tN = t0 + T régulière de [t0 , t0 + T ] d’un pas noté h, on cherche
à déterminer des valeurs approchées y0 , y1 , . . . , yN des valeurs y (tn ) prises par la solution exacte y. On appelle
méthode à un pas une méthode permettant de calculer yn+1 à partir de la seule valeur antérieure yn .

1 Méthodes à un pas :
1.1 Généralités :
Définition 1.1. Les méthodes à un pas sont les méthodes de résolution numérique qui peuvent s’écrire sous la
forme
yn+1 = yn + hn Φ (tn , yn , h) , 0 ≤ n < N,
où Φ : [t0 , t0 + T ] × R × R → R est une fonction que l’on supposera continue.
Exemple. La méthode d’Euler est la méthode à un pas associée à la fonction Φ(t, y, h) = f (t, y), et définie par
la formule de récurrence yn+1 = yn + hn f (tn , yn ).
Définition 1.2. L’erreur de consistance en relative à une solution exacte y est l’erreur

en = z (tn+1 ) − yn+1 , 0≤n<N

1
produite par application de l’algorithme yn+1 = yn + hn Φ (tn , yn , h) à partir de la valeur yn = y (tn ). Autrement
dit, cette erreur mesure l’écart entre la valeur exacte y (tn+1 ) au temps tn+1 , et la valeur approchée yn+1 issue
de la valeur yn = y (tn ) prise comme valeur initiale au temps tn (une seule étape de l’algorithme est donc mise
en jeu). En termes de la fonction Φ, on a

en = y (tn+1 ) − y (tn ) − hΦ (tn , y (tn ) , h) .

1.2 Consistance, stabilité, convergence et ordre d’une méthode

Définition 1.3. On dit qu une methode a un pas est consistante avec problème 1 si pour toute solution y de
celui-ci, on a
NX−1
lim |en | = 0
N →+∞
n=0
où
en = y (tn+1 ) − y (tn ) − ∆T Φ (tn , y (tn ) , ∆T )
est l’erreur de consistance.
Remarque. La consistance consiste en réalité à vérifier que la méthode employée est en accord (ou consistante)
avec le problème initial.
Théorème 1.4. une méthode explicite à un pas est consistante si et seulement si

∀(t, y) ∈ [t0 , t0 + T ] × R, Φ(t, y, 0) = f (t, y).

Démonstration. Soit εn = y (tn+1 ) − y (tn ) − ∆T Φ (tn , y (tn ) , ∆T ) (0 ≤ n ≤ N − 1). Par le théorème des
acroissements finis, on sait qu’il existe cn ∈] tn , tn+1 [ tel que

εn = ∆T f (cn , y (cn )) − ∆T Φ (tn , y (tn ) , ∆T ) = ∆T (αn + βn )

avec 
αn = f (cn , y (cn )) − Φ (cn , y (cn ) , 0) ,
βn = Φ (cn , y (cn ) , 0) − Φ (tn , y (tn ) , ∆T ) .
Comme la fonction Φ̃ : [t0 , t0 + T ] × [0, T ] → R, (t, h) 7→ Φ(t, y(t), h) est continue sur l’ensemble compact
[t0 , t0 + T ] × [0, T ], donc uniformément continue, on peut affirmer que

∀ε > 0, ∃N0 ∈ N∗ , ∀N ≥ N0 , ∀n ∈ {0, . . . , N − 1}, |βn | ≤ ε.

Ainsi, pour N ≥ N0 ,
N
X −1 N
X −1 N
X −1
|en | − h|αn | ≤ h |βn | ≤ εT
n=0 n=0 n=0
De plus, par définition de l’intégrale de Riemman,
N
X −1 Z t0 +T
lim h |αn | = |f (t, y(t)) − Φ(t, y(t), 0)|dt.
N →+∞ t0
n=0

Au vu de la définition et des remarques précédentes, une condition nécessaire et suffisante de consistance de la

méthode étudiée est donc que pour toute solution y de (1) on ait

∀t ∈ [t0 , t0 + T ] , f (t, y(t)) = Φ(t, y(t), 0).

Or, pour tout couple (t∗ , y ∗ ) ∈ [t0 , t0 + T ] × R, il existe par le théorème de Cauchy-Lipschitz une (unique)
solution y définie sur [t0 , t0 + T ] du problème de Cauchy suivant :
 0
y (t) = f (t, y(t)), t ∈ [t0 , t0 + T ] ,
y (t∗ ) = y ∗ .

En écrivant la relation précédente en t∗ pour cette fonction, également solution d’un problème de type (1),
on obtient :
Φ (t∗ , y ∗ , 0) = f (t∗ , y ∗ )
ce qui est bien le résultat recherché.

2
Définition 1.5. On dit qu’une méthode à un pas est stable si pour tout entier N ∈ N∗ et pour toute famille
(zn )0≤n≤N solution de la relation de récurrence perturbée par le réel εn :

z0 ∈ R fixé,
zn+1 = zn + hΦ (tn , zn , h) + εn 0 ≤ n ≤ N − 1,
on a la relation
N
X −1
max |zn − yn | ≤ M |z0 − y0 | + M 0 |εn |
0≤n≤N −1
n=0
où M et M 0 sont des constantes indépendantes de y0 , z0 et N .
Le théorème suivant fournit une condition suffisante de stabilité d’un schéma à un pas.
Théorème 1.6. Pour qu’une méthode explicite à un pas soit stable, il suffit que la fonction Φ soit Lipschitzienne
en y, à savoir qu’il existe Λ > 0 tel que
∀ (t, y1 , y2 , h) ∈ [t0 , t0 + T ] × R2 × [0, T ], |Φ (t, y2 , h) − φ (t, y1 , h)| ≤ Λ |y2 − y1 |
On peut alors prendre comme constantes de stabilité M = M 0 = eΛT .
Avant de passer à la preuve du théorème, citons en premier le lemme technique :
2 N


Lemme 1.7. Soit (zn )n∈N une suite réelle et (hn , αn )n∈N ∈ R+ tels que :

z0 ≥ 0
zn+1 ≤ (1 + hn ) zn + αn
Démonstration. En conservant les notations de la définition de la stabilité, on peut écrire sous les hypothèses
du théorème :
|yn+1 − zn+1 | = |yn − zn + h (Φ (tn , yn , h) − Φ (tn , zn , h)) − εn |
≤ (1 + Λh) |yn − zn | + |εn |
Grâce au Lemme 1.7, on a alors
n−1
X
|yn − zn | ≤ enΛh |y0 − z0 | + e(n−1−i)Λh |εi |
i=0
N
X −1
ΛT ΛT
≤e |y0 − z0 | + e |εi |
i=0

ce qui clôt la démonstration.

Définition 1.8. On dit qu’une méthode à un pas est convergente si pour toute solution y du problème (1)
lim max |y (tn ) − yn | = 0.
N →+∞ 0≤n≤N
y0 →y(t0 )

Théorème 1.9. Si une méthode explicite à un pas est stable et consistante, alors elle est convergente.
Démonstration. On pose εn = y (tn+1 ) − y (tn ) − hΦ (tn , y (tn ) , h) (0 ≤ n ≤ N − 1). La famille (zn )0≤n≤N −1
solution de la relation de récurrence perturbée

z0 = y (t0 ) ,
zn+1 = zn + hΦ (tn , zn , h) + εn (0 ≤ n ≤ N − 1)
vérifie facilement par construction même :
∀n ∈ {0, . . . , N }, zn = y (tn ) .
La méthode considérée étant stable et consistante, il vient respectivement
N
X −1
max |yn − y (tn )| ≤ M |y0 − y (t0 )| + M 0 |εn |
0≤n≤N −1
n=0

et
N
X −1
lim |εn | = 0.
N →+∞
n=0
En regroupant ces deux résultats, on a bien démontré la convergence de la méthode.

3
Définition 1.10. On dit qu’une méthode à un pas est d’ordre p ∈ N∗ si Φ et f sont p-fois continûment
différentiables sur leur ensemble de définition et si pour toute solution y du problème 1, il existe une constante
C > 0 indépendante de N telle que
N
X −1
∀N ∈ N∗ , |en | ≤ C(h)p
n=0

où
en = y (tn+1 ) − y (tn ) − hΦ (tn , y (tn ) , h) (0 ≤ n ≤ N − 1).
Théorème 1.11. Une méthode explicite à un pas est d’ordre p ∈ N∗ si et seulement si Φ et f sont p-fois
continûment différentiables sur leur ensemble de définition et si
∂l 1 [l]
∀l ∈ {0, . . . , p − 1}, ∀(t, y) ∈ [t0 , t0 + T ] × R, l
Φ(t, y, 0) = f (t, y), (2)
∂h l+1

où f [l] (t, y) désigne la l-ième dérivée totale de f suivant les caractéristiques du problème (1) :
 [0]
f (t, y) = f (t, y),
f [k+1] (t, y) = ∂t f [k] (t, y) + ∇y f [k] (t, y) · f (t, y) (0 ≤ k ≤ p − 1).

Démonstration. Condition suffisante : On remarque tout d’abord que si f ∈ C p ([t0 , t0 + T ] × R, R), alors
toute solution y de (1) appartient à C p+1 ([t0 , t0 + T ] , R) et vérifie

∀l ∈ {0, . . . , p}, ∀t ∈ [t0 , t0 + T ] , y (l+1) (t) = f [l] (t, y(t)).

On écrit alors la relation de Taylor-Lagrange à l’ordre p pour la fonction h 7→ Φ (tn , y (tn ) , h) entre 0 et h :
p−1
X (h)l ∂ l (h)p ∂ p
Φ (tn , y (tn ) , h) = l
Φ (tn , y (tn ) , 0) + Φ (tn , y (tn ) , λn )
l! ∂h p! ∂hp
l=0

(λn ∈] 0, h[) et à l’ordre p + 1 pour la fonction y entre tn et tn+1 :

p
X (h)k (h)p+1 (p+1)
y (tn+1 ) − y (tn ) = y (k) (tn ) + y (cn ) (cn ∈] tn , tn+1 [)
k! (p + 1)!
k=1
p
X (h)k [k−1] (h)p+1 (p+1)
= f (tn , y (tn )) + y (cn ) .
k! (p + 1)!
k=1

En soustrayant ces deux égalités après un changement d’indices (k = l + 1), il vient

εn = y (tn+1 ) − y (tn ) − hΦ (tn , y (tn ) , h)

p−1
(h)l+1 f [l] (tn , y (tn )) ∂ l Φ (tn , y (tn ) , 0)
X  
= −
l! l+1 ∂hl
l=0
(h)p+1 y (p+1) (cn ) ∂p
 
+ − Φ (tn , y (tn ) , λn )
p! p+1 ∂hp
(h)p+1 y (p+1) (cn ) ∂p
 
= − Φ (tn , y (tn ) , λn )
p! p+1 ∂hp

grâce à la relation 2. Ainsi, avec la régularité supposée de Φ et f

∃C > 0, ∀n ∈ {0, . . . , N − 1}, |εn | ≤ C(h)p+1 .

En sommant ces N inégalités, on trouve bien que la méthode est d’ordre p.

Condition nécessaire : on raisonne par l’absurde en supposant que la méthode est d’ordre p et que la
relation (4) est violée à partir d’un certain l ∈ {0, . . . , p − 1}. En écrivant les mêmes égalités de Taylor-Lagrange
que précédemment mais cette fois à l’ordre l + 1 et l + 2 respectivement, il vient

(h)l+1 f [l] (tn , y (tn )) ∂l

 
Φ (tn , y (tn ) , 0) + O (h)l+2


εn = − l
l! l+1 ∂h

4
puis en sommant
N −1 N −1
X εn 1 X
= h |ψ (tn )| + O(h)
n=0
hl l! n=0
[l] l
où ψ désigne la fonction continue ψ : [t0 , t0 + T ] → R, t 7→ f (t,y(t))
l+1
∂
− ∂hl Φ(t, y(t), 0). On déduit alors par

passage à la limite dans les deux membres de la dernière égalité lorsque h tend vers 0 que
Z t0 +T
1
0= |ψ(s)|ds
l! t0

ce qui implique immédiatement que pour toute fonction y solution de (1)

f [l] (t, y(t)) ∂l

∀t ∈ [t0 , t0 + T ] , ψ(t) = − Φ(t, y(t), 0) = 0.
l+1 ∂hl
En raisonnant comme dans la fin de la démonstration du Théorème 1.4 , on en déduit facilement que

∂l f [l] (t∗ , y ∗ )
∀ (t∗ , y ∗ ) ∈ [t0 , t0 + T ] × R Φ (t ∗ ∗
, y , 0) =
∂hl l+1
ce qui achève la démonstration par absurde.

2 Les méthodes d’Euler

2.1 Méthode d’Euler explicite
La méthode d’Euler explicite est un cas particulier de méthode explicite à un pas correspondant au choix
de Φ suivant :
∀(t, y, h) ∈ [t0 , t0 + T ] × Rm × [0, T ], Φ(t, y, h) = f (t, y).
Définition 2.1. On appelle méthode de résolution approchée d’Euler explicite (à pas de temps constant) du
problème (1), la construction pour tout N ∈ N∗ de la suite finie (yn )0≤n≤N telle que

y0 ∈ Rm fixé,

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

Dans toute la suite, on désignera par en l’erreur commise en remplaçant y (tn ) par yn :

en = yn − y (tn )
Rt
Remarque. La méthode d’Euler revient en fait à approcher l’intégrale tnn+1 y 0 (t)dt par la valeur (tn+1 − tn ) f (tn , yn ),
c’est à dire à utiliser la méthode de quadrature élémentaire des rectangles à gauche.
Théorème 2.2. La méthode d’Euler explicite est consistante, stable, convergente et d’ordre 1 (si f est C 1 ).
De plus, on a l’inégalité suivante valable pour tout n ∈ {0, . . . , N } et toute solution y du problème (1) :

enLh − 1
ken k ≤ ω (h, y 0 ) + enLh ke0 k
L
où ω désigne le module de continuité d’une fonction :

ω(δ, g) = max kg (t0 ) − g(t)k

{(t,t0 )∈[t0 ,t0 +T ]2 /|t0 −t|≤δ}

et L la constante de Lipschitz de f par rapport à la seconde variable.

Avant de passer à la preuve rappelons en premier le lemme suivant :
Lemme 2.3. Soit (zn )n∈N une suite réelle telle que

z0 ≥ 0
zn+1 ≤ Azn + B (n ∈ N)

5
où A et B sont deux réels tels que A ≥ 1 et B ≥ 0. On a l’inégalité suivante valable pour tout n ∈ N :

enQ − 1
zn ≤ enQ z0 + B
Q
avec Q = A − 1 et la convention
enQ − 1
=n lorsque Q = 0.
Q
Démonstration. Les propriétés de consistance, de stabilité, de convergence et l’ordre de la méthode d’Euler se
déduisent immédiatement de l’étude générale du paragraphe précédent.
Pour démontrer l’estimation de l’erreur dans la seconde partie du théorème, on considère une fonction y
solution du problème (1). On a :
Z tn+1
y (tn+1 ) = y (tn ) + y 0 (t)dt
tn
Z tn+1
= y (tn ) + hf (tn , y (tn )) + (y 0 (t) − y 0 (tn )) dt
tn
= y (tn ) + hf (tn , y (tn )) + εn

(5) de la suite des approximations :

yn+1 = yn + hf (tn , yn )
En soustrayant ces deux égalités et en utilisant la définition du module de continuité, on obtient

ky (tn+1 ) − yn+1 k ≤ ky (tn ) − yn k + Lh ky (tn ) − yn k + hω (h, y 0 ) .

On se trouve alors dans les conditions d’application du Lemme 2.3. On conclut alors que :

enLh − 1
ky (tn ) − yn k = ken k ≤ enLh ke0 k + hω (h, y 0 ) .
Lh
.
Corollaire 2.4. Soit la suite de fonctions (uN )N ∈N définies sur [t0 , t0 + T ] par interpolation affine entre les
points (tn )0≤n≤N pour lesquels uN (tn ) = yn où la famille (yn )0≤n≤N est construite suivant la méthode d’Euler
avec y0 ∈ Rm fixé. Alors uN converge dans C 0 ([t0 , t0 + T ] , Rm ) vers l’unique solution y de (1) telle que
y (t0 ) = y0 .

Démonstration. Soit N fixé et t ∈ [t0 , t0 + T [ . On note n l’unique entier compris entre 0 et N − 1 tel que
t ∈ [tn , tn+1 [.On a alors l’estimation suivante :

kuN (t) − y(t)k ≤ kuN (t) − yn k + kyn − y (tn )k + ky (tn ) − y(t)k

≤ kuN (t) − yn k + ken k + ω(h, y).

uN étant affine sur [tn , tn+1 ], on a

kuN (t) − yn k ≤ kyn+1 − yn k = h kf (tn , yn )k

Il suffit donc pour conclure à l’uniforme convergence de uN vers y de montrer que yn reste borné indépen-
dament de n et N . Pour cela, on remarque d’abord que

kyn+1 k ≤ kyn k + h kf (tn , yn )k ≤ kyn k + h (kf (tn , 0)k + L kyn k)

≤ (1 + hL) kyn k + h max kf (t, 0)k
t∈[t0 ,t0 +T ]

et on utilise ensuite le Lemme 1 :

eT L − 1
kyn k ≤ eT L ky0 k + max kf (t, 0)k
L t∈[t0 ,t0 +T ]

6
2.2 Méthode d’Euler implicite
Il existe une version implicite de la méthode d’Euler qui peut s’avérer à l’usage mieux conditionnée que la
version explicite (voir exemple ultérieur).
Définition 2.5. On appelle méthode de résolution approchée d’Euler implicite de l’équation (1), la construction
valable pour tout N ∈ N∗ tel que hL < 1 d’une suite finie (yn )0≤n≤N vérifiant

y0 ∈ Rm fixé,

yn+1 = yn + hf (tn , yn+1 ) 0≤n≤N −1

Théorème 2.6. La méthode d’Euler implicite est un cas particulier de méthode explicite à un pas. Elle est
L
consistante, stable, convergente et d’ordre 1 (si f est C 1 ). De plus, en notant L1 = 1−Lh , on a une estimation
de l’erreur commise :
eL1 (tn −t0 ) − 1
∀n ∈ {0, . . . , N }, ken k ≤ ω (δN , y 0 ) + eL1 (tn −t0 ) ke0 k
L1
si f est seulement continue. Si f est continûment dérivable, on a de plus
Z tn
δN
ken k ≤ eL1 (tn −s) ky 00 (s)k ds + eL1 (tn −t0 ) ke0 k .
1 − LδN t0

Démonstration. On se limite ici à la démonstration des deux estimations d’erreur. Pour cela, on écrit

en+1 = en + h (f (tn+1 , yn+1 ) − f (tn+1 , y (tn+1 )) + εn , n ∈ {0, . . . , N − 1}

avec Z tn+1 
εn = y (tn ) − y (tn+1 ) + hf (tn+1 , y (tn+1 )) = (y 0 t) − y 0 (tn+1 ) dt.
tn

En particulier
ken+1 k ≤ ken k + hL ken+1 k + kεn k .
On remarque alors que  
−1 L
0 < (1 − hL) ≤1+h = 1 + hL1
1 − Lh
et en utilisant à nouveau le Lemme 1.7, il vient
n−1
X
L1 (tn −t0 )
ken k ≤ e ke0 k + eL1 (tn −ti+1 ) kεi k
i=0

ce qui permet ensuite de conclure.

3 Les méthodes de Runge-Kutta

La famille des méthodes de Runge-Kutta explicites construites dans ce paragraphe (et à laquelle appartient
la méthode d’Euler explicite) est un exemple important de méthodes explicites à un pas.
On suppose pour simplifier que m = 1 et que le pas de temps est constant :
T

h= N ,
tn = t0 + nh (0 ≤ n ≤ N ).

Définition 3.1. Soit q ∈ N∗ et (ai,j )1≤j<i≤q , (bi )1≤i≤q et (ci )1≤i≤q trois familles de coefficients réels (avec
ci ∈ [0, 1] ). On appelle méthode de résolution approchée de Runge-Kutta explicite à un pas (et à pas de temps
constant) du problème (1), la construction pour tout N ∈ N∗ de la suite finie (yn )0≤n≤N telle que

y0 ∈ Rm fixé,


 t = t + c ∆T, 1 ≤ j ≤ q,

n,j n j
Pi−1
y n,i = yn + ∆T ai,j f (tn,j , yn,j ) , 1 ≤ i ≤ q,

 Pj=1
q
yn+1 = yn + ∆T j=1 bj f (tn,j , yn,j ) , 0 ≤ n ≤ N − 1.


7
 On représente conventionnellement une méthode de Runge-Kutta par le tableau de ses coefficients rangés
de la manière suivante :

c1 0
c2 a2,1 0
..
c3 a3,1 a3,2 .
.. .. .. .. ..
. . . . .
cq aq,1 aq,2 ... aq,q−1 0
b1 b2 ... bq−1 bq
Remarque. Pour mieux comprendre la définition des méthodes de Runge-Kutta, les formules donner sont à
rapprocher de celles valables pour toute solution exacte y du problème (1) :
( Rc
y (tn,i ) = y (tn ) + ∆T 0 i f (tn + uh, y (tn + uh)) du, 1 ≤ i ≤ q,
R1
y (tn+1 ) = y (tn ) + h 0 f (tn + uh, y (tn + uh)) du 0 ≤ n ≤ N − 1.

Le principe des méthodes de Runge-Kutta consiste donc à approcher successivement par une méthode de
quadrature, y (tn,i ) pour tout i ∈ {1, . . . , q}, puis y (tn+1 ) à l’aide des précédentes valeurs calculées.
Ainsi, les familles de coefficients (ai,j )1≤j<i≤q et (bi )1≤i≤q sont associées aux méthodes de quadrature :
( Rc Pi−1
i
g(t)dt ' j=1 ai,j g (cj ) , 1 ≤ i ≤ q,
R01 Pq
0
g(t)dt ' j=1 bj g (cj ) .

On obtient ainsi deux premières conditions (qui seront supposées vérifiées dans toutes la suite) sur les
familles de coefficients pour que les méthodes de quadrature soient au moins d’ordre 0 (c’est à dire exactes sur
les constantes) :
i−1
X
∀i ∈ {1, . . . , q}, ci = ai,j
j=1

et
q
X
1= bj .
j=1

En particulier, ceci impose que c1 = 0 (et aussi par conséquent tn,1 = tn et yn,1 = yn ).

3.1 Stabilité des méthodes de Runge-Kutta explicites

Les méthodes de Runge-Kutta explicites sont un cas particulier de méthodes explicites à un pas correspondant
à une fonction Φ égale à :
Xq
Φ(t, y, h) = bj f (t + cj h, yj )
j=1

où la famille (yi )1≤i≤q est définie pout tout (t, y, h) par

i−1
X
yi = y + h ai,j f (t + cj h, yj ) , 1 ≤ i ≤ q.
j=1

On a alors le résultat suivant :

Proposition 3.2. Les méthodes de Runge-Kutta explicites sont stables.
Démonstration. Il suffit de montrer que Φ est Lipschitzienne par rapport à la variable y. Pour cela, soit
i−1
X
α = max |ai,j | .
1≤i≤q
j=1

Si (yi )1≤i≤q et (zi )1≤i≤q sont deux familles construites suivant la formule (11) (à partir de y et z respecti-
vement), on montre aisément par récurrence l’inégalité
|yi − zi | ≤ 1 + (αLh) + . . . + (αLh)i−1 |y − z|

8
où L désigne la constante de Lispschitz de f par rapport à sa deuxième variable. On peut alors écrire pour tout
(t, h) ∈ [t0 , t0 + T ] × [0, ∆T ] :
q
X
|Φ(t, y, h) − Φ(t, z, h)| ≤ |bj | L |yj − zj | ≤ Λ|y − z|
j=1

où
q
X
|bj | 1 + (αLh) + . . . + (αLh)j−1 ,


Λ=L
j=1

ce qui achève la démonstration.

Théorème 3.3. Toute méthode de Runge-Kutta explicite vérifiant :

i−1
X q
X
∀i ∈ {1, . . . , q}, ci = ai,j et 1 = bj , (3)
j=1 j=1

est convergente et d’ordre 1 (si f est C 1 ). Une condition nécessaire et suffisante pour qu’elle soit d’ordre 2 (si
f est C 2 ) est que ses coefficients vérifient en outre :
q
X 1
bj cj = . (4)
j=1
2

Démonstration. En reprenant l’expression de la fonction Φ donnée au sous-paragraphe précédent, on s’aperçoit

aisément que
Xq
Φ(t, y, 0) = bj f (t, y) = f (t, y)
j=1

Grâce à l’étude générale précédemment effectuée, on en déduit que toute méthode de Runge-Kutta vérifiant (3)
est convergente car consistante et stable. De plus son ordre est au moins égal à 1 (si f est C 1 ). Afin d’améliorer
éventuellement cette valeur (si f est C 2 ), on calcule
q  
∂Φ X ∂yi
(t, y, h) = bj cj ∂t f (t + cj h, yj ) + ∂y f (t + cj h, yj )
∂h j=1
∂h

avec
j−1 j−1  
∂yi X X ∂yk
= aj,k f (t + ck h, yk ) + h aj,k ck ∂t f (t + ck h, yk ) + ∂y f (t + ck h, yk ) .
∂h ∂h
k=1 k=1

En particulier,
j−1
∂yi X
(t, y, 0) = aj,k f (t, y) = cj f (t, y)
∂h
k=1

et  
q q
∂Φ X X
(t, y, 0) = bj cj (∂t f (t, y) + ∂y f (t, y)f (t, y)) =  bj cj  f [1] (t, y).
∂h j=1 j=1

On obtient donc bien la condition 4 nécessaire et suffisante pour qu’une méthode de Runge-Kutta soit d’ordre
2.

3.2 Exemples
(i) q = 1 : on a nécessairement b1 = 1 et la méthode se réduit à :

yn+1 = yn + (tn+1 − tn ) f (tn , yn ) .

9
 On retouve la méthode d’Euler explicite étudiée au paragraphe précédent.
(ii) q = 2 : Pour tout α ∈ [0, 1], on construit des méthodes de Runge-Kutta d’ordre 2 avec le tableau

0 0
α α 0
1 1
1− 2α 2α

En particulier, lorsque α = 1, la méthode (dite de Heun ou d’Euler Cauchy) s’écrit

∆T
yn+1 = yn + [f (tn , yn ) + f (tn + ∆T, yn + ∆T f (tn , yn ))]
2
et lorsque α = 12 , la méthode (dite d’Euler modifié ou du point milieu) s’écrit
 
∆T ∆T
yn+1 = yn + ∆T f tn + , yn + f (tn , yn ) .
2 2

(iii) q = 4 : un exemple de méthode de Runge-Kutta fréquemment utilisé dans la pratique est le suivant :

0
1 1
2 2
1 1
2 0 2
1 0 0 1
1 2 2 1
6 6 6 6

On peut montrer que cette méthode est d’ordre 4.

4 Mise en œuvre et difficultés

Beaucoup de problèmes applicatifs se ramènent au terme de leur modélisation à la résolution d’un système
d’équations différentielles. Si on désire simuler celui-ci numériquement avec les méthodes exposées précédem-
ment, il est indispensable au préalable de s’assurer que le problème est bien posé.

4.1 Problème bien posé mathématiquement et numériquement

La notion de problème bien (ou mal) posé peut s’entendre à deux niveaux. Le premier niveau est purement
théorique : on dit qu’un problème de Cauchy est mathématiquement bien posé si la solution est unique et
dépend continûment de la donnée initiale. Par exemple, le problème
 0 p
y (t) = 2 |y(t)|, t ∈ R+ ,
y(0) = 0

est mathématiquement mal posé car il existe une infinité de solutions différentes de la solution nulle : en effet,
pour tout a ≥ 0 (
0 si t ∈ [0, a]
y(t) =
(t − a)2 si t > a
est solution. Le second niveau est d’ordre numérique : on dit qu’un problème de
Cauchy est numériquement bien posé s’il est bien posé mathématiquement et si la dépendance en la donnée
initiale est numériquement contrôlable : par exemple, le problème
 0 p
y (t) = 2 |y(t)|, t ∈ R+
y(0) = 0

est mathématiquement mal posé car il existe une infinité de solutions différentes de la solution nulle : en effet,
pour tout a ≥ 0 
0 si t ∈ [0, a],
y(t) =
(t − a)2 si t > a
est solution. Le second niveau est d’ordre numérique : on dit qu’un problème de

10
 Cauchy est numériquement bien posé s’il est bien posé mathématiquement et si la dépendance en la donnée
initiale est numériquement contrôlable : par exemple, le problème
 0
y (t) = 3y(t) − 1, t ∈ [0, 30]
y(0) = 31

est mathématiquement bien posé mais numériquement mal posé : en effet, il possède bien une unique solution
y(t) = 13 , mais la solution ỹ(t) obtenue en perturbant légèrement la donnée initiale, ỹ(0) = 31 + ε, est égale à

1
ỹ(t) = + εe3t
3
Par exemple, si ε = 10−17 (précision machine), on a ỹ(30) − y(30) ' 1022 .
De manière manifeste, aucune simulation numérique sur cette équation ne pourra donc être accomplie. Une
façon de s’assurer qu’un problème est numériquement bien posé pourra consister à déterminer une estimation
a priori de la solution en fonction des données initiales.

4.2 Conditionnement
Lorsqu’un problème de Cauchy est numériquement bien posé, il est encore essentiel avant toute simulation de
s’assurer qu’il est bien conditionné pour la méthode numérique envisagée, c’està dire que celle-ci peut approcher
de manière satisfaisante la solution exacte avec un pas de temps raisonnable. Soit par exemple (voir [Dem]) le
problème de Cauchy  0
y (t) = −150y(t) − 30, t ∈ [0, 1],
y(0) = 51
Il est mathématiquement bien posé (solution y(t) = 15 ) et aussi numériquement car la solution ỹ(t) corres-
pondant à une donnée initiale ỹ(0) = 15 + ε est égale à :

1
ỹ(t) = + εe−150t .
5
La méthode d’Euler explicite à pas constant appliquée à ce problème donne :
 
1 150 30
yn+1 = yn + (−150yn + 30) = 1 − yn + , n ∈ {0, . . . , N − 1}
N N N
soit   
1 150 1
yn+1 − = 1− yn −
5 N 5
puis enfin  n  
1 150 1
yn − = 1− y0 − .
5 N 5
En prenant N = 50 et y0 = 15 + ε, on obtient ainsi yn = 15 + (−2)n ε, approximation totalement aberrante
(par exemple, y50 ' 51 + 1015 ε ).
La condition sur le pas de temps pour obtenir une bonne approximation de la solution exacte est en fait
N > 75. On peut donc légitimement considérer que le problème est mal conditionné pour la méthode d’Euler
explicite. Par contre, le lecteur pourra vérifier qu’il est bien conditionné pour la méthode d’Euler implicite.

4.3 Problèmes raides

Il existe certains problèmes numériquement bien posés mais mal conditionnés pour toutes les méthodes
usuelles. On parle alors de problèmes raides. Ceux- ci interviennent dans de nombreux domaines applicatifs et
doivent faire l’objet d’études préalables complexes : mécanique (problème à trois corps ), cinétique chimique
(réactions très rapides, réactions nucléaires, photodissociation), etc... Une simulation numérique peut éventuel-
lement s’adapter à ces circonstances particulières : localement, un pas de temps très faible devra être adopté. Les
solveurs incorporés dans les logiciels Scilab (fonction ode) ou Matlab ont recours à des stratégies d’adaptation
de pas de temps permettant, dans une certaine mesure, de parvenir à un résultat satisfaisant.

11