Equations Différentielles II

STEP, MINES ParisTech∗

5 mars 2021 (#72befad)

Table des matières

Introduction 2

Objectifs du cours 3

Limites du schéma d’Euler 4

Systèmes raides . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
Systèmes hamiltoniens . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

Méthodes à un pas 6
Principe . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
Exemples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
Définition implicite de Φ . . . . . . . . . . . . . . . . . . . . . . . . . . 8

Analyse d’erreur 9
Erreur de troncature locale . . . . . . . . . . . . . . . . . . . . . . . . 9
Consistance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
Condition nécessaire et suffisante de consistance . . . . . . . . . . 10
Condition nécessaire et suffisante de consistance d’ordre ≥ p . . . 12
Stabilité . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
Stabilité . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
Condition suffisante de stabilité . . . . . . . . . . . . . . . . . . . 14
Convergence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
Convergence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
Théorème de Lax . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
Condition suffisante de convergence . . . . . . . . . . . . . . . . . 15
Erreurs d’arrondi et pas optimal . . . . . . . . . . . . . . . . . . . . . 16
∗ Ce document est un des produits du projet m boisgera/CDIS, initié par la collaboration de

(S)ébastien Boisgérault (CAOR), (T)homas Romary et (E)milie Chautru (GEOSCIENCES),

(P)auline Bernard (CAS), avec la contribution de Gabriel Stoltz (Ecole des Ponts ParisTech,
CERMICS). Il est mis à disposition selon les termes de la licence Creative Commons “attribution
– pas d’utilisation commerciale – partage dans les mêmes conditions” 4.0 internationale.

1
Annexe – Choix du pas de temps 16
Pas fixe . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
Adaptation du pas de temps . . . . . . . . . . . . . . . . . . . . . . . . 17

Projet numérique 19

Exercices 21
Consistance et ordre de schémas . . . . . . . . . . . . . . . . . . . . . 21
Convergence de schémas . . . . . . . . . . . . . . . . . . . . . . . . . . 21
Explicite ou implicite ? . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
Euler symplectique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

Corrections 23
Consistance et ordre de schémas . . . . . . . . . . . . . . . . . . . . . 23
Convergence de schémas . . . . . . . . . . . . . . . . . . . . . . . . . . 26
Explicite ou implicite ? . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
Euler symplectique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

Références 29

Introduction
Ce chapitre est consacré à la résolution numérique d’équations différentielles

ẋ = f (t, x) , x(t0 ) = x0 .

La nécessité de développer des méthodes d’intégration numériques vient du

constat que seule une infime partie des équations différentielles sont résolubles
exactement. Or, on a parfois besoin de connaître le plus précisément possible le
comportement futur d’un système dynamique :
— soit en temps fini, par exemple pour déterminer la trajectoire d’une fusée
pour la mise en orbite d’un satellite;
— soit en temps long, par exemple pour déterminer un cycle limite asympto-
tique (dynamique de population) ou bien se prononcer sur la stabilité de
notre système solaire.
La méthode la plus connue est la méthode d’Euler datant de 1768, qui consiste
à implémenter
xj+1 = xj + ∆t f (tj , xj ) x0 = x0
pour un pas de temps ∆t suffisamment petit. Cette méthode appartient à
la famille des méthodes explicites, c’est-à-dire que xj+1 est directement et
explicitement défini en fonction de xj . En 1824, Cauchy montre la convergence
de cette méthode lorsque le pas de temps ∆t tend vers 0, et prouve ainsi
l’existence et l’unicité des solutions (en fait, il utilise plutôt la version implicite
de la méthode d’Euler).

2
Même si la méthode d’Euler suffit dans les cas simples, elle exige parfois de
recourir à des pas très faibles pour obtenir une précision acceptable sur des
temps longs (voir Systèmes raides (p. 4)). Parfois, le compromis entre précision à
chaque itération et accumulation des erreurs d’arrondis devient même impossible.
De plus, cette méthode n’est pas adaptée à la simulation de certains systèmes
dont certaines propriétés cruciales (comme la conservation de l’énergie) ne sont
pas préservées (voir Systèmes Hamiltoniens (p. 5)). Au cours des derniers siècles,
les scientifiques ont donc progressivement développé des méthodes de plus en
plus complexes et performantes : schémas multi-pas d’ordre supérieur, méthodes
implicites, variation du pas, schémas symplectiques etc.
En fait, dans l’histoire des équations différentielles, c’est souvent la mécanique
céleste qui a été motrice des plus grandes avancées. Au milieu du XIXe siècle,
les astronomes Adams et Le Verrier prédisent mathématiquement l’existence
et la position de la planète Neptune et l’on entend parler pour la première
fois de méthodes multi-pas. Ensuite, les progrès se sont enchaînés au rythme
des modèles physiques. La première tendance a été de rechercher des schémas
permettant toujours plus de précision à pas plus grand. Parmi les dates clés, on
peut citer la publication en 1895 de la première méthode de Runge-Kutta par
Runge, puis en 1901, de la populaire méthode de Runge-Kutta d’ordre 4 par
Kutta, et ensuite en 1910, de l’extrapolation de Richardson permettant la montée
en ordre et donc le recours à des pas plus grand pour une même précision. Mais
au milieu du XXe siècle, on découvre des systèmes, dits raides (Hirschfelder,
1952), pour lesquels cette montée en ordre ne suffit pas et pour lesquels il faut
repenser de nouveaux schémas (Dalquist, 1968). Enfin, à partir des années 80,
les scientifiques développent l’intégration numérique géométrique, c’est-à-dire
qui préservent les propriétés structurelles du système (symétrie, conservation
d’énergie etc.), utile en particulier pour la simulation des systèmes hamiltoniens.

Objectifs du cours
Ce cours a pour but de sensibiliser aux problèmes apparaissant lors de la
simulation numérique des solutions d’équations différentielles, et de donner les
bases d’analyse d’erreur numérique. Pour un exposé plus approfondi, on pourra
par exemple se référer à (Demailly 2006).
En première lecture :
— comprendre les limites d’un schéma d’Euler.
— comprendre qu’un schéma numérique à un pas consiste à discrétiser une
intégrale, en connaître quelques-uns autres que le schéma d’Euler.
— comprendre les notions de consistance/convergence d’un schéma et leur
ordre. Savoir montrer que le schéma d’Euler explicite est convergent
d’ordre 1.
En deuxième lecture :

3
 — comprendre comment fonctionne un schéma implicite et comment l’implé-
menter.
— comprendre que la convergence est la combinaison de deux concepts : la
consistance et la stabilité ; savoir utiliser ces notions pour évaluer l’impact
des erreurs d’arrondi.
— savoir calculer l’ordre de consistance et montrer la convergence de schémas
de base.
— comprendre l’apport de schémas symplectiques pour les systèmes hamil-
toniens.

Limites du schéma d’Euler

La première limite du schéma d’Euler est qu’il est d’ordre 1, c’est-à-dire qu’il
produit une erreur en ∆t2 à chaque pas. Nous verrons dans la suite des algorithmes
d’ordre supérieur qui permettent d’utiliser un pas plus grand pour une précision
donnée. Mais au delà de cette problématique, il existe des systèmes pour lesquels
de telles méthodes (même d’ordre supérieur) échouent. En voici deux exemples
célèbres.

Systèmes raides

La dénomination systèmes raides a été introduite en 1952 par Hirschfelder pour

désigner des systèmes comprenant des dynamiques aux constantes de temps
très différentes. Dans ce cas, le pas nécessaire pour simuler avec précision les
dynamiques très rapides est si petit, qu’il est alors impossible de simuler assez
longtemps pour observer les parties lentes. La particularité de ces systèmes est
que cette décroissance du pas apparaît alors que la solution est parfaitement
régulière, et non pas proche de singularités. C’est le cas des systèmes linéaires
ẋ = Ax + b
avec A de Hurwitz quand le rapport entre les parties réelles maximales et
minimales des valeurs propres devient très grand. Ce phénomène peut notamment
apparaître dans un simple système masse/ressort
mÿ = −ρẏ − ky
 0 1 
qui se met sous la forme précédente avec x = (y, ẏ) ∈ R2 et A = − k − ρ .
√ m m
Lorsque les valeurs propres sont réelles (i.e. ρ > 2 mk), leur rapport est donné
par r
mk
1+ 1−4 2
ρ
r
mk
1− 1−4 2
ρ

4
qui explose lorsque mk
ρ2 tend vers 0. Par exemple, lorsque les frottements sont
très grands par rapport à la raideur du ressort, ou bien lorsque ρ et k sont du
même ordre de grandeur et très grands.
Plus généralement, la coexistence de dynamiques très lentes à très rapides
apparaît en cinétique chimique ou en biologie. La réaction de Robertson (1966)
0.04
A −→ B (lente)
3×107
B + B −→ B + C (très rapide)
4
10
B + C −→ A + C (rapide)

modélisée par

ẋa = −0.04xb + 104 xb xc

ẋb = 0.04xa − 104 xb xc − 3 × 107 x2b
ẋc = 3 × 107 x2b

en est un exemple classique, souvent utilisée pour tester les schémas numériques.
Il s’avère que pour ces systèmes, des schémas dits implicites performent beaucoup
mieux car ils autorisent l’utilisation de pas plus grands pour une même précision
et plus de stabilité (voir l’exercice Explicite ou Implicite? (p. 21)). Pour plus
de détails voir (Hairer and Wanner 1996). L’impossibilité de trouver un pas
approprié avec un schéma d’Euler explicite pour ce système est aussi illustré
dans le notebook Equations Differentielles [Link].

Systèmes hamiltoniens

La mécanique hamiltonienne permet typiquement de modéliser le comportement

de systèmes dont une certaine énergie est conservée au cours du temps. Il peut
s’agir par exemple de planètes en interaction gravitationnelle, de particules en
interaction électromagnétique, etc.
Par exemple, dans un problème à N corps en interaction gravitationnelle, l’ha-
miltonien s’écrit 1
N
X 1 > X mi mj
H(q, p) = pi pi − G
i=1
2m i kqi − qk k
1≤i<k≤N

1. Pour obtenir l’hamiltonien, on commence par définir le lagrangien L(t, q, q̇), puis la
quantité de mouvement p = ∇q̇ L(t, q, q̇), et enfin l’hamiltonien H(t, q, p) est obtenu par
transformée de Legendre. Notons que dans ce cas général où H peut dépendre explicitement
du temps (par exemple si de l’énergie est injectée ou prélevée par une action extérieure au
d
système), on a dt H(t, q(t), p(t)) = ∇t H(t, q(t), p(t)), donc l’hamiltonien varie selon cet effet
extérieur, et n’est plus constant.

5
où qi ∈ R3 désigne la position de chaque corps, mi sa masse, et pi = mi q̇i ∈ R3
sa quantité de mouvement. Le comportement de chaque corps est alors régi par
la dynamique hamiltonienne 2
1
q̇i = ∇pi H(q, p) = pi
mi
X mi mj
ṗi = −∇qi H(q, p) = −G (qi − qk )
kqi − qk k3
k6=i

On a alors le long des trajectoires

d
H(q(t), p(t)) = h∇q H(t, q(t), p(t)), q̇i + h∇p H(t, q(t), p(t)), ṗi = 0
dt
et l’énergie H(q, p) est donc conservée.
Or, lorsqu’on essaye de simuler le système solaire avec un schéma d’Euler
(explicite), l’énergie augmente peu à peu à chaque révolution et les trajectoires
sont des spirales divergentes. Avec un schéma d’Euler implicite, Jupiter et Saturne
s’effondre vers le soleil et sont éjectées du système solaire ! Même des schémas
d’ordre supérieur ne permettent pas de simuler correctement ce système sur
des temps “courts” sur l’échelle de temps astronomique (à moins de prendre
des pas déraisonnablement petits). En fait, le problème c’est que ces méthodes
d’intégration ne préservent pas les propriétés structurelles des solutions telles
que la conservation de l’énergie. Il faut donc développer des schémas particuliers,
appelés symplectiques, comme illustré sur un simple oscillateur dans l’exercice
Schéma symplectique (p. 22). Pour aller plus loin sur ces méthodes, voir (Hairer,
Lubich, and Wanner 2010). Un exemple simple de système à deux corps est aussi
donné dans le notebook Equations Differentielles [Link].

Méthodes à un pas

Principe

Pour approximer les solutions d’une équation différentielle sur un intervalle [0, T ],
les méthodes numériques à un pas se basent sur la représentation intégrale
Z t J−1
X Z tj+1
x(t) = x0 + f (s, x(s))ds = x0 + f (s, x(s))ds
t0 j=0 tj

2. L’application des lois de Newton donnerait directement

X X mi mj (qi − qk )
mi ai = mi q¨i = Fk = −G
kqi − qk k2 kqi − qk k
k6=i k6=i

où Fk sont les forces de gravitation exercées par chaque corps k sur le corps i.

6
où t0 < t1 < . . . < tJ avec tJ = T . L’idée est d’approximer les intégrales
R tj+1
tj
f (s, x(s)) sur des intervalles [tj , tj+1 ] suffisamment petits.

Dans la suite, on note xj l’approximation au temps tj de la valeur exacte x(tj )

et ∆tj = tj+1 − tj le jème pas de temps. L’idée est de calculer récursivement
xj+1 = xj + ∆tj Φ(tj , xj , ∆tj )
où Φ(tj , xj , ∆tj ) doit donc approximer
Z tj+1
1
f (s, x(s))ds .
tj+1 − tj tj
Les différentes méthodes de quadrature, i.e. d’approximation de l’intégrale,
peuvent donc être mises à profit. La difficulté ici est que seule la valeur initiale
f (tj , x(tj )) de f est connue (ou du moins estimée) à l’itération j, par f (tj , xj ).
On distingue donc les méthodes explicites où Φ(tj , xj , ∆tj ) est écrite directement
explicitement en fonction de la valeur initiale xj , et les méthodes implicites où
cette expression n’est connue qu’implicitement et des étapes intermédiaires de
calcul sont nécessaires.

Exemples

1. Méthodes explicites:
— Euler explicite: l’intégrale est approximée par l’aire d’un rectangle
déterminé par la valeur initiale de f à gauche de l’intervalle, i.e.

xj+1 = xj + ∆tj f (tj , xj ) .

— méthode de Heun : l’intégrale est approximée par l’aire d’un trapèze
déterminé par la valeur initiale de f et une approximation de sa valeur
finale, i.e.,
∆tj 

xj+1 = xj + f (tj , xj ) + f tj+1 , xj + ∆tj f (tj , xj ) .
2
— schéma de Runge–Kutta d’ordre 4:
j

 F1 = f (t

 
j, x )

∆tj j ∆tj


F f t , x F


 2
 = j + + 1
2 2
 
 ∆tj j ∆tj

 F 3 = f tj + , x + F 2
2 2




 j
F4 = f (tj + ∆tj , x + ∆tj F3 ),


et on pose
F1 + 2F2 + 2F3 + F4
xj+1 = xj + ∆t .
6

7
 2. Méthodes implicites:
— Euler implicite : l’intégrale est approximée par l’aire d’un rectangle
déterminé par la valeur finale de f à droite de l’intervalle, i.e.

xj+1 = xj + ∆tj f (tj+1 , xj+1 ) .

— méthode des trapèzes (ou Crank–Nicolson) : l’intégrale est approximée
par l’aire du trapèze déterminé par les valeurs initiales et finales de f ,
i.e.
∆tj  
xj+1 = xj + f (tj , xj ) + f (tj+1 , xj+1 ) .
2
— méthode du point milieu : l’intégrale est approximée par l’aire d’un
rectangle déterminé par une approximation de la valeur de f au milieu
de l’intervale, i.e.

tj + tj+1 xj + xj+1
 
xj+1 = xj + ∆tj f , .
2 2
On peut bien sûr construire des méthodes plus compliquées et plus précises pour
des méthodes de Runge–Kutta d’ordre supérieur (explicites ou implicites).

Définition implicite de Φ

Dans les schémas implicites, l’application Φ est définie de manière implicite. Par
exemple, pour le schéma d’Euler, on a :
 
Φ(tj , xj , ∆tj ) = f tj + ∆tj , xj + ∆tj Φ(tj , xj , ∆tj ) .

Il faut donc s’assurer que Φ est bien définie, c’est-à-dire qu’il existe bien xj+1
tel que
xj+1 = xj + ∆tj f (tj+1 , xj+1 ) .
Pour cela, nous pouvons voir xj+1 comme le point fixe de l’application Fj définie
par
Fj (x) = xj + ∆tj f (tj+1 , x) .
à xj , ∆tj , tj+1 fixés. L’existence (et l’unicité) de ce point fixe peut alors être
démontrée par le théorème de point fixe de Banach. Si x 7→ f (tj+1 , x) est
Lipschitzienne, c’est-à-dire s’il existe Lj tel que

kf (tj+1 , xa ) − f (tj+1 , xb )k ≤ Lj kxa − xb k ∀(xa , xb ) ∈ Rn × Rn ,

alors Fj : Rn → Rn est contractante pour un pas de temps ∆tj suffisamment

petit puisque
kFj (xa ) − Fj (xb )k ≤ ∆tj Lj kxa − xb k .
Puisque Rn est complet, on déduit par le théorème du point fixe que xj+1 existe
bien.

8
En pratique, on peut utiliser la méthode itérative de construction de ce point
fixe donnée par la preuve du théorème pour approcher xj+1 . Une stratégie est
de partir de la valeur donnée par le schéma d’Euler explicite

xj,0 = xj + ∆tj f (tj , xj )

et affiner ensuite par l’algorithme du point fixe en itérant

xj,k+1 = F (xj,k )

jusqu’à ce que l’évolution relative

xj,k+1 − xj,k
xj,0
devienne inférieure à un seuil choisi par l’utilisateur. Puisque la suite (xj,k )k∈N
est de Cauchy, on sait que cette algorithme s’arrête en un nombre fini d’itérations.
Un tel schéma est plus lourd en terme de calculs qu’un algorithme explicite mais
il apporte en général plus de stabilité et permet souvent d’utiliser un pas plus
grand. C’est en particulier utile pour les systèmes raides, comme illustré dans
l’exercice Explicite ou Implicite? (p. 21).

Analyse d’erreur
L’objectif de l’analyse d’erreur a priori est de donner une estimation de l’erreur
commise par la méthode numérique en fonction des paramètres du problème
(temps d’intégration, pas de temps, propriétés de f ). L’idée générale est de
remarquer qu’à chaque pas de temps, on commet une erreur d’intégration locale
(erreur de troncature dans la discrétisation de l’intégrale, à laquelle s’ajoutent
souvent des erreurs d’arrondi), et que ces erreurs locales s’accumulent au fil
des pas. Le contrôle de cette accumulation demande l’introduction d’une notion
de stabilité adéquate, alors que les erreurs locales sont liées à une notion de
consistance. L’alliance de stabilité et de consistante donne une propriété de
convergence qui est souhaitée lors de l’implémentation de méthodes numériques.

Erreur de troncature locale

L’erreur de troncature locale à l’itération j est l’erreur que l’on commet en une
seule itération lors de l’approximation de l’intégrale pour passer de xj à xj+1 .
C’est donc l’erreur théorique que l’on obtiendrait si l’on appliquait le schéma
numérique à la solution exacte x(tj ). Elle est ainsi définie comme

x(tj+1 ) − x(tj ) − ∆tj Φ(tj , x(tj ), ∆tj )

η j+1 := .
∆tj

9
Définition – Consistance

On note ∆t = max0≤j≤J−1 ∆tj le pas de temps maximal. On dit qu’une méthode

numérique est consistante si
 
lim max kη j k = 0 ,
∆t→0 1≤j≤J

et qu’elle est consistante d’ordre ≥ p s’il existe une constante cs telle que, pour
tout 0 ≤ j ≤ J − 1,
kη j+1 k ≤ cs (∆tj )p .
Une méthode est donc consistante d’ordre p si elle est consistante d’ordre ≥ p,
mais pas ≥ p + 1.

Théorème – Condition nécessaire et suffisante de consistance

Si Φ est continue, alors le schéma est consistant si et seulement si

Φ(t, x, 0) = f (t, x) ∀(t, x) ∈ [0, T ] × Rn .

Démonstration Soit C un ensemble fermé et borné tel que x(t) ∈ C pour

tout t ∈ [0, T ]. On note toujours ∆t = max0≤j≤J−1 ∆tj . Par la représentation
intégrale des solutions,
Z tj+1
x(tj+1 ) = x(tj ) + f (x(s), s)ds
tj

l’erreur de troncature locale s’écrit

Z tj+1 
j+1 1 
η = f (s, x(s)) − Φ(tj , x(tj ), ∆tj ) ds .
∆tj tj

Si le schéma est consistant alors cette erreur tend vers 0 lorsque ∆t tend vers 0
(pour n’importe quel système et n’importe quelle trajectoire). Or,
Z tj +∆t
1  
lim f (s, x(s))−Φ(tj , x(tj ), ∆tj ) ds = f (tj , x(tj ))−Φ(tj , x(tj ), 0)
∆t→0 ∆tj tj

qui doit donc être nul.

Réciproquement, supposons Φ(·, ·, 0) = f . On doit montrer que l’erreur de
consistance tend vers 0 lorsque ∆t tend vers 0 (uniformément en j = 1, . . . , N ).
Soit ε > 0. Par la continuité de Φ et f et puisque Φ(·, ·, 0) = f , il existe ∆1 > 0
tel que si ∆t ≤ ∆1 , alors

kΦ(t, x, ∆tj ) − f (t, x)k ≤ ε ∀(t, x) ∈ [0, T ] × C ∀j = 1, . . . , N .

10
Donc Z tj+1
1
kη j+1 k ≤ ε + f (s, x(s)) − f (tj , x(tj )) ds .
∆tj tj

Puisque s 7→ f (s, x(s)) est continue sur l’intervalle fermé et borné [0, T ], elle y
est uniformément continue, donc il existe ∆2 > 0 tel que si ∆t ≤ ∆2 ,

f (s, x(s)) − f (tj , x(tj )) ≤ ε ∀s ∈ [tj , tj+1 ] ∀j = 1, . . . , N

et donc kη j+1 k ≤ 2ε pour tout j. Le schéma est donc bien consistant. 

Exemple – Schémas consistants Reprenons les exemples donnés plus haut.

— Euler explicite : Φ(t, x, ∆t) = f (t, x) indépendamment de ∆t donc la
condition est trivialement satisfaite.
— Méthode de Heun : Φ(t, x, ∆t) = f (t,x)+f (t,x+∆tf
2
(t,x))
donne bien f (t, x)
si ∆t = 0.
— Runge Kutta d’ordre 4 : lorsque ∆t = 0, F1 = F2 = F3 = F4 = f (t, x)
donc Φ(t, x, 0) = F1 +2F2 +2F
6
3 +F4
= f (t, x).
De même, la consistance des méthodes implicites s’obtiennent en remarquant
que xj+1 = xj lorsque ∆t = 0.
Cette condition suffisante permet donc de prouver facilement le caractère consis-
tant d’un schéma. Cependant, en pratique, on s’intéresse surtout à son ordre de
consistance. Pour cela, l’erreur de consistance se calcule souvent par des dévelop-
pements de Taylor des solutions lorsque celles-ci sont suffisamment régulières, et
la constante cs s’exprime alors comme une borne sur les dérivées des solutions.
En fait, on remarque que lorsque f est continue, la solution est C 1 (par définition
de nos solutions). Mais puisque ẋ(t) = f (t, x(t)), ẋ hérite de la régularité de f :
si f est C k alors les solutions x sont C k+1 . Le calcul de l’ordre de consistance
dans le cas dus schéma d’Euler explicite est donné ci-dessous. Pour les autres
schémas, voir l’exercice Consistance de schémas (p. 21)

Exemple – Ordre de consistance du schéma d’Euler explicite L’erreur

de troncature s’écrit
 
x(tj + ∆tj ) − x(tj ) + ∆tj f (tj , x(tj ))
η j+1 = .
∆tj

Or, si f est C 1 , alors x est C 2 et par application la formule de Taylor avec reste
intégral, on a
Z 1
x(tj + ∆tj ) = x(tj ) + ∆tf (tj , x(tj )) + ∆t2j ẍ(tj + s∆tj )(1 − s)ds,
0

11
en utilisant ẋ(tj ) = f (tj , x(tj )). Ceci donne donc
Z 1
∆tj ∆tj
kη j+1 k ≤ ∆tj ẍ(tj +s∆tj )(1−s)ds ≤ max kẍ(t)k ≤ max kẍ(t)k .
0 2 t∈[tj ,tj+1 ] 2 t∈[0,T ]
Le schéma d’Euler explicite est donc consistant d’ordre ≥ 1 avec
maxt∈[0,T ] kẍ(t)k
cs = .
2

Notons qu’en utilisant ẋ(t) = f (t, x(t)),

ẍ(t) = ∂t f (t, x(t)) + ∂x f (t, x(t)) · f (t, x(t)),
et on peut exprimer cs en fonction de bornes sur x et sur les dérivées de
f . Plus généralement, en dérivant successivement lorsque f est C k , notons
f [k] : R × Rn → Rn la fonction dépendant des dérivées successives de f telle que
x(k+1) (t) = f [k] (t, x(t)) .
On a alors le théorème suivant, généralisant les calculs précédents.

Théorème – Condition nécessaire et suffisante de consistance d’ordre

≥p

Si Φ et f sont C p alors le schéma est consistant d’ordre ≥ p si et seulement si

∂kΦ 1
k
(t, x, 0) = f [k] (t, x) ∀0 ≤ k ≤ p − 1 , ∀(t, x) ∈ [0, T ] × Rn .
∂∆t k+1

Démonstration L’erreur de troncature s’écrit

 
x(tj + ∆tj ) − x(tj ) + ∆tj Φ(tj , x(tj ), ∆tj )
η j+1 = .
∆tj
Or, si f est C p , alors x est C p+1 et par application la formule de Taylor avec
reste intégral, on a
p
∆tkj (k) ∆tp+1
Z 1
j
X
x(tj + ∆tj ) = x(tj ) + x (tj ) + x(p+1) (tj + s∆tj )(1 − s)p ds
k! p! 0
k=1
p−1
∆tk+1 ∆tp+1
Z 1
j j
X
= x(tj ) + f [k] (tj , x(tj )) + x(p+1) (tj + s∆tj )(1 − s)p ds
(k + 1)! p! 0
k=0

Par ailleurs, puisque Φ est C p ,

p−1
∆tk ∂ k Φ ∆tp
Z 1 p
X ∂ Φ
Φ(t, x, ∆t) = (t, x, 0)+ (t, x, tj +s∆tj )(1−s)p−1 ds
k! ∂∆tk (p − 1)! 0 ∂∆tp
k=0

12
Il s’ensuit que
p−1
∆tkj ∂kΦ
 
X 1
η j+1 = f [k] (tj , x(tj )) − (t j , x(tj ), 0) + ∆tpj Rj
k=0
k! k+1 ∂∆tk

∂pΦ
où Rj est borné (uniformément en ∆tj ) par continuité de ∂∆tp et x(p+1) .
η j+1
Maintenant, si le schéma est d’ordre ≥ p, alors ∆tp
doit rester borné lorsque
j
∆tj tend vers 0 et on obtient bien la condition du théorème. Réciproquement, si
la condition du théorème est vérifiée, on obtient directement que η j+1 est borné
en ∆tpj . 

Exemple – Ordre de consistance de schémas On a vu que si f est C 1 ,

le schéma d’Euler explicite est consistant d’ordre ≥ 1. On peut le retrouver ici
en appliquant le critère précédent pour p = 1 puisque Φ(·, ·, 0) = f . Ensuite, si f
est C 2 , on constate que
∂Φ
(t, x, 0) = 0 6= f [1] (t, x) = ∂t f (t, x) + ∂x f (t, x) · f (t, x)
∂∆t
donc le schéma d’Euler explicite est consistant d’ordre exactement 1.
Par contre, toujours si f est C 2 , on constate que l’on a bien pour le schéma de
Heun
∂Φ
(t, x, 0) = f [1] (t, x) = ∂t f (t, x) + ∂x f (t, x) · f (t, x)
∂∆t
donc ce schéma est d’ordre ≥ 2. On peut vérifier que si f est C 3 , le critère n’est
pas vérifié à l’ordre supérieur donc il est consistant d’ordre égal à 2.

Stabilité

Une fois que l’on a étudié l’erreur locale commise en une itération, on s’intéresse
à la manière dont elle va se propager au fur et à mesure des itérations. Pour cela,
la notion de stabilité quantifie la robustesse de l’approximation numérique par
rapport à l’accumulation des erreurs locales et perturbations.

Définition – Stabilité

On dit qu’une méthode numérique est stable s’il existe une constante S(T ) > 0
(indépendente des ∆tj ) telle que, pour toutes suites x = {xj }1≤j≤J et z =
{z j }1≤j≤J vérifiant
( j+1
x = xj + ∆tj Φ(tj , xj , ∆tj ),
z j+1 = z j + ∆tj Φ(tj , z j , ∆tj ) + δ j+1 ,

13
on ait
 J
X 
max kxj − z j k ≤ S(T ) kx0 − z 0 k + kδ j k .
1≤j≤J
j=1

Théorème – Condition suffisante de stabilité

Si Φ sont Lipschitziennes en x, c’est-à-dire il existe L > 0 tel que pour tout

0 ≤ j ≤ J,

kΦ(tj , xa , ∆tj ) − Φ(tj , xb , ∆tj )k ≤ Lkxa − xb k ∀(xa , xb ) ∈ Rn × Rn

alors le schéma est stable avec S(T ) = eLT .

Démonstration On a alors

kxj+1 − z j+1 k ≤ kδ j+1 k + (1 + ∆tj L)kxj − z j k ≤ kδ j+1 k + e∆tj L kxj − z j k

puisque 1 + x ≤ ex pour tout x ∈ R. Par récurrence, on montre alors que pour

tout 1 ≤ j ≤ J,
j
X
kxj − z j k ≤ e(tj −t0 )L kx0 − z 0 k + e(tj −tk )L kδ k k .
k=1

Il s’ensuit que
j
!
X
j j TL 0 0 k
kx − z k ≤ e kx − z k + kδ k ,
k=1

ce qui donne le résultat. 

Convergence

La combinaison de consistance et de stabilité donne une propriété dite de

convergence qui dit que l’erreur commise par le schéma par rapport à la vraie
solution converge vers 0 lorsque le pas de temps converge vers 0. C’est une
propriété cruciale pour un schéma numérique.

Définition – Convergence

Soit ∆t = max0≤j≤J−1 ∆tj . Un schéma numérique est convergent si

lim max kxj − x(tj )k = 0

∆t→0 1≤j≤J

14
lorsque x0 = x(t0 ). S’il existe p ∈ N>0 et cv > 0 (indépendent de ∆t) tel que
max kxj − x(tj )k ≤ cv (∆t)p
1≤j≤J

on dit que le schéma est convergent à l’ordre p.

Théorème – Théorème de Lax

Une méthode stable et consistante (à l’ordre p) est convergente (à l’ordre p).

Démonstration Notons z j = x(tj ). On remarque que

z j+1 = z j + ∆tj Φ(tj , z j , ∆tj ) + ∆tj η j+1 ,
où η est l’erreur de consistance. D’après la propriété de stabilité, on a donc
J
X
kxj − x(tj )k ≤ S(T ) ∆tj−1 kη j k ,
j=1

et par consistance
J
X
kxj − x(tj )k ≤ S(T ) cs ∆tj−1 (∆tj−1 )p ≤ cs S(T ) T (∆t)p .
j=1

Théorème – Condition suffisante de convergence

L’inconvénient du théorème de Lax est qu’il faut prouver la stabilité pour obtenir
la convergence. Or la seule condition suffisante dont nous disposions à cet effet,
est le caractère globalement Lipschitzien de x 7→ Φ(t, x, ∆t). Mais il s’agit d’une
condition très forte. En fait, il est possible de prouver la convergence sous la
condition plus faible que x 7→ Φ(t, x, ∆t) est “localement Lipschitzienne” :
Si

1. le schéma est consistant d’ordre p,

2. pour tout boule fermée B de Rn , il existe L > 0, ∆tm > 0 tels que pour
tout t ∈ [0, T ] et pour tout ∆t ∈ [0, ∆tm ],
kΦ(tj , xa , ∆tj ) − Φ(tj , xb , ∆tj )k ≤ Lkxa − xb k ∀(xa , xb ) ∈ B × B

Alors il existe un pas de temps maximal ∆tmax > 0 tel que le schéma est
convergent d’ordre p.
L’hypothèse 2. est en particulier vérifiée si x 7→ Φ(t, x, ∆t) est C 1 d’après une
version un peu plus générale du théorème des accroissement finis.

15
Exemple – Convergence du schéma d’Euler explicite On a déjà montré
que le schéma d’Euler explicite est consistant d’ordre 1. Par ailleurs, si f est C 1
par rapport à x, alors x 7→ Φ(t, x) = f (t, x, ∆t) est C 1 . Donc d’après le théorème
précédent, le schéma est convergent d’ordre 1.

Erreurs d’arrondi et pas optimal

A chaque itération, lorsque la machine calcule xj+1 , elle commet en plus de

l’erreur de troncature de l’intégrale des erreurs d’arrondi de l’ordre de la précision
machine. La solution obtenue est donc en fait donnée par

x̂j+1 = x̂j + ∆tj Φ(tj , x̂j , ∆tj ) + ρj+1 + εj+1

au lieu de
xj+1 = xj + ∆tj Φ(tj , xj , ∆tj ) ,
où ρ modélise l’erreur commise sur le calcul de Φ et  l’erreur sur l’addition
finale. La stabilité nous donne alors l’écart
J
X
j j
max kx − x̂ k ≤ S(T ) ∆tj−1 kρj k + kεj k .
0≤j≤J
j=1

En considérant une borne ε des εj et ρ des ρj , on obtient

 
j j ε
max kx − x̂ k ≤ S(T )(T ρ + Jε) ≤ S(T )T ρ + ,
0≤j≤J min ∆tj
et donc finalement, en supposant l’algorithme convergent d’ordre p,

max kx(tj ) − x̂j k ≤ max kx(tj ) − xj k + kxj − x̂j k

0≤j≤J 0≤j≤J
 
p ε
≤ cv (max ∆tj ) + S(T )T ρ + .
j minj ∆tj
Les paramètres ε et ρ sont typiquement petits de l’ordre d’un facteur de la
précision machine. Cependant, on voit que plus le pas de temps décroit, plus il y
a d’itérations et plus les erreurs d’arrondi se propagent. D’un autre côté, plus
il augmente, plus les erreurs de quadrature augmentent. En supposant le pas
constant, il y a donc un pas “optimal’ ’ donné par
  1
S(T )T ε p+1
∆topt = .
cv p

Annexe – Choix du pas de temps

Jusqu’à présent, on a présenté des schémas dépendant de pas de temps ∆tj , sans
jamais dire comment les choisir. Le plus simple est de choisir un pas ∆t fixe

16
mais il est difficile de savoir à l’avance quel pas est nécessaire. En particulier,
comment savoir si la solution obtenue est suffisamment précise, sans connaître la
vraie ?

Pas fixe

Une voie empirique est de fixer un pas, lancer la simulation, puis fixer un pas
plus petit, relancer la simulation, jusqu’à ce que les résultats ne semble plus
changer (au sens de ce qui nous intéresse d’observer). Notons que la connaissance
des constantes de temps présentes dans le système peut aider à fixer un premier
ordre de grandeur du pas. On pourrait aussi directement choisir le pas ∆topt
obtenu plus haut en prenant en compte les erreurs d’arrondis. Mais les constantes
cv et S(T ) sont souvent mal connues et conservatives.

Adaptation du pas de temps

Les méthodes à pas fixe exploitent la convergence des schémas, mais

— on ne peut pas prendre un pas de temps arbitrairement petit car on est
contraint par le temps de simulation.
— on n’a aucune idée de l’erreur commise et on n’est jamais sûr d’avoir la
bonne solution.
— l’utilisation d’un pas très petit peut n’être nécessaire qu’autour de certains
points sensibles (proches de singularités par exemple) et consomme des
ressources inutiles ailleurs.
L’idée serait donc plutôt d’adapter la valeur du pas ∆tj à chaque itération. En
d’autres termes, on se fixe une tolérance d’erreur que l’on juge acceptable et on
modifie le pas de temps en ligne, selon si l’on estime être au-dessus ou en-dessous
du seuil d’erreur. Mais cela suppose d’avoir une idée de l’erreur commise. . . Il
existe justement des moyens de l’estimer.
Tout d’abord, de quelle erreur parle-t-on ?
— erreur globale ? L’idéal serait de contrôler max0≤j≤N kxj − x(tj )k. Or la
stabilité nous dit que
J
X
max kxj − x(tj )k ≤ S(T ) ∆tj−1 kη j k
0≤j≤N
j=1

avec η j les erreurs de consistances locales. Donc si on se fixe une tolérance

sur l’erreur globale Tolg , on a
Tolg
kη j k ≤ =⇒ max kxj − x(tj )k ≤ Tolg .
T S(T ) 0≤j≤N

17
 En d’autre termes, Tolg nous fixe une erreur maximale locale sur η j ,
à chaque itération. Notons cependant que cette borne ne prend pas en
compte la propagation des erreurs d’arrondis : plus ∆t diminue, plus
l’erreur globale risque d’augmenter. Ce phénomène devrait donc en toute
rigueur aussi nous donner un pas de temps minimal ∆tmin . Notons que
tous ces calculs dépendent des constantes cv et S(T ) qui sont souvent
mal connues ou très conservatives.
— erreur (absolue) locale ? A chaque itération, une erreur locale est commise
dûe à l’approximation de l’intégrale. Cette erreur est donnée par
!
Z tj+1
ej+1 = x̃(tj+1 ) − xj+1 = xj + f (s, x̃(s))ds
tj

où x̃ est la solution de ẋ = f (t, x) qui serait initialisée à xj au temps tj .

Notons que si on avait xj = x(tj ), on aurait exactement ej+1 = ∆tj η j+1 ,
où η j+1 est l’erreur de consistance. On se donne donc une tolérance
d’erreur locale
kej+1 k ≤ Tolabs .
— erreur relative ? Fixer une erreur absolue est parfois trop contraignant et
n’a de sens que si les solutions gardent un certain ordre de grandeur. En
effet, l’erreur acceptable quand la solution vaut 1000 n’est peut-être pas
la même que lorsqu’elle vaut 1. On peut donc plutôt exiger une certaine
erreur relative Tolrel , i.e.,

kej+1 k
≤ Tolrel .
kxj k

En général, les solvers assurent (approxivement)

kej+1 k ≤ Tolabs + Tolrel kxj k .

Par défaut, dans les solvers de Numpy, Tolabs = 10−6 et Tolrel = 10−3 .
Mais pour cela nous devons trouver un moyen d’estimer l’erreur locale. C’est
souvent fait en utilisant une même méthode à deux pas différents (par exemple
∆tj et ∆tj /2), ou bien en imbriquant des schémas de Runge-Kutta d’ordres
différents.
En fait, il est possible de montrer (exercice) que si f est C 1 , on a pour un schéma
d’Euler explicite

f (tj+1 , xj+1 ) − f (tj , xj )

kej+1 k = ∆tj + o(∆t2j )
2

On peut donc estimer à chaque itération l’erreur commise ej+1 et adapter

le pas selon si celle-ci est inférieure ou supérieure au seuil de tolérance. En
effet, puisque l’on sait par ailleurs que ej+1 = O(∆t2j ), une possible stratégie

18
d’adaptation est de prendre
r
Tolabs
∆tnew = ∆tj
kej+1 k
(éventuellement avec une marge de sécurité)
La fonction correspondante
def solve_euler_explicit_variable_step(f, x0, t0, tf, dtmin, dtmax, atol):
...
return t, x
est fournie dans le notebook Equations Differentielles [Link].

Projet numérique
Les équations de Lotka-Volterra, ou “modèle proie-prédateur”, sont couramment
utilisées pour décrire la dynamique de systèmes biologiques dans lesquels un
prédateur et sa proie interagissent dans un milieu commun. Elles ont été proposées
indépendamment par A. J. Lotka en 1925 et V. Volterra en 1926 et s’écrivent de
la manière suivante :

ẋ1 = x1 (α − βx2 )
ẋ2 = −x2 (γ − δx1 )

où x1 et x2 désignent le nombre (positif) de proies et de prédateurs respectivement

et α, β, γ, δ sont des paramètres strictement positifs.

1. Donner une interprétation physique à chaque terme de la dynamique.

Montrer qu’il existe deux points d’équilibre (0, 0) et x̄ ∈ R>0 × R>0 . Que
peut-on dire de leur stabilité à ce stade ?
2. A l’aide des fonctions meshgrid et quiver, visualiser graphiquement le
champ de vecteurs. Intuiter le comportement des solutions. On pourra
aussi utiliser streamplot pour visualiser le portrait de phase.
3. Par le théorème de Cauchy-Lipschitz, démontrer que toute solution initia-
lisée dans R>0 × R>0 reste dans R>0 × R>0 sur son ensemble de définition.
4. On considère la fonction

H(x1 , x2 ) = δx1 − γ ln x1 + βx2 − α ln x2

définie sur R>0 × R>0 . Calculer la dérivée de H le long des solutions

initialisées dans R>0 × R>0 . En déduire que toute solution maximale
initialisée dans R>0 × R>0 est définie sur R.
5. Représenter les courbes de niveau de H. Où se trouve x̄ ? Qu’en conclut-on
sur le comportement des solutions ? En déduire (graphiquement) que x̄
est stable, au sens de la définition de stabilité.

19
On souhaite maintenant simuler numériquement les trajectoires.

6. Coder une fonction du type

def solve_euler_explicit(f, x0, dt, t0, tf):
...
return t, x
prenant en entrée une fonction f : R×Rn → Rn quelconque, une condition
initiale x0 , un pas de temps dt, les temps initiaux et finaux, et renvoyant
le vecteur des temps tj et de la solution xj du schéma d’Euler explicite
appliqué à ẋ = f (t, x). La tester sur une équation différentielle aux
solutions exactes connues. Vérifier la convergence du schéma lorsque dt
tend vers 0. Comment visualiser graphiquement l’ordre de convergence ?
7. Utiliser le schéma d’Euler explicite pour simuler les équations de Lotka-
Volterra. Que constate-t-on en temps long ? Cette résolution vous semble-
t-elle fidèle à la réalité ? On pourra tracer l’évolution de la fonction H.
8. Coder maintenant une fonction du type
def solve_euler_implicit(f, x0, dt, t0, tf, itermax = 100):
...
return t, x
donnant la solution d’un schéma d’Euler implicite appliqué à ẋ = f (t, x)
selon la méthode présentée dans le cours. Vérifier de nouveau sa conver-
gence sur des solutions connues. Que se passe-t-il cette fois-ci sur les
équations de Lotka-Volterra ?

On propose maintenant de modifier ces schémas de façon à stabiliser H et assurer

sa conservation le long des solutions numériques.

9. Expliquer pourquoi les solutions de

ẋ1 = x1 (α − βx2 ) − u1 (x1 , x2 )(H(x1 , x2 ) − H0 )
ẋ2 = −x2 (γ − δx1 ) − u2 (x1 , x2 )(H(x1 , x2 ) − H0 )
sont identiques à celles de Lotka-Volterra si H0 = H(x(0)) pour tout
choix de u : R2 → R2 continûment différentiable.
10. Soit H0 ∈ R. Calculer la dérivée de H − H0 le long des solutions de ce
nouveau système. Montrer que l’on peut choisir u tel que
d
(H(x(t)) − H0 ) = −kk∇H(x(t))k2 (H(x(t)) − H0 ) .
dt
En déduire qu’alors H(x(t)) converge exponentiellement vers H0 lorsque
t tend vers l’infini si x reste à une distance strictement positive de x̄.
11. En déduire comment modifier l’implémentation du schéma d’Euler pour
assurer la stabilité de H. Quel est le rôle de k ? Peut-il être choisi arbi-
trairement grand ? Pourquoi ? On pourra exprimer H(xj+1 ) − H(x0 ) en
fonction de H(xj ) − H(x0 ) au premier ordre en dt.

20
Exercices

Consistance et ordre de schémas

Supposons f de classe C 2 . Montrer que :

Question 1 le schéma de Heun est consistant d’ordre ≥ 2 et égal à 2 si f est

C 3 . (Solution p. 23.)

Question 2 le schéma d’Euler implicite est consistant d’ordre ≥ 1 (Solution p.

24.)

Question 3 la méthode des trapèzes est consistante d’ordre ≥ 2. (Solution p.

25.)

Question 4 le schéma du point milieu est consistant d’ordre ≥ 2. (Solution p.

25.)

Question 5 le schéma de Runge-Kutta d’ordre 4 est bien consistant d’ordre 4

si f est C 5 . (Solution p. 26.)
On supposera le pas suffisamment petit pour que les schémas implicites soient
définis.

Convergence de schémas

Sous l’hypothèse que f est C 1 , montrer que les schémas de Heun et d’Euler
implicite sont convergents. (Solution p. 26.)

Explicite ou implicite ?

Question 1 Comparer les performances des schémas d’Euler implicites et

explicites à pas fixe dans le cas de ẋ = −λx, x(0) = 1, et ẋ = λx, x(0) = 1, sur
un horizon de temps T donné. (Solution p. 27.)

21
Question 2 Lorsqu’on modélise des systèmes chimiques ou biologiques, on
obtient souvent des réactions aux constantes de temps très différentes. Vaut-il
mieux utiliser un schéma d’Euler implicite ou explicite pour simuler
 
−1 0
ẋ = x
0 −µ

avec µ >> 1 ? (Solution p. 28.)

Euler symplectique

Pour ω > 0 donné, considérons le système

ẋ1 = x2 , ẋ2 (t) = −ω 2 x1

de condition initiale x(0) = (1, 0). On rappelle que pour une suite de la forme
xj+1 = Axj converge vers 0 si les valeurs propres de A sont à l’intérieur du cercle
unité et diverge si au moins une valeur propre est à l’extérieur.

Question 1 Montrer que pour n’importe quel pas ∆t fixé, un schéma d’Euler
explicite donne une solution divergente, et un schéma d’Euler implicite donne
une solution qui converge vers 0. Lequel a raison ? (Solution p. 28.)
On définit maintenant le schéma suivant qui “mélange’ ’ les schémas d’Euler
implicites et explicites :

xj+1
1 = xj1 + ∆t xj2
xj+1
2 = xj2 − ∆t ω 2 xj+1
1

Question 2 Montrer que la quantité ω 2 x21 + x22 + ∆t ω 2 x1 x2 est conservée.

Quelle est alors la forme des solutions obtenues dans le plan de phase si ω∆t < 2 ?
En déduire la pertinence de ce schéma. On parle de schéma symplectique, car il
conserve les volumes. (Solution p. 28.)

Question 3 En écrivant le schéma sous la forme xj+1 = Axj , montrer qu’il

diverge par contre si ∆t ω > 2. (Solution p. 29.)

Question 4 Plus généralement, proposer un schéma pour simuler un système

Hamiltonien du type

q̇ = ∇p H(q, p)
ṗ = −∇q H(q, p)

22
où (q, p) ∈ RN × RN sont les positions généralisées et quantités de mouvement,
H est le Hamiltonien que l’on pourra vérifier être conservé le long des trajectoires.
(Solution p. 29.)
A noter que les conclusions de cet exercice sont les mêmes si l’on utilise un schéma
d’Euler implicite sur la première composante et un schéma d’Euler explicite sur
la deuxième. Ces deux schémas s’appellent respectivement Euler symplectique A
et B.

Corrections

Consistance et ordre de schémas

Vu que f est C 2 , la dérivée seconde (en temps) des solutions s’écrit

d 
ẍ(t) = f (t, x(t)) = ∂t f (t, x(t)) + ∂x f (t, x(t))ẋ(t)
dt
= ∂t f (t, x(t)) + ∂x f (t, x(t))f (t, x(t)).

On a donc
f [1] (t, x) = ∂t f (t, x) + ∂x f (t, x)f (t, x)
et la formule de Taylor le long des solutions donne

x(tj+1 ) = x(tj ) + ∆tf (tj , x(tj ))

∆t2  
+ ∂t f (tj , x(tj )) + ∂x f (tj , x(tj ))f (tj , x(tj )) + O(∆t3 ).
2
(1)
Pour les calculs de consistance, deux options possibles:
— soit on compare à la main (1) aux développements en puissances de ∆t
de la solution numérique xj+1 , partant de x(tj ).
— soit on utilise la condition nécessaire et suffisante de consistance (p. 12).

Question 1 A la main, pour la méthode de Heun,

 i
j+1 1 ∆t h 
η = x(tj+1 ) − x(tj ) − f (tj , x(tj )) + f tj+1 , x(tj ) + ∆tf (tj , x(tj )) .
∆t 2

Or,
 
f tj+1 ,x(tj ) + ∆tf (tj , x(tj ))
= f (tj , x(tj )) + ∂t f (tj , x(tj ))(tj+1 − tj ) + ∂x f (tj , x(tj ))∆tf (tj , x(tj )) + O(∆t2 )
 
= f (tj , x(tj )) + ∆t ∂t f (tj , x(tj )) + ∂x f (tj , x(tj ))f (tj , x(tj )) + O(∆t2 ).

23
On en déduit que η j+1 = O(∆t2 ) en utilisant (1). Donc on a une consistance
d’ordre ≥ 2.
Sinon il suffit de constater que
∂Φ 1
Φ(t, x, 0) = f (t, x) , (t, x, 0) = f [1] (t, x) .
∂∆t 2
2
∂ Φ 1 [2]
Par contre, si f est C 3 , ∂∆t 2 (t, x, 0) 6= 3 f (t, x) donc le schéma n’est pas d’ordre
3 et donc il est d’ordre égal exactement à 2.

Question 2 Pour le schéma d’Euler implicite

x(tj+1 ) − x(tj ) − ∆tf (tj+1 , xj+1 )

η j+1 = ,
∆t
où xj+1 est solution de

xj+1 = x(tj ) + ∆tf tj+1 , xj+1 .

Comme vu dans la section Définition implicite de Φ (p. 8), cette définition

implicite de xj+1 admet une unique solution pour ∆t suffisamment petit si f est
Lipschitzienne par rapport à x. On aimerait dire ici que xj+1 = x(tj ) + O(∆t)
pour ∆t suffisamment petit. Pour faire ça proprement, fixons (tj , x(tj )) et posons


F (x, ∆t) = x − x(tj ) + ∆tf (tj + ∆t, x)

qui est de classe C 1 . On a alors

F (x(tj ), 0) = 0 , F (xj+1 , ∆t) = 0 .

Puisque ∂x F (x(tj ), 0) = Id est inversible, le théorème des fonctions implicites

nous dit que pour ∆t suffisamment petit, il existe une fonction ψ de classe C 1
telle que
F (xj+1 , ∆t) = 0 , xj+1 = ψ(∆t)
au voisinage de (x(tj ), 0). Puisque ψ est continue, il s’ensuit donc bien que
xj+1 = x(tj ) + O(∆t). On a donc

f (tj+1 , xj+1 ) = f (tj , x(tj )) + ∂t f (tj , x(tj ))(tj+1 − tj ) + ∂x f (tj , x(tj ))(xj+1 − x(tj )) + O(khk2 )
= f (tj , x(tj )) + O(∆t)

avec l’incrément h = (tj+1 − tj , xj+1 − x(tj )) qui vérifie khk = O(∆t).

Ainsi, toujours au vu de (1),
1 h  i
η j+1 = x(tj+1 ) − x(tj ) + ∆tf (tj , x(tj )) + O(∆t2 ) = O(∆t).
∆t

24
Question 3 Pour la méthode des trapèzes, on a l’erreur
∆t


x(tj+1 ) − x(tj ) − f (tj , x(tj )) + f (tj+1 , xj+1 )
η j+1 = 2
,
∆t
où xj+1 est solution de
∆t
xj+1 = x(tj ) + f (tj , x(tj )) + f (tj+1 , xj+1 ) .


2

Cette fois-ci on voudrait montrer que xj+1 = x(tj ) + ∆tf (tj , x(tj )) + O(∆t2 ).
Pour cela, on redéfinit

∆t
f (tj , x(tj )) + f (tj + ∆t, xj+1 )



F (x, ∆t) = x − x(tj ) +
2
qui est de classe C 1 . On a alors

F (x(tj ), 0) = 0 , F (xj+1 , ∆t) = 0 .

Puisque ∂x F (x(tj ), 0) = Id est inversible, le théorème des fonctions implicites

nous dit que pour ∆t suffisamment petit, il existe une fonction ψ de classe C 1
telle que
F (xj+1 , ∆t) = 0 , xj+1 = ψ(∆t)
au voisinage de (x(tj ), 0) et de plus,

ψ 0 (0) = Id−1 · ∂∆t F (x(tj ), 0) = f (tj , x(tj )) .

Puisque ψ est de classe C 1 , on a donc bien

xj+1 = x(tj ) + ∆tf (tj , x(tj )) + O(∆t2 ) .

Il s’ensuit que

f (tj+1 , xj+1 ) = f (tj , x(tj ))+∆t (∂t f (tj , x(tj )) + ∂x f (tj , x(tj ))f (tj , x(tj )))+O(∆t2 )

soit

∆t2
xj+1 = x(tj )+∆tf (tj , x(tj )+ (∂t f (tj , x(tj )) + ∂x f (tj , x(tj ))f (tj , x(tj )))+O(∆t3 )
2

On en déduit donc η j+1 = O(∆t2 ).

Question 4 Enfin, pour la méthode du point milieu, on a l’erreur

x(tj )+xj+1
 
x(tj+1 ) − x(tj ) − ∆tf tj + ∆t
2 , 2
η j+1 = ,
∆t

25
où xj+1 est solution de
∆t x(tj ) + xj+1
 
j+1
x = x(tj ) + ∆tf tj + , .
2 2
On montre de la même façon qu’à la question précédente que pour ∆t suffisam-
ment petit, xj+1 = x(tj ) + ∆tf (tj , x(tj )) + O(∆t2 ) et donc
∆t x(tj ) + xj+1
 
f tj + ,
2 2
x(tj ) + xj+1
 
∆t
= f (tj , x(tj )) + ∂t f (tj , x(tj )) + ∂x f (tj , x(tj )) − x(tj ) + O(khk2 )
2 2
∆t
= f (tj , x(tj )) + (∂t f (tj , x(tj )) + ∂x f (tj , x(tj ))f (tj , x(tj ))) + O(∆t2 )
2
x(tj )+xj+1
 
avec l’incrément h = ∆t 2 , 2 − x(t j ) , soit

∆t2
xj+1 = x(tj )+∆tf (tj , x(tj )+ (∂t f (tj , x(tj ))+∂x f (tj , x(tj ))f (tj , x(tj )))+O(∆t3 )
2
On en déduit donc η j+1 = O(∆t2 ).

Question 5 On vérifie que le critère est vérifié pour 0 ≤ k ≤ 3 mais pas pour
k = 4!

Convergence de schémas

Tout d’abord, dans l’exercice précédent, nous avons montré que les schémas de
Heun et d’Euler implicite étaient consistants d’ordre 2 et 1 respectivement. Il ne
nous reste donc plus qu’à montrer que Φ est localement lipschitzienne (ou C 1 )
par rapport à x pour ∆t suffisamment petit, pour en déduire la convergence à
l’ordre 2 et 1 respectivement.
Pour le schéma de Heun,
1

Φ(t, x, ∆t) = f (t, x) + f t + ∆t, x + ∆tf (t, x)
2
donc Φ est C 1 par rapport à x si f l’est.
Prenons maintenant le schéma d’Euler implicite. Pour ∆t ≤ ∆tm , Φ est définie
par  
Φ(t, x, ∆t) = f t + ∆t, x + ∆tΦ(t, x, ∆t) .
Soit B un compact de Rn . Soit B 0 un compact tel que x + ∆tΦ(t, x, ∆t) ∈ B 0
pour tout x ∈ B, tout t ∈ [0, T ] et tout ∆t ∈ [0, ∆tm ]. Puisque f est continue,
et C 1 par rapport à x, il existe Lf > 0 tel que
kf (t, xa ) − f (t, xb )k ≤ Lf kxa − xb k ∀(xa , xb , t) ∈ B 0 × B 0 × [0, T + ∆tm ] .

26
Ceci est vrai par le théorème des accroissements finis appliqué à x 7→ f (t, x) et
pour t dans un intervalle fermé et borné. Prenons maintenant (xa , xb ) ∈ B × B,
t ∈ [0, T ] et ∆t ∈ [0, ∆tm ], alors

kΦ(t, xa , ∆t) − Φ(t, xb , ∆t)k

   
= kf t + ∆t, xa + ∆tΦ(t, xa , ∆t) − f t + ∆t, xb + ∆tΦ(t, xb , ∆t) k
≤ Lf (kxa − xb k + ∆tkΦ(t, xa , ∆t) − Φ(t, xb ∆t)k)

soit
Lf
kΦ(t, xa , ∆t) − Φ(t, xb , ∆t)k ≤ kxa − xb k
1 − Lf ∆t
si ∆t < 1/Lf . Donc Φ est bien localement lipschitzienne par rapport à x.

Explicite ou implicite ?

Question 1 Prenons d’abord ẋ = −λx, x(0) = 1, dont la solution exacte est

x(t) = e−λt .
Le schéma d’Euler explicite donne

xj+1 = xj − λ∆txj = (1 − λ∆t)j

soit  λ
T T
xJ = (1 − λ∆t)J = (1 − λ∆t) ∆t = (1 − λ∆t) λ∆t .

On a bien
lim xJ = e−λT .
∆t→0

Cependant, il faut |1 − λ∆t| < 1 pour que la solution converge au moins vers
0. Sinon, pour λ∆t = 2, xJ = (−1)J , qui n’a rien à voir avec la solution. Pire,
pour λ∆t = 2, l’algorithme diverge. Il faut donc adapter ∆t à la constante de
temps λ du système. Ceci peut poser problème lorsque l’on simule des systèmes
sur des temps longs (par rapport à λ)
De l’autre côté, le schéma d’Euler implicite donne

xj+1 = xj − λ∆txj+1

soit
1 1
xJ = = T
(1 + λ∆t)J (1 + λ∆t) ∆t
qui tend vers 0 quelque soit le pas ∆t ! On parle de stabilité inconditionnelle.
Ceci est très pratique pour des simulations sur temps longs, où la condition
λ∆t < 1 est trop contraignante.

27
Prenons maintenant ẋ = λx, x(0) = 1, dont la solution exacte est x(t) = eλt .
Cette fois-ci, Euler explicite donne
T
xJ = (1 + λ∆t) ∆t

qui fait maintenant sens même pour des pas grands. Par contre, Euler implicite
donne
1
xJ = T
(1 − λ∆t) ∆t
qui n’est pas défini pour λ∆t = 1 et qui explose pour des valeurs proche de 1.

Question 2 Lorsque l’on a deux dynamiques asymptotiquement stables aux

constantes de temps très différentes la condition de stabilité de Euler explicite
exige de choisir un pas câlé sur la plus petite constante de temps, i.e. il faut
∆t < µ1 . Ceci est très exigeant car il faut attendre un nombre d’itérations de
l’ordre de µ pour voir l’évolution du système lent. Par contre, une méthode
implicite permet de choisir librement le pas de temps en fonction des performances
souhaitées.

Euler symplectique

Question 1 Dans le cas d’Euler explicite, xj+1 = Axj avec

 
1 ∆t
A=
−∆t ω 2 1
p
dont les valeur propres sont 1 ± iω∆t de norme 1 + ∆t2 ω 2 > 1. Donc les
solutions divergent.
Dans le cas d’Euler implicite, xj+1 = Axj avec
 
1 1 ∆t
A= 2 2
1 + ∆t ω −∆t ω2 1
p
dont les valeurs propres sont 1/(1 ± iω∆t) de norme 1/ 1 + ∆t2 ω 2 < 1. Donc
les solutions convergent vers 0.
Or on peut vérifier que le long des vraies solutions, l’énergie ω 2 x21 + x22 est
constante donc les trajectoires sont bornées et ne peuvent pas converger vers
zéro. Aucun des deux schémas n’approxime les solutions correctement sur le
long-terme.

Question 2 On vérifie par le calcul que

!
∆t w2
> ω2
ω 2 x21 + x22 2
+ ∆t ω x1 x2 = x ∆t w2
2 x
2 1

28
 2 4
est constante. Pour ω 2 − ∆t4ω > 0, soit ω∆t < 2, cette matrice est définie
positive, donc les solutions restent sur une ellipse. Cette ellipse se rapproche
de la vraie solution lorsque ∆t tend vers 0. Ce schéma est donc approprié pour
simuler les trajectoires sur un temps long.

Question 3 L’algorithme symplectique est décrit par xj+1 = Axj avec

 
1 ∆t
A=
−∆t ω 2 1 − ∆t2 ω 2
dont le polynôme caractéristique s’écrit
s2 − (2 − ∆t2 ω 2 )s + 1
On a les cas suivants :
— si (1−∆t2 ω 2 )2 −4 < 0, i.e., si ω∆t < 2, les valeurs propres sont imaginaires
conjuguées et de module 1.
— si ω∆t > 2, les valeurs propres sont réelles de produit 1, donc l’une est
supérieure à 1 est le schéma diverge.
— dans le cas extrême où ω∆t = 2, il y a une valeur propre double en -1.

Question 4 Pour un système hamiltonien, on peut donc proposer

q j+1 = q j + ∆t ∇p H(q j , pj )
pj+1 = pj − ∆t ∇q H(q j+1 , pj+1 )
ou bien
q j+1 = q j + ∆t ∇p H(q j+1 , pj+1 )
pj+1 = pj − ∆t ∇q H(q j , pj )
pour ∆t suffisamment petit.

Références
Demailly, J.-P. 2006. Analyse Numérique et équations Différentielles. EDP
Sciences. Grenoble Sciences.
Hairer, E., C. Lubich, and G. Wanner. 2010. Geometric Numerical Integration :
Structure-Preserving Algorithms for Ordinary Differential Equations. Edited
by Springer Series in Computational Mathematics. 2nd ed. Springer-Verlag,
Berlin.
Hairer, E., and G. Wanner. 1996. Solving Ordinary Differential Equations II.
Stiff and Differential-Algebraic Problems. Edited by Springer Series in Com-
putational Mathematics. 2nd ed. Springer-Verlag, Berlin.

29