UPMC Master P&A/SDUEE

UE MU4PY209
Méthodes Numériques et Calcul Scientifique

Résolution numérique des

équations différentielles ordinaires (EDO)

2019–2020 [email protected]
EDO TABLE DES MATIÈRES TABLE DES MATIÈRES

Table des matières

1 Introduction 6

1.1 Problème différentiel . . . . . . . . . . . . . . . . . . . . . . . . . 6

1.2 Deux types de problèmes différentiels à résoudre . . . . . . . . . . . 7

1.3 Équations différentielles scalaires du 1er ordre . . . . . . . . . . . . 8

1.4 Unicité et problème bien posé : conditions suffisantes . . . . . . . . . 9

1.5 Méthodes de résolution numérique et notations . . . . . . . . . . . . 10

2 Méthodes à un pas 12

2.1 Méthodes du premier ordre . . . . . . . . . . . . . . . . . . . . . . 13

2.1.1 Méthode d’Euler progressive (explicite) . . . . . . . . . . . . 13

2.1.2 Méthode d’Euler rétrograde (implicite) . . . . . . . . . . . . . 15

MNCS 1 2019-2020
EDO TABLE DES MATIÈRES TABLE DES MATIÈRES

2.2 Méthodes du deuxième ordre . . . . . . . . . . . . . . . . . . . . . 18

2.2.1 Méthode du point milieu . . . . . . . . . . . . . . . . . . . . 18

2.2.2 Méthode d’Euler modifiée . . . . . . . . . . . . . . . . . . . 20

2.2.3 Méthode de Heun . . . . . . . . . . . . . . . . . . . . . . . 23

2.3 Méthodes de Runge Kutta . . . . . . . . . . . . . . . . . . . . . . 24

2.3.1 Méthode de Runge Kutta d’ordre 3 . . . . . . . . . . . . . . 24

2.3.2 Méthode de Runge Kutta d’ordre 4 . . . . . . . . . . . . . . 26

2.4 Erreur absolue en fonction du pas et de l’ordre . . . . . . . . . . . . 27

2.5 Exemple de l’équation logistique . . . . . . . . . . . . . . . . . . . 28

2.5.1 Exemple d’erreur totale maximale en simple précision . . . . . 31

2.5.2 Exemple d’erreur totale maximale en double précision . . . . . 32

2.5.3 Comparaison des erreurs maximales simple/double précision . 33

MNCS 2 2019-2020
EDO TABLE DES MATIÈRES TABLE DES MATIÈRES

3 Méthodes à plusieurs pas 34

3.1 Méthodes d’Adams . . . . . . . . . . . . . . . . . . . . . . . . . . 34

3.1.1 Adams Bashforth : explicite, pas de terme en f (ti+1 , ui+1 ) . 35

3.1.2 Adams Moulton : implicite, terme en f (ti+1 , ui+1 ) . . . . . . 36

3.1.3 Comparaison méthodes à un pas et Adams explicite . . . . . . 38

3.1.4 Méthodes de prédicteur correcteur . . . . . . . . . . . . . . 39

3.2 Méthodes adaptatives . . . . . . . . . . . . . . . . . . . . . . . . 40

3.2.1 Exemple : méthode de Runge Kutta Fehlberg . . . . . . . . . 41

3.3 Méthodes d’extrapolation de Gragg . . . . . . . . . . . . . . . . . . 42

3.3.1 Principe de l’extrapolation . . . . . . . . . . . . . . . . . . . 42

3.3.2 Comparaison méthodes à un pas et extrapolation de Gragg . . 45

4 Les EDO du premier ordre en pratique 46

MNCS 3 2019-2020
EDO TABLE DES MATIÈRES TABLE DES MATIÈRES

4.1 Échelles de temps et problèmes raides . . . . . . . . . . . . . . . . 46

4.2 Validation des résultats . . . . . . . . . . . . . . . . . . . . . . . . 47

4.3 Structure des programmes de résolution d’EDO du 1er ordre . . . . . 48

5 Systèmes d’EDO du 1er ordre 49

5.1 Méthodes scalaires explicites . . . . . . . . . . . . . . . . . . . . . 49

5.2 Équations de Lotka-Volterra . . . . . . . . . . . . . . . . . . . . . 52

6 Équations différentielles d’ordre supérieur 58

6.1 Exemple . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59

6.2 Exemple d’EDO d’ordre 2 : le pendule . . . . . . . . . . . . . . . . 60

7 Implémentation vectorielle 69

7.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69

MNCS 4 2019-2020
EDO TABLE DES MATIÈRES TABLE DES MATIÈRES

7.2 En fortran (norme 2003) . . . . . . . . . . . . . . . . . . . . . . . 70

7.3 En C89 avec des tableaux dynamiques . . . . . . . . . . . . . . . . 73

7.4 En C99 avec des tableaux automatiques . . . . . . . . . . . . . . . 76

Bibliographie 79

MNCS 5 2019-2020
EDO 1 Introduction

1 Introduction

1.1 Problème différentiel

— équation différentielle scalaire d’ordre n

dn y dn−1 y
 
dy
=f t, y, , ...,
dtn dt dtn−1

où f est la fonction second membre donnée

⇒ famille de solutions y(t) à n paramètres
— ensemble de n conditions imposées
⇒ choix d’ une solution dans la famille

MNCS 6 2019-2020
EDO 1 Introduction 1.2 Deux types de problèmes différentiels à résoudre

1.2 Deux types de problèmes différentiels à résoudre

— Conditions initiales données pour une seule valeur t0 de t, par exemple

(n−1)
y(t0 ) = y0 , y 0 (t0 ) = y00 , . . . , y (n−1) (t0 ) = y0

Problème de conditions initiales ou de Cauchy

— Conditions données pour des valeurs distinctes de la variable indépendante t,

par exemple :

y(t0 ) = y0 , y(t1 ) = y1 , . . . , y(tn−1 ) = yn−1

Problème de conditions aux limites (non traité, sauf problème de tir).

MNCS 7 2019-2020
EDO 1 Introduction 1.3 Équations différentielles scalaires du 1er ordre

1.3 Équations différentielles scalaires du 1er ordre

Étudier d’abord les équations différentielles scalaires du premier ordre.

⇒ famille de solutions y(t) à un paramètre (y0 )

dy
= f (t, y(t)) avec y(t0 ) = y0 condition initiale
dt

Les EDO d’ordre supérieur se ramènent à des systèmes différentiels couplés du

premier ordre (EDO vectorielles du premier ordre).

MNCS 8 2019-2020
EDO 1 Introduction 1.4 Unicité et problème bien posé : conditions suffisantes

1.4 Unicité et problème bien posé : conditions suffisantes

La condition de Lipschitz

|f (t, y2 ) − f (t, y1 )| 6 K |y2 − y1 |

assure l’unicité de la solution.

 
 ∂f 
 (t, y) 6 K dans un domaine convexe ⇒ condition de Lipschitz vérifiée.
 ∂y 
Les erreurs d’arrondi amènent à toujours résoudre un problème perturbé.

Problème bien posé si : le problème faiblement perturbé (second membre ou

condition initiale) possède une solution proche de celle du problème original.
La condition de Lipschitz assure que le problème est bien posé.

MNCS 9 2019-2020
EDO 1 Introduction 1.5 Méthodes de résolution numérique et notations

1.5 Méthodes de résolution numérique et notations

Résolution numérique approchée sur l’intervalle [t0 , t0 + L] de longueur L

Discrétisation par découpage de l’intervalle de longueur L selon un pas constant h

Échantillonnage de la solution aux instantsti = t0 + ih pour 1 6 i 6 n.

Solution numérique : ui = approximation de y(ti )

À partir de la condition initiale u0 = y(t0 ) imposée,

faire une boucle sur les abscisses ti pour calculer l’approximation ui+1 à ti+1
→ approximer ainsi de proche en proche la solution sur l’intervalle L.
⇒ accumulation des erreurs dans la boucle
À chaque pas de la boucle, pour calculer ui+1 , on peut s’appuyer :
— sur la dernière valeur calculée ui : méthodes à un pas
— sur plusieurs valeurs ui−k (k > 0) antérieurement calculées :
méthodes à plusieurs pas (initialisation nécessaire par méthode à un pas)

MNCS 10 2019-2020
EDO 1 Introduction 1.5 Méthodes de résolution numérique et notations

Méthode à pas
y exacte constant
famille de solutions exactes approx.
dépendant de y0 Découpage de
l’intervalle de
yi+1 erreur
longueur L selon
ui+1 erreur locale cumulée
un pas fixe
h = L/n.
ui la solution
passant par (ti , ui ) ui = approximat.
de y(ti )
h
y0
Un pas :
ti → ti+1
t0 ti ti+1 t ui → ui+1
MNCS 11 2019-2020
EDO 2 Méthodes à un pas

2 Méthodes à un pas
Constituent l’algorithme de base qui permet d’estimer la valeur de la solution à
l’instant ti+1 = ti + h, connaissant seulement ui , celle à ti .
La valeur à estimer peut être approchée par un développement limité de Taylor :

dy h2 d2 y
y(ti + h) = y(ti ) + h (ti ) + 2
(ti ) + · · · (1)
dt 2 dt
Ordre n de la méthode = plus grande puissance de h prise en compte dans
l’approximation.
— Somme des termes négligés = erreur de troncature locale ∝ hn+1
déterministe, augmente si le pas h augmente et si l’ordre de la méthode diminue
— Précision finie des opérations sur les réels ⇒ erreur d’arrondi aléatoire
augmente lorsque les calculs se compliquent, en particulier si le pas h diminue.
Indépendamment du coût (en temps de calcul) des opérations, et des cas où la
fonction est tabulée, ne pas croire que diminuer le pas améliore toujours la
qualité du résultat : un compromis doit être trouvé entre ces deux types d’erreurs.

MNCS 12 2019-2020
EDO 2 Méthodes à un pas 2.1 Méthodes du premier ordre

2.1 Méthodes du premier ordre

2.1.1 Méthode d’Euler progressive (explicite)

Méthode du premier ordre d’intérêt pédagogique, à éviter en pratique

ui+1 = ui + hf (ti , ui ) (2)

Exemple : stabilité

dy y
=− ⇒ solution analytique y = y0 e−t/τ ⇒ yn = y0 (e−h/τ )n
dt τ
h n
ui+1 = ui − ui ⇒ solution numérique un = y0 (1 − h/τ )
τ
Si τ> 0, la solution exacte vérifie y(∞) = 0,
Mais pour l’approximation, un → 0 ⇐⇒ |1 − h/τ | < 1 ⇐⇒ 0 < h < 2τ .
Condition de stabilité : h < 2τ (pas h petit)
Mais, si h > τ , alors (1 − h/τ ) < 0 : alternance de signe de la solution un .

MNCS 13 2019-2020
EDO 2 Méthodes à un pas 2.1 Méthodes du premier ordre

y Méthode d’Euler

Méthode explicite qui ne

nécessite qu’une
seule évaluation
ui+1
de la fonction second
k1 h membre f par pas :

k1 k1 = f (ti , ui )
ui facilement instable
h
ui+1 − ui
= f (ti , ui )
ti ti+1 t h
voir dérivée avant

MNCS 14 2019-2020
EDO 2 Méthodes à un pas 2.1 Méthodes du premier ordre

2.1.2 Méthode d’Euler rétrograde (implicite)

ui+1 = ui + hf (ti+1 , ui+1 ) (3)

Méthode implicite : résolution itérative , plus difficile à mettre en œuvre, sauf si la

forme de f (t, u) permet le calcul analytique de ui+1 à partir de l’équation (3).
Avantage : meilleure stabilité que la méthode progressive explicite.

Exemple : stabilité

dy y
=− ⇒ solution analytique y = y0 e−t/τ ⇒ yn = y0 (e−h/τ )n
dt τ
h ui
ui+1 = ui − ui+1 ⇒ solution numérique ui+1 =
τ 1 + h/τ
y0
un = n
(1 + h/τ )
Si τ > 0, y(∞) = 0, et aussi un → 0 ∀τ > 0, ∀h > 0 solution stable

MNCS 15 2019-2020
EDO 2 Méthodes à un pas 2.1 Méthodes du premier ordre

Mise en œuvre de la méthode d’Euler rétrograde : résolution de l’équation

implicite par itération

ui+1 = ui + hf (ti+1 , ui+1 )

Itérer l’application g pour rechercher son point fixe

v20 = g(v2 ) = ui + hf (t2 , v2 )

Ce point fixe est la solution de l’équation implicite.

— Utilise plusieurs évaluations du second membre, sans calcul de ses dérivées.

— Très peu d’itérations nécessaires

Initialisation par le prédicteur avec Euler progressif

t2 = ti + h
k1 = f (ti , ui )
v2 = ui + hk1

MNCS 16 2019-2020
EDO 2 Méthodes à un pas 2.1 Méthodes du premier ordre

Boucle pour recherche du point fixe de g(v2 ) = v20 = ui + hf (t2 , v2 )

k2 = f (t2 , v2 )
v20 = ui + hk2
δv2 = v20 − v2
arrêt si |δv2 |2 6 α2 |v2 |2
v2 = v20
 
∂f 
La fonction g est contractante si g 0 (v2 ) = h  ∂v 61

2
ce qui est vérifié si le pas est assez faible.
 
 ∂f 
Rappel : la condition de Lipschitz est  ∂v
2
6K
Critère d’arrêt : choisir α faible, mais α > ε.

MNCS 17 2019-2020
EDO 2 Méthodes à un pas 2.2 Méthodes du deuxième ordre

2.2 Méthodes du deuxième ordre

Première idée : augmenter le nombre de termes du développement de Taylor :

rarement utilisé, car nécessite l’évaluation des dérivées partielles de f .

dy d2 y ∂f ∂f dy ∂f ∂f
= f (t, y(t)) ⇒ = + = +f (6)
dt dt2 ∂t ∂y dt ∂t ∂y
Préférer utiliser plusieurs évaluations du second membre f en des points adaptés.

Centrer l’évaluation de la dérivée au point milieu tm = (ti + ti+1 )/2.

h dy 1 h2 d2 y 3
y(ti + h) = y(tm ) + (tm ) + 2
(t m ) + O(h ) (7a)
2 dt 2 4 dt
h dy 1 h2 d2 y 3
y(ti ) = y(tm ) − (tm ) + (t m ) + O(h ) (7b)
2 dt 2 4 dt2
Par différence, (approximation locale parabolique, voir aussi dérivée centrée à 2 termes)
dy
y(ti + h) − y(ti ) = h (tm ) + O(h3 )
dt

MNCS 18 2019-2020
EDO 2 Méthodes à un pas 2.2 Méthodes du deuxième ordre

2.2.1 Méthode du point milieu

Nécessite l’évaluation du second membre f en 2 points :

en (ti , ui ) et au milieu (ti+1/2 = ti + h/2, ui+1/2 ) d’un pas (hors grille).
 
h h
ui+1 = ui + hf ti + , ui + f (ti , ui )
2 2

k1 = f (ti , ui ) (8a)

h h
(ui+1/2 calculé via Euler) k2 = f (ti + , u i + k1 ) (8b)
2 2
ui+1 = ui + hk2 (8c)

MNCS 19 2019-2020
EDO 2 Méthodes à un pas 2.2 Méthodes du deuxième ordre

ui+1 Méthode du
k2 point milieu

k2 h Méthode explicite

k2 qui nécessite deux

évaluations du second
k1 h/2
k1 membre par pas dont
ui
une hors grille.

ti ti + h/2 ti+1 t

MNCS 20 2019-2020
EDO 2 Méthodes à un pas 2.2 Méthodes du deuxième ordre

2.2.2 Méthode d’Euler modifiée

En appliquant 7a et 7b à la dérivée et en faisant la somme, on peut remplacer la

dérivée au milieu par la moyenne des dérivées aux extrémités de l’intervalle
(voir méthode de quadrature dite des trapèzes) :
dy dy dy
(ti ) + (ti+1 ) = 2 (tm ) + O(h2 )
dt dt dt
D’où une approximation n’utilisant pas la valeur de f au point milieu tm :
h
ui+1 = ui + [f (ti , ui ) + f (ti+1 , ui+1 )]
2
De nouveau, méthode a priori implicite, plus stable, mais plus lourde.
⇒ Contournement du problème en utilisant l’approximation d’Euler explicite (voir 2)
pour évaluer ui+1 intervenant dans f .

h
ui+1 = ui + [f (ti , ui ) + f (ti+1 , ui + hf (ti , ui ))]
2

MNCS 21 2019-2020
EDO 2 Méthodes à un pas 2.2 Méthodes du deuxième ordre

Bilan : méthode de type prédicteur-correcteur équivalent à

— un demi-pas avec la pente initiale k1

— et un demi-pas avec la pente k2 du point prédit par Euler progressif.

k1 = f (ti , ui ) (9a)

k2 = f (ti+1 , ui + k1 h) (9b)

h
ui+1 = ui + [k1 + k2 ] (9c)
2
Remarques

— deuxième ordre comme point milieu mais sans évaluation hors grille
— la résolution de l’équation implicite peut se faire en itérant la correction jusqu’à
ce qu’elle devienne négligeable.

MNCS 22 2019-2020
EDO 2 Méthodes à un pas 2.2 Méthodes du deuxième ordre

corrigé Méthode d’Euler

ui+1 modifiée
k2 h/2 k2 k2
Méthode explicite qui
prédit
nécessite deux éva-
k1 h/2 luations de la fonction
k1
ui par pas en des points
de la grille.
ti ti + h/2 ti+1 t

MNCS 23 2019-2020
EDO 2 Méthodes à un pas 2.3 Méthodes de Runge Kutta

2.2.3 Méthode de Heun

k1 = f (ti , ui ) (10a)
2 2
k2 = f (ti + h, ui + k1 h) (10b)
3 3
h
ui+1 = ui + [k1 + 3k2 ] (10c)
4

2.3 Méthodes de Runge Kutta

Plus généralement, avec r évaluations de f , on peut atteindre une méthode d’ordre

r si r 6 4. Pour atteindre l’ordre 5, six évaluations sont nécessaires.
=⇒ la méthode de Runge Kutta d’ordre 4 est très utilisée.
Les méthodes de Runge-Kutta sont stables.

MNCS 24 2019-2020
EDO 2 Méthodes à un pas 2.3 Méthodes de Runge Kutta

2.3.1 Méthode de Runge Kutta d’ordre 3

k1 = f (ti , ui ) (11a)

h h
k2 = f (ti + , ui + k1 ) (11b)
2 2
k3 = f (ti + h, ui + (2k2 − k1 )h) (11c)

h
ui+1 = ui + (k1 + 4k2 + k3 ) (11d)
6

MNCS 25 2019-2020
EDO 2 Méthodes à un pas 2.3 Méthodes de Runge Kutta

2.3.2 Méthode de Runge Kutta d’ordre 4

k1 = f (ti , ui ) (12a)

h h
k2 = f (ti + , ui + k1 ) (12b)
2 2
h h
k3 = f (ti + , ui + k2 ) (12c)
2 2
k4 = f (ti + h, ui + k3 h) (12d)

h
ui+1 = ui + (k1 + 2k2 + 2k3 + k4 ) (12e)
6

MNCS 26 2019-2020
EDO 2 Méthodes à un pas 2.4 Erreur absolue en fonction du pas et de l’ordre

2.4 Erreur absolue en fonction du pas et de l’ordre

nombre de pas = L/h =⇒ erreur globale ∼ erreur locale × L/h

TABLE 1 – Erreur de troncature seule

Méthode ordre erreur locale erreur globale

Euler explicite 1 ∝ h2 ∝h
Point milieu – Euler modifiée 2 ∝ h3 ∝ h2
Runge-Kutta 3 3 ∝ h4 ∝ h3
Runge-Kutta 4 4 ∝ h5 ∝ h4

Erreur d’arrondi locale indépendante de h ⇒ erreur d’arrondi globale ∝ 1/h

MNCS 27 2019-2020
EDO 2 Méthodes à un pas 2.5 Exemple de l’équation logistique

dy
2.5 Exemple de l’équation logistique = y(1 − y/2)
dt
sol_exacte_logistique
2.0 [ 0.0e+00 , 2.0e+01] h= 2.00e-02 y0=1.00e-01 n= 1001

1.5
logistique

1.0

0.5

0.0
0 5 10 15 20
temps
F IGURE 1 – Solution analytique de l’équation logistique pour t0 = 0, y(t0 ) = 0.1,
a = 1 et k = 2.
MNCS 28 2019-2020
EDO 2 Méthodes à un pas 2.5 Exemple de l’équation logistique

fichier euler.dat ordre 1 erreur max 0.830E−02 position de l’err max 0.296E+01

0.001

0.000

−0.001

Erreur = Estimé − Exact −0.002

−0.003

−0.004

−0.005

−0.006

−0.007

−0.008

−0.009
0 2 4 6 8 10 12 14 16 18 20
estimé − exact
temps

F IGURE 2 – Erreur dans l’intégration de l’équation logistique avec la méthode d’Euler

pour h = 0, 02. L’allure régulière montre que l’erreur de troncature domine.
Erreur de troncature locale liée à la courbure de la solution

MNCS 29 2019-2020
EDO 2 Méthodes à un pas 2.5 Exemple de l’équation logistique

fichier rk4.dat ordre 4 erreur max 0.954E−06 position de l’err max 0.190E+01

1e−06

8e−07

6e−07

Erreur = Exact −Estimé

4e−07

2e−07

0e+00

−2e−07

−4e−07

−6e−07

−8e−07

−1e−06
0 2 4 6 8 10 12 14 16 18 20
exact − estimé
temps

F IGURE 3 – Erreur dans l’intégration de l’équation logistique avec Runge Kutta

d’ordre 4 pour h = 0, 02. L’allure bruitée est caractéristique de l’erreur d’arrondi
et on retrouve les niveaux de quantification des réels sur 32 bits.

MNCS 30 2019-2020
EDO 2 Méthodes à un pas 2.5 Exemple de l’équation logistique

2.5.1 Exemple d’erreur totale maximale en simple précision (32 bits)

Valeur absolue de l’erreur absolue pour y’= y(1−y/2) sur [0,20] y(0)=0.1 simple précision
100
euler sp
milieu sp
rk3 sp
10
−1 rk4 sp pente +1

10−2
erreur maximale

pente +2
−3
10
pente −1

−4
10 pente +3

10−5

pente +4
−6
10

10−7 −4
10 10−3 10−2 10−1 100
pas
MNCS 31 2019-2020
EDO 2 Méthodes à un pas 2.5 Exemple de l’équation logistique

2.5.2 Exemple d’erreur totale maximale en double précision (64 bits)

Erreur absolue y’= y(1−y/2) sur [0,20] y(0)=0.1 double précision

100

10−2 pente +1

10−4

pente +2
10−6
erreur maximale

10−8
pente +3

10−10

10−12 pente +4

−14 euler dp
10
milieu dp
rk3 dp
rk4 dp
10−16 −4
10 10−3 10−2 10−1 100
pas
MNCS 32 2019-2020
EDO 2 Méthodes à un pas 2.5 Exemple de l’équation logistique

2.5.3 Comparaison des erreurs maximales simple/double précision

Erreur absolue y’= y(1−y/2) sur [0,20] y(0)=0.1 simple/double précision

100

10−2

10−4

10−6
erreur maximale

10−8

10−10

euler sp
10−12 milieu sp
rk3 sp
rk4 sp
−14 euler dp
10
milieu dp
rk3 dp
rk4 dp
10−16 −4
10 10−3 10−2 10−1 100
pas
MNCS 33 2019-2020
EDO 3 Méthodes à plusieurs pas

3 Méthodes à plusieurs pas

3.1 Méthodes d’Adams

Principe : les erreurs augmentent avec l’intégration, les points les plus proches de
la valeur initiale ont tendance à être plus fiables. Pour calculer ui+1 , on peut
s’appuyer non seulement sur la dernière valeur estimée ui , mais sur les m précédentes.

— si le calcul invoque la pente au point recherché f (ti+1 , ui+1 ), la méthode est

implicite (voir 2.1.2) : A DAMS -M OULTON

— sinon elle est explicite : A DAMS -B ASHFORTH .

Dans les deux cas, il faut initialiser le calcul par une méthode à un pas sur les m
premiers points.

Le calcul réutilise les évaluations antérieures du second membre.

=⇒ stocker ces valeurs pour économiser les calculs

MNCS 34 2019-2020
EDO 3 Méthodes à plusieurs pas 3.1 Méthodes d’Adams

3.1.1 Adams Bashforth : explicite, pas de terme en f (ti+1 , ui+1 )

Adams-Bashforth utilise m points à gauche de ti+1 :

m
X m
X
ui+1 = ui + βh αj f (ti−j+1 , ui−j+1 ) avec β αj = 1
j=1 j=1

méthode d’ordre m, mais une seule nouvelle évaluation du second membre à

chaque pas
m β α1 α2 α3 α4
1 1 1
2 1/2 3 −1
3 1/12 23 −16 5
4 1/24 55 −59 37 −9

MNCS 35 2019-2020
EDO 3 Méthodes à plusieurs pas 3.1 Méthodes d’Adams

Exemple : Adams-Bashforth à deux pas

h
ui+1 = ui + [3f (ti , ui ) − f (ti−1 , ui−1 )]
2

y Adams-Bashforth à
deux pas
ui+1
−k−1 h/2
Méthode explicite
3k0 h/2
h/2 d’ordre deux mais
k0 seulement une éva-
ui luation nouvelle du

k−1 second membre à

ui−1 chaque pas
3h/2 ⇒ mémoriser les
ti−1 ti ti+1 t seconds membres.

MNCS 36 2019-2020
EDO 3 Méthodes à plusieurs pas 3.1 Méthodes d’Adams

3.1.2 Adams Moulton : implicite, terme en f (ti+1 , ui+1 )

Adams-Moulton utilise m + 1 points dont ti+1 donc implicite :

m
X m
X
ui+1 = ui + βh αj f (ti−j+1 , ui−j+1 ) avec β αj = 1
j=0 j=0

méthode d’ordre m + 1 (nombre d’évaluations du second membre)

m β α0 α1 α2 α3
0 1 1
1 1/2 1 1
2 1/12 5 8 −1
3 1/24 9 19 −5 1

Éviter les difficultés de l’implicite en utilisant un prédicteur explicite de ui+1 ,

injecté ensuite dans l’expression d’Adams-Moulton vue comme correcteur .

MNCS 37 2019-2020
EDO 3 Méthodes à plusieurs pas 3.1 Méthodes d’Adams

3.1.3 Comparaison méthodes à un pas et Adams explicite

Erreur absolue y’= y(1−y/2) sur [0,20] y(0)=0.1 simple précision

0
10
euler sp
milieu sp
rk3 sp
10
−1 rk4 sp
adams−bashforth2 sp
adams−bashforth3 sp
adams−bashforth4 sp
−2
10
erreur maximale

−3
10

10−4

10−5

10−6

10−7 −4
10 10−3 10−2 10−1 100
pas
MNCS 38 2019-2020
EDO 3 Méthodes à plusieurs pas 3.1 Méthodes d’Adams

3.1.4 Méthodes de prédicteur correcteur

Principe : bénéficier des qualités d’une méthode implicite mais l’appliquer à une
estimation obtenue par une méthode explicite du même ordre (voir Euler modifiée).

— prédiction de ui+1 par une méthode explicite

— correction de ui+1 par une formule implicite où f (ti+1 , y(ti+1 )) a été

approximé par la prédiction f (ti+1 , ui+1 ).

Exemple : méthode d’Euler modifiée

Une itération de la partie correction est possible.
L’ordre est celui du correcteur, mais la stabilité dépend plus du prédicteur.
Ces méthodes permettent d’estimer l’erreur de troncature à partir de la différence
entre prédicteur et correcteur =⇒ adaptation du pas

MNCS 39 2019-2020
EDO 3 Méthodes à plusieurs pas 3.2 Méthodes adaptatives

3.2 Méthodes adaptatives

Principe : ajuster le pas localement pour obtenir une précision imposée.

Estimer l’erreur de troncature par l’écart entre deux solutions numériques :

y(ti+1 ) − ui+1
ui , d’ordre n ηi = ∝ hn (13a)
h
? y(ti+1 ) − u?i+1
u?i , d’ordre n + 1 ηi = ∝ hn+1 (13b)
h
u?i+1 − ui+1
or h  1 donc ηi ≈ ∝ hn (13c)
h
Modification du pas d’un facteur q

n qn ?
ηi (hq) ≈ q ηi (h) ≈ (ui+1 − ui+1 ) (14)
h
Pour obtenir une précision globale ∆y ≈ Lδ , imposer localement :

|ηi (hq)| ≈ δ (15)

MNCS 40 2019-2020
EDO 3 Méthodes à plusieurs pas 3.2 Méthodes adaptatives

d’où le facteur q à appliquer au pas :

 1/n
hδ
q≈ (16)
u?i+1 − ui+1
— Si q < 1, refuser ui+1 et diminuer le pas
— Si q > 1, accepter ui+1 et augmenter le prochain pas

B En pratique, limiter le pas à un intervalle raisonnable et éviter des variations brutales.

On suppose ici que le pas choisi permet de négliger l’erreur d’arrondi.

3.2.1 Exemple : méthode de Runge Kutta Fehlberg

Une méthode de Runge Kutta d’ordre 5 (6 évaluations de f )

permet de contrôler la précision obtenue par un Runge Kutta d’ordre 4
utilisant les évaluations de f aux mêmes points que celle d’ordre 5
(les poids ne sont pas ceux de la méthode d’ordre 4 classique).

MNCS 41 2019-2020
EDO 3 Méthodes à plusieurs pas 3.3 Méthodes d’extrapolation de Gragg

3.3 Méthodes d’extrapolation de Gragg

3.3.1 Principe de l’extrapolation

t0 h t0 + L
Subdivisions successives
k=1
h/p2
Découpage de l’intervalle de pas
grossier h en sous-intervalles
k=2
h/p3 de pas fin hk = h/pk de plus
en plus petits.
k=3 Développement polynomial de
l’erreur de troncature en fonction
pas fin : h/pkm
du pas pour extrapoler au pas nul
pkm pas k = km (hk → 0).

MNCS 42 2019-2020
EDO 3 Méthodes à plusieurs pas 3.3 Méthodes d’extrapolation de Gragg

Exemple de la méthode d’Euler : l’erreur de troncature globale sur un pas

grossier h est du premier ordre en fonction du pas fin h/pk .
Par exemple, pour les subdivisions p1 et p2 :
   2
h h
y(t + h) = v1 + a1 + a2 + ··· (17)
p1 p1
   2
h h
y(t + h) = v2 + a1 + a2 + ··· (18)
p2 p2
Combinaison linéaire de ces deux estimateurs v1 et v2
⇒ éliminer le terme d’ordre 1 de l’erreur (a1 inconnu)
⇒ nouvel estimateur w2,2 d’ordre 2 tel que :
 2
h
y(t + h) = w2,2 + b2 + ···
p2
(p2 /p1 ) w2,1 − w1,1 w2,1 − w1,1
w2,2 = = w2,1 + (19)
p2 /p1 − 1 p2 /p1 − 1

MNCS 43 2019-2020
EDO 3 Méthodes à plusieurs pas 3.3 Méthodes d’extrapolation de Gragg

Autre exemple avec la méthode du point milieu :

— l’erreur de troncature globale sur h est d’ordre 2 en fonction du pas h/pk ;

— pas de termes d’ordre impair dans le développement de l’erreur.

Combinaison linéaire de deux estimateurs avec des pas fins différents :

⇒ élimination du terme d’ordre 2 de l’erreur de troncature
⇒ erreur de troncature d’ordre 4
Itérer le processus avec une suite de subdivisions et de combinaisons linéaires
d’estimateurs pour augmenter l’ordre de l’erreur de troncature.
B Mais amélioration limitée par l’erreur d’arrondi...

L’écart wk+1,k+1 − wk,k donne une estimation de l’erreur de troncature si on

retient la solution wk,k . En ajustant, pour chaque intervalle de largeur h, le nombre
km de subdivisions pour respecter une erreur absolue imposée, on obtient une
version adaptative de la méthode de Gragg.

MNCS 44 2019-2020
EDO 3 Méthodes à plusieurs pas 3.3 Méthodes d’extrapolation de Gragg

3.3.2 Comparaison méthodes à un pas et extrapolation de Gragg

Erreur absolue y’= y(1−y/2) sur [0,20] y(0)=0.1 simple précision

0
10
euler sp
milieu sp
rk3 sp
10−1 rk4 sp
mil−mod sp
gragg2−4sp
gragg2−8sp
10−2
erreur maximale

10−3

10−4

10−5

10−6

10−7 −4
10 10−3 10−2 10−1 100
pas
MNCS 45 2019-2020
EDO 4 Les EDO du premier ordre en pratique

4 Les EDO du premier ordre en pratique

4.1 Échelles de temps et problèmes raides

Ne pas oublier que chaque problème différentiel possède une ou plusieurs échelles
de temps propres (périodes ou pseudo-périodes, constantes de temps).
La solution ne peut être représentée correctement qu’avec un pas assez inférieur
au plus petit de ces temps propres.
Cette analyse impose donc une valeur maximale pour le pas.

Certains problèmes différentiels qualifiés de raides comportent des échelles de

temps très différentes : leur intégration numérique s’avère délicate et coûteuse (pas
faible pour respecter le temps court, mais nombreux pour accéder au temps long).
Il existe des méthodes spécifiques des EDO raides qui ne sont pas présentées ici.

MNCS 46 2019-2020
EDO 4 Les EDO du premier ordre en pratique 4.2 Validation des résultats

4.2 Validation des résultats

Validation via une solution analytique d’un problème simplifié

Lorsqu’une solution analytique est disponible (par exemple pour certaines valeurs
de paramètres qui permettent de simplifier l’EDO), sa comparaison avec la solution
numérique permet de tester la méthode. Le calcul de l’erreur dans le domaine où la
troncature domine permet d’extrapoler l’effet d’un changement de pas connaissant
l’ordre de la méthode.

Validation sans solution analytique

Dans le cas où aucune solution analytique de référence n’est disponible, la

validation s’appuie sur les mêmes outils que les méthodes adaptatives :

— diminution du pas (division par 2)

— augmentation de l’ordre de la méthode
— extrapolation de Gragg
— calcul d’invariants (énergie par exemple)
MNCS 47 2019-2020
EDO 4 Les EDO du premier ordre 4.3 Structure des programmes de résolution d’EDO du 1er ordre
en pratique

4.3 Structure des programmes de résolution d’EDO du 1er ordre

1. un algorithme de base (appliquant une méthode d’ordre 1, 2, 3 ou 4 à la

fonction second membre f passée en argument) permettant d’avancer d’un
pas dans l’intégration de l’équation différentielle
2. éventuellement une procédure qui choisit le pas le plus grand possible
compatible avec la précision attendue et contrôle la progression de l’intégration
(elle pourrait comporter un algorithme adaptatif)
3. un programme d’interface avec l’utilisateur qui choisit la méthode, le second
membre, lit les paramètres (conditions initiales par ex.), déclenche et arrête la
boucle d’intégration et stocke les résultats.
4. un module comportant les fonctions seconds membres de l’équation
différentielle et les éventuelles solutions analytiques exactes ou approchées
5. un module d’utilitaires notamment pour écrire les résultats dans un fichier pour
visualisation ultérieure.

MNCS 48 2019-2020
EDO 5 Systèmes d’EDO du 1er ordre

5 Systèmes d’équations différentielles du 1er ordre

5.1 Extension des méthodes scalaires explicites aux vecteurs

Système de n équations différentielles couplées du considérer les

premier ordre associées à n conditions initiales vecteurs
#– et #–
y f.

dy1
   
= f1 (t, y1 , y2 , . . . , yn ) y1 f1
dt  
y 
 
f 
dy2  2  2
  et  
= f2 (t, y1 , y2 , . . . , yn )
dt . . .
 
. . .
 
... = ...
dyn yn fn
= fn (t, y1 , y2 , . . . , yn )
dt

MNCS 49 2019-2020
EDO 5 Systèmes d’EDO du 1er ordre 5.1 Méthodes scalaires explicites

Les méthodes explicites de résolution des équations différentielles scalaires du

premier ordre s’appliquent aux systèmes.
#–
dy #– #–)
= f (t, y
dt
À chaque étape, effectuer les calculs sur chaque composante
avant de passer à l’étape suivante : exemple avec point milieu

Étape 1 : vecteur des pentes au bord gauche de l’intervalle

# – #–
k1 = f (t1 , y# –1 )

k1,1 = f1 (t1 , y1,1 , y1,2 , . . . , y1,n )

k1,2 = f2 (t1 , y1,1 , y1,2 , . . . , y1,n )
... = ...
k1,n = fn (t1 , y1,1 , y1,2 , . . . , y1,n )

MNCS 50 2019-2020
EDO 5 Systèmes d’EDO du 1er ordre 5.1 Méthodes scalaires explicites

avant de calculer...

Étape 2 : vecteur des pentes au point milieu prédit

# – #– #–
k2 = f (t1 + h/2, y# –1 + k1 h/2)

k2,1 = f1 (t1 + h/2, y1,1 + k1,1 h/2, y1,2 + k1,2 h/2, . . . , y1,n + k1,n h/2)
k2,2 = f2 (t1 + h/2, y1,1 + k1,1 h/2, y1,2 + k1,2 h/2, . . . , y1,n + k1,n h/2)
... = ...
k2,n = fn (t1 + h/2, y1,1 + k1,1 h/2, y1,2 + k1,2 h/2, . . . , y1,n + k1,n h/2)

Étape 3 : vecteur résultat au bord droit de l’intervalle

#–
u = #– + h k# –
u
i+1 i 2

MNCS 51 2019-2020
EDO 5 Systèmes d’EDO du 1er ordre 5.2 Équations de Lotka-Volterra

5.2 Exemple de système non-linéaire couplé du premier ordre :

équations de Lotka-Volterra
Deux populations en conflit : modèle proies (y1 ) – prédateurs (y2 )
a1 = 1/τ1 = taux de croissance de y1 (proies) en l’absence de y2 (prédateurs)
a2 = 1/τ2 = taux de décroissance de y2 (prédateurs) en l’absence de y1 (proies)

Termes de couplage non-linéaires en y1 y2 (rencontre des 2 espèces)

a1
y2 = taux de destruction des proies par les prédateurs
k2
a2
y1 = taux de croissance des prédateurs au détriment des proies
k1
 
dy1 y2
= +a1 y1 1 − (20a)
dt k2
 
dy2 y1
= −a2 y2 1 − (20b)
dt k1
Solutions périodiques
MNCS 52 2019-2020
EDO 5 Systèmes d’EDO du 1er ordre 5.2 Équations de Lotka-Volterra

Lotka-Volterra : cycle dans le plan de phase

En éliminant le temps, on obtient un invariant , donc des solutions périodiques :

dy2 a2 y2 1 − y1 /k1
=− ⇒ y1a2 y2a1 e−a1 y2 /k2 −a2 y1 /k1 = Cte
dy1 a1 y1 1 − y2 /k2
Tangentes horizontales pour y1 = k1 (ou y2 = 0) : équilibre des prédateurs
Tangentes verticales pour y2 = k2 (ou y1 = 0) : équilibre des proies
k1 = k2 = 1, a1 = 1, a2 = 0, 2
Plan de phase de Lotka Volterra Plan de phase de Lotka Volterra
valeurs initiales : proies(0) = 10 prédateurs(0) = 1 valeurs initiales : proies(0) = 10 prédateurs(0) = 1
pas = 0,1 ordre de la méthode : 1 pas = 0,1 ordre de la méthode : 4
8.0 8.0
prédateurs= f(proies) prédateurs= f(proies)
7.0 7.0

6.0 6.0

5.0 5.0
prédateurs

prédateurs
4.0 4.0

3.0 3.0

2.0 2.0

1.0 1.0

0.0 0.0
0 5 10 15 20 0 5 10 15 20
proies proies

Méthode d’Euler : non périodique Runge Kutta d’ordre 4 : cycle correct

MNCS 53 2019-2020
EDO 5 Systèmes d’EDO du 1er ordre 5.2 Équations de Lotka-Volterra

Résolution numérique de Lotka-Volterra : k1 = k2 = 1, a1 = 1, a2 = 0, 2, h = 0, 1

Échelles linéaires

Lotka Volterra valeurs initiales : proies(0) = 10 prédateurs(0) = 1 Lotka Volterra valeurs initiales : proies(0) = 10 prédateurs(0) = 1
pas = 0,1 ordre de la méthode : 1 pas = 0,1 ordre de la méthode : 4
20 20
proies proies
prédateurs prédateurs

15 15
Populations

Populations
10 10

5 5

0 0
0 20 40 60 80 100 0 20 40 60 80 100
temps temps

Méthode d’Euler progressive : Méthode de Runge Kutta d’ordre 4 :

Les solutions divergent Cycle stable

MNCS 54 2019-2020
EDO 5 Systèmes d’EDO du 1er ordre 5.2 Équations de Lotka-Volterra

Résolution numérique de Lotka-Volterra : k1 = k2 = 1, a1 = 1 et a2 = 0, 2.

Échelle log en ordonnée

Lotka Volterra valeurs initiales : proies(0) = 10 prédateurs(0) = 1 Lotka Volterra valeurs initiales : proies(0) = 10 prédateurs(0) = 1
pas = 0,1 ordre de la méthode : 1 pas = 0,1 ordre de la méthode : 4
2 2
10 10
proies proies
prédateurs prédateurs
100 100

10−2 10−2
Populations

Populations
10−4 10−4

10−6 10−6

10−8 10−8

−10 −10
10 10
0 20 40 60 80 100 0 20 40 60 80 100
temps temps

Méthode d’Euler progressive : divergente Méthode de Runge Kutta d’ordre 4 : stable

d ln y2
Pente avec peu de proies : ≈ −1/τ2 d’où facteur e−4 ≈ 1/54 sur une durée de 20 = 4τ2 .
dt
d ln y1
Pente avec peu de prédateurs : ≈ 1/τ1 d’où facteur 100 sur durée de 4, 6 = 4, 6τ1
dt
MNCS 55 2019-2020
EDO 5 Systèmes d’EDO du 1er ordre 5.2 Équations de Lotka-Volterra

Changement de variables : passage au logarithme des populations

En normalisant les populations par leur valeurs à l’équilibre et en prenant le

logarithme, on introduit les variables z1 = ln(y1 /k1 ) et z2 = ln(y2 /k2 ).
Le système différentiel prend alors la forme séparable :
dz1 1
= (1 − exp(z2 )) (21)
dt τ1
dz2 1
= − (1 − exp(z1 )) (22)
dt τ2
où la dérivée de z1 ne dépend que de z2 et réciproquement.

Le système transformé possède un invariant Λ

Λ(z1 , z2 ) = τ1 (exp(z1 ) − z1 ) + τ2 (exp(z2 ) − z2 ) (23)

qui ne dépend donc que des conditions initiales. Il est constant sur une orbite
(trajectoire fermée) dans le plan de phase.

MNCS 56 2019-2020
EDO 5 Systèmes d’EDO du 1er ordre 5.2 Équations de Lotka-Volterra

Exemple de résolution avec Euler progressive :

Etu.jal | Lotka-Voltera: LV-euler_h0.05_50.dat | dt= 0.05 , z1(0)= 2 , z2(0)= -0.2
1.5
Diagramme de Phase
18

22
18
14
1.0

0.5
Composante Z2(t)

0.0

0.5

1.0 10

18
1.5
14

2.0
18

22
6 4 2 0 2
Composante Z1(t)

La trajectoire de la solution coupe des iso-Λ. Elle ne respecte pas l’invariant.

MNCS 57 2019-2020
EDO 6 Équations différentielles d’ordre supérieur

6 Équations différentielles d’ordre supérieur

n n−1
 
d y dy d y
= f t, y, , . . . ,
dtn dt dtn−1
Une EDO scalaire d’ordre n se ramène à un système de n équations différentielles
du premier ordre couplées en posant :

       
y1 y y10 y2
       
 y   y0   y0   y3 
 2   2 
= =⇒  =
 
    
. . .   . . .  . . .   ... 
       
yn y (n−1) yn0 f (t, y1 , y2 , . . . , yn )

MNCS 58 2019-2020
EDO 6 Équations différentielles d’ordre supérieur 6.1 Exemple

6.1 Exemple

Système linéaire du second ordre avec excitation h(t)

d2 y dy
2
= a + by + h(t)
dt dt
Poser
       
y y y10 y2
 1 =   ⇒  = 
y2 y0 y20 ay2 + by1 + h(t)

Condition initiale vectorielle : position y(t0 ) et vitesse y 0 (t0 )

Remarque Système différentiel d’ordre p de dimension n

⇒ système différentiel couplé du premier ordre à np dimensions.

MNCS 59 2019-2020
EDO 6 Équations différentielles d’ordre supérieur 6.2 Exemple d’EDO d’ordre 2 : le pendule

6.2 Exemple d’EDO d’ordre 2 : le pendule

Pendule non linéaire (y = position angulaire)

d2 y 2 2 l
= −k sin(y) où k = g/l (24)
dt2
y
Pendule linéarisé (cas des petites amplitudes) : sin(y) ≈ y
mg
d2 y 2
2
= −k y (25)
dt
l’équation linéarisée admet une solution analytique en A cos(kt) + B sin(kt).

Exprimer cette EDO non linéaire du second ordre sous la forme d’un système
différentiel couplé de dimension 2 mais du premier ordre.
       
y y y10 y2
 1 =   =⇒  = 
y2 y0 y20 −k 2 sin(y1 )

MNCS 60 2019-2020
EDO 6 Équations différentielles d’ordre supérieur 6.2 Exemple d’EDO d’ordre 2 : le pendule

Résolution système non-linéaire, avec le vecteur des valeurs initiales :

     
y1 (0) y(0) position angulaire
  =  dy = 
y2 (0) (0) = a vitesse angulaire
dt

Énergie mécanique conservée (après division par ml2 ) :

 2
1 dy
+ k 2 (1 − cos y) = constante
2 dt
Cas où y(0) = 0 (départ en position d’équilibre stable)
 2  2
1 dy 1 dy
+ k 2 (1 − cos y) = (0)
2 dt 2 dt
Vitesse angulaire minimale pour y = π (position d’équilibre instable si atteinte).
dy
Si a = (0) > 2k (seuil) ⇒ la vitesse angulaire ne s’annule pas (apériodique).
dt
Étude de la transition périodique–apériodique selon a dans le cas où k = 1

MNCS 61 2019-2020
EDO 6 Équations différentielles d’ordre supérieur 6.2 Exemple d’EDO d’ordre 2 : le pendule

Comparaisons non-linéaire (Runge-Kutta 4)–analytique linéarisé : y(t)

pendule + sol. linéarisée: méth. ordre 4 vitesse ang. init=0.2 pendule + sol. linéarisée: méth. ordre 4 vitesse ang. init=1
méthode d’ordre 4 méthode d’ordre 4
de 0.000000 à 50.000000 avec un pas de 0.01 de 0.000000 à 50.000000 avec un pas de 0.01
valeurs initiales 0.0000000e+00 2.0000000e−01 valeurs initiales 0.0000000e+00 1.0000000e+00
0.25 1.5

0.20
1.0
0.15

0.10
0.5
0.05
theta

theta
0.00 0.0

−0.05
−0.05

−0.5
−0.5
−0.10
−0.10

−0.15
−0.15
−1.0
−1.0

−0.20
−0.20

−0.25
−0.25 −1.5
−1.5
0 5 10 15 20 25 30 35 40 45 50 0 5 10 15 20 25 30 35 40 45 50
estimé estimé
linéaire temps linéaire temps

a = 0.2  1 linéarisable a = 1 périodique non sinusoïdal

pendule + sol. linéarisée: méth. ordre 4 vitesse ang. init=1.98 pendule + sol. linéarisée: méth. ordre 4 vitesse ang. init=2.02
méthode d’ordre 4 méthode d’ordre 4
de 0.000000 à 50.000000 avec un pas de 0.01 de 0.000000 à 50.000000 avec un pas de 0.01
valeurs initiales 0.0000000e+00 1.9800000e+00 valeurs initiales 0.0000000e+00 2.0200000e+00
3 50

2 40

1 30
theta

theta
0 20

−1
−1 10

−2
−2 0

−3
−3 −10
−10
0 5 10 15 20 25 30 35 40 45 50 0 5 10 15 20 25 30 35 40 45 50
estimé estimé
linéaire temps linéaire temps

MNCS
a = 1.98 périodique non sinusoïdal 62
a = 2.02 apériodique 2019-2020
EDO 6 Équations différentielles d’ordre supérieur 6.2 Exemple d’EDO d’ordre 2 : le pendule

Comparaisons non-linéaire (RK 4)–analytique linéarisé : plan de phase y 0 (y)

plan de phase pendule: méth. ordre 4 vitesse ang. init.=0.2 plan de phase pendule: méth. ordre 4 vitesse ang. init.=1
méthode d’ordre 4 méthode d’ordre 4
de 0.000000 à 50.000000 avec un pas de 0.01 de 0.000000 à 50.000000 avec un pas de 0.01
valeurs initiales 0.0000000e+00 2.0000000e−01 valeurs initiales 0.0000000e+00 1.0000000e+00
0.25 1.0

0.20 0.8

0.15 0.6

0.10 0.4

0.05 0.2
dtheta/dt

dtheta/dt
0.00 0.0

−0.05 −0.2

−0.10 −0.4

−0.15 −0.6

−0.20 −0.8

−0.25 −1.0
−0.25 −0.20 −0.15 −0.10 −0.05 0.00 0.05 0.10 0.15 0.20 0.25 −1.5 −1.0 −0.5 0.0 0.5 1.0 1.5
estimé estimé
linéarisé theta linéarisé theta

a = 0.2  1 linéarisable a = 1 périodique non sinusoïdal

plan de phase pendule: méth. ordre 4 vitesse ang. init.=1.98 plan de phase pendule: méth. ordre 4 vitesse ang. init.=2.02
méthode d’ordre 4 méthode d’ordre 4
de 0.000000 à 50.000000 avec un pas de 0.01 de 0.000000 à 50.000000 avec un pas de 0.01
valeurs initiales 0.0000000e+00 1.9800000e+00 valeurs initiales 0.0000000e+00 2.0200000e+00
2.0 2.5

2.0
1.5
1.5
1.0
1.0
0.5
0.5
dtheta/dt

dtheta/dt
0.0 0.0

−0.5
−0.5
−1.0
−1.0
−1.5
−1.5
−2.0

−2.0 −2.5
−3 −2 −1 0 1 2 3 −10 0 10 20 30 40 50
estimé estimé
linéarisé theta linéarisé theta

MNCS
a = 1.98 périodique non sinusoïdal 63
a = 2.02 apériodique 2019-2020
EDO 6 Équations différentielles d’ordre supérieur 6.2 Exemple d’EDO d’ordre 2 : le pendule

Comparaisons non-linéaire (Euler)–analytique linéarisé y(t)

pendule + sol. linéarisée: méth. ordre 1 vitesse ang. init=0.2 pendule + sol. linéarisée: méth. ordre 1 vitesse ang. init=1
méthode d’ordre 1 méthode d’ordre 1
de 0.000000 à 50.000000 avec un pas de 0.01 de 0.000000 à 50.000000 avec un pas de 0.01
valeurs initiales 0.0000000e+00 2.0000000e−01 valeurs initiales 0.0000000e+00 1.0000000e+00
0.3 1.5

0.2 1.0

0.1 0.5
theta

theta
0.0 0.0

−0.1
−0.1 −0.5
−0.5

−0.2
−0.2 −1.0
−1.0

−0.3
−0.3 −1.5
−1.5
0 5 10 15 20 25 30 35 40 45 50 0 5 10 15 20 25 30 35 40 45 50
estimé estimé
linéaire temps linéaire temps

a = 0.2 a=1
pendule + sol. linéarisée: méth. ordre 1 vitesse ang. init=1.98 pendule + sol. linéarisée: méth. ordre 1 vitesse ang. init=2.02
méthode d’ordre 1 méthode d’ordre 1
de 0.000000 à 50.000000 avec un pas de 0.01 de 0.000000 à 50.000000 avec un pas de 0.01
valeurs initiales 0.0000000e+00 1.9800000e+00 valeurs initiales 0.0000000e+00 2.0200000e+00
5 60

0
50
−5
−5

−10
−10 40

−15
−15
30
theta

theta
−20
−20

20
−25
−25

−30
−30
10
−35
−35
0
−40
−40

−45
−45 −10
−10
0 5 10 15 20 25 30 35 40 45 50 0 5 10 15 20 25 30 35 40 45 50
estimé estimé
linéaire temps linéaire temps

MNCS
a = 1.98 apériodique selon Euler ! 64
a = 2.02 2019-2020
EDO 6 Équations différentielles d’ordre supérieur 6.2 Exemple d’EDO d’ordre 2 : le pendule

Comparaisons non-linéaire (Euler)–analytique linéarisé : plan de phase y 0 (y)

plan de phase pendule: méth. ordre 1 vitesse ang. init.=0.2 plan de phase pendule: méth. ordre 1 vitesse ang. init.=1
méthode d’ordre 1 méthode d’ordre 1
de 0.000000 à 50.000000 avec un pas de 0.01 de 0.000000 à 50.000000 avec un pas de 0.01
valeurs initiales 0.0000000e+00 2.0000000e−01 valeurs initiales 0.0000000e+00 1.0000000e+00
0.3 1.5

0.2 1.0

0.1 0.5
dtheta/dt

dtheta/dt
0.0 0.0

−0.1 −0.5

−0.2 −1.0

−0.3 −1.5
−0.3 −0.2 −0.1 0.0 0.1 0.2 0.3 −1.5 −1.0 −0.5 0.0 0.5 1.0 1.5
estimé estimé
linéarisé theta linéarisé theta

a = 0.2 a=1
plan de phase pendule: méth. ordre 1 vitesse ang. init.=1.98 plan de phase pendule: méth. ordre 1 vitesse ang. init.=2.02
méthode d’ordre 1 méthode d’ordre 1
de 0.000000 à 50.000000 avec un pas de 0.01 de 0.000000 à 50.000000 avec un pas de 0.01
valeurs initiales 0.0000000e+00 1.9800000e+00 valeurs initiales 0.0000000e+00 2.0200000e+00
2.0 2.5

1.5 2.0

1.5
1.0
1.0
0.5
0.5
0.0
dtheta/dt

dtheta/dt
0.0
−0.5
−0.5
−1.0
−1.0
−1.5
−1.5

−2.0 −2.0

−2.5 −2.5
−45 −40 −35 −30 −25 −20 −15 −10 −5 0 5 −10 0 10 20 30 40 50 60
estimé estimé
linéarisé theta linéarisé theta

MNCS
a = 1.98 apériodique selon Euler ! 65
a = 2.02 2019-2020
EDO 6 Équations différentielles d’ordre supérieur 6.2 Exemple d’EDO d’ordre 2 : le pendule

Stabilité à long terme avec Euler progressive et rétrograde

Pendule linéarisé sans frottement représenté dans l’espace des phases :

comportement à long terme d’un système non dissipatif

Même méthode sur les 2 composantes (position et vitesse) h = 0.025

pendulelin-Euler-esp-phases pendulelin-EulerBack-esp-phases
[ 0.0e+00 , 5.0e+01] h= 2.50e-02 y0[0]=0.00e+00 y0[1]=1.00e+00 n= 2001 [ 0.0e+00 , 5.0e+01] h= 2.50e-02 y0[0]=0.00e+00 y0[1]=1.00e+00 n= 2001
2 2

1 1
vitesse

vitesse
0 0

1 1

2 2
3 2 1 0 1 2 3 3 2 1 0 1 2 3
position position

F IGURE 4 – Euler progressive : in- F IGURE 5 – Euler rétrograde : stable,

stable, amplification contraction

MNCS 66 2019-2020
EDO 6 Équations différentielles d’ordre supérieur 6.2 Exemple d’EDO d’ordre 2 : le pendule

Méthode d’Euler symplectique Méthodes progressives et méthodes rétrogrades

présentent chacune des inconvénients antagonistes : dans le cas d’un système à
deux composantes, on peut espérer les compenser en appliquant une méthode
progressive sur z1 et rétrograde sur z2 .

dz1
z1 (t + h) = z1 (t) + h (t, z1 (t), z2 (t)) (26)
dt
dz2
z2 (t + h) = z2 (t) + h (t + h, z1 (t + h), z2 (t + h) ) (27)
dt
Dans le cas d’un système séparable , la méthode rétrograde ne nécessite plus
dz
d’itération car dt2 ne dépend pas de z2 et z1 (t + h) a été calculé avant.

dz1
z1 (t + h) = z1 (t) + h (t, z2 (t)) (28)
dt
dz2
z2 (t + h) = z2 (t) + h (t + h, z1 (t + h)) (29)
dt

MNCS 67 2019-2020
EDO 6 Équations différentielles d’ordre supérieur 6.2 Exemple d’EDO d’ordre 2 : le pendule

Exemple de méthode symplectique : alternance des méthodes progressive et

rétrograde entre position et vitesse

pendulelin-EulerForBack-esp-phases pendulelin-EulerBackFor-esp-phases
[ 0.0e+00 , 5.0e+01] h= 2.50e-02 y0[0]=0.00e+00 y0[1]=1.00e+00 n= 2001 [ 0.0e+00 , 5.0e+01] h= 2.50e-02 y0[0]=0.00e+00 y0[1]=1.00e+00 n= 2001
1.5 1.5

1.0 1.0

0.5 0.5
vitesse

vitesse
0.0 0.0

0.5 0.5

1.0 1.0

1.5 1.5
2.0 1.5 1.0 0.5 0.0 0.5 1.0 1.5 2.0 2.0 1.5 1.0 0.5 0.0 0.5 1.0 1.5 2.0
position position

F IGURE 6 – Euler mixte F IGURE 7 – Euler mixte

— progressive sur position, — rétrograde sur position,

— rétrograde sur vitesse. — progressive sur vitesse.

MNCS 68 2019-2020
EDO 7 Implémentation vectorielle

7 Mise en œuvre vectorielle des méthodes à un pas

7.1 Introduction

— Les méthodes d’intégration doivent fonctionner quelle que soit la taille p des
vecteurs qui représentent la solution
#– et le second membre #–
y f de l’EDO.
— Il en est de même pour l’interface formelle de la fonction second membre en
fortran ou le pointeur de fonction second membre en C.
— C’est le programme principal qui fixera cette taille.
Il devra donc choisir un second membre de la même dimension.
— Les tailles des tableaux des seconds membres effectifs seront héritées du
programme principal et non déclarées explicitement.
#–
Mais seules les p composantes effectives de f (2 pour le pendule : dérivée et
#–
dérivée seconde) seront calculées à partir des p composantes de y .

MNCS 69 2019-2020
EDO 7 Implémentation vectorielle 7.2 En fortran (norme 2003)

7.2 En fortran (norme 2003)

Utiliser des fonctions à argument tableau de rang 1
#–
y
d’étendue p déterminée à l’exécution (nombre p d’EDO scalaires d’ordre 1)
et à résultat tableau de même étendue que
#– pour :
y
1. le second membre de l’équation différentielle :

MODULE abstrait
ABSTRACT INTERFACE ! de la fct générique IInd mb de l’EDO
FUNCTION fty(t, y) ! dy/dt
REAL, DIMENSION(:),INTENT(in) :: y ! variable vecteur
REAL, INTENT(in) :: t
REAL, DIMENSION(SIZE(y)) :: fty ! vecteur résultat
END FUNCTION fty
END INTERFACE
END MODULE abstrait
#–
L’étendue p du vecteur résultat f effectif sera donc fixée par le programme
principal via
#– et non par la fonction second membre.
y
MNCS 70 2019-2020
EDO 7 Implémentation vectorielle 7.2 En fortran (norme 2003)

2. chacune des méthodes à un pas (Euler, point milieu et Runge Kutta) :

#–
les pentes locales ki seront des tableaux locaux, par exemple automatiques.

FUNCTION u2_rk4(u1, t1, h, f)

USE abstrait ! où est définie l’interface abstraite fty
REAL, DIMENSION(:), INTENT(IN) :: u1 ! valeur initiale
REAL, INTENT(IN) :: t1 ! instant initial
REAL, INTENT(IN) :: h ! pas
PROCEDURE(fty) :: f ! déclaration de l’interface de f
REAL, DIMENSION(SIZE(u1)) :: u2_rk4 ! valeur estimée à t1+h
! variables locales de même étendue que u1
REAL, DIMENSION(SIZE(u1)) :: k1, k2, k3, k4 ! pentes locales
...
u2_rk4 = u1 + ...
END FUNCTION u2_rk4

MNCS 71 2019-2020
EDO 7 Implémentation vectorielle 7.2 En fortran (norme 2003)

3. Dans le programme principal (et dans la procédure d’écriture sur fichier), les
solutions vectorielles (analytique et par intégration) sont représentées par des
tableaux 2D alloués dynamiquement (n instants, p composantes).
Les étendues n et p sont donc choisies à l’exécution, sachant que p doit être
correspondre au nombre effectif de composantes du second membre étudié.
Mais la dimension temporelle n’est pas « vue » par les méthodes : elles
travaillent sur des vecteurs (d’étendue p) dans un intervalle [ti , ti+1 ], à i fixé.

MNCS 72 2019-2020
EDO 7 Implémentation vectorielle 7.3 En C89 avec des tableaux dynamiques

7.3 En C89 avec des tableaux dynamiques sur le tas

Utiliser des fonctions à « argument tableau 1D »

#– de taille déterminée à l’exécution
y
#–
et rendant un pointeur vers un tableau alloué sur le tas de même taille que y .

1. le second membre de l’équation différentielle sera alloué par la fonction

#– #–
f ( y , t) qui rend le pointeur vers ce tableau, dont la fonction appelante (la
méthode d’intégration) devra prendre en charge la libération ;

float * pendule(float t, float *u, int p){

float *second_membre= NULL; /* tableau 1D */
second_membre = float1d(p); /* allocation sur le tas */
second_membre[0] = u[1];
second_membre[1] = -sin(u[0]);
return second_membre; /* valeur de retour = pointeur */
}

MNCS 73 2019-2020
EDO 7 Implémentation vectorielle 7.3 En C89 avec des tableaux dynamiques

2. pour chacune des méthodes à un pas (Euler, point milieu et Runge Kutta) les
#–
pentes locales ki seront des tableaux alloués et libérés localement, car leur
nombre dépend de la méthode ; en revanche, le résultat
#–
u i+1 qui est aussi
vectoriel sera passé en argument (sous forme pointeur plus nombre
d’éléments), son allocation et libération prises en charge par l’appelant.

void u2_milieu(int p, float *u1, float t1, float h,

float* (* ptr_f) (float, float*, int),
float * u2){
/* permettant d’avancer d’un pas en temps*/
float *k1 = NULL; /* vecteur pente */
float *k2 = NULL; /* vecteur pente */
float *u12 = NULL; /* vecteur intermédiaire */
k1 = (*ptr_f)(t1, u1, p); /* allocation par la fct IInd membre */
u12 = float1d (p); /* allocation locale */
/* ... */
float1d_libere(k1); /* libération des vecteurs des pentes */
/* ... */

MNCS 74 2019-2020
EDO 7 Implémentation vectorielle 7.3 En C89 avec des tableaux dynamiques

3. Dans le programme principal (et dans la procédure d’écriture sur fichier), les
solutions vectorielles (analytique et par intégration) sont représentées par des
tableaux 2D. Mais la dimension temporelle n’est pas « vue » par les méthodes :
elles travaillent sur des vecteurs de taille p dans un intervalle [ti , ti+1 ]. Cela
impose que les composantes des vecteurs soient contigües en mémoire, donc
le deuxième indice est celui des composantes, le premier celui du temps.

/* allocation dans le main des tableaux 2D */

u = float2d(n, p); /* n instants et p equations */
/* appel de la méthode du point milieu par exemple */
u2_milieu(p, u[i], t[i], h, &pendule, u[i+1]);
/* donc u[i] est un tableau 1D = vecteur des composantes de u_i *

MNCS 75 2019-2020
EDO 7 Implémentation vectorielle 7.4 En C99 avec des tableaux automatiques

7.4 En C99 avec des tableaux automatiques

#–
Fonctions à « argument tableau 1D » y de taille p déterminée à l’exécution
Déclaration tardive des tableaux automatiques ⇒ éviter les tableaux dynamiques
Mais une fonction ne peut pas rendre un tableau ⇒ fonctions à résultat void
⇒ déclaration du tableau argument par l’appelant
et remplissage par la fonction appelée
1. la fonction second membre de l’équation différentielle remplit le tableau
#– #–
second_mb de taille p représentant f ( y , t)
qui a été déclaré par l’appelant (la méthode) avec la taille fixée par le main

// version C99 avec tableaux automatiques

void pendule(float t, int p, float u[p], float second_mb[p]){
// p = 2 ici = dimension des vecteurs u et second_membre
second_mb[0] = u[1];
second_mb[1] = -sin(u[0]);
return;
}
MNCS 76 2019-2020
EDO 7 Implémentation vectorielle 7.4 En C99 avec des tableaux automatiques

2. pour chacune des méthodes à un pas (Euler, point milieu et Runge Kutta) les
#–
pentes locales ki seront des tableaux locaux automatiques, car leur nombre
#–
dépend de la méthode ; en revanche, le résultat u i+1 qui est aussi vectoriel
sera passé en argument, sa déclaration étant prise en charge par l’appelant.

// C99 avec tableaux automatiques

void u2_milieu(int p, float u1[p], float t1, float h,
void (* ptr_f) (float, int, float[], float[]),
float u2[p]) {
// methode permettant d’avancer d’un pas en temps
float k1[p] ; // vecteur pente à gauche
float k2[p] ; // vecteur pente au milieu
float u12[p]; // vecteur intermédiaire
(*ptr_f)(t1, p, u1, k1); // évaluation du second membre en u1
// résultat dans le vecteur local k1
...
// tableau u2[p] déclaré dans l’appelant et rempli ici (boucle)

MNCS 77 2019-2020
EDO 7 Implémentation vectorielle 7.4 En C99 avec des tableaux automatiques

3. Dans le programme principal (et dans la procédure d’écriture sur fichier), les
solutions vectorielles (analytique et par intégration) sont représentées par des
tableaux 2D. Mais la dimension temporelle n’est pas « vue » par les méthodes :
elles travaillent sur des vecteurs de taille p dans un intervalle [ti , ti+1 ]. Cela
impose que les composantes des vecteurs soient contigües en mémoire, donc
le deuxième indice est celui des composantes, le premier celui du temps.

// tableaux automatiques 2D C99 déclarés dans le main

// p choisi selon la dimension du second membre
float u[n][p]; // n instants et p equations
// dans la boucle sur les instants i :
// appel de la méthode du point milieu par ex.
u2_milieu(p, u[i], t[i], h, &pendule, u[i+1]);
// u[i] et u[i+1] : vecteurs à p composantes
// t[i] : scalaire

MNCS 78 2019-2020
EDO RÉFÉRENCES RÉFÉRENCES

Références

A KAI , T ERRENCE J., Applied Numerical Methods for Engineers, 410 pages (Wiley,
1994), ISBN 0-471-57523-2.

B URDEN , R ICHARD L. et J. D OUGLAS FAIRES, Numerical Analysis, 847 pages

(Thompson, Brooks/Cole, 2005), huitième édition, ISBN 0-534-40499-5.

D EMAILLY, J.-P., Analyse numérique et équations différentielles, 350 pages (EDP

Sciences, 2006), troisième édition, ISBN 978-2-86883-891-9.

G UILPIN , C H ., Manuel de calcul numérique appliqué, 577 pages (EDP Sciences,

1999), ISBN 2-86883-406-X.

R APPAZ , J ACQUES et M ARCO P ICASSO, Introduction à l’analyse numérique, 268

pages (Presses polytechniques et universitaires romandes, 2010), ISBN
978-2-88074-851-7.

MNCS 79 2019-2020