Univ.
Paris Cité 2024-2025 M1 MMA, Modélisation Déterministe, Projet
Modélisation Déterministe
Projet - Informations Complémentaires
Ce document présente un certain nombre de notions utiles pour réaliser le projet.
Ce document ne comporte aucune question qui sera évaluée.
Table des matières
1 Introduction 2
2 Stabilité, consistance, et convergence d’un schéma 3
2.1 Définition de la convergence d’un schéma numérique . . . . . . . . . . . . . . 3
2.2 Stabilité . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
2.3 Consistance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.4 Convergence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.5 Choix du schéma en pratique . . . . . . . . . . . . . . . . . . . . . . . . . . 6
3 Quelques schémas numériques 7
3.1 Méthode d’Euler explicite . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
3.2 Méthode d’Euler implicite . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
3.3 Méthodes de Crank-Nicolson et du point milieu . . . . . . . . . . . . . . . . 10
1
Univ. Paris Cité 2024-2025 M1 MMA, Modélisation Déterministe, Projet
1 Introduction
Considérons le problème de Cauchy
(
u0 (t) = f (u(t)), t ∈ [0, T ],
(PdC)
u(0) = u0 ,
où T est un temps final fixé (T = +∞ si la solution est globale, T < +∞ sinon 1 ), et f est
une fonction L-lipschitzienne sur son domaine de définition.
Si on ne sait pas résoudre le problème (PdC), pour déterminer ce que “vaut la solution” il
faudrait être capable de donner les valeurs u(t) de la solution pour n’importe quel temps
t ∈ [0, T ]. On rencontre alors deux problèmes :
(P1) il n’est pas possible d’accéder à l’ensemble non dénombrable des valeurs u(t), t ∈ [0, T ] ;
(P2) en général, on ne dispose pas d’une formule explicite définissant u.
Pour résoudre le problème (P1), on restreint l’estimation à un sous-ensemble dénombrable
de [0, T ]. On définit un nombre fini N + 1 de temps, notés tn , régulièrement 2 espacés d’un
pas de temps h, avec h = T /N , et on essaie d’estimer u en chaque temps tn . On cherche
donc à estimer les grandeurs
u(tn ), tn = nh, n ∈ {0, . . . , N }.
L’ensemble des valeurs tn est appelé maillage de l’intervalle [0, T ], les temps tn sont les
nœuds du maillage. L’intervalle [0, T ] est “découpé” en N sous-intervalles [tn , tn+1 ] grâce
aux N + 1 nœuds du maillage. La valeur de la solution étant donnée au premier nœud (t0 )
par la condition initiale u0 , on cherche à estimer N valeurs de la solution : u(t1 ),. . ., u(tN ).
Le problème (P2) est – malheureusement – classique, au sens où on dispose rarement d’une
formule explicite des solutions. On essaye alors de proposer une méthode d’estimation des
grandeurs u(tn ), et de déterminer la fiabilité de cette méthode afin de savoir si on peut avoir
confiance en elle (ou plus souvent, à quel point on peut avoir confiance en elle).
Un schéma numérique va consister à construire des grandeurs notées un censées appro-
cher les valeurs u(tn ). Idéalement, on souhaiterait que cette approximation soit bonne et
s’améliore lorsque le pas de temps h est réduit.
1. dans le cas où la solution maximale n’est pas globale, on a vu que le domaine de définition de la
solution est nécessairement ouvert à droite, on devrait donc écrire [0, T [ ; on utilisera toutefois un intervalle
fermé dans ce document car T correspondra toujours à la plus grande valeur pour laquelle on calculera une
approximation numérique de la solution de l’EDO.
2. On peut être amené à considérer des temps qui ne seraient pas régulièrement espacés, mais on ne
considèrera ici que ce cas là.
2
Univ. Paris Cité 2024-2025 M1 MMA, Modélisation Déterministe, Projet
2 Stabilité, consistance, et convergence d’un schéma
On suppose dans le reste de ce document que la fonction f dans (PdC) est continue et
globalement lipschitzienne, c’est-à-dire il existe une constante L > 0 telle que
∀(u, v) ∈ (Rd )2 , kf (u) − f (v)k ≤ Lku − vk.
On se place donc dans les conditions du théorème de Cauchy-Lipschitz global, et on note
dans toute la suite u l’unique solution globale de (PdC) définie sur [0, T ].
Considérons une définition abstraite d’un schéma numérique, permettant de décrire une
grande variété de schémas (en particulier, les exemples de schémas donnés dans la section 3).
Ce cadre est celui des schémas à un pas.
Définition 2.1 (Schéma à un pas). On dit qu’un schéma est à un pas s’il s’écrit sous la
forme
un+1 = un + h Φ(un+1 , un , h), (2.1)
et s’il existe hmax > 0 tel que cette équation définisse un+1 de manière unique pour tout
h ≤ hmax .
Dans cette définition, Φ est une fonction de Rd × Rd × [0, hmax ] à valeurs dans Rd .
En pratique : Avant d’utiliser un schéma numérique, il faut démontrer que ce schéma
est bien défini, c’est-à-dire qu’il permet de définir une unique solution approchée. Cette
question ne se pose réellement que lorsque le schéma est implicite.
2.1 Définition de la convergence d’un schéma numérique
Définition 2.2 (Convergence d’un schéma). On dit qu’un schéma est convergent si
lim max kun − u(tn )k = 0.
h→0 0≤n≤N
Le schéma est dit convergent d’ordre p si la vitesse de convergence est en O(hp ), c’est-à-dire
s’il existe hmax et une constante C > 0 indépendante de h ∈ [0, hmax ] telle que
∀h ∈ [0, hmax ], max kun − u(tn )k ≤ Chp .
0≤n≤N
On va à présent essayer de répondre à la question : comment montrer la convergence d’un
schéma numérique à un pas ?
3
Univ. Paris Cité 2024-2025 M1 MMA, Modélisation Déterministe, Projet
2.2 Stabilité
Les deux notions centrales permettant de conclure à la convergence d’un schéma sont la
stabilité et la consistance.
Définition 2.3 (Stabilité). On dit qu’un schéma à un pas (2.1) est stable s’il existe une
constante M > 0 telle que, pour tout h suffisamment petit, toute suite (εn )0≤n≤N −1 , toutes
suites un , ũn définies pour 0 ≤ n ≤ N − 1 par les récurrences
un+1 = un + h Φ(un+1 , un , h),
ũn+1 = ũn + h Φ(ũn+1 , ũn , h) + εn ,
u0 = ũ0 ,
on a nécessairement X
max kun − ũn k ≤ M kεn k.
0≤n≤N
0≤n≤N −1
Autrement dit, un schéma est dit stable s’il n’amplifie pas outre mesure les erreurs commises
en approchant les solutions. Le théorème suivant permet de conclure à la stabilité d’un
schéma numérique.
Theorème 2.1. On suppose que Φ est globalement lipschitzienne par rapport à ses deux
premières variables, uniformément par rapport à la troisième, c’est-à-dire qu’il existe une
constante K > 0 telle que ∀h ∈ [0, hmax ] et ∀(u1 , u2 , v1 , v2 ) ∈ R4d
kΦ(v1 , u1 , h) − Φ(v2 , u2 , h)k ≤ K (ku1 − u2 k + kv1 − v2 k) .
Alors le schéma (2.1) est bien défini, stable pour tout h ≤ hmax , et hmax < 1/K.
En quoi la notion de stabilité va-t-elle permettre d’établir la convergence ? Elle est absolu-
ment centrale, car pour que les écarts entre la valeur approchée un et la valeur exacte u(tn )
tendent vers 0, il faut que les cumuls d’erreurs effectuées au fur et à mesure des itérations
soient contrôlés. Il est donc nécessaire qu’un schéma soit stable pour qu’il converge.
Remarquons toutefois que la notion de stabilité est totalement indépendante du problème
de Cauchy : elle dit juste qu’on a défini une récurrence qui n’amplifie pas les erreurs, et rien
de plus. Par exemple, le schéma un+1 = un (avec Φ identiquement nulle) est parfaitement
stable, mais il n’a aucun lien avec le problème qu’on veut résoudre.
En pratique : La stabilité se démontre en général facilement à condition de connaître
quelques résultats de base sur les fonctions lipschitziennes.
4
Univ. Paris Cité 2024-2025 M1 MMA, Modélisation Déterministe, Projet
2.3 Consistance
Considérons à présent la notion de consistance. Il s’agit de la propriété qui va faire le lien
entre le schéma numérique stable et la résolution de l’EDO.
Définition 2.4 (Erreur de consistance). On appelle erreur de consistance la suite de Rd
définie à partir de la solution de l’EDO u par
u(tn+1 ) − u(tn )
rn = − Φ(u(tn+1 ), u(tn ), h), 0 ≤ n ≤ N − 1.
h
Définition 2.5 (Consistance). On dit que le schéma est consistant si
lim max krn k = 0.
h→0 0≤n≤N −1
Le schéma est dit consistant d’ordre p si la vitesse de convergence est en O(hp ) : autrement
dit, s’il existe une constante C > 0 indépendante de h telle que
max krn k ≤ Chp .
0≤n≤N −1
La proposition suivante permet de conclure à la consistance d’un schéma explicite, unique-
ment, et sans donner l’ordre de convergence. Ca reste un résultat intéressant et pratique.
Proposition 2.1. Supposons que le schéma à un pas est explicite, c’est-à-dire s’écrit
un+1 = un + hΦ(un , h),
avec Φ : Rd × [0, hmax ] 7→ Rd . Alors il est consistant si et seulement si, pour tout u ∈ Rd ,
Φ(u, 0) = f (u).
En pratique : On peut estimer l’erreur de consistance grâce à des développements de
Taylor. Ceux-ci requièrent suffisamment de régularité sur la fonction f , qui se traduisent en
de la régularité sur la fonction u. En effet, on a le lemme suivant :
Lemme 2.1. Soit k ∈ N. Si f ∈ C k (Rd ), alors la solution u de (PdC) est C k+1 ([0, T ]).
Démonstration. La preuve se fait par récurrence : pour k = 0, f étant continue, la dérivée
de u est égale à la composée de f , qui est continue, et de t 7→ (t, u(t)), tout autant continue
(u est au moins dérivable pour être solution de l’EDO), donc u0 est continue, ce qui implique
que u est de classe C 1 . La récurrence se poursuit alors avec des arguments similaires.
Pour le calcul pratique de la consistance, l’idée est de tout ramener à un même temps
(tn ou tn+1 , par exemple), et de penser à utiliser que u satisfait l’équation différentielle
u0 (t) = f (u(t)).
5
Univ. Paris Cité 2024-2025 M1 MMA, Modélisation Déterministe, Projet
2.4 Convergence
Le théorème essentiel permettant de montrer la convergence d’un schéma, appelé parfois
Théorème de Lax 3 , peut à présent être énoncé et démontré.
Theorème 2.2 (Théorème de Lax). Un schéma (2.1) stable et consistant est convergent.
S’il est de plus consistant d’ordre p, alors il est convergent d’ordre p.
Démonstration. Supposons que le schéma (2.1) est stable et consistant. Pour rappel (Defi-
nition 2.2), on veut montrer que
lim max kun − u(tn )k = 0.
h→0 0≤n≤N
Bien que la suite (u(tn ))0≤n≤N ne vérifie pas la relation de récurrence définissant le schéma,
on peut tout à fait écrire
u(tn+1 ) = u(tn ) + h Φ(u(tn+1 ), u(tn ), h)
u(tn+1 ) − u(tn )
+h − Φ(u(tn+1 ), u(tn ), h)
h
= u(tn ) + h Φ(u(tn+1 ), u(tn ), h) + h rn .
Autrement dit, cette suite vérifie la récurrence à une erreur près, qui est donnée par la suite
(h rn )n∈{0,...,N } . Notons
εn := h rn , n ∈ {0, . . . , N }.
La stabilité du schéma numérique implique (on utilise u(tn ) = ũn dans la définition de la
stabilité)
X
max kun − u(tn )k ≤ M kεn k
0≤n≤N
0≤n≤N −1
X
= Mh krn k
0≤n≤N −1
≤ M h N max krn k
0≤n≤N
= M T max krn k.
0≤n≤N
Le schéma étant, par hypothèse, consistant, on obtient bien par passage à la limite la
convergence. On note que l’ordre de la consistance définit l’ordre de la convergence.
2.5 Choix du schéma en pratique
En pratique, si on veut résoudre numériquement une EDO, on va se demander quel schéma
choisir.
3. Peter Lax, mathématicien américain, d’origine hongroise, toujours en vie en 2024 (98 ans).
6
Univ. Paris Cité 2024-2025 M1 MMA, Modélisation Déterministe, Projet
Un premier critère pertinent pour choisir un schéma est bien évidemment son ordre de
convergence : plus l’ordre de convergence est élevé, meilleure sera a priori l’approximation.
Il est donc crucial de savoir démontrer qu’un schéma est convergent et à quel ordre.
Toutefois, si le seul critère d’utilisation des schémas est celui de leur ordre de convergence, il
semble qu’il faille systématiquement choisir un schéma d’ordre élevé, voire le schéma d’ordre
le plus élevé possible s’il en existe un. Quelle est alors l’utilité de schémas d’ordre peu élevé ?
Et s’il y a plusieurs schémas d’ordre de convergence “le plus élevé possible”, comment choisir
entre deux schémas du même ordre ?
Notons qu’en plus d’avoir un ordre de convergence élevé, on attend d’un schéma numérique
que son implémentation soit économe numériquement, c’est-à-dire que les calculs permettant
de construire une solution approchée ne soient pas trop longs. Il s’agit du deuxième critère
pour sélectionner un schéma : son coût de calcul.
En effet, il se peut que l’ensemble des itérations d’un schéma d’ordre 1 soit réalisé en moins
de temps que les itérations d’un schéma d’ordre 2. Il existe donc en général un compromis à
trouver entre la précision du schéma et sa rapidité d’exécution. Selon l’application, cela peut
conduire à préférer un schéma d’ordre moindre selon les ressources de calcul disponibles.
On sait construire des schémas d’ordre arbitrairement grand, mais ceux-ci sont de plus en
plus coûteux en temps de calcul.
3 Quelques schémas numériques
Remarquons que la solution exacte u de l’EDO (PdC) vérifie en chaque noeud tn d’un
maillage donné de l’intervalle [0, T ],
u0 (tn ) = f (u(tn )) (3.1)
mais également, par intégration de (PdC) sur l’intervalle [tn , tn+1 ],
Z tn+1
u(tn+1 ) = u(tn ) + f (u(s))ds. (3.2)
tn
Il ne s’agit pas d’approximations : ces égalités sont vraies car u est la solution exacte de
l’EDO et que tn ∈ [0, T ] pour tout n.
Schématiquement, une méthode numérique pour construire une solution approchée de (PdC)
va consister
1) soit à construire une formule d’approximation de la dérivée dans (3.1),
2) soit à construire une formule d’approximation de l’intégrale dans (3.2).
La mise en œuvre des méthodes se fait alors systématiquement en deux temps :
Etape 1 : on écrit l’approximation de l’intégrale ou de la dérivée,
Etape 2 : on “remplace” les valeurs u(tn ) (inconnues) par un , ce qui mène à une équation
de récurrence pour calculer les valeurs un .
7
Univ. Paris Cité 2024-2025 M1 MMA, Modélisation Déterministe, Projet
3.1 Méthode d’Euler explicite
Considérons le schéma numérique s’appuyant sur une méthode de quadrature standard, la
méthode des rectangles à gauche (Figure 1.A).
Cette méthode correspond à l’approximation
Z b
F (s)ds ≈ (b − a)F (a).
a
Pour un intervalle de taille b − a = h, on obtient l’approximation
Z tn+1
f (u(s))ds ≈ (tn+1 − tn )f (u(tn )) = hf (u(tn )).
tn
En remplaçant les valeurs u(tn ) par leurs valeurs approchées un dans l’expression (3.2), on
obtient alors
un+1 = un + hf (un ), n ∈ {0, . . . , N − 1}.
Il s’agit du schéma d’Euler explicite. Ce schéma est dit explicite car on peut explicitement
calculer un+1 à partir de un grâce à la formule ci-dessus.
On peut également aboutir à ce schéma à partir de l’expression (3.1) en utilisant l’approxi-
mation de la dérivée u0 (tn ) suivante (Figure 2.A)
u(tn+1 ) − u(tn )
u0 (tn ) ≈ .
h
3.2 Méthode d’Euler implicite
En approchant cette fois l’intégrale dans (3.2) par la méthode des rectangles à droite (Figure
1.B), c’est-à-dire en utilisant l’approximation
Z b
F (s)ds ≈ (b − a)F (b),
a
on obtient l’approximation
Z tn+1
f (u(s))ds ≈ (tn+1 − tn )f (u(tn+1 )) = hf (u(tn+1 )).
tn
En remplaçant les valeurs u(tn ) par leurs valeurs approchées un dans l’expression (3.2), on
obtient alors
un+1 = un + hf (un+1 ), n ∈ {0, . . . , N − 1}.
Il s’agit du schéma d’Euler implicite. Ce schéma est dit implicite car on ne peut pas calculer
un+1 à partir de un à l’aide de la formule ci-dessus. Si un est connu, trouver un+1 revient
à résoudre une équation d’inconnue un+1 , a priori non-linéaire, ce qui induit des difficultés
théoriques et numériques supplémentaires : il faut à la fois s’assurer que cette équation
admet une unique solution, et par ailleurs avoir recours à une méthode numérique adaptée
pour la recherche de zéros d’une fonction. En l’occurrence, il s’agit de trouver (l’unique ?)
zéro de la fonction z 7→ z − un − hf (z).
8
Univ. Paris Cité 2024-2025 M1 MMA, Modélisation Déterministe, Projet
A. B.
C. D.
Figure 1: Méthodes de quadrature : (A) méthode des rectangles à gauche, (B) méthode
des rectangles à droite, (C) méthode des trapèzes, (D) méthode du point milieu.
A. B.
Figure 2: Illustrations d’approximations de la dérivée en un point tn . (A) L’approximation
s’appuie sur la valeur en tn et en tn+1 . (B) L’approximation s’appuie sur la valeur en tn et
en tn−1 . Dans les deux cas la dérivée en tn est illustrée par la droite en noir et pointillés.
9
Univ. Paris Cité 2024-2025 M1 MMA, Modélisation Déterministe, Projet
3.3 Méthodes de Crank-Nicolson et du point milieu
La méthode des trapèzes (Figure 1.C) est définie par
Z b
F (b) + F (a)
F (s)ds ≈ (b − a) .
a 2
Si on applique cette formule à (3.2), on obtient l’approximation
h
un+1 = un + (f (u(tn+1 )) + f (u(tn ))) , n ∈ {0, . . . , N − 1}.
2
Il s’agit de la méthode de Crank-Nicolson, nommée après le physicien anglais Crank et la
mathématicienne anglaise Nicolson. Comme la méthode d’Euler implicite, cette méthode
est implicite car on ne peut pas exprimer un+1 en fonction de un uniquement. Il faut donc
résoudre un problème non-linéaire pour l’utiliser.
Si on utilise la méthode de quadrature dite du point milieu (Figure 1.D), définie par
Z b
b−a
F (s)ds ≈ (b − a)F a + .
a 2
alors on obtient l’approximation
Z tn+1
h
f (u(s))ds ≈ hf u tn + .
tn 2
A ce stade, on se rend compte que cette formule fait apparaître un temps, tn + h/2, qui
n’est pas inclus dans le maillage. On ne peut donc pas a priori utiliser l’approximation de
la solution u en ce temps là. L’astuce (qui est en soi le cœur de la méthode) consiste à
appliquer le schéma d’Euler explicite au terme u(tn + h/2), comme si le temps tn + h/2 était
un nœud du maillage, c’est-à-dire
h h
u tn + = un + f (un )
2 2
pour obtenir le schéma numérique dit du point milieu,
h
un+1 = un + hf un + f (un ) , n ∈ {0, . . . , N − 1}.
2
Ce schéma est explicite, mais contrairement au schéma d’Euler explicite il nécessite l’éva-
luation de la fonction f deux fois dans la même itération. Cet “inconvénient” est compensé
par l’efficacité de la méthode (cf. Section 2).
Ce schéma numérique a été proposé par Runge, mathématicien allemand, à la fin du XIXème
siècle. Son originalité (c’est aussi le cas pour la méthode de Crank-Nicolson) tient à une éva-
luation de la fonction f deux fois durant la même itération, et ceci a mené au développement
de méthodes dites de Runge-Kutta s’appuyant sur ce principe.
10