Résolution numérique des équations

différentielles

H. Allouche

Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche

Rappels

Soit I0 un intevalle de R et soit

F : I0 × R→R
(t, y ) 7→ F (t, y )

une fonction à deux variables

On appelle équation différentielle de 1er ordre tout
equation de la forme:

y 0 (t) = F (t, y (t)), t ∈ I0 . (1)

La résolution consiste à touver une fonction y (t) continue

et dérivable sur I0 vérifiant (1).
Lorsque la fonction F est linéaire l’équation est dite
linéaire (que l’on sait résoudre analytiquement).
Lorque la fonction F est non linéaire l’équation est dite
non linéaire (Il y a très peu de type d’équation non
linéaire que l’on sait résoudre analytiquement).
Il s’agit de trouver une solution approchée de cette
fonction "inconnue".
Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche
On appelle problème de Cauchy :
 0
y (t) = F (t, y (t)), t ∈ I0 = [t0 , t0 + T ]
(2)
y (t0 ) = µ donnée

(appelé aussi problème avec conditions initiales)

Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche

Existence et unicité

Défintion
Soit U un ouvert de R R
× m . F une fonction définie de U
R
dans , elle est dite localement lipshitizienne par
rapport 2ème variable si ∀K = [a, b] × B, où B est une boule
R
fermée de m avec K ⊂ U alors il existe c > 0 tel que:

kF (t, x) − F (t, y )k ≤ ckx − y k, ∀(x, y ) ∈ K .

Exemple
Si F (t, x) est de classe C 1 sur U alors F (t, x) est
localement lipshitzienne par rapport à x ∈ m en effet: R
Si K = [a, b] × B ⊂ U, on pose M = supK kDx F (t, x)k (Dx : la
differentielle par rapport à x).
Z 1
F (t, x) − F (t, y ) = Dx F (t, y + θ(x − y ))(x − y )dθ
0
ce qui donne
kF (t, x) − F (t, y )k ≤ Mkx − y k, ∀(x, y ) ∈ K .
Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche
Théorème
Si F (t, x) est une fonction localement lipshitzienne par
R
rapport à x (x ∈ m ) alors il existe une solution et une
seule y (t) du problème de Cauchy (1) définie de I0 dans m R
et qui vérifie la propriété suivante:
Si y1 (t) est une solution du même problème sur un
intervalle I1 contenant t0 alors I1 ⊂ I0 et y1 (t) = y (t) pour
tout t ∈ I1 .

Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche

Méthode d’Euler

On se donne une subdivision de l’intervalle I0 = [t0 , t0 + T ]

t0 < t1 < t2 < · · · < tn < tn+1 < · · · < tN = t0 + T

on pose hn = tn+1 − tn et = max0≤n≤N−1 hn

sur le "petit" intervalle [tn , tn+1 ] la solution vérifie :
Z t
y (t) = y (tn ) + F (s, y (s))ds
tn

ce qui donne
Z tn+1
y (tn+1 ) = y (tn ) + F (s, y (s))ds
tn

si yn est une approximation de y (tn ) alors sur l’intervalle

]tn , tn+1 [ on approche F (s, y (s)) par F (tn , yn )

Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche

On obtient ainsi la méthode récursive

yn+1 = yn + hn F (tn , yn )
y0 = µ

appelée méthode d’Euler.

Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche

Interprétation graphique

Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche

Estimation d’erreur de la méthode d’Euler

L’erreur entre la solution exacte et la solution approchée

au point tn est
en = y (tn ) − yn
On pose

εn = y (tn+1 ) − y (tn ) − hn F (tn , y (tn ))

où y est la solution exacte du problème de Cauchy.

Théorème
n−1
X
∀n ≥ 0, |en | ≤ ec(tn −t0 ) |e0 | + ec(tn −ti+1 ) |εi |
i=0

où c est la constante: |F (t, x) − F (t, y )| ≤ c|x − y |.

Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche

Lemme de Grönwall discret

Soient (un ), (vn ) et (tn ) trois suites. On pose hn = tn+1 − tn

Si ∀n ≥ 0, un+1 ≤ (1 + hn λ)un + vn

alors
n−1
X
∀n ≥ 0, un ≤ eλ(tn −t0 ) u0 + eλ(tn −ti+1 ) vi
i=0

Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche

Preuve

par récurrence sur n

pour n = 0, |e0 | ≤ |e0 |, vraie
On constate que

en+1 = en + hn (F (tn , y (tn )) − F (tn , yn )) + εn

ce qui donne

|en+1 | ≤ |en | + hn c|y (tn ) − yn | + |εn | = (1 + hn c)|en | + |εn |

on utilise le fait que: 1 + x ≤ ex ∀x ∈ R et on aura:

|en+1 | ≤ ehn c |en | + |εn |

et l’hypothèse de récurrence appliquée à |en | donne le
résultat pour n + 1.

Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche

Corollaire
Si c 6= 0

ec(tn −t0 ) − 1
|en | ≤ ec(tn −t0 ) |e0 | + ω(y 0 , h)
c
où ω(y 0 , h) est le module de continuité de y 0 de pas
h.
Si y ∈ C 2 (I0 )
Z tn
|en | ≤ ec(tn −t0 ) |e0 | + h ec(tn −s) |y 00 (s)|ds
t0

Si la solution est de calsse C 1 ou C 2 la méthode est

convergente.
La méthode est dite convergente si
lim |e0 | = 0 =⇒ lim max |en | = 0.
h→0 h→0 n

Seulement cette convegence est très lente, elle est

d’ordre h. Cette lenteur laisse la méthode très instable.
Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche
Preuve

La preuve du premier point à faire comme exercice.

Preuve du second point :
Si y ∈ C 2 (I0 ) alors deux intégrations par partie donne
Z tn+1 Z tn+1
εn = (tn −s)y 00 (s)ds ce qui donne |εn | ≤ h |y 00 (s)|ds
tn tn
en remplaçant dans l’expression précedente de l’estimation
d’erreur de |en | on obtient
n−1
X
|en | ≤ ec(tn −t0 ) |e0 | + ec(tn −ti+1 ) |εi |
i=0
n−1
X Z ti+1
≤ ec(tn −t0 ) |e0 | + ec(tn −ti+1 ) h |y 00 (s)|ds
i=0 ti
Z tn
≤ ec(tn −t0 ) |e0 | + h ec(tn −s) |y 00 (s)|ds
t0

Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche

On pose M = maxs∈I0 |y 00 (s)|, et on aura
Z tn
|en | ≤ ec(tn −t0 ) |e0 | + hM ec(tn −s) ds
t0

ec(tn −t0 ) − 1
≤ ec(tn −t0 ) |e0 | + hM
c

et par consequent

ecT − 1
max |en | ≤ ecT |e0 | + hM
0≤n≤N c

Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche

Exemple

Considérons le problème de Cauchy

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

La solution exacte de ce problème est y (t) = tg(t).

La résolution approchée par la méthode d’Euler avec la pas
h = 0.01 donne

Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche

 n tn yn |y (tn ) − yn |
0 0.000000 0.000000 0.000000
1 0.010000 0.010000 0.000000
2 0.020000 0.020001 0.000002
3 0.030000 0.030005 0.000004
4 0.040000 0.040014 0.000007
5 0.050000 0.050030 0.000012
6 0.060000 0.060055 0.000017
7 0.070000 0.070091 0.000023
8 0.080000 0.080140 0.000031
9 0.090000 0.090204 0.000039
10 0.100000 0.100286 0.000049
.
. .
. .
. .
.
. . . .
71 0.710000 0.854804 0.004725
72 0.720000 0.872111 0.004957
73 0.730000 0.889717 0.005201
74 0.740000 0.907633 0.005457
75 0.750000 0.925870 0.005726
76 0.760000 0.944443 0.006009
77 0.770000 0.963363 0.006306
78 0.780000 0.982643 0.006618
79 0.790000 1.002299 0.006947
80 0.800000 1.022345 0.007293

Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche

 n tn yn |y (tn ) − yn |
100 1.000000 1.536977 0.020431
101 1.010000 1.570600 0.021606
102 1.020000 1.605268 0.022863
103 1.030000 1.641037 0.024207
104 1.040000 1.677967 0.025648
105 1.050000 1.716122 0.027193
106 1.060000 1.755573 0.028852
107 1.070000 1.796393 0.030635
108 1.080000 1.838664 0.032553
109 1.090000 1.882470 0.034621
110 1.100000 1.927907 0.036852
.
. .
. .
. .
.
. . . .
141 1.410000 5.561946 0.603408
142 1.420000 5.881299 0.699819
143 1.430000 6.237196 0.818265
144 1.440000 6.636222 0.965601
145 1.450000 7.086616 1.151472
146 1.460000 7.598817 1.389784
147 1.470000 8.186237 1.701129
148 1.480000 8.866383 2.116999
149 1.490000 9.662510 2.687348
150 1.500000 10.606151 3.495269

Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche

 15
la solution exacte y(x)=tg(x)
Solution approchée par la méthode d’Euler
h=0.01

10

0
0 0.5 1 1.5

Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche

 3.5

3 Erreur de la méthode d’Euler avec le pas h=0.01

2.5

1.5

0.5

0
0 0.5 1 1.5

Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche

Schéma d’Euler perturbé

En tenant compte des erreurs d’arrondis le schéma d’Euler

perturbé sera:

ȳn+1 = ȳn + hn F (tn , ȳn ) + αn

et on suppose que les erreurs d’arrondis αn sont majorées

par une constante |αn | ≤ α.
Montrer que
Z tn n−1
X
|ēn | ≤ ec(tn −t0 ) |ē0 | + h ec(tn −s) |y 00 (s)|ds + α ec(tn −ti+1 ) .
t0 i=0

les effets des erreurs d’arrondis sont

N−1 N−1
X X α
α ec(tN −ti+1 ) ≥ +α e0 = Nα =
h
i=0 i=0

ces effets tendent vers +∞ (comme h1 ) lorsque h → 0.

Même si la solution y est dans C 2 (I), on ne peut pas
conclure que la méthode d’Euler converge.
Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche
Notons ϕ(h) l’estimation d’erreur globale (tenant compte
des erreurs d’arrondis)
Nous pouvons conclure que si e0 = 0 alors

1
ϕ(h) = c1 h + c2
h
dont le graphe est:

Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche

Problème
Pour assurer la stabilité de la méthode d’Euler comment
faut-il choisir les hn ?
∂F
Montrer que si ∂y < 0 alors pour tout h tel que

2 ∂F
h< où λ(t) = (t, y (t))
maxt0 ≤t≤t0 +T |λ(t)| ∂y
la méthode d’Euler est stable.

Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche

Considérons le problème de Cauchy :

y 0 (t) = 1 − y 2 , t ∈ [0, +∞[ et y (0) = 0 (P1 )

Déterminer le seuil du choix du pas h qui assure la

stabilité.
Même question pour le problème de Cauchy suivant :

y 0 (t) = −2y (t), t ∈ [0, +∞[ et y (0) = 1 (P2 )

Même question pour le problème de Cauchy suivant :

y 0 (t) = arctg(3y ) − 3y + t, t ∈ [0, +∞[ et y (0) = 1

(P3 )

Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche

Résolution numérique du problème (P1 ) par la méthode d’Euler

solution exacte

solution approchØe solution exacte

solution approchØe

le pas h = 0.95 le pas h = 1.05

le pas optimal qui assure la stablité est h = 1,

e2t − 1
la solution exacte est : y (t) = 2t
e +1

Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche

Résolution numérique du problème (P2 ) par la méthode d’Euler

Solution exacte
Solution approchØe

solution exacte
solution approchØe

le pas h = 0.9 le pas h = 1.1

le pas optimal qui assure la stablité est h = 1,

la solution exacte est : y (t) = e−2t

Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche

Résolution numérique du problème (P3 ) par la méthode d’Euler

h = 0.75 et h = 0.05
2
le pas optimal qui assure la stablité est h = 3

Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche

Méthodes à un pas

L’amélioration de la méthode d’Euler passe par l’étude des

méthodes à un pas.
Une méthode à un pas est de la forme :

yn+1 = yn + hn ψ(tn , yn , hn )
y0 = µ (donnée)
R
où ψ est une fonction continue de I0 × × [0, h0 ] dans R.
la méthode d’Euler correspond à ψ(t, y , h) = F (t, y ).

Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche

Consistance, stabilité et convergence

Consistance :
Soit y (t) la solution exacte du problème de Cauchy, on pose

y (tn+1 ) − y (tn )
εn = − ψ(tn , y (tn ), hn )
hn
c’est l’erreur obtenue en injectant la solution exacte dans la
méthode
La méthode est dite consistante si

lim εn = 0
h→0

où
h= max hn
0≤n≤N−1

Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche

Stabilité :
Soit (yn ) la suite obtenue par la méthode à un pas.
et soit (zn ) la suite obtenue par la même méthode mais
perturbée.

zn+1 = zn + hn (ψ(tn , yn , hn ) + αn ))

où les perturbations αn sont bornées.

On dit que la méthode est stable ssi il existe deux constantes

M1 et M2 tels que :

max |yn − zn | ≤ M1 |y0 − z0 | + M2 max |αn |

n n

Il est nécessaire d’étudier la stabilité, une méthode instable ne

présente aucun intérêt.

Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche

Convergence

Définition
Soit l’erreur à l’instant tn , en = y (tn ) − yn , où y (t) est la solution
exacte. Nous dirons que la méthode est convergente ssi

lim |e0 | = 0 =⇒ lim max |en | = 0.

h→0 h→0 n

Théorème
Si une méthode est stable et consistante alors elle est
convergente.

Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche

Preuve :
On a :
∆y (tn )
= ψ(tn , yn , hn ) + εn
hn
La méthode est consistante alors on a :
∆y (tn )
lim max |εn | = lim max − ψ(tn , yn , hn ) = 0
h→0 n h→0 n hn

la méthode est stable alors il existe deux constantes M1 et M2

telles que

max |en | ≤ M1 |e0 | + M2 max |εn |

0<n≤N−1 n

et donc
lim |e0 | = 0 =⇒ lim max |en | = 0.
h→0 h→0 n

Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche

lemme 1
Une condition nécessaire et suffisante pour qu’une méthode à
un pas soit consistante est que

ψ(t, y , 0) = F (t, y ), ∀t ∈ I0 , ∀y ∈ R.
lemme 2
Si ψ(t, y , h) vérifie la condition de Lipschitz par rapport à la
variable y càd il existe une constante c indépendante de h telle
que
|ψ(t, y , h) − ψ(t, ȳ , h)| ≤ c|y − ȳ |
et ceci ∀t ∈ I0 , ∀y , ȳ ∈ R
et ∀h ∈ [0, h0 ]
alors la méthode est stable.

Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche

Ordre de convergence

Définition
Une méthode à un pas est d’ordre p (p > 0) si

∆y (tn )
max − ψ(tn , y (tn ), hn ) ≤ Khp
0≤n≤N−1 hn

où K est une constante indépendante de h.

Soit f (t, y ) une fonction à deux variables, on note :

f (0) (t, y ) = f (t, y )
∂f ∂f
f (1) (t, y ) = + f
∂t ∂y
.. ..
. .
∂f (k −1) ∂f (k −1)
f (k ) (t, y ) = + f
∂t ∂y

Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche

Théorème
Si F (t, y ) est p fois continuememnt différentiable sur I0 × R et si
∂ψ ∂pψ
ψ, ,··· , p existent et sont continues
∂h ∂h
sur I0 × R × [0, h0] alors la méthode est d’ordre p ssi

ψ(t, y , 0) = F (t, y )
∂ψ 1
(t, y , 0) = F (1) (t, y )
∂h 2
.. ..
. .
∂ p−1 ψ 1 (p−1)
(t, y , 0) = F (t, y )
∂hp−1 p

Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche

Méthode de Taylor

On considère le problème de Cauchy (2).

Supposons que F est p fois continuement différentiable sur
R
I0 × . La méthode de Taylor est définie par le schéma récursif
suivant : 
yn+1 = yn + hn ψ(tn , yn , hn )
y0 = µ (donnée)
avec
h hp−1 (p−1)
ψ(t, y , h) = F (t, y ) + F (1) (t, y ) + · · · + F (t, y )
2 p!
Il est évident que ψ vérifie :
ψ(t, y , 0) = F (t, y )
∂ψ 1
(t, y , 0) = F (1) (t, y )
∂h 2
.. ..
. .
∂ p−1 ψ 1 (p−1)
(t, y , 0) = F (t, y )
∂hp−1 p
Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche
on en déduit que la méthode est consistante et si elle est
convergente, elle sera d’ordre p.
Si de plus toutes les fonctions F (k ) (t, y ), k = 0, · · · , p − 1 sont
lipschitziennes par rappot à y càd

∃ck , |F (k ) (t, y ) − F (k ) (t, ȳ )| ≤ ck |y − ȳ |, ∀y , ȳ ∈ R

alors la méthode de Taylor est stable et par conséquent, elle
est convergente.

hp−1
 
h
|ψ(t, y , h) − ψ(t, ȳ , h)| ≤ c0 + c1 + · · · + cp−1 |y − ȳ |
2! p!
h0p−1
!
h0
≤ c0 + c1 + · · · + cp−1 |y − ȳ |
2! p!
≤ c|y − ȳ |

Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche

Méthode de Heun

On considère le problème de Cauchy (2).

On se place dans les conditions d’existence de la solution , F
est lipschitzienne sur . R
|F (t, y ) − F (t, ȳ )| ≤ c|y − ȳ |, ∀y , ȳ ∈ R
La méthode de Heun est définie par le schéma récursif suivant :

yn+1 = yn + hn ψ(tn , yn , hn )
y0 = µ (donnée)

avec
 
1 3 2 2
ψ(t, y , h) = F (t, y ) + F t + h, y + hF (t, y )
4 4 3 3

Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche

La méthode de Heun est stable en effet
1
|ψ(t, y , h) − ψ(t, ȳ , h)| ≤ |F (t, y ) − F (t, ȳ ))|
 4   
3 2 2 2 2
+ F t + h, y + hF (t, y ) − F t + h, ȳ + hF (t, ȳ )
4 3 3 3 3
   
1 3 2 2
≤ c|y − ȳ | + c y + hF (t, y ) − ȳ + hF (t, ȳ )
4 4 3 3
1
≤ c|y − ȳ | + ch|F (t, y ) − F (t, ȳ )|
2
1 2
≤ (c + c h0 )|y − ȳ |.
2
On a ψ(t, y , 0) = F (t, y ),
la méthode de Heun est consistante et par consequent, elle est
convergente.

Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche

Quel est son ordre ?

 
∂ψ 1 ∂F 2 2
(t, y , h) = t + h, y + hF (t, y )
∂h 2 ∂t 3 3
 
1 ∂F 2 2
+ t + h, y + hF (t, y )
2 ∂y 3 3

∂ψ
On a ψ(t, y , 0) = F (t, y ), ∂h (t, y , 0) = 12 F (1) (t, y )
∂2ψ
et après calcul de ∂h2
(t, y , h), on trouve

∂2ψ 1
2
(t, y , 0) 6= F (2) (t, y )
∂h 3
et donc la méthode est d’ordre 2.

Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche

Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche
Les méthodes Runge-Kutta

On considère le problème de Cauchy (2).

On se place dans les conditions d’existence de la solution
, F est lipschitzienne sur . R
|F (t, y ) − F (t, ȳ )| ≤ c|y − ȳ |, ∀y , ȳ ∈ R
Les méthodes Runge-Kutta sont des méthodes d’approximation
à un pas, elles sont définies par:
 Pr
 yi = y + h j=1 aij F (t + cj h, yj ) i = 1, . . . , r
ki = F (t +P
ci h, yi ) i = 1, . . . , r
r
ψ(t, y , h) = j=1 bj kj


Le schéma récursif est alors:

 Pr
 yn,i = yn + h j=1 aij F (tn + cj h, yn,j ) i = 1, . . . , r
kn,i = F (tn + P
ci h, yn,i ) i = 1, . . . , r
r
yn+1 = yn + hn j=1 bj kn,j


Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche

Elles sont représentées par

c1 a11 ··· a1r

c2 a21 ··· a2r
.
. .
. .. .
.
. . . .
cr ar 1 ··· arr
b1 ··· br

Elles sont dites explicites si c1 = 0 et aij = 0 pour

i ≤ j. Elles sont implicites dans le cas contraire.

Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche

Exemples

Méthode d’Euler explicite Méthode d’Euler implicite

0 0 1 0
(r = 1) (r = 1)
1 1

yn1 = yn et kn,1 = F (tn , yn )

yn1 = yn+1 ' y (tn1 ) = y (tn+1 ) et
ce qui donne kn,1 = F (tn+1 , yn+1 ) ce qui donne
yn+1 = yn + hn F (tn , yn ) yn+1 = yn + hn F (tn+1 , yn+1 )
Autre exemple:
0 0 0
α α 0 (r = 2)
1 1
1 − 2α 2α

Pour α = 12 , on retrouve la méthode de la tangente

améliorée.
Pour α = 23 , on retrouve la méthode de Heun.

Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche

Autre exemple :

1 1 1 1 1
2 − √
2 3 4 4 − √
2 3

1 1 1 1 1
2 + √
4 + √
4
2 3 2 3 (r = 2)

1 1
2 2

est une méthode Runge-Kutta implicite de rang 2 et d’ordre

4.

Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche

Idée derrière les methodes Runge-Kutta

Nous avons vu que la solution de l’équation differentielle

sur l’intevalle [tn , tn+1 ] est
Z tn+1
y (tn+1 ) = y (tn ) + F (s, y (s))ds
tn

Avant de passer de tn à tn+1 , on cherche des instants

intermédiaires tni = tn + ci h, i = 1, · · · , r de sorte à pouvoir
corriger les pentes avant d’arriver à tn+1 , ceci peut se
réaliser de la façon suivante, on a:
Z tni
y (tni ) = y (tn ) + F (s, y (s))ds
tn

R tni
On approche l’intégrale tn F (s, y (s))ds par une méthode de quadrature de
la forme Z tn i r
X
F (s, y (s))ds ' h aij F (tnj , y (tnj ))
tn j=1

Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche

et la valeur approchée intermédiaire yni sera
r
X
yni = yn + h aij F (tnj , y (tnj ))
j=1

R tn+1
et l’intégrale tn
F (s, y (s))ds par une autre méthode de
qudrature
Z tn+1 r
X
F (s, y (s))ds ' h bj F (tnj , ynj )
tn j=1

et la valeur approchée yn+1 à l’instant tn+1 sera

r
X
yn+1 = yn + h bj F (tnj , ynj )
j=1

mais comment déterminer les coefficients ci , bi et aij ,

i = 1, · · · , r et j = 1, · · · , r ?

Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche

Par application du théorème page 34

On pose

c1 0 ··· 0
 
b1 1
     
a11 ··· a1r  0 c2 
A =  ··· ··· ···  C =  . et e =  .
 b=
     
. .. . .
.  .
. 
. . .
ar 1 ··· arr  . . 
br 1
0 ··· 0 cr

ordre p Conditions

1 bT .e = 1

1
2 bT .(Ce) = bT .(Ae) =
2
1
3 bT .(C 2 e) = bT .(CAe) = bT .(Ae)2 =
3
1
bT .(ACe) = bT .(A2 e) =
6
.
.
. ···
N.B: Le carré d’un vecteur désigne le carré des ses composantes.

Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche

Méthode de Runge-Kutta classique RK4

C’est une méthode explicite, elle est définie par le schéma

récursif suivant :

yn+1 = yn + hn ψ(tn , yn , hn )
y0 = µ (donnée)

où ψ(t, y , h) est définie par :

Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche

 1
ψ(t, y , h) = (k1 + 2k2 + 2k3 + k4 )
6
et

k1 = F (t, y )
 
h h
k2 = F t + , y + k1
2 2
 
h h
k3 = F t + , y + k2
2 2
k4 = F (t + h, y + hk3 )

0 0 0 0 0
1 1
2 2 0 0 0
1 1
elle est représenté par : 2 0 2 0 0
1 0 0 1 0
1 2 2 1
6 6 6 6
Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche
On montre que cette méthode est stable et que

ψ(t, y , 0) = F (t, y )
∂ψ 1
(t, y , 0) = F (1) (t, y )
∂h 2
∂2ψ 1 (2)
(t, y , 0) = F (t, y )
∂h2 3
∂3ψ 1 (3)
(t, y , 0) = F (t, y )
∂h3 4

càd la méthode est convergente et d’ordre 4.

C’est pour cette raison cette méthode est appelée Runge-Kutta
4 (RK4)

Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche

 15

la solution exacte y(t)=tg(t)

Solution approchée par la méthode Runge−Kutta

avec pas 0.01
10

0
0 0.5 1 1.5

Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche

 −4
x 10
6
Erreur de la méthode Runge−Kutta avec pas 0.01

0
0 0.5 1 1.5

Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche

Résolution numérique du problème (P1 ) par la méthode
Runge-Kutta RK4
1 1

0.8 0.8

Solution exacte
Solution exacte
0.6 0.6
Solution Approchee
Solution approchee

0.4 0.4

0.2 0.2

0 0
0 5 10 15 20 0 5 10 15 20

le pas h = 0.95 le pas h = 1.05

e2t − 1
la solution exacte est : y (t) =
e2t + 1

Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche

Résolution numérique du problème (P2 ) par la méthode
Runge-Kutta RK4
0.9 0.9

0.8 0.8

0.7 0.7
Solution exacte
0.6 0.6
Solution approchee Solution exacte
0.5 0.5

0.4 0.4 Solution approchee

0.3 0.3

0.2 0.2

0.1 0.1

0 0

−0.1 −0.1
0 2 4 6 8 10 0 2 4 6 8 10

le pas h = 0.9 le pas h = 1.1

la solution exacte est : y (t) = e−2t

Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche

illustration de la performance des méthodes implicites

Résolution numérique du problème (P1 ) par les deux

méthodes : Euler explicite et Euler implicite

le pas h = 0.95 le pas h = 1.05

Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche

Représentation des erreurs absolues des deux méthodes :
Euler explicite (blue) et Euler implicite (rouge)

le pas h = 0.95 le pas h = 1.05

Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche

Equations différentielles d’ordre élevé

Soit I0 un intevalle de . R
Considérons l’équation différentielle d’ordre m

y (m) (t) = ϕ(t, y (t), y 0 (t), · · · , y (m−1) ), t ∈ I0 . (Eq)

où ϕ est une fonction de I0 × m dans R R.

On pose
y1 (t) = y (t)
y2 (t) = y10 (t)
..
.
0
ym (t) = ym−1 (t)
et Y (t) = (y1 (t), · · · , ym (t)), Y 0 (t) = (y10 (t), · · · , ym
0 (t))

L’équation précédente devienne :

Y 0 (t) = (y2 , y3 , · · · , ym , ϕ(t, y (t), y 0 (t), · · · , y (m−1) )), t ∈ I0

Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche

que l’on peut écrire sous la forme :

Y 0 (t) = (y2 , y3 , · · · , ym−1 , ϕ(t, y1 , y2 , · · · , ym )), t ∈ I0

Pour tout Y ∈ Rm , Y = (y1 , y2 , · · · , ym ), on pose

F (t, Y ) = (y2 , y3 , · · · , ym , ϕ(t, y1 , y2 , · · · , ym ))

et F sera une fonction vectorielle définie de I0 × m dans m . R R

Ainsi nous avons transformé l’équation diffrentielle (Eq) d’ordre
m en un système d’équations différentielles de taille m et
d’ordre 1 qui est

Y 0 (t) = F (t, Y (t)), t ∈ I0 .

Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche

La fonction prédéfinie ode45

Soit le système d’équations différentielles:

 0
 y1 (t) = F1 (t, y1 (t), y2 (t))
t ∈ [t0 , t0 + T ]
 0
y2 (t) = F2 (t, y1 (t), y2 (t))

avec les conditions initiales suivantes:

y1 (t0 ) = α et y2 (t0 ) = β

Exercice: Ecrire les méthodes déjà vues adaptées à la

résolution approchée de ce système.
Tout d’abord, on définit la fonction à 3 variables
 
F1 (t, x, y )
F (t, x, y ) =
F2 (t, x, y )
dans un fichier qu’on appélera: F.m "de façon
vectorielle"

Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche

de la manière suivante:
function d=F(t,v)
d=[F1 (t, v (1), v (2));F1 (t, v (1), v (2))];
Le programme sera écrit dans un autre fichier.m comme
suit:
[t,y]=ode45(’F’,[t0 , t0 + T ],[α;β]);
la solution sera récupérée dans le tableau à deux
colonnes: y. Le tableau t contiendra les points de la
subdivision de l’intervalle [t0 , t0 + T ].
Si on veut traçer la courbe solution (t, y1 (t)) pour
t ∈ [t0 , t0 + T ] alors on trace la première composante du
vecteur y en fonction de t comme suit: plot(t,y(:,1));
Si on veut traçer la courbe solution (t, y2 (t)) pour
t ∈ [t0 , t0 + T ] alors on trace la deuxième composante du
vecteur y en fonction de t comme suit: plot(t,y(:,2));
Si on veut traçer la courbe paramétrique solution
(y1 (t), y2 (t)) pour t ∈ [t0 , t0 + T ] alors on trace le vecteur y
comme suit: plot(y(:,1),y(:,2));

Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche

Comparaison entre Euler et ode45

dx


 (t) = y 2 (t)
dt






dy t ∈ [0, 10]
q
Considérons le problème
(t) = 1 + x 2 (t)
dt







x(0) = 1 y (0) = 0.4


14
y (t) ode45
1
12 y2(t) ode45
y1(t) Euler h=0.2
10 y2(t) Euler h=0.2

−2
0 2 4 6 8 10
Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche
Cet exemple met en défaut la méthode d’Euler

0.6
y (t) ode45
2
0.4 y2(t) Euler h=0.05

0.2

−0.2

−0.4

−0.6

−0.8

−1

−1.2

−1.4
0 1 2 3 4 5 6 7 8 9 10

Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche

Problème des trois corps

Considérons les trois corps suivants : La terre, la lune

et un satellite. Et on suppose que la terre et la lune
soient en rotation circulaire dans un plan (vitesse
constante) autour du centre de gravité du système. Le
mouvement du satellite est dû à la force d’attraction des
deux corps, celle-ci est proportionnelle à m1r 2m2 où m1 , m2
sont les masses des corps et r la distance entre eux.

1 − µ : la masse de la terre µ = 0.012277471 : la masse de

la lune. ε : la masse du satellite ε → 0 (négligeable)
L’origine est au centre de gravité du système,
Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche
on a A = −µeit (pour la terre) et B = (1 − µ)eit (pour la
lune) En appliquant la loi de Newton
(masse×accélération=force), on obtient pour le mouvement
du satellite (Y ∈ ) C
ε(1 − µ) A − Y εµ B−Y
εY 00 = + (1)
kA − Y k2 kA − Y k kB − Y k2 kB − Y k
Introduisons la variable y = e−it Y = y1 + iy2 , et on montre
que
−µ − y 1−µ−y
y 00 + 2iy 0 − y = (1 − µ) 3
+µ (2)
kµ + y k k1 − µ − y k3
La transformation de cette équation différentielle du
second ordre non linéaire donne le système d’équation
différentiel de première ordre suivant:
 0
 y1 (t) = y3 (t)
0


 y (t) = y4 (t)
 2

1−µ

0 µ
y3 (t) = y1 (t) + 2y4 (t) − 3
(y1 (t) + µ) − 3 (y1 (t) − 1 + µ)

 r 1 r 2

 0 1−µ µ
 y4 (t) = y2 (t) − 2y3 (t) − y2 (t) − 3 y2 (t)


r13 r2

Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche

où
q q
r1 = (y1 (t) + µ)2 + y22 (t) et r2 = (y1 (t) − 1 + µ)2 + y22 (t)

Système que l’on peut écrire sous la forme

Y 0 (t) = F (t, Y (t)), t ∈ I0 .

où F est la fonction vectorielle définie de I0 × 4 dans R

R 4
par
 
  y3
y1  y4 
1−µ µ
 
 y2   
F (t, 
 ) =  y1 + 2y4 −

3
(y1 + µ) − 3 (y1 − 1 + µ) 
y3   r 1 r2


y4  (1 − µ) µ 
y2 − 2y3 − 3
y2 − 3 y2
r1 r2

Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche

 Programme "code Matlab" pour la résolution de ce système:
La fonction F
1 function dy = F (t , y )
2 mu = 0 . 0 1 2 3 2 0 8 7 8 5 3 2 2 0 8 4 ;
3 v (1) = y (3) ;
4 v (2) = y (4) ;
5 r1 = sqrt (( y (1) + mu ) ^2+ y (2) ^2) ;
6 r2 = sqrt (( y (1) -1+ mu ) ^2+ y (2) ^2) ;
7 v (3) = y (1) +2* y (4) -((1 - mu ) *( y (1) + mu ) /( r1 ^3) ) - mu *( y (1) -1+ mu ) /(
r2 ^3) ;
8 v (4) = y (2) -2* y (3) -((1 - mu ) * y (2) /( r1 ^3) ) - mu * y (2) /( r2 ^3) ;
9 dy =[ v (1) ; v (2) ; v (3) ; v (4) ];
Programme d’appel
1 [T , Y ]= ode45 ( ’F ’ ,[0
30] ,[0.994;0;0; -2.00158510637908252240537862224]) ;
2 N = length ( Y (: ,1) ) ;
3 plot (0 ,0 , ’o ’) ;
4 text (0.08 ,0 , ’ Centre ’ , ’ FontSize ’ ,8) ;
5 hold on ;
6 for i =1: N
7 plot ( Y (i ,1) ,Y (i ,2) , ’. ’) ;
8 pause (0.1) ;
9 end

Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche

 4

0 Terre

−1

−2

−3

−4
−3 −2 −1 0 1 2 3 4

Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche

Réaliser les fonctions Euler.m et RK44.m avec code Matlab
utilisant respectivement les méthodes d’Euler et de
Range-Kutta classique et résoudre le même problème en
replaçant la fonction prédéfinie ode45 par ces deux
fonctions et représenter dans chaque cas la solution
obtenue.

Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche

Système d’équations différentielles de Lotka-Volterra

qui modélise l’intéraction et l’évolution de deux

populations y1 et y2 où y1 sont les proies et y2 sont les
prédateurs.

dy1
 dt (t) = c1 y1 (t)(1 − b1 y1 (t) − d2 y2 (t))




 dy2 (t) = −c2 y2 (t)(1 − d1 y1 (t))




dt
où c1 et c2 sont les facteurs de croissance des deux
populations, d1 et d2 tiennent compte de l’interaction
entre les deux populations, tandis que b1 est lié à la
quantité de nourriture disponible pour la population des
proies y1 . avec les données suivantes:
c1 = 0.08, b1 = 0, d2 = 0.05, c2 = 0.06 et d1 = 0.002/c2 , la
population initiale des proies est 40 et la population
initiale des prédateurs est 20, tracer les courbes
d’évolution des deux espèces sur une durée de 120
(unité=ans). Donner une estimation des deux espèces après
50 ans.
Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche
programmation des methodes implicites

Ecrire une fonction "code Matlab"

function s=Newton(f,x0,eps) qui prends en entrée
l’équation, la valeur de démarrage x0 et la precision
eps avec laquelle sera calculée la solution et qui
retourne la solution approchée de l’éqution f (x) = 0
par la méthode de Newton.
Ecrire un programme en code Matlab ou un pseudo code
pour chacune des fonctions suivantes:
function y=EulerImpl(F,y0,N,a,b) et
function y=CrankNicholson(F,y0,N,a,b)
qui prennent en entrée la fonction F (t, y ) définissant
l’équation différentielle, la valeur initiale y 0, le
nombre N définissant le pas de discrétisation et
l’intervalle de résolution a et b et qui retourne la
solution approchée dans le tableau y respectivement
par chacune des méthodes implicites Euler retrograde
(ou implicite) et Crank-Nicholson.
Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche
Appliquer ces fonctions avec le même pas de
discrétisation h pour résoudre de façon approchée le
problème suivant

1
y 0 (t) = sin(2t) − y cos(t); t ∈ [0, 2π]; y (0) = 0; (1)
2
dont la solution exacte est
y (t) = −1 + exp(− sin(t)) + sin(t)
Représenter dans la même figure la courbe de la
solution exacte et les deux autres courbes
"approchées" obtenues par les deux méthodes.
Représenter la courbe erreur de chacune des deux
méthodes et interpréter les résultats obtenus.
Donner une estimation numérique de l’ordre de
convergence de chacune des deux méthodes. (Y-a-il
confirmation de résultats théoriques)

Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche

Oscillateur de Van der Pol
On considère l’équation différentielle y 00 (t) = λ(1 − y 2 (t))y 0 (t) − y (t)
où t varie dans l’intervalle
I = [0, 6.663], avec les conditions initiales suivantes: y (0) = 2 et
y 0 (0) = 0.
(Voir [Link]
Transformer cette équation différentielle du second ordre avec
conditions initiales en système différentiel de premier ordre avec
conditions initiales:
 0
 Y (t) = F (t, Y(t)),  y (t) 
y (0)
1 où Y(t) =
 Y(0) = y 0 (t)
y2 (0)

où F est une fonction vectorielle définie de I × 2 dans 2 . R R

Ecrire une fonction "code matlab" dont les paramètres
  d’entrée

y1
sont: t, y 1 et y 2 et qui retourne le vecteur F t, .
y2
Ecrire une fonction "code matlab" dont le prototype est
function Y= RK44(F,lambda,Y0,I,N)
R
où Y est un vecteur de d , I: l’intervalle de résolution, Y0: le
vecteur initial et F: la fonction définissant le système
différentiel, définie de I × d . R
Pour différentes valeurs de λ tracer (y1 (t), y2 (t)) t ∈ I.
Comparer ces trajectoires avec celles obtenues en utilisant ode45.
Cours: Méthodes Numériques Résolution numérique des E.D.O H. Allouche