Analyse numérique

Ibrahim ALAME

ESTP

20/01/2020

Ibrahim ALAME (ESTP) Analyse numérique 20/01/2020 1 / 34

Equations différentielles
Problème de Cauchy
y0 (t) = f(t; y(t)) dans I

(?)
y(t0 ) = y0

I un intervalle de R
f une fonction à deux variables de I × Rn dans Rn .
y : I → Rn la fonction inconnue.
t0 (= 0) est un point de I .

 0
y (t) = f1 (t, y1 (t), y2 (t), ...yn (t))
 10


y2 (t) = f2 (t, y1 (t), y2 (t), ...yn (t))

 ······
 0
yn (t) = fn (t, y1 (t), y2 (t), ...yn (t))

Ibrahim ALAME (ESTP) Analyse numérique 20/01/2020 2 / 34

Une équation différentielle d’ordre p de la forme
y (p) = f (t, y , y 0 , · · · , y (p−1) )
où p > 1 se ramène à un système équivalent en posant :
 0
z (t) = z2 (t)
 10


z2 (t) = z3 (t)

 ······
 0
zp (t) = f (t, z1 (t), z2 (t), ...zp−1 (t))
autrement dit :
z0 (t) = F(t; z(t))

Exemple
y 00 (t) − ay 0 (t) − by (t) = f (t)
On pose z1 = y et z2 = y 0 d’où le système
 0
z1 (t) = z2 (t)
z20 (t) = bz1 (t) + az2 (t)
Ibrahim ALAME (ESTP) Analyse numérique 20/01/2020 3 / 34
Theorem (de Cauchy-Lipschitz)
Soit f une application continue :

f : (t, z) ∈ (0, T ) × Rn 7→ f (t, z) ∈ Rn

Si f vérifie la condition de Lipschitz,

∀t ∈ [0, T ], ∀|z1 |, |z2 | ∈ Rn , |f (t, z1 ) − f (t, z2 )| ≤ L|z1 − z2 |

alors il y a existence et unicité locale de la solution du problème de Cauchy.

Le théorème est admis dans ce cours, il s’obtient via un argument de point

fixe pour les applications contractantes dans l’espace de Banach...

Ibrahim ALAME (ESTP) Analyse numérique 20/01/2020 4 / 34

Nous considérons maintenant trois exemples typiques :
y 0 (t) = ay (t)

f (t, y ) = ay =⇒ |f (t, y1 ) − f (t, y2 )| = a|y1 − y2 |

y 0 (t) = [y (t)]2 .
Nous vérifions facilement que f (t, y ) = y 2 est localement
lipschitzienne car f est de classe C 1 donc localement lipschitzienne
d’après le Théorème des accroissements fnis.
y 0 (t) = 3 y (t).
p
p
La fonction f (t, y ) = 3 y (t) n’est pas lipschitzienne car elle n’est pas
dérivable en 0 et donc le théorème ne s’applique pas !

Ibrahim ALAME (ESTP) Analyse numérique 20/01/2020 5 / 34

La méthode d’Euler
y

ei
εi

h= Tn
x
0=x0 xi T =xn

On cherche à approcher y (ti ) ' yi . Taylor à l’ordre 2 :

T
y (ti + h) = y (ti ) + hy 0 (ti ) + o(h) où h = ti+1 − ti =
n
or y 0 (ti ) = f (ti , y (ti )) d’où "l’idée" de considérer le schéma suivant

yi+1 = yi + hf (ti , yi )
y0 = y (0)
Ibrahim ALAME (ESTP) Analyse numérique 20/01/2020 6 / 34
Etude de l’erreur
ei = y (ti ) − yi l’erreur commise au point ti
εi+1 = y (ti+1 ) − y (ti ) − hf (ti , y (ti )) l’erreur de consistance, c’est
l’erreur systématique commise sur yi+1 : c’est l’erreur de passage de la
valeur à l’instant ti supposée exacte à la valeur approchée en ti+1 .

|ei+1 | = |y (ti+1 ) − yi+1 |

= |y (ti+1 ) − y (ti ) − hf (ti , y (ti ))
+y (ti ) + hf (ti , y (ti ))
−yi + yi
+hf (ti , yi ) − hf (ti , yi )
−yi+1 |
≤ |εi+1 | + |ei | + h |f (ti , y (ti )) − f (ti , yi )|
Si f est lipschitzienne en y :

|ei+1 | ≤ |εi+1 | + (1 + Lh) |ei |

Ibrahim ALAME (ESTP) Analyse numérique 20/01/2020 7 / 34

d’autre part

|εi+1 | = |y (ti+1 ) − y (ti ) − hf (ti , y (ti ))| 

R t Rt
=  ti i+1 y 0 (s)ds − ti i+1 f (ti , y (ti ))ds 

R 
 t
=  ti i+1 (y 0 (s) − y 0 (ti ))ds 


≤ h · ω(h, y 0 )

où ω(h, y 0 ) est le module de continuité de la fonction y 0 sur [0, T ].

Finalement
|ei+1 | ≤ (1 + Lh) |ei | + h · ω(h, y 0 )

Ibrahim ALAME (ESTP) Analyse numérique 20/01/2020 8 / 34

Lemme
Soit (θi )i une suite de nombres positifs satisfaisant :

θi+1 ≤ (1 + A)θi + B

où A et B sont des constantes strictement positives alors :

exp(iA) − 1
θi ≤ exp(iA)θ0 + B
A
Démonstration par récurrence. Grâce au lemme, on a :
exp(iLh) − 1
|ei | ≤ exp(iLh)|e0 | + h · ω(h, y 0 )
Lh
or e0 = 0 d’où le théorème :
Theorem
Si f est L-lipschitzienne alors

exp(LT ) − 1
max |ei | ≤ · ω(h, y 0 )
0≤i≤N L
Ibrahim ALAME (ESTP) Analyse numérique 20/01/2020 9 / 34
à partir de ce théorème on peut démontrer le résultat suivant :
Proposition
Si yh est la fonction affine par morceaux telle que yh (ti ) = yi alors

ky − yh k∞ −→ 0 quand h −→ 0

Ibrahim ALAME (ESTP) Analyse numérique 20/01/2020 10 / 34

Etude générale des méthodes à un pas

Les méthodes à un pas sont des méthodes de la forme :


yi+1 = yi + hΦ(t, yi , h) pour i = 0, 1, · · · , N − 1
y0 = y0,h

Où Φ est une fonction continue sur [0, T ] × R × [0, H], H désignant un pas
de discrétisation maximal.

Ibrahim ALAME (ESTP) Analyse numérique 20/01/2020 11 / 34

La méthode est dite consistante si Σh −→ 0 quand h −→ 0
N−1
X
Σh = |y (ti+1 ) − y (ti ) − hΦ(ti , y (ti ), h)|
i=0

Theorem
La méthode à un pas est consistante ssi pour tout t ∈ [0, T ] et y ∈ R :

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

εi+1 = y (ti+1 ) − y (ti ) − hΦ(ti , y (ti ), h)

Rt Rt
= ti i+1 y 0 (s)ds − ti i+1 Φ(ti , y (ti ), h)ds
R ti+1
= ti (f (s, y (s)) − Φ(ti , y (ti ), h))ds
|f (s, y (s)) − Φ(ti , y (ti ), h)| ≤ |f (s, y (s)) − f (ti , y (ti ))|
+ |f (ti , y (ti )) − Φ(ti , y (ti ), 0)|
+ |Φ(ti , y (ti ), 0) − Φ(ti , y (ti ), h)|
Un argument de continuité uniforme de f ,y et Φ permet de conclure.
Ibrahim ALAME (ESTP) Analyse numérique 20/01/2020 12 / 34
Stabilité
On perturbe le schéma de la façon suivante

ỹi+1 = ỹi + hΦ(t, ỹi , h) + εi pour i = 0, 1, · · · , N − 1
ỹ0 = ỹ0,h

La méthode est dite stable s’il existe deux constantes S1 et S2 telles que :
N−1
X
max ≤ S1 |y0,h − ỹ0,h | + S2 |εi |
i
i=0

Theorem (Condition suffisante de stabilité)

La méthode à un pas est stable si Φ est L-lipschitzienne par rapport à la
deuxième variable y avec L indépendant de deux autres variables t et h.

On note θi la différence entre la solution yi et sa perturbation ỹi . A l’aide

de l’inégalité triangulaire on obtient
θi+1 ≤ (1 + Lh)θi + |εi |
Ibrahim ALAME (ESTP) Analyse numérique 20/01/2020 13 / 34
 θi+1 θi |εi |
i+1
≤ i
+
(1 + Lh) (1 + Lh) (1 + Lh)i+1
Par sommation (télescopique) :
i−1
θi θ0 X |εk |
i
≤ 0
+
(1 + Lh) (1 + Lh) (1 + Lh)k+1
k=0

i−1
X
θi ≤ (1 + Lh)i θ0 + (1 + Lh)i−k |εi |
k=0
N−1
X
max θi ≤ exp(LT )θ0 + exp(LT ) |εi |
i
k=0

D’où le résultat avec S1 = S2 = exp(LT ).

Theorem
Tout méthode stable et consistante converge.
Ibrahim ALAME (ESTP) Analyse numérique 20/01/2020 14 / 34
Définition : La méthode est d’ordre p ≥ 1 s’il existe C > 0 telle que
|εi | = |y (ti+1 ) − y (ti ) − hΦ(ti , y (ti ), h)| ≤ Chp+1

Théorème : La méthode est d’ordre p si pour tout t, z et k ≤ p − 1 :

Φ(t, z, 0) = f (t, z)
∂Φ 1 [1]
(t, z, 0) = f (t, z)
∂h 2
..
.
∂k Φ 1
k
(t, z, 0) = f [k] (t, z)
∂h k +1
Dans ce cas
∂p Φ
 
1 1
|εi | = hp+1 [p]
f (ti , y (ti )) − (ti , y (ti ), 0) + o(hp+1 )
i! p+1 ∂hp

∂f [k] (t,z) ∂f [k] (t,z)

où f [k] est défini par f [0] = f et f [k+1] (t, z) = ∂t + ∂z × f (t, z)

Ibrahim ALAME (ESTP) Analyse numérique 20/01/2020 15 / 34

exercise
Soit la méthode à un pas définie par :

Φ(t, z, h) = α1 f (t, z) + α2 f (t + β1 h, z + β2 hf (t, z))

Déterminer les conditions sur les coefficients pour que la méthode soit
d’ordre le plus élevé possible.

La méthode est consistante d’ordre 1 si : Φ(t, z, 0) = f (t, z) cela

implique α1 + α2 = 1, car

Φ(t, z, 0) = α1 f (t, z) + α2 f (t + 0, z + 0)

La méthode est d’ordre 2 si

∂Φ 1
(t, z, 0) = f [1] (t, z)
∂h 2
∂f (t,z) ∂f (t,z)
où f [1] (t, z) = ∂t + ∂z × f (t, z).
Ibrahim ALAME (ESTP) Analyse numérique 20/01/2020 16 / 34
On a
∂Φ
∂h (t, z, h) = α2 β1 ∂f
∂t (t + β1 h, z + β2 hf (t, z))
∂f
+α2 β2 ∂z (t + β1 h, z + β2 hf (t, z)) · f (t, z)

Nous avons l’égalité en h = 0 si α2 β1 = α2 β2 = 12 . On vérifie facilement

que l’on ne peut pas aller plus loin. Finalement, si on pose α2 = α :

h h
Φ(t, z, h) = (1 − α)f (t, z) + αf (t + ,z + f (t, z))
2α 2α

On retrouve
La méthode de la tangente améliorée si α = 1.
La méthode d’Euler modifiée si α = 21 .
La méthode de Heun si α = 43 .

Ibrahim ALAME (ESTP) Analyse numérique 20/01/2020 17 / 34

La méthode d’Euler améliorée (α = 1)

h h
Φ(t, z, h) = f (t + , z + f (t, z))
2 2
y

h= Tn
x
0=x0 xi T =xn

yi+1 = yi + hf (ti + h2 , yi + h2 f (ti , yi ))



y0 = y (0)

Ibrahim ALAME (ESTP) Analyse numérique 20/01/2020 18 / 34

La méthode d’Euler modifiée (α = 1/2)

1
Φ(t, z, h) = [f (t, z) + f (t + h, z + hf (t, z))]
2
y

h= Tn
x
0=x0 xi T =xn

h

yi+1 = yi + 2 [f (ti , yi ) + f (ti + h, yi + hf (ti , yi ))]
y0 = y (0)

Ibrahim ALAME (ESTP) Analyse numérique 20/01/2020 19 / 34

La méthode de Heun si (α = 34 ).

 
1 2 2h
Φ(t, z, h) = f (t, z) + 3f (t + h, z + f (t, z))
4 3 3

h= Tn
x
0=x0 xi T =xn

h 2h 2h
  
yi+1 = yi + 4 f (ti , yi ) + 3f (ti + 3 , yi + 3 f (ti , yi ))
y0 = y (0)

Ibrahim ALAME (ESTP) Analyse numérique 20/01/2020 20 / 34

Méthode de Runge-Kutta

L’idée commence dans l’écriture intégrale de l’équation différentielle :

Z ti+1 Z ti+1
0
y (s) = f (s, y (s)) =⇒ y (s)ds = f (s, y (s))ds
ti ti
Z ti+1
y (ti+1 ) − y (ti ) − f (s, y (s))ds
ti

Puis on utilise une formule de quadrature pour approcher l’intégrale

q
X
y (ti+1 ) = y (ti ) + bj f (ti,j , y (ti,j ))
j=0

Les ti,j sont des points intermédiaire entre ti et ti + h.

Ibrahim ALAME (ESTP) Analyse numérique 20/01/2020 21 / 34

Méthode de Runge-Kutta
On écrit ti,j = ti + cj h. Une nouvelle formule de quadrature à j points
permet d’approcher yi,j avec des coefficients dependant de j :

ti,j = ti + cj h
= y (ti ) + j−1
P
y (ti,j ) k=0 ai,k f (ti,k , y (ti,k ))

En résumé, on calcul les valeurs de y aux points intermédiaires par des

formules de quadrature faisant intervenir les points précédents et les valeurs
déjà calculées.
On schématise souvent la méthode de Runge et Kutta par le tableau

0 0
c2 a2,1
c3 a3,1 a3,2
.. .. ..
. . .
cq aq,1 aq,2 · · · aq,q−1
b1 b2 · · · bq−1 bq
Ibrahim ALAME (ESTP) Analyse numérique 20/01/2020 22 / 34
Méthode de Runge-Kutta
Ce tableau représente le schéma suivant :

 
 ti,1 = ti + 0 × h = ti

 yi,1 = yi + 0 × f (·, ·) = yi







 ti,2 = ti + c2 h
 yi,2 = yi + a2,1 f (ti,1 , yi,1 )




ti,3 = ti + c3 h




yi,3 = yi + a3,1 f (ti,1 , yi,1 ) + a3,2 f (ti,2 , yi,2 )
 .
.
 .


 
ti,q = ti + cq h




yi,q = yi + aq,1 f (ti,1 , yi,1 ) + aq,2 f (ti,2 , yi,2 ) + · · · + aq,q−1 f (ti,q−1 , yi,q−1








yi+1 = yi + b1 f (ti,1 , yi,1 ) + b2 f (ti,2 , yi,2 ) + · · · + bq−1 f (ti,q−1 , yi,q−1 ) + bq


Ibrahim ALAME (ESTP) Analyse numérique 20/01/2020 23 / 34

La méthode de Runge-Kutta classique q = 4 (RK4) :
0 0
1 1
2 2
1 1
2 0 2
1 0 0 1
1 2 2 1
6 6 6 6
 
 ti,1 = ti

 yi,1 = yi 1






 ti,2 = ti + 2 h
  yi,2 = yi + 21 f (ti,1 , yi,1 )



ti,3 = ti + 12 h
1
  yi,3 = yi + 2 f (ti,2 , yi,2 )



ti,4 = ti + h







 yi,4 = yi + f (ti,3 , yi,3 )
yi+1 = yi + 16 f (ti,1 , yi,1 ) + 26 f (ti,2 , yi,2 ) + 26 f (ti,3 , yi,3 ) + 16 f (ti,4 , yi,4 )


Nous admettons que l’ordre de cette méthode est 4.

Ibrahim ALAME (ESTP) Analyse numérique 20/01/2020 24 / 34
Cette méthode s’écrit aussi plus couramment :

 k1 = f (ti , yi )
 2 = f (ti + h2 , yi + h2 k1 )

k


k3 = f (ti + h2 , yi + h2 k2 )
k = f (xi + h, yi + hk3 )

 4



yi+1 = yi + (k1 + 2k2 + 2k3 + k4 )/6

Ibrahim ALAME (ESTP) Analyse numérique 20/01/2020 25 / 34

exercise
Estimer l’erreur en t = 1 commise en appliquant la méthode de Runge et
Kutta à 4 points (RK4) à l’équation différentielle
 0
y (t) = −y (t) t ∈ [0, 1]
y (0) = 1

Ibrahim ALAME (ESTP) Analyse numérique 20/01/2020 26 / 34

Nous avons


 yi,1 = yi
yi,2 = yi + h2 f (ti,1 , yi,1 ) = (1 − h2 )yi



 2
yi,3 = yi + 12 f (ti,2 , yi,2 ) = yi − h2 yi,2 = (1 − h2 + h4 )yi




 y = y + f (t , y ) = y − h y = (1 − h + h2 − h3 )y


i,4 i i,3 i,3 i 2 i,3 2 4 i

yi+1 = yi + 16 f (ti,1 , yi,1 ) + 26 f (ti,2 , yi,2 ) + 26 f (ti,3 , yi,3 ) + 61 f (ti,4 , yi,4 )




yi − 61 (yi,1 + 2yi,2 + 2yi,3 + yi,4 )

= 




2 3 4
1 − h + h − h + h yi

=


2 6 24

Le schéma se réduit à
(  
h2 h3 h4
yi+1 = 1 − h + 2 − 6 + 24 yi
y0 = 1
4 n
 2 3

D’où yn = 1 − h + h2 − h6 + h24 or h = n1 donc
 n  
1 1 1 1 −1 1 1
yn = 1 − + 2 − 3 + =e 1− + o( 4 )
n 2n 6n 24n4 120n4 n
Ibrahim ALAME (ESTP) Analyse numérique 20/01/2020 27 / 34
d’où l’erreur
1
e = y (tn ) − yn ∼ −
120 e n4
Pour n = 10, nous avons y10 = 0, 3678798, la solution exacte étant
y (t) = e−t et e−1 = 0, 3678794. y10 est donc une approximation de e−1 à
mieux que 10−6 près.

Ibrahim ALAME (ESTP) Analyse numérique 20/01/2020 28 / 34

Soit l’équadiff
(x + 1)y 0 + 2y = sin(x)e −x/5


y (0) = −1/5
dont la solution exacte est

1912/5 − e −x/5 ((325x − 450) cos x + (65x − 235) sin x)

y (x) =
338x 2 + 676x + 338

Ibrahim ALAME (ESTP) Analyse numérique 20/01/2020 29 / 34

Le problème de Cauchy est
( −x/5
y 0 = sin(x)ex+1 −2y = F (x, y )
y (0) = −1/5

Ibrahim ALAME (ESTP) Analyse numérique 20/01/2020 30 / 34

Ibrahim ALAME (ESTP) Analyse numérique 20/01/2020 31 / 34
La méthode d’Euler améliorée

h h
Φ(t, z, h) = f (t + , z + f (t, z))
2 2

Ibrahim ALAME (ESTP) Analyse numérique 20/01/2020 32 / 34

Ibrahim ALAME (ESTP) Analyse numérique 20/01/2020 33 / 34
La méthode de Runge-Kutta

Ibrahim ALAME (ESTP) Analyse numérique 20/01/2020 34 / 34

La méthode de Runge-Kutta

Ibrahim ALAME (ESTP) Analyse numérique 20/01/2020 35 / 34