Équations différentielles

c.i.
c.l.

Géraldine Faure
[email protected]

Université Clermont Auvergne

2020

L3PhysNum (EUPI) Équations différentielles 2020 1 / 30

Equations différentielles à conditions initiales

Méthodes
Théorème de Cauchy
Méthodes Numériques
Euler explicite, implicite
Heun
RK2
RK4
Adams

L3PhysNum (EUPI) Équations différentielles 2020 2 / 30

Introduction

Deux catégories :
• problèmes à conditions initiales (C.I.)
• problèmes à conditions aux limites (C.L.)

F y , y 0 , y 00 , . . . , y (p) = g(t) (1)
On cherche la fonction y (t) sur un intervalle, g(t) est le second membre de l’équation.
Si g(t) = 0 on dit que l’équation différentielle est homogène.
Si p = 1 l’équation différentielle est d’ordre 1, et si p > 1, on peut toujours se ramener à un
système d’équations différentielles du 1er ordre.

L3PhysNum (EUPI) Équations différentielles 2020 3 / 30

 Équations différentielles C.I.

Problème de Cauchy
Soit : (
dy (t)
= f (y (t), t) (2)
dt
y (t0 ) = y0

où y : [a, b] → Rn et f : R × Rn → Rn et t0 ∈ [a, b]

Équation différentielle du second ordre


dy (t)

 = z(t)
 2
d y (t) dy (t)

 dt
 = f (t, y (t), )  dz(t)


 2 dt dt = f (t, y (t), z(t))
dt (4)
y(t0 ) = y0 (3)  y(t0 ) = y0
dy 
dy
y00
 

 = 
 = z(t0 ) = y00
dy t=t0 
dy t=t0

d
   
y (t) z(t)
= (5)
dt z(t) f (t, y , z)

L3PhysNum (EUPI) Équations différentielles 2020 4 / 30

 Équations différentielles C.I.

Généralisation
 dy1
 dt
= f1 (y1 , y2 , . . . , yn , t)
 dy2
= f2 (y1 , y2 , . . . , yn , t)
dt (6)
 .dy. n.

dt
= fn (y1 , y2 , . . . , yn , t)
avec les conditions initiales : 
 y1 (0) = y10
y2 (0) = y20
(7)
 ...
yn (0) = yn0

On peut écrire ce système sous la forme

d→
− −
→ −

y
dt
= f (→
y , t) (8)
−
→
y (0) = −
→
y 0

avec −
→y = (y1 , y2 , . . . , yn )
−
→
f = (f1 , f2 , . . . , fn )
et −
→
y0 = (y10 , y20 , . . . , y0n )

L3PhysNum (EUPI) Équations différentielles 2020 5 / 30

Théorème de Cauchy-Lipschitz

Théorème
Si
1 f : [t0 , t0 + a] × Rn → Rn est continue
2 Il existe L tel que pour tout (x , y , z) on ait :

||f (x , y ) − f (x , z)|| ≤ L||y − z||

alors le problème de Cauchy (

dy (t)
= f (y (t), t)
dt
y (t0 ) = y0
a une solution unique.
∂f
La condition 2 est remplie si ∂y (t, y ) est bornée uniformément en t.
Dans la suite on cherchera à résoudre des problèmes de Cauchy (fonction f convenable,
Lipschitzienne)

L3PhysNum (EUPI) Équations différentielles 2020 6 / 30

Méthodes Numériques : Principe

principe
On subdivise l’intervalle I en N sous-intervalles de longueur h, h est donc le pas de discrétisation.
On intègre l’équation différentielle
y 0 (t) = f (t, y (t))
entre tn et tn+1 (avec h = tn+1 − tn )
Z tn+1
y (tn+1 ) − y (tn ) = f (t, y (t))dt (9)
tn

Pour chaque nœud tn = t0 + nh avec 1 ≤ n ≤ N, on cherche la valeur inconnue un qui approche

y (tn )
L’ensemble des valeurs (u0 = y0 , u1 , u2 , . . . , uN ) représente la solution numérique.

L3PhysNum (EUPI) Équations différentielles 2020 7 / 30

Suivant la méthode de calcul du membre de droite de l’équation (9) on obtient différents schémas
de calcul.
Méthode d’Euler explicite : la plus simple à mettre en œuvre
Méthode d’Euler implicite
Méthode d’Euler modifiée
Méthode de Heun
Méthode de Runge-Kutta à l’ordre 2 et à l’ordre 4
Méthode d’Adams (pas multiples)

L3PhysNum (EUPI) Équations différentielles 2020 8 / 30

Méthode Euler progressif

schéma explicite d’Euler

rectangle à gauche Z tn+1
f (t, y (t))dt ≈ hf (tn , y (tn )) (10)
tn

On obtient à partir de (9) et (10) :

n
u0 = y (t0 ) = y0
(11)
un+1 = un + hf (tn , un ) ∀n = 0, 1, 2, . . . , N − 1

Ce schéma est dit explicite puisque le terme un+1 est fonction de un . A partir de la connaissance
de la c.i. y0 , u0 = y0 , le processus peut être exécuté.

L3PhysNum (EUPI) Équations différentielles 2020 9 / 30

Méthode Euler progressif

Exercice


dy
= f (t, y ) = t(1 − y )
dt (12)
y (0) = 0

Appliquer le schéma de Euler progressif (explicite) sur l’intervalle [0; 3], avec un pas h = 0, 5.
Pour cela, faites un tableau dont la dernière colonne correspondra à la valeur absolue de l’erreur,
t2
sachant que la solution exacte est y (t) = 1 − e − 2

A vos calculettes ! !

L3PhysNum (EUPI) Équations différentielles 2020 10 / 30

Solution

n tn yn yexacte |erreur |
0 0,000 0,000 0,0 0
1 0,500 0,000 0,117 0,117
2 1,000 0,250 0,393 0,143
3 1,500 0,625 0,675 0,05
4 2,000 0,906 0,865 0,041
5 2,500 0.9997 0,956 0,0437
6 3,000 1,000 0,989 0,011

Table – Résolution de (12), par Euler explicite avec h = 0, 5 et t0 = 0. On a les récurrences

suivantes :tn+1 = tn + h = t0 + nh = 0, 5 + tn et yn+1 = yn + 0, 5tn (1 − yn )

L3PhysNum (EUPI) Équations différentielles 2020 11 / 30

Méthode Euler rétrograde

schéma implicite d’Euler

rectangle à droite Z tn+1
f (t, y (t))dt ≈ hf (tn+1 , y (tn+1 )) (13)
tn

A partir des équations (9) et (13), on obtient :

n
u0 = y (t0 ) = y0
(14)
un+1 − hf (tn+1 , un+1 ) = un ∀n = 0, 1, 2, . . . , N − 1

Ce schéma est dit implicite puisque le terme un+1 est fonction de un et un+1 . Il n’est pas possible
d’obtenir directement un+1 .
deux possibilités :
résoudre l’équation (14) un+1 − hf (tn+1 , un+1 ) = un pour déterminer un+1
estimer un+1 par la méthode d’Euler progressif

L3PhysNum (EUPI) Équations différentielles 2020 12 / 30

Méthode Euler rétrograde

schéma implicite d’Euler

Soit ug
n+1 une estimation de un+1 , on obtient alors :
(
u0 = y0
ug
n+1 = un + hf (tn , un ) (15)
un+1 = un + hf (tn+1 , un+1
˜ )

L3PhysNum (EUPI) Équations différentielles 2020 13 / 30

Méthode Euler modifiée

schéma d’Euler modifiée

quadrature du point milieu
Z tn+1
h h
f (t, y (t))dt ≈ hf (tn + , y (tn + )) (16)
tn
2 2

A partir des équations (9) et (16), on obtient :


u0 = y (t0 ) = y0
(17)
un+1 = un + hf (tn + h2 , un+1/2 ) ∀n = 0, 1, 2, . . . , N − 1

où un+1/2 est une approximation de y (tn + h2 ) = y (tn+1/2 ).

De même que précédemment (Euler rétrograde), on utilise Euler progressif pour prédire un+1/2
par u^
n+1/2
On obtient le schéma d’Euler modifié suivant :
(
u0 = y0
u^
n+1/2 = un + h2 f (tn , un ) (18)
un+1 = un + hf (tn+1/2 , u^ n+1/2 ) ∀n = 0, 1, 2, . . . , N − 1

L3PhysNum (EUPI) Équations différentielles 2020 14 / 30

Méthode de Heun

schéma de Heun
utilisation de la formule des Trapèzes, on obtient le schéma de Heun suivant :

u0 = y0
(19)
un+1 − h2 f (tn+1 , un+1 ) = un + h2 f (tn , un ) ∀n = 0, 1, 2, . . . , N − 1

En prédisant un+1 on obtient :

(
u0 = y0
ug
n+1 = un + hf(tn , un )  (20)
un+1 = un + h2 f (tn , un ) + f (tn+1 , ug
n+1 )

L3PhysNum (EUPI) Équations différentielles 2020 15 / 30

A partir de Simpson

Z tn+1
h 
f (t, y )dt ≈ f (tn , yn ) + 4f (tn+1/2 , yn+1/2 ) + f (tn+1 , yn+1 ) (21)
tn
6

en utilisant le même raisonnement que précédemment on obtient le schéma RK2 suivant :


 u0 = y0
un + h2 f (tn , un )
 u^ =
n+1/2
ug = un + hf (tn , un )
 n+1
 h
 
un+1 = un + 6
f (tn , un ) + 4f (tn+1/2 , u^
n+1/2 ) + f (tn+1 , u
g n+1 )
(22)

L3PhysNum (EUPI) Équations différentielles 2020 16 / 30

Méthode Runge Kutta ordre 2, RK2

RK2
On utilise les différences centrées au 1er ordre

yi+1 − yi = hf (yi+1/2 , ti+1/2 ) (23)

h
Estimation de yi+1/2 par y^
i+1/2 = yi + 2 f (yi , ti )
Le schéma Runge Kutta d’ordre 2 est alors :
(
u0 = y0
u^
i+1/2 = ui + h2 f (ui , ti ) (24)
ui+1 = ui + hf (u^ i+1/2 , ti+1/2 )

L3PhysNum (EUPI) Équations différentielles 2020 17 / 30

Méthode Runge Kutta ordre 4, RK4

Z tn+1
h 
f (t, y )dt ≈ f (tn , yn ) + 4f (tn+1/2 , yn+1/2 ) + f (tn+1 , yn+1 ) (25)
tn
6

On sépare le point milieu en 2 :

Z tn+1
h 
f (t, y )dt ≈ f (tn , yn ) + 2f (tn+1/2 , yn+1/2 ) + 2f (tn+1/2 , yn+1/2 ) + f (tn+1 , yn+1 ) (26)
tn
6

en utilisant le même raisonnement que précédemment on obtient le schéma RK4 suivant :


 u0 = y0
 u^ = un + h2 f (tn , un )

 n+1/2
u\
n+1/2 = un + h2 f (tn+1/2 , u^
n+1/2 )


 ug
n+1 = un + hf(tn+1/2 , u\ n+1/2 ) 
un+1 = un + h6 f (tn , un ) + 2f (tn+1/2 , u^
n+1/2 ) + 2f (tn+1/2 , u\
n+1/2 ) + f (tn+1 , u n+1 )
 g
(27)

L3PhysNum (EUPI) Équations différentielles 2020 18 / 30

Méthode Runge Kutta ordre 4, RK4

Schéma RK4
Souvent on écrit la méthode RK4 en utilisant les coefficients k1 , k2 , k3 et k4 suivants :

k = hf (tn , un )
 1

k2 = hf (tn+1/2 , un + k21 )
(28)
 k3 = hf (tn+1/2 , un + k22 )

k4 = hf (tn+1 , un + k3 )
1
un+1 = un + (k1 + 2k2 + 2k3 + k4 ) (29)
6

L3PhysNum (EUPI) Équations différentielles 2020 19 / 30

RK4

Exercice

On reprend l’exercice précédent :


dy
= f (t, y ) = t(1 − y )
dt (30)
y (0) = 0

Appliquer le schéma RK4 sur l’intervalle [0; 3], avec un pas h = 0, 5. Pour cela, faites un tableau
dont la dernière colonne correspondra à la valeur absolue de l’erreur, sachant que la solution
t2
exacte est y (t) = 1 − e − 2
Compléter le tableau suivant.

n tn k1 k2 k3 k4 un yexacte |erreur |
0 0
1 0,5
2 1,0
3 1,5
4 2,0
5 2,5
6 3,0

L3PhysNum (EUPI) Équations différentielles 2020 20 / 30

RK4 - correction

n tn k1 k2 k3 k4 un yexacte |erreur |
0 0 0 0
1 0,5 0,000 0,125 0,117 0,221 0,1174 0,117
2 1,0 0,221 0,289 0,277 0,303 0,393 0,393
3 1,5 0,3035 0,284 0,290 0,237 0,675 0,675
4 2,0 0,244 0,178 0,206 0,118 0,863 0,865
5 2,5 0,136 0,077 0,111 0,033 0,954 0,956
6 3,0 0,058 0,023 0,047 -1,769e-3 0,986 0,989
l’erreur est inférieure à 10−3 pour n < 4 et de l’ordre de 10−3 pour n ≥ 4.

L3PhysNum (EUPI) Équations différentielles 2020 21 / 30

Algorithme - Programmation
Reprenons le schéma de Euler progressif (11) :

A chaque itération, on ajoute à la liste

Exemple Euler explicite X et Y les valeurs de t et y
respectivement.
Les listes X et Y sont retournées.
Variables :
entrées :
Possibilité de les convertir en nd.array.
y0 , t0 réel, C.I. Plusieurs façons d’écrire en python cet
Tf : borne supérieure de l’intervalle algo.
N : nombre de subdivision de l’intervalle (ou h)
f : la fonction
locales : t, y , h réels
sorties : X et Y listes des points solutions sur l’intervalle
Début :
X ← [t0 ]
Y ← [y0 ]
t ← t0
y ← y0
Tf −t 0
h← N
tant que t ≤ Tf Faire
y ← y + hf (t, y )
t ←t+h
X ← ajout de x
Y ← ajout de y
fin tant que
retourner X et Y
Fin
Algorithme 1 : Euler progressif

Faites de même avec les autres schémas numériques.

L3PhysNum (EUPI) Équations différentielles 2020 22 / 30
Algorithme - Programmation

Pour utiliser les équations (6) et (7), c’est-à-dire pour des équations différentielles d’ordre p ou
pour des systèmes de p équations différentielles du 1er ordre, les mêmes schémas numériques
précédents sont utilisés sous forme vectorielle. Pour cela les propriétés du module numpy sont
nécessaires.

Équation différentielle d’ordre p

Rappel :  −
→0 −
→ →
y = F (t, −
y)
−−→ − (31)
y (a) = →
α
avec −
→
y = (y0 , y1 , . . . , yp−1 )  
y1
y2
−
→ −
F (t, →
 
y)= ...  (32)
yp−1
 
f (t, −
→y)

L3PhysNum (EUPI) Équations différentielles 2020 23 / 30

Algorithme - Programmation

Équation différentielle d’ordre p -suite-

Le système est alors :
   
y0 y1
y1 y2
d    
 ... = ...  (33)
dt  yp−2
  yp−1

yp−1 f (y0 , y1 , . . . , yp−1 , t)

avec p conditions initiales −

→
α = (y (a), y 0 (a), . . . , y (p−1) (a)) = (y0 (a), y1 (a), . . . , yp−1 (a))

algo Euler
 →
− −
→ −→
dY
= F ( Y , t)
−
→dt (34)
Y0 (C.I.)

L3PhysNum (EUPI) Équations différentielles 2020 24 / 30

Algorithme - Programmation

Exemple Euler explicite vectorielle

 −
→ −−− →
Y0 = y (t0 )
−−→ − → −
→ − → (35)
Yn+1 = Yn + h F (Yn , tn )
Variables :
entrées→−:
F : ≡ F tableau (nd.array) (y1 , y2 , yp−1 , f (y0 , . . . , t))
x et y C.I. →
−
xf : borne supérieure de l’intervalle algo F
h : pas de calcul
sorties : X et Y listes des points solutions sur l’intervalle Variables :
entrées :
Début : t et y
X ← [] sortie :
Y ← [] g : vecteur (nd.array)
X ←ajout x Début :
Y ← ajout y g[0] ← y [1]
tant que x ≤ xf Faire g[1] ← y [2]
...
y ← y + hF (t, y )
g[p − 1] ← f (y0 , y1 , . . . , yp−1 , t)
x ←x +h retourner g
X ← ajout de x Fin
Y ← ajout de y
fin tant que
retourner X et Y
Fin
Algorithme 2 : Euler progressif vectorielle
Appel Euler
A, B = euler (F , t, y , tfin, h)
avec t et y les c.i.
A et B solution (A échantillonnage sur t
B = (B[0] ≡ y la solution, B[1] ≡ y 0 , . . .)
L3PhysNum (EUPI) Équations différentielles 2020 25 / 30
Méthode d’Adams - Basforth

Adams-Bashforth
A partir du système 
dy
= f (y , t)
dt (36)
y (t0 ) = y0

Utilisons le développement en série de Taylor au point t = ti :

h2 00 h3 000
yi+1 = yi + hyi0 + y + y + ... (37)
2 i 3! i
avec y 0 (ti ) = yi0 = f (yi , ti ) et y 00 (ti ) = yi00 = f 0 (yi0 , ti ) etc

h2 0 h3 00
yi+1 = yi + hfi + f + f + ... (38)
2 i 3! i

Ordre 1
yi+1 = yi + hfi (39)
On retrouve la méthode d’Euler

L3PhysNum (EUPI) Équations différentielles 2020 26 / 30

Méthode d’Adams - Basforth

Ordre 2
h2 0
yi+1 = yi + hfi + f + O(h)3 (40)
2 i
f −f
Utilisons les différences à gauche pour la dérivée :fi0 = i hi−1 + O(h)
On obtient alors
h2 fi − fi−1
 
yi+1 = yi + hfi + + O(h) + O(h)3
2 h
h
yi+1 = yi + (3fi − fi−1 ) + O(h)3 (41)
2
Formule d’ordre 2 puisque l’on néglige les termes en h3
Pour calculer yi+1 il faut connaître fi et fi−1 .
Pour lancer le processus il faut connaître y0 et y1 .
y0 Ok c.i.
y1 on le calcule par une méthode vue précédemment du même ordre, ici on prendra RK2

L3PhysNum (EUPI) Équations différentielles 2020 27 / 30

Méthode d’Adams - Basforth

Ordre 3
h2 0 h3 00
yi+1 = yi + hfi + fi + f + O(h)4
2 3! i
On utilise les différences à gauche pour fi0 et fi00

fi − fi−1 h
fi0 = + fi00 + O(h)2
h 2

fi − 2fi−1 + fi−2
fi00 = + O(h)
h2
En reportant dans la relation précédente on trouve finalement :

h
yi+1 = yi + (23fi − 16fi−1 + 5fi−2 ) + O(h)4 (42)
12
Le calcul de y1 et y2 par RK4 (ordre 4), puis le processus de l’équation (42) peut être réalisé.

L3PhysNum (EUPI) Équations différentielles 2020 28 / 30

Méthode d’Adams - Basforth

Ordre n > 3
On calcul de la même façon et on obtient :
n
X
yi+1 = yi + h βnk fi−k + O(h)n+2 (43)
k=0

La formule la plus couramment utilisée est celle d’ordre 4

h
yi+1 = yi + (55fi − 59fi−1 + 37fi−2 − 9fi−3 ) + O(h)5 (44)
24

L3PhysNum (EUPI) Équations différentielles 2020 29 / 30

Méthode d’Adams - Basforth

k
0 1 2 3 4 5
n

0 1

3 −1
1 2 2

23 −16 5
2 12 12 12

55 −59 37 −9
3 24 24 24 24

1901 −2774 2616 −1274 251

4 720 720 720 720 720

4277 −7923 9982 −7298 2877 −475

5 1440 1440 1440 1440 1440 1440

Table – Coefficients βnk de l’équation (43) jusqu’à l’ordre 6 pour la méthode Adams-Bashforth

L3PhysNum (EUPI) Équations différentielles 2020 30 / 30