Introduction à la modélisation

G. Barles

23 novembre 2015

Table des matières

1 Introduction 2

2 Modélisation et équations différentielles ordinaires 6

2.1 Un peu de physique du point matériel . . . . . . . . . . . . . 6
2.2 L’exemple du pendule . . . . . . . . . . . . . . . . . . . . . . 7
2.3 Étude mathématique : problèmes et solutions ? . . . . . . . . 8
2.4 Rappels sur les schémas numériques pour les EDO . . . . . . 9
2.4.1 Étude de l’erreur (I) . . . . . . . . . . . . . . . . . . . 10
2.4.2 Étude de l’erreur (II) . . . . . . . . . . . . . . . . . . . 13
2.4.3 Quelques élément sur les méthodes de Runge-Kutta . 14
2.5 À vous : quelques modèles à traiter . . . . . . . . . . . . . . . 16

3 Modélisation et équations aux dérivées partielles 18

3.1 Évolution de la température le long d’une barre : modélisation
et premières remarques . . . . . . . . . . . . . . . . . . . . . . 18
3.2 Évolution de la température le long d’une barre : étude mathématique 19
3.3 Évolution de la température le long d’une barre : étude numérique 23
3.4 Mise en place des principaux schémas . . . . . . . . . . . . . 23
3.5 Convergence des schémas . . . . . . . . . . . . . . . . . . . . 25
3.6 Aproximation par domaine de calcul borné . . . . . . . . . . 26
3.7 Équation de la chaleur en domaines bornés . . . . . . . . . . 27
3.8 Quelques éléments de modélisation du trafic routier . . . . . . 28

4 Modélisation et optimisation 31

5 Appendice : Compléments 34
5.1 Étude générale des méthodes à un pas . . . . . . . . . . . . . 34
5.1.1 Propriétés importantes d’une méthode à un pas . . . . 34
5.1.2 Condition nécessaire et suffisante de consistance . . . 35
5.1.3 Condition suffisante de stabilité . . . . . . . . . . . . . 36
5.1.4 Ordre d’un schéma . . . . . . . . . . . . . . . . . . . . 36
5.1.5 Exemples . . . . . . . . . . . . . . . . . . . . . . . . . 38

1
1 Introduction
Le but de ce cours est de décrire quelques aspects de la modélisation de
phénomènes essentiellement déterministes mais avant cela, on peut légitimemement
se poser la question :

Qu’est-ce que la modélisation ? Quel est son but ?

Donner une définition précise de la modélisation n’est pas chose aisée

car elle peut apparaı̂tre sous différentes formes suivant les domaines ou les
phénomènes à modéliser. Nous en verrons plus loin quelques exemples assez
différents. Une tentative de définition (à affiner !) pourrait être : l’objectif de
la modélisation, c’est d’écrire un certain nombre d’équations mathématiques
dont les solutions décrivent “assez bien” le phénomène auquel on s’intéresse,
dans le but de prédire l’évolution d’un système à partir d’une situation
donnée ou de mieux comprendre les raisons de cette évolution.

En général, sur une problématique donnée, la première modélisation vise

à comprendre les causes du phénomène : le modèle doit donc être suffi-
samment simple pour pouvoir être analysé complètement et savoir quel(s)
terme(s) dans l’équation provoque tel ou tel effet. Ensuite, une fois convaincu
de la pertinence de cette première modélisation, on peut la compliquer pour
prendre en compte des aspects de plus en plus sophistiqués.

Ceci crée des différences entre les disciplines : dans le cas des sciences
physiques ou des sciences pour l’ingénieur, la modélisation s’appuie sur des
siècles de pratique et d’expérience ainsi que sur l’existence d’un nombre im-
portant de lois physiques bien établies ; on s’intéresse aussi à des phénomènes
clairement déterministes, les mêmes causes induisants les mêmes effets, iné-
luctablement.

À cause de ces caractéristiques, dans de nombreuses situations voire la

plupart, la modélisation est vraiment utilisée dans le but de prédire le com-
portement de tel ou tel système (écoulement autour d’une aile d’avion, tran-
sition de phase, résistance des matériaux, trajectoire de fusée...). Et il est
clair aussi que la modélisation se révèle être un outil précieux sinon indis-
pensable pour les industriels.

Dans le cas d’autres sciences (biologie, économie,...), l’intervention des

êtres vivants perturbe un peu le déterminisme, même si on espère conser-
ver via une compréhension plus fine des phénomènes, un certain aspect
prévisible, au moins statistiquement. En tout cas, l’aspect historique est
complètement différent car le souci de modélisation (en tout cas à un cer-
tain niveau de complexité) est assez récent dans le temps. Dans ces disci-

2
plines, on doit mettre en évidence des “lois de comportement” qui s’appuie
souvent sur des méthodes de modélisation “traditionnelles” mais on observe
un peu plus l’intervention d’aspects aléatoires, déjà présent en mécanique
statistique. L’exemple de la mécanique statistique montre bien comment
l’aléatoire peut conduire à des phénomènes déterministes si on a un effet de
moyennisation.

Ces sciences où l’intervention de la modélisation est plus récentes (ou du

moins avec une approche scientifique plus récente = moins d’un siècle) nous
ramènent loin en arrière, à ce qu’ont certainement été les premiers pas en
physique mais aussi dans des contextes assez différents : par exemple, dans
le cas des sciences humaines, on peut parfois se demander si une quelconque
modélisation est possible tant l’intervention irraisonnée de la psychologie
humaine peut avoir des conséquences tout à fait imprévisibles. Pensons, par
exemple, aux situations de panique à la bourse qui sont souvent dispropor-
tionnées par rapport à la situation économique réelle. Parmi ces sciences plus
récentes où la modélisation a quelquefois été une conséquence de l’apparition
et des progrès de l’informatique, on peut évoquer les sciences économiques,
parmi lesquelles la finance, les sciences d’aides à la décision (théorie des jeux,
apprentissage...), l’informatique et le traitement d’images.

Ce dernier exemple très à la mode recouvre beaucoup de problématiques

différentes mais avec un point commun important : il ne s’agit pas ici de
prédire mais de comprendre comment structurer et/ou traiter une image.
Une image ou son traitement doit ressembler à quel phénomène physique ?
quels sont ses constituants essentiels ? A cause de cela, l’analyse d’images
fait appel à des modèles physiques et mathématiques très sophistiqués aussi
bien pour le débruitage (amélioration de la qualité d’une image) que pour
la simplification d’image (détection de contour...).

Comment modéliser ? Quelles sont les étapes principales ?

Malheureusement (ou heureusement ?), la diversité des situations décrites

ci-dessus empêche l’existence de règles uniformes. Mais il existe un certain
nombre de passages obligés.

La première difficulté, c’est d’identifier clairement les variables et les

paramètres importants du système considéré ; mathématiquement cela si-
gnifie : quelles sont les inconnues pertinentes que l’on va utiliser dans nos
équations et de quelles variables dépendent-elles ? Quels sont les paramètres
additionnels ? En d’autres termes, quels sont les facteurs explicatifs du phénomène
considéré et quelles sont les grandeurs que l’on estime essentielles ?

La deuxième étape consiste à trouver des relations entre ces inconnues,

3
ces variables et ces paramètres (lois de comportement) ainsi que des
équations régissant l’évolution des inconnues.

Cette deuxième étape conduit très souvent à un nombre d’équations

trop important ou à des équations trop complexes. Il est donc nécessaire de
simplifier ces équations en analysant l’importance des différents termes, ce
qui amène à supprimer des inconnues, des variables, des paramètres et des
équations. La question est : peut-on négliger tel ou tel phénomène ?
sait-on a priori que telle ou telle quantité est “petite” et donc on
peut la supposer nulle ? Bien entendu, on peut se tromper et revenir en
arrière. C’est en grande partie cela le travail de modélisation car plus un
modèle est simple, plus il est utilisable et par là même utile.

Ensuite, quand c’est possible, on procède à une étude mathématique

des équations restantes. Evidemment, une situation idéale est celle où ces
équations ont des solutions totalement explicites, mais c’est rare : dans
ce cas, on a un accès direct aux solutions mais aussi (et surtout) à leurs
dépendances en les différentes variables et les divers paramètres. Ceci per-
met des études simples de sensibilités et donne une compréhension plus aisée
des phénomènes.

Malheureusement, dans beaucoup de cas, on ne sait pas dire grand chose

mais cette étude doit au moins analyser le type des équations considérées et
d’en déduire les méthodes numériques adaptées pour “calculer numériquement”
la (ou les) solution(s) avec efficacité. Car on ne sait bien calculer numériquement
que les solutions d’équations pour lesquelles on a un peu de théorie (au moins
qualitative) et encore...

En général, on procède soit avant l’étude mathématique, soit au moment

de l’implémentation numérique à l’adimensionalisation des équations. Les
ordinateurs ne gérant vraiment bien que les quantités entre 1 et 10, on essaie
de faire en sorte que les quantités considérées soient dans cet intervalle et
pour cela, on les normalise par une quantité de référence. Exemple : si on
doit calculer des longueurs l de l’ordre de 104 mètres, on introduira ˜l = 10−4 l
avec des chances pour que ˜l soit de l’ordre de 1. Si on s’intéresse à des fusées
ou des escargots, l’adimensionalisation de la vitesse ne se fera pas sur les
mêmes bases ! De manière générale, si x est une variable ou une fonction à
calculer, on s’intéressera plutôt à :
x
x̃ = ,
x̄
où x̄ est soit une valeur moyenne de x, soit une valeur caractéristique de x.

Ensuite, il faut être prudent quant à l’interprétation des simulations

4
numériques : les méthodes numériques peuvent générer des arte fact de cal-
cul n’ayant rien à voir avec la réalité... cela a conduit pendant quelques
années des patients sur le billard parce que l’image ”montrait” une région
suspecte (avec possibilité de tumeur cancéreuse) . Une remarque essentielle :
en général, on teste les méthodes numériques sur des exemples simples puis
on les utilise dans des cas plus complexes. On ne fait jamais confiance à une
méthode numérique les yeux fermés.

Validation et limites de la modélisation

Une fois la simulation réalisée, on compare (en général) ses résultats

avec l’observation ; cette deuxième étape permet de juger la qualité de la
modélisation et elle conduit, éventuellement, à la modifier : le processus de
modélisation est, le plus souvent, un processus itératif (100 fois sur le métier
remettez votre ouvrage...).

La question des limites des modèles est trop souvent négligée par les
utilisateurs et même par les scientifiques eux-mêmes : quelles ont été les
hypothèses utilisées pour élaborer le modèle ? Sont-elles satisfaites dans la
situation considérée ?

Il faut aussi faire attention au fait qu’un bon ajustement à l’expérience

ne prouve rien, par exemple quand le nombre de mesures est insuffisant. De
plus, lors de la comparaison modèle-expérience, il faut prendre en compte les
erreurs de mesure. Dans le cadre de phénomènes aléatoires, un grand nombre
d’expérience est nécessaire pour échapper aux cas exceptionnels. Dans tous
les cas, il faut être conscient que seule la répétition d’un nombre important
d’expériences aux conclusions concordantes valide de façon pragmatique (et
temporaire) un modèle.

Dans certains domaines d’observation, il n’ y a pas d’expérience di-

recte possible, il faut alors observer non le phénomène lui-même mais ses
conséquences . C’est le cas notamment dans le domaine de la relativité et
de l’astronomie (effet Doppler pour constater l’expansion de l’univers).

D’autres difficultés peuvent apparaı̂tre, à commencer par la sensibilité

aux conditions initiales ; par exemple lorsque l’on veut étudier la trajectoire
d’une feuille de papier. Le fameux effet “aile de papillon” ! Cela se produit
aussi quand le modèle n’assure pas une solution mathématique unique : le
fait de simuler le processus avec une donnée précise ou même avec plusieurs
données initiales ne renseigne absolument pas sur ce qui va se passer avec
une donnée initiale voisine !

Pour d’autres raisons, certains phénomènes ne sont pas non plus prédictibles

5
comme en sciences économiques, où les modèles complexes sont (ou de-
vraient) être utilisés pour mieux comprendre les phénomènes, pour donner
une tendance et pas de manière prédictive ou alors en étant bien conscient
des limites d’un tel modèle.

Enfin, des problèmes purement techniques peuvent compliquer les simu-

lations numériques : puissance de calcul trop faible (météo) ou systèmes
linéaires mal déterminés (ils n’ont de solutions que sous certaines condi-
tions).

Dans ce cours, nous allons examiner trois types de situations en fonctions

du type de problématiques mathématiques à laquelle on aboutit après la
modélisation : une première partie où les modèles obtenus sont des équations
différentielles ordinaires, une seconde où l’on obtient des équations aux
dérivées partielles et une dernière où le thème central sera l’optimisation.

2 Modélisation et équations différentielles ordinaires

2.1 Un peu de physique du point matériel
Dans beaucoup de situations physiques, l’objet d’étude peut être assi-
milé à un point matériel : planète dans le système solaire, bille sur un plan
incliné, l’hypothèse “point matériel” provient soit de la taille “petite” (re-
lativement), soit du fait que seul le comportement du centre de gravité du
corps considéré joue un rôle.

En Mécanique du point matériel, une loi fondamentale est le principe

fondamental de la dynamique que l’on peut résumer par la formule :

F = mγ ,

où F désigne la “résultante” (la somme) des forces extérieures, m la masse

du point matériel et γ son accélération. Si on note x(t) la position du point
matériel à l’instant t, alors γ = ẍ(t) où “ ˙ ” désigne la dérivée par rapport
à t et donc ici “ ¨ ” est une dérivée seconde.
Les forces les plus communes sont la gravitation, les forces électro-statiques
et on dit que F dérive d’un potentiel si :

F = −∇U ,

où U est le potentiel. Ceci nous amène à parler d’énergie car si on est dans
le cas d’une force dérivant d’un potentiel, le principe fondamental de la
dynamique nous dit que :

mẍ(t) = −∇U (x(t)) ,

6
et on voit facilement (en dérivant) que :
1
E(t) = m(ẋ(t))2 + U (x(t)) = constante .
2
La quantité constante (“conservé”) E(t) est l’énergie du système, le premier
terme désignant l’énergie cinétique et le second l’énergie potentielle.

2.2 L’exemple du pendule

Pour mettre en application la physique du point matériel, on s’intéresse
au pendule simple qui est une masse ponctuelle fixée à l’extrémité d’un fil
sans masse, inextensible et sans raideur, lui même ayant son autre extrémité
fixée en un point O. Il oscille sous l’effet de la pesanteur. Si sa masse est
normalisée à m = 1, l’équation s’écrit :
→
−
ẍ(t) = →
−
g +R ,
→
−
où →
−
g est la force de pesanteur, R est la tension du fil (une force qui exprime
que la masse ponctuelle est attachée au fil et donc ne tombe pas) et x(t) est
la position de la masse ponctuelle à l’instant t.

On suppose que le modèle est écrit dans R2 et il s’agit d’écrire les

équations de manière plus utilisable. Pour cela, on prend O comme origine
et on utilise des coordonnées polaires : si r est la longueur du fil alors :

x(t) = r exp(iθ(t)) .

On peut faire les calculs soit en repère mobile soit via les nombres complexes :
dans ce dernier cas, on a :

ẋ(t) = iθ̇(t)r exp(iθ(t)) et ẍ(t) = iθ̈(t)r exp(iθ(t)) − (θ̇(t))2 r exp(iθ(t)) .

Le premier terme dans le calcul de ẍ(t) est l’accellération tangentielle (dans

la direction orthogonale au fil i exp(iθ(t))) et le second terme est la force
centrifuge (dans la direction radiale exp(iθ(t)))).
Comme il n’y a claiment pas de mouvement dans la direction radiale (le
poids, la tension du fil et la force centrifuge s’équilibrant parfaitement), il
reste la direction tangentielle où seule la résultante de la gravité intervient :
si l’axe θ = 0 est la verticale dirigée vers le centre de la terre, on a :
→
−
g = g(∈ R) = g cos(θ(t)) exp(iθ(t)) − g sin(θ(t))i exp(iθ(t)) (1) ,

et donc :
θ̈(t)r = −g sin(θ(t)) .
(1). Il est à noter ici que l’on écrit tout dans le repère mobile formé des vecteurs exp(iθ(t))
et i exp(iθ(t)).

7
 En normalisant les constantes à 1 (r = g = 1), on a donc une équation
différentielle ordinaire simple :

θ̈(t) = − sin(θ(t)) ,

qui nécessite la connaissance de (par exemple) θ(0), θ̇(0) pour être résolue.
Quelques remarques :
1. Il s’agit d’une équation non linéaire qui n’est pas évidente à résoudre...
2. En multipliant par θ̇(t) et en intégrant, on voit que :
1
(θ̇(t))2 − cos(θ(t)) = constante .
2
On a affaire à un système “conservatif” où l’énergie est préservée, le
premier terme étant, comme ci-dessus, un terme d’énergie cinétique
et le second un terme d’énergie potentielle.
3. Nous avons fait beaucoup d’hypothèses simplificatrices pour arriver
à ce modèle élémentaire : point matériel, absence de frottements (en
pratique, le pendule s’arrête... ici, non), la terre est plate,...etc.

2.3 Étude mathématique : problèmes et solutions ?

On va dans cette section regarder plusieurs propriétés des solutions de
cette équation différentielle.
Théorème 1. (Existence, unicité, stabilité (I)) Pour toute donnée ini-
tiale θ(0) = a, θ̇(0) = b, l’équation du pendule a une unique solution définie
pour tous temps. De plus cette solution dépend continûment de a et b.
Démonstration. Un bon cours d’edo et aucune difficulté...

Une attention est toujours portée sur les points d’équilibre, c’est-à-dire
ceux pour lesquels la solution ne bouge pas. Pour le pendule, cela semble
assez évident mais en résolvant θ̇(t) = θ̈(t) = 0 pour tout t, on trouve deux
équilibre : (0, 0) et (π, 0) (modulo 2π).

L’étude de la stabilité de ces points d’équilibre est une première occasion

de faire intervenir une méthode très utilisée : la linéarisation. Si on connait
une trajectoire θ(·) (par exemple explicitement), on cherche des solutions
proches en considérant θ(·) + h(·) où h(·) est supposée “petite” (par exemple
pour la norme du sup ou la norme C 1 ). Ceci conduit à écrire :

θ̈(t) + ḧ(t) = − sin(θ(t) + h(t)) = − sin(θ(t)) − cos(θ(t))h(t) + · · · .

Si on utilise l’équation satisfaite par θ et si on néglige les · · · en ne conservant

que le terme “linéaire”, on aboutit à :

ḧ(t) = − cos(θ(t))h(t) .

8
Si on utilise la solution d’équilibre θ(t) ≡ 0, on est amené à considérer l’edo
ḧ(t) = −h(t) dont les solutions sont des sin et cos et on conclut à la stabilité.
Par contre, pour la solution d’équilibre θ(t) ≡ π, on est amené à considérer
l’edo ḧ(t) = h(t) dont les solutions sont exp(t) et exp(−t), la première étant
non bornée pour les t > 0.
Un modèle plus général prend en compte un terme de frottement qui
s’oppose au mouvement, par exemple en −αẋ(t). On vérifie facilement que
la nouvelle équation est de la forme :

θ̈(t) = −αθ̇(t) − sin(θ(t)) ,

puis on prouve de la même manière l’existence, l’unicité et la stabilité des

solutions : on a même des résultats de stabilité meilleurs, ce qui ne doit pas
surprendre du point de vue physique ?

Il reste une question à aborder : si on ne connait pas les solutions expli-

citement, on fait comment ? On construit des schémas numériques.
Pour ce faire, on écrit l’équation sous la forme :

ẏ(t) = f (t, y(t)) ,

où, ici, y(t) = (θ(t), θ̇(t)) et f (t, (y1 , y2 )) := (y2 , − sin(y1 )).
Pour simplifier la présentation, on utilisera l’hypothèse simplifiée sui-
vante (f est “globalement lipschitzienne” en y) :

(GL) La fonction f : [0, T ] × Rn → Rn est continue et il existe une

constante L telle que :

|f (t, y1 ) − f (t, y2 )| ≤ L|y1 − y2 | ,

pour tous t ∈ [0, T ] et y1 , y2 ∈ Rn .

2.4 Rappels sur les schémas numériques pour les EDO

On va d’abord présenter la méthode d’Euler.

Pour résoudre numériquement l’EDO, on va se donner une grille, c’est-

à-dire des points t0 = 0 < t1 < t2 < t3 < · · · < tN = T , et on va essayer de
calculer une “bonne” approximation des valeurs de la solution en tous ces
points, c’est-à-dire des valeurs yi ' y(ti ).
Numériquement le sens de “bonne approximation” n’est pas absolu car,
d’une part, l’ordre de grandeur joue un rôle (une approximation à 100 000
ans près peut être excellente si on est géologue...) et, d’autre part, le temps
mis pour calculer la solution peut être un facteur important (quel est l’intérêt
d’avoir un résultat très précis s’il faut 10 ans pour l’obtenir ?). Il y a toujours,

9
dans les méthodes numériques un “rapport qualité - prix ” : précision vs
temps de calcul ou complexité.
Essentiellement il y a deux approches pour calculer les yi qui, dans le
cas de la méthode d’Euler, vont aboutir au même résultat : soit on approche
directement l’EDO par “différences finies” en utilisant une approximation
de la dérivée y 0 (t), soit on intègre l’EDO, se rapprochant ainsi de la preuve
d’existence.
L’approche par “différences finies” consiste à approcher y 0 (ti ) par différences
finies ; par exemple :
y(ti+1 ) − y(ti )
y 0 (ti ) ' .
ti+1 − ti
L’EDO se réécrit alors sous la forme :
y(ti+1 ) − y(ti )
' f (ti , y(ti )) ,
ti+1 − ti
et donc :
y(ti+1 ) ' y(ti ) + (ti+1 − ti )f (ti , y(ti )) .
Ceci suggère que l’on peut calculer les yk via la relation de récurrence :

yi+1 = yi + (ti+1 − ti )f (ti , yi ) ,

le terme y0 étant connu (donnée initiale). C’est la méthode d’Euler.

La deuxième approche, qui va nous conduire au même résultat mais avec
une philosophie très différente, consiste à intégrer l’équation de ti à ti+1 :
Z ti+1
y(ti+1 ) = y(ti ) + f (s, y(s))ds .
ti

Si on applique la méthode des rectangles à l’intégrale de la manière suivante :

Z ti+1
f (s, y(s))ds ' (ti+1 − ti )f (ti , y(ti )) ,
ti

on retrouve les calculs ci-dessus et la méthode d’Euler.

Remarque : En pratique, il peut être intéressant d’utiliser une subdivision
adaptée avec plus de points aux endroits où la fonction y varie beaucoup et
moins de points aux endroits où les variations sont faibles. Mais choisir une
telle subdivision (ou faire en sorte que l’ordinateur choisisse automatique-
ment cette subdivision = schémas adaptatifs) n’est pas toujours simple.

2.4.1 Étude de l’erreur (I)

Désormais nous nous plaçons dans le cadre d’une grille uniforme :
T
ti+1 − ti = h = .
N

10
La méthode d’Euler s’écrit alors :

yi+1 = yi + hf (ti , yi ) .

On notera :
ei = y(ti ) − yi ,
l’erreur commise au point ti et :

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

l’erreur de “consistance” ; c’est l’erreur systématique commise sur yi+1 :

même si on a calculé exactement la valeur à l’instant ti (yi ' y(ti )), on a
une erreur sur y(ti+1 ) qui est εi .
Pour évaluer l’erreur, on procède comme suit :

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

= [y(ti+1 ) − y(ti ) − hf (ti , y(ti ))] + [y(ti ) − yi ] + [yi + hf (ti , yi )] +
h[f (ti , y(ti )) − hf (ti , yi )] − yi+1
= εi+1 + ei + h[f (ti , y(ti )) − hf (ti , yi )] .

D’où :

|ei+1 | ≤ |εi+1 | + |ei | + h|f (ti , y(ti )) − hf (ti , yi )|

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

en utilisant le caractère lipschitz de f .

Il reste à estimer |εi+1 |. En utilisant l’équation, on a :

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

Z ti+1 Z ti+1
= | y 0 (s)ds − y 0 (ti )ds|
ti ti
Z ti+1
= | [y 0 (s) − y 0 (ti )]ds|
ti
≤ h sup |y 0 (s) − y 0 (ti )|
s∈[ti ,ti+1 ]
0
≤ h.ω(h, y ) ,

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

Finalement on a l’estimation suivante de l’erreur ei+1 :

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

Pour conclure on utilise le :

11
Lemme 1. (Lemme de Gronwall discret)
Si (θi )i est une suite de réels positifs qui satisfait :

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

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

exp(iA) − 1
θi ≤ exp(iA)θ0 + B.
A

On utilise d’abord le lemme avec A = Lh et B = h.ω(h, y 0 ) :

exp(Lih) − 1
|ei | ≤ exp(Lih)|e0 | + h.ω(h, y 0 ) .
Lh
Mais ih = ti et (a priori) e0 = 0 [ou du moins, e0 est très petit] donc :

exp(Lti ) − 1 exp(LT ) − 1
|ei | ≤ ω(h, y 0 ) ≤ ω(h, y 0 ) .
L L
On vient donc de prouver le :

Théorème 2. Sous l’hypothèse (GL) :

exp(LT ) − 1
max |ei | ≤ ω(h, y 0 ) .
0≤i≤N L

En particulier, max |ei | → 0 quand N → +∞.

0≤i≤N

On donne maintenant la :
Preuve du Lemme de Gronwall discret : On note (ui )i la suite définie
par :
θi
ui = .
(1 + A)i
En divisant la propriété satisfait par la suite (θi )i par (1 + A)i+1 , on voit
que :
B
ui+1 ≤ ui + .
(1 + A)i+1
Une récurrence immédiate montre que :
i−1
X B B 1 − ai
ui ≤ u0 + k+1
= u0 + ,
(1 + A) 1+A 1−a
k=0

où a = 1/(1 + A). Comme :

1 1+A
= ,
1−a A

12
il en résulte :
1 − ai
ui ≤ u0 + B ,
A
et donc :
(1 + A)i − 1
θi ≤ (1 + A)i u0 + B .
A
Il reste à remarquer que 1 + A ≤ exp(A), ce qui est clair puisque :
A2 An
exp(A) = 1 + A + + ··· + + ··· .
2! n!


À titre d’exercice, on pourra démontrer la :

Proposition 1. Si yh est la fonction affine par morceaux telle que yh (ti ) =
yi pour tout i, on a :

||yh − y||∞ → 0 quand h → 0 .

2.4.2 Étude de l’erreur (II)

Le section précédente donne une estimation de convergence qui dépend
du module de continuité de y 0 . Mais la fonction y est inconnue donc ce
résultat n’est pas satisfaisant car il n’est pas explicite. Nous allons mainte-
nant donner une autre estimation qui ne dépend que des données, c’est-à-dire
de f .
Pour cela, on étudie le module de continuité de y 0 :

|y 0 (t) − y 0 (s)| = |f (t, y(t)) − f (s, y(s))|

= |f (t, y(t)) − f (s, y(t)) + f (s, y(t)) − f (s, y(s))|
≤ |f (t, y(t)) − f (s, y(t))| + |f (s, y(t)) − f (s, y(s))|
≤ |f (t, y(t)) − f (s, y(t))| + L|y(t) − y(s)|

D’après les résultats sur les équations différentielles, |y(t)| ≤ D pour une
certaine constante D sur l’intervalle [0, T ] et le premier terme est estimé par
le module de continuité de f sur [0, T ] × B(0, D), noté ωD (·, f ).
Quant au second, par le Théorème des Accroissements Finis :

|y(t) − y(s)| ≤ Mf |t − s| ,

où :
Mf = max |f (t, y)| .
[0,T ]×B(0,D)

Finalement :
ω(h, y 0 ) ≤ ωD (h, f ) + LMf h ,
ce qui donne le résultat suivant qui était notre objectif :

13
Théorème 3. Sous l’hypothèse (GL) :
exp(LT ) − 1
max |ei | ≤ (ωD (h, f ) + LMf h) .
0≤i≤N L
On renvoit à l’appendice pour l’étude générale des méthodes à un pas
qui s’écrivent 
yi+1 = yi + hΦ(ti , yi , h)
y0 = y0,h
où Φ est une fonction continue sur [0, T ] × Rn × [0, H], H désignant un pas
de discrétisation maximal.

2.4.3 Quelques élément sur les méthodes de Runge-Kutta

Ces méthodes sont les plus utilisées : elles sont rentrées “en standard”
dans la plupart des logiciels. Comment marchent-elles ?
On repart de l’idée fondamentale qui consiste à écrire :
Z ti+1
y(ti+1 ) = y(ti ) + f (s, y(s))ds .
ti

Nous avons vu dans la deuxième partie que pour calculer une intégrale, on
utilise une formule de quadrature :
Z 1 q
X
ψ(s)ds ' bj ψ(cj ) ,
0 j=0

où c0 < c1 < · · · < cq , ce qui, en prenant ψ(s) = f (ti + sh, y(ti + sh)), nous
donne :
q
X
y(ti+1 ) ' y(ti ) + h bj f (ti,j , y(ti,j )) ,
j=0

où ti,j = ti + cj h. Ceci suggère une méthode que l’on peut écrire :
q
X
yi+1 = yi + h bj ki,j ,
j=0

où ki,j est une approximation de f (ti,j , y(ti,j ))

Le problème, c’est qu’il faut encore calculer des approximations yi,j des
y(ti,j ) pour avoir celle de ki,j et ceci est fait via :
q
X
yi,j = yi + h aj,k f (ti,k , y(ti,k )) .
k=0

Cette procédure a l’air d’induire des équations non-linéaires couplées diffi-

ciles à résoudre et pour que ce ne soit pas le cas, on suppose que les yi,j ne

14
dépendent que des points déjà calculés, c’est-à-dire des yi,k pour k < j. On
a donc :
j−1
X
yi,j = yi + h aj,k f (ti,k , y(ti,k )) ,
k=0

et les yi,j , ainsi que les ki,j , sont calculés de proche en proche.
On résume souvent une méthode de Runge-Kutta grâce à un tableau de
la forme :

c1 a1,1 · · · a1,q
.. .. .
. . · · · ..
cq aq,1 · · · aq,q
b1 · · · bq

Exemples :
• q = 1 : C’est la méthode d’Euler basée sur la méthode des rectangles.
• q = 2 C’est l’exemple de la section précédente, basée sur la méthode des
trapèzes, avec le tableau :

0 0 0
β β 0
1 1
1 − 2β 2β

NB : β = (2α)−1 .

• q = 4 : la méthode de Runge-Kutta “classique” (la plus utilisée) basée

sur la formule de Simpson :

0 0 0 0 0
1/2 1/2 0 0 0
1/2 0 1/2 0 0
1 0 0 1 0
1/6 2/6 2/6 1/6
La méthode s’écrit aussi :
1
Φ(t, y, h) = [k1 + 2k2 + 2k3 + k4 ]
6

15
avec :

k1 = f (t, y)
h h
k2 = f (t + , y + k1 )
2 2
h h
k3 = f (t + , y + k2 )
2 2
k4 = f (t + h, y + hk3 )

La méthode est d’ordre 4 mais il vaut mieux avoir Maple pour le vérifier !

2.5 À vous : quelques modèles à traiter

Exercice 1. Décrire le mouvement d’une bille roulant sans frottement sur
un plan incliné.

Exercice 2. Mouvements des planètes, équations de Képler

On considère une planète assimilée à un point matériel qui tourne autour
du soleil (dans R3 ). Si x(t) désigne sa position à l’instant t, le soleil étant à
l’origine, on a, d’après la loi de la gravitation universelle :
c
ẍ(t) = − x(t) ,
|x(t)|3

où c est une constante qui prend en compte les différentes masses et la
constante de gravitation g.
(i) Montrer que, si ∧ est le produit vectoriel dans R3 , x(t)∧ ẋ(t) est constant.
En déduire que le mouvement est plan et suit la “loi des aires” (l’aire ba-
layée par le vecteur x(t) pendant l’intervalle de temps [t1 , t2 ] ne dépend que
de t2 − t1 .
(ii) Comme le mouvement est plan, on peut écrire x(t) = r(t) exp(iθ(t)).
Donner les équations satisfaites par r(t) et θ(t) et retrouver la loi des aires.
(iii) On suppose que θ est une fonction strictement croissante de t [on pourra
se demander pourquoi on peut supposer cela]. Si ψ(θ) = t (au moins locale-
ment) et on s’intéresse à u(θ) = [r(ψ(θ))]−1 . Prouver que :

u00 (θ) = −u(θ) + l ,

où l est une constante et en déduire que :

p
r(ψ(θ)) = ,
1 + e cos(θ)

et que la planète décrit une ellipse.

(iv) La loi de la gravitation universelle dérive-t-elle d’un potentiel ? Si oui,
écrire l’énergie conservée et dire quand la vitesse de la planète est minimale
et maximale sur sa trajectoire.

16
Exercice 3. Équation logistique :
Quand on s’intéresse à l’évolution d’une population de taille N (t) à l’instant
t, on doit prendre en compte un taux de natalité n et de mortalité d et
l’équation la plus simple que l’on peut écrire est :

Ṅ (t) = nN (t) − dN (t) = aN (t) ,

où a = n − d. Comme échauffement, on pourra calculer les solutions de

cette équation pour a constant ou dépendant du temps t et interpréter les
résultats.
Mais ces modèles sont trop naı̈fs et les biologistes préfèrent des lois “logisti-
ques” :
Ṅ (t) = γN (t)(N ∗ − N (t)) ,
où γ et N ∗ sont des constantes positives. L’interprétation de cette équation,
c’est que le milieu où vit cette population ne peut nourrir qu’un nombre limité
d’individu (essentiellement N ∗ ) et, arrivé à ce seuil, la population ne peut
plus croı̂tre. Résoudre explicitement cette edo et interpréter les résultats.
Même question pour la dynamique de type Michaelis-Menten :
aN (t)
Ṅ (t) = ,
b + N (t)
où a, b > 0, équation utilisée en chimie et en biologie.

Exercice 4. Epidémiologie, évolution de population

La biologie fournit un nombre de modèles très variés : un des plus connus
est le modèle proies-prédateurs qui conduit au système différentiel de Volterra-
Lotka

dL(t)
= aL − bLR ,
dt
dR(t)
= −pR + qLR .
dt
où L(t) est le nombre de proies à l’instant t, R(t) le nombre de prédateurs
et a, b, p, q sont des constantes positives.
Les différents termes s’interprètent comme suit :
– AL est la croissance du nombre des proies grâce aux naissances.
– −bLR est sa diminution due aux repas des prédateurs !
– −pR est la diminution du nombre de prédateurs en l’absence de proies.
– Dès qu’il y a présence de proies, le nombre de prédateurs augmente et
ceci d’autant plus vite qu’il y a de proies d’où le terme qLR.
Il s’agit là encore d’un cas modèle car l’étude mathématique de ce système
d’équations n’est pas trop difficile. De plus, on simule facilement (en utili-
sant par exemple la méthode de Runge-Kutta) les évolutions des populations
régies par ces équations à partir d’un quelconque état initial.

17
 Étudier ce système en montrant qu’il existe une intégrale première de la
forme :
c1 L(t) + c2 R(t) + c3 log(L(t)) + c4 log(R(t)) ,
en caculant le ou les points d’équilibre et en discutant sa (ou leur) stabilité.
Programmer un schéma d’Euler et un schémas RK4 : que constate-t-on ?

3 Modélisation et équations aux dérivées partielles

3.1 Évolution de la température le long d’une barre : modélisation
et premières remarques
On considère une barre infiniment longue et on note x la coordonnée
le long de la barre. A l’instant t = 0, on observe une température T0 (x)
au point x. Cette chaleur va se diffuser le long de la barre et on voudrait
connaı̂tre son évolution au cours du temps.
On note T (x, t) la température au point x à l’intant t. La variation de
la quantité de chaleur dans un intervalle [x − ∆x, x + ∆x] est égale au flux
de chaleur à travers les bord (les points x − ∆x, x + ∆x) et la loi de Fourier
indique que ce flux de chaleur est proportionnel au gradient de T . D’où :
Z x+∆x 
d ∂T ∂T
T (y, t)dy = −λ (x − ∆x, t) + λ (x + ∆x, t) ,
dt x−∆x ∂x ∂x

où λ > 0 est la conductivité thermique du milieu. Le changement de signe

s’explique par le changement de signe du flux. On en déduit :
x+∆x x+∆x
∂2T
Z Z
∂T
dy = λ (y, t)dt ,
x−∆x ∂t x−∆x ∂x2

Comme ∆x est arbitraire, la fonction T satisfait l’équation de la chaleur :

∂T ∂2T
(1) −λ 2 =0 dans R × (0, tmax ) .
∂t ∂x
(on prendra désormais λ = 1/2 pour simplifier).
La première remarque c’est que la modélisation s’appuie effectivement
ici sur des lois de la physique qui existent et sont bien établies.
L’équation ci-dessus revient à peu près à dire que si on connait la température
T à l’instant t, on la calcule à l’instant t + ∆t en moyennant la température
autour de x de la manière suivante :
1h √ √ i
(2) T (x, t + ∆t) ' T (x + ∆t, t) + T (x − ∆t, t) .
2
Ce cas est un exemple modèle car :

18
 1. On connait explicitement la solution :
Z +∞
|x − y|2
 
1
(3) T (x, t) = √ T0 (y) exp − dy .
2π −∞ 2t

Et, grâce à cela, on a une étude mathématique complète. Rarissime !

2. Du point de vue physique, ce modèle a un problème majeur : les
informations se propagent à vitesse infinie ce qui n’est pas physique-
ment acceptable. En effet, si T0 > 0 sur un petit intervalle centré
en 0 et nul ailleurs, alors T (x, t) > 0 pour tout temps t > 0 et
pour tout x. On n’a pas ici de propriété de vitesse finie de propaga-
tion : c’est une différence fondamentale entre équations paraboliques
et équations hyperboliques que nous verrons par la suite.
3. Pourtant, la confrontation avec l’expérience est très satisfaisante : ce
qui prouve que la modélisation même très simplifiée peut donner de
bons résultats en dépit de problèmes apparents. Et justement l’avan-
tage de ce modèle, c’est sa simplicité.
4. Du point de vue simulation, on ne calcule JAMAIS les valeurs de
la fonction T avec la formule explicite (3) (sauf évidemment quand
l’intégrale se calcule explicitement...) mais plutôt avec des schémas
numériques de type (2), en général plus efficaces que celui-là.

3.2 Évolution de la température le long d’une barre : étude

mathématique
On commence par le :
Théorème 4. Si T0 est une fonction continue qui satisfait la propriété :
2
(4) lim T0 (x)ε−A|x| = 0
|x|→+∞

pour une certaine constante A > 0, alors la fonction T définie par (3) est
une solution de (1) avec la donnée initiale :

(5) T (x, 0) = T0 (x) dans R ,

et cette solution est définie sur un intervalle de temps [0, tmax [ avec tmax =
1/(2A). De plus, cette solution est unique dans la classe des solutions qui
croissent au plus en exp(Cx2 ) à l’infini.
Remarque : si t ∈]0, tmax [ avec tmax = 1/(2A), l’intégrale écrite dans la
proposition a bien un sens, d’où le choix de tmax qui est optimal.
Preuve : On va calculer la solution fondamentale de l’équation (1) i.e. la
solution ρ de (1) associée à la donnée initiale suivante :

(6) ρ(x, 0) = δ0 (x) dans R ,

19
où δ0 est la masse de Dirac centrée en x = 0. On notera ρ(x, t) cette solution.
La solution T sera alors donnée à partir de ρ via la formule :
 
T (x, t) = ρ(·, t) ∗ T0 (x)

où la convolution agit dans la variable d’espace seulement. En effet, d’après

les propriétés de la convolution, si ρ est une fonction régulière,
∂T 1 ∂2T  ∂ρ 1 ∂ 2 ρ 
− = − ∗ T0 = 0
∂t 2 ∂x2 ∂t 2 ∂x2
et :  
ρ(·, t) ∗ T0 (x) → T0 (x) quand t → 0

(ces propriétés deviendront plus claires après le calcul explicite de ρ).

Maintenant, pour calculer ρ, on raisonne formellement en introduisant
la transformée de Fourier en x :
Z
ρ̂(ξ, t) = ρ(x, t)e−iξx dx
R

et on applique cette transformée de Fourier à l’équation. Il est bien connu

que :

∂ρ
c ∂ ρ̂
(ξ, t) = (ξ, t)
∂t ∂t
2ρ
∂d
(ξ, t) = −ξ 2 ρ̂(ξ, t)
∂x2
d’où :  ∂ ρ̂ 
+ ξ 2 ρ̂ (ξ, t) = 0 dans R × (0, tmax )
∂t
avec :
ρ̂(ξ, 0) = 1 dans R .
Il en résulte de l’intégration de cette équation différentielle en t que :
 tξ 2 
ρ̂(ξ, t) = exp − dans R × (0, tmax ) .
2
x2

r
2 π
Or la transformée de Fourier de exp(−aξ ) est : exp − , de sorte
a 4a
que l’on obtient une formule explicite pour ρ :
1  x2 
ρ(x, t) = √ exp − ,
2πt 2t
puisqu’en effet la transformée de Fourier inverse est :
Z
1
ρ(x, t) = ρ̂(ξ, t)εiξx dξ .
2π R

20
On peut aussi vérifier directement et de manière élémentaire (mais cela de-
mande un peu de travail autour du Théorème de Lebesgue de dérivation sous
le signe somme) que la fonction T donnée dans l’énoncé de la proposition
est bien solution de l’équation (exo !).

Pour l’unicité, on va prouver que l’on a une propriété un peu plus forte :
le Principe du Maximum.
Lemme 2. Si T1 et T2 sont des fonctions régulières (C 2 en x et C 1 en t)
sur R × [0, tmax [ qui satisfont la condition de croissance :

|T1 (x, t)|, |T2 (x, t)| ≤ C exp(Cx2 ) dans R × (0, tmax ) ,

pour une certainne constante C > 0 et :

∂T1 1 ∂ 2 T1
− ≤0 dans R × (0, tmax ) (T1 est sous-solution de l’équation),
∂t 2 ∂x2
∂T2 1 ∂ 2 T2
− ≥0 dans R × (0, tmax ) (T2 est sursolution de l’équation),
∂t 2 ∂x2
T1 (x, 0) ≤ T2 (x, 0) dans R ,
alors :
T1 (x, t) ≤ T2 (x, t) dans R × (0, tmax ) .
Évidemment ce résultat donne l’unicicité car si T1 et T2 sont deux solu-
tions qui satisfont la condition de croissance alors on peut faire jouer à T1 le
rôle de sous-solution et à T2 le rôle de sursolution et on conclut que T1 ≤ T2
dans R × (0, tmax ). Mais, en inversant les rôles, on a l’inégalité opposée donc
l’égalité.

Pour prouver le lemme, on note T = T1 − T2 qui satisfait :

|T (x, t)| ≤ 2C exp(Cx2 ) dans R × (0, tmax ) ,

pour une certainne constante C > 0 et :

∂T 1 ∂2T
− ≤0 dans R × (0, tmax ) ,
∂t 2 ∂x2
T (x, 0) ≤ 0 dans R .
On considère maintenant :

χ(x, t) := exp((K1 t + K2 )Cx2 ) ,

où K1 , K2 > 1 sont des constantes à fixer. Des calculs élémentaires montrent
que, si on choisit K2 > 1 quelconque et K1 assez grand alors :
∂χ 1 ∂ 2 χ
− >0 dans R × (0, τ ] ,
∂t 2 ∂x2

21
pour τ assez petit (typiquement K1 τ = K2 ).

On raisonne alors de la manière suivante : pour 0 < α  1, on considère :

sup (T (x, t) − αχ(x, t)) .

R×[0,τ ]

Ce “sup” est en fait un “max” car T (x, t)−αχ(x, t) → −∞ quand |x| → +∞.
Il existe donc (x0 , t0 ) ∈ R × [0, τ ] qui est un point de maximum de T − χ
sur R × [0, τ ].
Se présente alors deux cas :

– t0 = 0 alors le maximum est négatifs puisque T (x, 0) ≤ 0 pour tout x

et χ ≥ 0. On a donc, si tel est le cas, T (x, t) − αχ(x, t) ≤ 0 dans R × [0, τ ].

– t0 > 0 : comme on est en un point de maximum,

∂T ∂χ ∂2T ∂2χ
(x0 , t0 ) = α (x0 , t0 ) et (x 0 , t0 ) ≤ α (x0 , t0 ) ,
∂t ∂t ∂x2 ∂x2
mais on aurait alors :
∂T 1 ∂2T ∂χ ∂2χ
0≥ (x0 , t0 ) − (x 0 , t0 ) ≥ α( (x0 , t0 ) − (x0 , t0 )) > 0 ,
∂t 2 ∂x2 ∂t ∂x2
une contradiction.
On est donc toujours dans le premier cas, T (x, t) − αχ(x, t) ≤ 0 dans
R × [0, τ ] pour tout α > 0 et on conclut en faisant tendre α vers 0.
NB : une petite erreur s’est glissée dans le raisonnement ci-dessus : la trou-
ver ! 

Exercice 5. L’équation de Black & Scholes en Finance :

Elle s’écrit :
∂U ∂U 1 ∂2U
− + rU − rP − σ2P 2 =0 pour P > 0, t ≤ T ,
∂t ∂P 2 ∂P 2
où U est le prix d’une option d’achat qui rapporte (P − K)+ si on attend
la maturité T [on a donc U (P, T ) = (P − K)+ ], P désigne le prix d’une
action, σ sa volatilité (taux de fluctuation aléatoire de l’action), r est le
taux court (taux de caisse d’épargne). Le but de l’exercice est de calculer U
“explicitement” en se ramenant à (3) via divers changement de variables.
On commencera par poser P = exp(x) et :

u(x, t) = U (P, T − t) ,

puis :
v(x, t) = exp(Kt)u(x − ct, t) .

22
3.3 Évolution de la température le long d’une barre : étude
numérique
On va s’intéresser à la résolution numérique de :
∂u 1 ∂ 2 u
(7) − =0 dans R × (0, T ) .
∂t 2 ∂x2

(8) u(x, 0) = u0 (x) dans R .

Dans un premier temps, on ne va pas se soucier d’une difficulté essentielle

de la résolution numérique de (7)-(8) : le fait que l’équation soit posée dans
un domaine non borné...
Pour calculer numériquement la solution de (7)-(8), on introduit une
grille en x et en t, en posant :

xj = j∆x, j ∈ Z,
tn = n∆t, n ∈ {0, 1, · · · , N }, N ∆t = T,

et on notera unj une approximation de u(xj , tn ). Il est à remarquer cette fois

que j ∈ Z, ce qui n’est évidemment pas très compatible avec un calcul sur
ordinateur ! Pour simplifier la présentation, on va donc d’abord présenter
des schémas numériques définis pour j ∈ Z puis on expliquera comment se
ramener à un nombre fini de valeurs de j.

3.4 Mise en place des principaux schémas

Il y a trois schémas numériques “classiques” pour calculer la solution
de (7)-(8) qui sont les suivants : le schéma explicite standard (SES), le
schéma implicite standard (SIS) et les θ-schémas de Crank-Nickolson ; nous
n’étudierons ici en détail que les deux premiers.

Schéma Explicite Standard :

un+1
j − unj 1 unj+1 + unj−1 − 2unj
(SES) − = 0 , n ∈ N, j ∈ Z .
∆t 2 (∆x)2
Comme dans le cas des équations hyperboliques (déjà vu ?), on étudie la
stabilité par une analyse de Fourier ; si unj = exp(ikj∆x), alors :

∆t
un+1
j = g(λ, k∆x)unj , avec λ= .
2(∆x)2
Le calcul de g donne :


g(λ, k∆x) = 1 + 2λ cos(k∆x) − 1
 k∆x 
= 1 − 4λ sin2 .
2

23
Tout d’abord, comme λ ≥ 0, on a bien g(λ, k∆x) ≤ 1. Pour avoir maintenant
g(λ, k∆x) ≥ −1, il est nécessaire d’avoir λ ≤ 1/2, i.e.

∆t
≤ 1,
(∆x)2

c’est la condition de CFL. Il est à noter que cette condition est équivalente à
la monotonie du schéma ; en effet, on peut le réécrire sous la forme suivante :

un+1
j = λunj+1 + λunj−1 + (1 − 2λ)unj

de telle sorte que, si λ ≤ 1/2, on a la propriété de monotonie :

unj ≥ 0 pour tout j ∈ Z =⇒ un+1

j ≥ 0 pour tout j ∈ Z ,

ou encore :

unj ≥ vjn pour tout j ∈ Z =⇒ un+1

j ≥ vjn+1 pour tout j ∈ Z .

On retrouve ainsi une version “discrète” de la propriété de principe de maxi-

mum (et de comparaison) satisfaite par l’équation.
Le schéma (SES) est consistant, un calcul simple montre qu’il est d’ordre
2 en x et d’ordre 1 en t.

Remarque : La condition de CFL est désagréable car elle implique que

l’on doit prendre ∆t de l’ordre de (∆x)2 , c’est-à-dire extrêmement petit. On
“avance” ainsi très peu vite en temps et pour calculer u(x, t) il faut beaucoup
d’itérations. En revanche, à chaque pas d’itération, le schéma calcule vite la
solution puisqu’il est explicite (2) .

Schéma Implicite Standard :

un+1
j − unj n+1 n+1 n+1
1 uj+1 + uj−1 − 2uj
(SIS) − = 0 , n ∈ N, j ∈ Z .
∆t 2 (∆x)2

Pour étudier la stabilité de ce schéma, on peut adopter deux stratégies : une

par Fourier, l’autre par principe de maximum.
Par Fourier, on a :
h
i−1
g(λ, k∆x) = 1 − 2λ cos(k∆x) − 1 ≤ 1,

car cos(k∆x) − 1 ≤ 0, donc le schéma implicite est inconditionnellement

stable.
(2). ceci est à comparer avec le schéma implicite standard pour lequel on est pas obligé de
prendre ∆t petit, mais en revanche il faut résoudre un système linéaire à chaque itération
- le choix d’un “bon” schéma n’est donc pas évident a priori.

24
 Cette propriété peut être obtenue par un argument heuristique de type
principe de maximum. Supposons (unj )j∈Z borné et montrons que :

max |un+1
j | ≤ max |unj | .
j∈Z j∈Z

On ne va démontrer que l’inégalité :

max un+1
j ≤ max unj ,
j∈Z j∈Z

celle avec les valeurs absolues s’obtenant ou bien en raisonnant de la même

façon avec le min, ou bien en changeant les un+1
j , unj en −un+1
j , −unj et en
utilisant le résultat pour le max pour les −un+1
j , −unj .
Supposons donc que le premier maximum est atteint en (n + 1, j0 ). Cela
entraine, en utilisant le schéma :
1 ∆t
un+1 un+1 n+1 n+1
= unj0 .


j0 − 2 j0 +1 + uj0 −1 − 2uj0
2 (∆x)

Mais un+1 n+1

j0 +1 ≤ uj0 et un+1 n+1
j0 −1 ≤ uj0 , donc :

un+1 n+1 n+1

j0 +1 + uj0 −1 − 2uj0 ≤ 0,

(ce qui n’est rien d’autre qu’une version “discrète” du fait que la dérivée
seconde est négative au sens large en un point de maximum !) et ainsi :

max un+1
j = un+1
j0 ≤ unj0 ≤ max unj .
j∈Z j∈Z

Il est à noter que cette propriété signifie que le schéma est incondition-
nellement monotone, i.e.

unj ≤ 0 pour tout j ∈ Z =⇒ un+1

j ≤0 pour tout j ∈ Z .

Comme le schéma (SES), le schéma (SIS) est consistant et d’ordre 1 en

temps, 2 en espace.

3.5 Convergence des schémas

On a le résultat suivant :

Théorème 5. On suppose que u0 est bornée, uniformément continue sur

R. Alors la solution (unj )n,j des schémas (SES) quand la condition de CFL
est satisfaite et (SIS) (inconditionnellement) converge uniformément vers u
sur R × [0, T ] pour tout T > 0, i.e.

max unj − u(j∆x, n∆t) → 0 quand ∆x + ∆t → 0 .

j∈Z
n∈N

25
Preuve :
(k)
Etape 1 : on traite le cas où u0 est de classe C 4 avec u0 bornées pour
0 ≤ k ≤ 4. On injecte les valeurs de u(j∆x, n∆t) dans le schéma et on utilise
la monotonie pour montrer que :

unj − u(j∆x, n∆t) ≤ Cn∆t (∆x)2 + ∆t ,

 

où C[(∆x)2 + ∆t] est le terme d’erreur que l’on commet à chaque étape en
remplaçant unj par u(j∆x, n∆t), et :

Cn∆t (∆x)2 + ∆t
 

et une “sur-solution” qui permet de prendre en compte cette erreur.

Etape 2 : on approche u0 qui est bornée et uniformément continue par
des fonctions (uε0 )ε qui sont de classe C 4 avec les quatre permières dérivées
bornées. Cette approximation uniforme (obtenue par convolution) est rendue
possible par l’uniforme continuité de u0 . On raisonne ensuite comme dans
le Chapitre 1, en utilisant le formule explicite qui donne u, uε ,... etc. 
(k)
Remarque : Dans le cas où u0 est C 4 , avec u0 borné pour 0 ≤ k ≤ 4, on a
bien une convergence en (∆t)1 + (∆x)2 comme le prévoit l’ordre du schéma.

3.6 Aproximation par domaine de calcul borné

A cause de la vitesse infinie de propagation des informations, ce problème
apparait comme non trivial. Pourtant, si l’on considère par exemple le problème
suivant :
∂uR 1 ∂ 2 uR
(1’) − =0 dans ] − R, R[×(0, ∞)
∂t 2 ∂x2
avec la condition initiale suivante :

(5’) uR (x, 0) = u0 (x) dans ] − R, R[

et la condition de bord :

(4’) uR (±R, t) = u0 (±R) pour t > 0 ,

on peut montrer que uR → u localement uniformément quand R → ∞. La

semi-discrétisation en temps que nous avons évoqué montre que ce problème
se ramène à ce que nous avons vu au Chapitre 1.
N.B. : la condition aux limites artificielle (4’) peut (en fait, doit !) être
remplacée la plupart du temps par une condition plus adaptée et plus astu-
cieuse ! ! !
Exercice 6. Programmer au moins l’un des schémas et examiner l’effet des
conditions aux limites artificielles.

26
3.7 Équation de la chaleur en domaines bornés
On peut s’intéresser à l’équation de la chaleur dans un domaine borné
(par exemple [0, 1]) :

∂u 1 ∂ 2 u
(9) − =0 dans [0, 1] × (0, T ) .
∂t 2 ∂x2
avec une donnée initiale :

(10) u(x, 0) = u0 (x) dans [0, 1] .

Mais il faut ajouter des conditions aux limites sur les bords de l’intervalle
[0, 1] (ou sur les ensembles {0} × (0, T ) et {1} × (0, T )). Parmi les diverses
conditions existentes, deux jouent un rôle plus importants :
1. Conditions de Dirichlet : u(a, t) = g(t) pour a = 0 et/ou 1. On fixe
donc la température à une extrémité de l’intervalle ou aux deux.
∂u
2. Conditions de Neumann : = h(t). Ici on fixe le flux de chaleur.
∂x
Par exemple, si h ≡ 0, on a une condition d’adiabaticité : aucun flux
de chaleur ne traverse la paroi (isolation parfaite).
Évidemment on peut combiner ces deux types de conditions au limites
en imposant Dirichlet à une extremité et Neumann à l’autre.

Pour traiter ce problème, nous allons disposer d’un outil différent de

type série de Fourier. Pour justifier son introduction, on va considérer des
conditions de Dirichlet homogène :

(11) u(0, t) = u(1, t) = 0 pour tous t > 0 ,

et on va chercher des solutions à variables séparées donc de la forme :

u(x, t) = ψ(t)φ(x) .

On a donc :
ψ 0 (t)φ(x) − ψ(t)φ00 (x) = 0 ,
pour tous x et t. Seule possibilité pour obtenir des solutions non triviales :

ψ 0 (t) = λψ(t) pour tous t > 0,

λφ(x) − φ00 (x) = 0 dans (0, 1) ,

sans oublier la condition aux limites :

φ(0) = φ(1) = 0 .

27
Pour ψ, on a ψ(t) = c exp(λt) où c ∈ R et pour φ, on a un problème de type
valeur propre qui conduit à :

λk = k 2 π 2 et φk (x) = sin(kπx) ,

pour k ∈ N∗ .
Phrase compliquée : la théorie des opérateurs compacts dans L2 ([0, 1])
implique que les φk forment une base hilbertienne de L2 ([0, 1]) (l’opérateur
compact est T : L2 ([0, 1]) → L2 ([0, 1]) tel que T (f ) = u où u est l’unique
solution de :
−u00 (x) = f dans (0, 1) ,
avec la condition aux limites :

u(0) = u(1) = 0 .)

Quoiqu’il en soit, si on peut écrire que :

+∞
X
u0 (x) = ak φk (x) ,
k=1

alors la solution est donnée par :

+∞
X
u(x, t) = ak exp(−k 2 π 2 t)φk (x) ,
k=1

et on voit facilement que même si u0 ∈ L2 ([0, 1]) alors u est C ∞ (car ? ? ? ?).

Exercice 7. Faire le même travail pour l’équation des cordes vibrantes

(fixées aux extrémités) :

∂2u 2
2∂ u
(12) − c =0 dans R × (0, T ) ,
∂t2 ∂x2
où c > 0. Quelles sont les différences avec l’équation de la chaleur ? Com-
ment s’interprète c ?

3.8 Quelques éléments de modélisation du trafic routier

On va reprendre quelques éléments de la modélisation conduisant à
l’équation de la chaleur mais dans un contexte différent.
Cette fois la droite réelle est une route à sens unique (de la gauche vers
la droite) où se déplacent des véhicules.
On note u(x, t) la densité de véhicules au point x à l’intant t. La variation
de la quantité de véhicules dans un intervalle [x − ∆x, x + ∆x] est égale au

28
flux entrant (par le point x − ∆x) et sortant (par le point x + ∆x). Si φ est
ce flux, on a :
Z x+∆x 
d
u(y, t)dy = φ(x − ∆x, t) − φ(x + ∆x, t) ,
dt x−∆x

Le changement de signe du flux est ici évident. On en déduit :

Z x+∆x Z x+∆x
∂u ∂φ(x, t)
dy = − (y, t)dt ,
x−∆x ∂t x−∆x ∂x

et comme ∆x est arbitraire, les fonctions u et φ sont liées par l’équation :

∂u ∂φ
(13) + =0 dans R × (0, +∞) .
∂t ∂x
Cette équation ayant deux fonctions inconnues, il n’est pas possible de
la résoudre sans lier ces deux fonctions par une loi de comportement liant le
flux de véhicule et sa densité. Si on normalise u de telle sorte que 0 ≤ u ≤ 1
alors on peut penser a une loi du type :
1
φ = u(1 − u) .
2
Le flux est petit si u est proche de 0 car il y a peu de véhicules, il est
petit quand il y a trop de véhicules (i.e. pour u proche de 1) à cause de
l’engorgement et u = 1/2 donne le flux maximal (trafic ni trop important,
ni trop faible).
En posant v = 1/2−u, comme φ = −v 2 /2+1/8, on se ramène à l’équation
de Burgers :
∂v ∂v
(14) +v =0 dans R × (0, +∞) .
∂t ∂x
Le cours sur les équation de transport indique que, pour résoudre de
telles équations, on doit chercher des courbes caractéristiques (x(t), t) le
long desquelles v est constante. En dérivant v(x(t), t) = constante, il vient :
∂v ∂v
+ ẋ(t) =0
∂t ∂x
et en identifiant avec l’équation, on voit que ẋ(t) = v(x(t), t) = v(x(0), 0).
Les courbes caractéristiques sont des droites et l’on voit facilement que ces
droites peuvent se couper : ceci signifie que les solutions de ces équations
ne sont pas régulières en général (sinon le calculs ci-dessus serait justifié et
on aurait une contradiction...) et il faut employer la notion de solution au
sens des distributions (ce qui ne suffit toujours pas...). Il y a, en général,
plusieurs solutions au sens des distributions et il faut employer la notion de
solutions entropiques qui ne rentre pas dans le cadre de ce cours.

29
 Le problème de Riemann est un cas particulièrement intéressant car on
sait calculer l’unique solution entropique pour l’équation de Burgers. C’est
le cas où la donnée initiale est donnée par :
(
v− si x < 0 ,
v(x, 0) =
v + si x ≥ 0 .

On a deux cas :
Si v− > v+ , la solution présente une discontinuité et elle est donnée par :

v+ + v−
t x= t
6

2












v(x, t) = v−
v(x, t) = v+











-
x
Si v− < v+ , la solution est continue et elle est donnée par :

t x = v− t x = v+ t
6  
 
 
 
 

 v(x, t) = x 
 t 
 

v(x, t) = v− 

 v(x, t) = v+

 

 
 
 
 -
x

Dans le cas du trafic routier, on a une interprétation assez simple de

ces solutions. Rappelons d’abord que la densité de véhicule u vaut 1 − v.

30
Si v− > v+ , comme le trafic a lieu vers la droite, il y a plus de véhicule en
aval qu’en amont : les véhicules arrivant de l’amont se retrouve donc face
à une concentration de véhicules plus grande devant eux et donc face à un
embouteillage (un “choc”). Au contraire, v− < v+ , il y a plus de véhicule en
amont qu’en aval : les véhicules arrivant de l’amont se retrouve donc face à
une concentration de véhicules plus faible devant eux. La route se dégage,
c’est une “onde de détente”.

Exercice 8. Vérifier que la première solution est TOUJOURS solution au

sens des distributions (que l’on ait v− > v+ ou v− < v+ ) et que la deuxième
l’est dans le cas où elle est bien définie, i.e. v− < v+ .

Reste à calculer numériquement ces solutions : les trois schémas les plus
classiques se mettent sous forme conservative :

vjn+1 = vjn − λ g(vjn , vj+1

n n
, vjn ) ,


) − g(vj−1
où g est une fonction telle que g(v, v) = v 2 /2 et λ = ∆t/∆x.
Le schéma le plus spécifique dans le cas non linéaire est celui de Godounov
où :
g(v− , v+ ) = w(0, v− , v+ ) ,
où w(0, v− , v+ ) est la solution du problème de Riemann associée aux données
v− , v+ qui ne dépend pas de t vu la forme des solutions (w = v+ , v− , 0 suivant
divers cas).
Le schémas de Lax-Friedrichs est celui où :
 
1 1 2
g(v− , v+ ) = (v+ + v−2
) − λ−1 (v+ − v− ) ,
2 2

et celui de Lax-Wendroff est donné par :

 
1 1 2 v+ + v− 2
g(v− , v+ ) = 2
(v+ + v− ) − λ−1 ( 2
)(v+ − v− ) .
2 2 2

Exercice 9. Vérifier la consistance de ces schémas, programmer les et com-

parer les.

4 Modélisation et optimisation
On se propose ici de traiter un problème particulier sur la composition
chimique d’un mélange de gaz à pression et température données (texte
proposé à l’épreuve de Calcul scientifique de l’agrégation).
Après réaction chimique, un mélange de différents gaz à température T et
pression p fixées doit normalement se trouver à l’équilibre thermodynamique.

31
On suppose que ce mélange est composé de ns espèces chimiques contenant
au total ne types d’atomes différents, que l’on nomme les éléments.
Une mole de l’espèce k est la quantité de matière correspondant à un
nombre de molécules de cette espèce égal au nombre d’Avogadro. On note Nk
le nombre de moles de l’espèce k et on introduit le vecteur N = (N1 , · · · , Nns )
de Rns et si N̄ = N1 + · · · + Nns , on appelle fraction molaire de l’espèce k,
la quantité :
Nk
Xk = .
N̄
Et X = (X1 , · · · , Xns ). À l’équilibre chimique, les Nk et les Xk sont des
inconnues.

Le nombre d’atomes de l’élément numéro j dans une molécule de l’espèce

k sera noté Ekj et on note E la matrice (Ekj )kj . C’est une matrice dont les
éléments sont des entiers positifs ou nuls donnés. Chaque ligne et chaque
colonne de E admet au moins un élément non nul.
Le nombre de moles d’atome j dans le mélange est noté Nje et on a :
ns
X
Nje = Nk .Ekj .
k=1

On regroupe également ces valeurs en un vecteur N e

On considère un système fermé, donc le nombre de moles de chaque type
d’atomes se conserve au cours de la réaction. Ceci signifie que le vecteur N e
est une donnée du problème.
Le vecteur des nombres de moles N est donc soumis à ne contraintes
d’égalité, linéaires, indépendantes si l’on fait l’hypothèse que la matrice E
est de rang ne et qui s’écrivent :

(15) N e = N.E = N0 .E ,

si N0 est la composition initiale du mélange de gaz (avant les réactions

conduisant à l’équilibre thermodynamique).
À p et T fixées, une composition chimique d’équilibre, c’est-à-dire le
vecteur N que l’on cherche à déterminer, est celle qui réalise un minimum
d’une fonction G (appelée enthalpie libre ou énergie de Gibbs), sous les
contraintes ci-dessus. Ici, G peut s’expliciter sous la forme :
ns
X ns
X Xns
(16) G(N ) = Nk (gk + log(Xk )) = Nk (gk + log(Nk ) − log( Nj )) ,
k=1 k=1 j=1

où gk est une constante à p et T fixées appelée enthalpie libre molaire de

l’espèce k.
Le problème mathématique à résoudre est donc de minimiser la fonction
G donnée par (16) sous la contrainte (15).

32
 Il s’agit d’abord de considérer l’existence et l’unicité de la solution. Pour
cela, on se ramène à un théorème classique puisque les Nk satisfont :

0 ≤ Nk ≤ Nksup := min(Nje /Ekj ) .

j

(pourquoi ?) Ceci (avec (15)) définit un compact de Rns et on vérifie facile-

ment que G est continue sur ce compact.

Pour l’unicité, on procède en deux temps : on calcule le gradient de G

dont les composantes valent :
∂G
= gk + log(Xk ) .
∂Nk
On vérifie d’abord que le minimum n’est pas atteint sur un point de la
frontière du domaine de contraintes (c’est-à-dire pour Nk = 0 ou Nk =
Nksup := minj (Nje /Ekj ). Pour cela, on montre que la dérivée normale tend
vers −∞ en un point où Nk = 0 et un tel point existe dans les deux cas
(même le second).
Puisque le point de minimum est atteint à l’intérieur, il s’agit de voir
que G est strictement convexe sur la contrainte. Or on a :
∂2G
= δk,j (Nk )−1 − (N̄ )−1 ,
∂Nk ∂Nj
et il faut tester D2 G suivant les directions possibles ∆N le long de la
contrainte, pour lesquelles on a ∆N.E = 0 et donc ∆N change de signe.
On a :
ns ns
!2
2
X [(∆N )k ]2 X
D G(N )∆N · ∆N = − (∆N )k /N̄ .
Nk
k=1 k=1

Par Cauchy-Schwarz, cette quantité est strictement positive puisque ∆N ne

peut être colinéaire au vecteur N .
Pour calculer numériquement l’équilibre, il est plus simple d’utiliser que
la contrainte est affine et d’introduire une base (fj )j de l’espace vectoriel
associé (donc une base du noyau de la transposée de E). En écrivant :

N = N0 + x1 .f1 + ... + xl .fl ,

(que vaut l ?), on se ramène à un problème d’optimisation sans contrainte sur

(x1 , · · · , xl ) pour lequel on peut appliquer des méthodes classiques (gradient
conjugué, Newton,...).
Exercice 10. Traiter complètement le cas du mélange H2 -O2 -H2 O (en pre-
nant des gk égaux respectivement à −18, −27, −45).
Puis celui faisant aussi intervenir les espèces OH, O et H (en prenant pour
ces espèces des gk égaux respectivement à −21, −2, +2).

33
5 Appendice : Compléments
5.1 Étude générale des méthodes à un pas
T
On conserve, dans cette section, une grille uniforme de pas h = N. Les
méthodes à un pas sont des méthodes de la forme :

yi+1 = yi + hΦ(ti , yi , h)
y0 = y0,h

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

5.1.1 Propriétés importantes d’une méthode à un pas

• CONSISTANCE

Définition 1. On appelle erreur de consistance de la méthode à un pas, la

quantité :
N
X −1
Σh = |y(ti+1 ) − y(ti ) − hΦ(ti , y(ti ), h)| .
i=0

La méthode est dite consistante si Σh → 0 quand h → 0.

La quantité y(ti+1 ) − y(ti ) − hΦ(ti , y(ti ), h) est l’analogue de ce que nous

avons noté εi ci-dessus ; nous conserverons au besoin cette notation.

• STABILITÉ

Définition 2. La méthode à un pas est dite stable s’il existe deux constantes
S1 , S2 telles que, si (ỹi )i est défini par :

ỹi+1 = ỹi + hΦ(ti , ỹi , h) + εi
,
ỹ0 = ỹ0,h

alors :
N
X −1
max |yi − ỹi | ≤ S1 |y0,h − ỹ0,h | + S2 |εi | .
i
i=0

Bien entendu, la méthode est dite convergente si max |yi − y(ti )| → 0

i
quand h → 0.

Théorème 6. Toute méthode stable et consistante converge à condition que

y0,h → y(0) quand h → 0.

34
 La dernière condition étant toujours satisfaite en pratique, l’étude de la
convergence des méthodes à un pas se réduit à l’étude de leur consistance
et de leur stabilité, ce qui est plus simple comme on va le voir.
Preuve : Si on note ỹi = y(ti ), on a par définition de εi (juste après la
définition de la consistance) :

ỹi+1 = ỹi + hΦ(ti , ỹi , h) + εi .

Et ỹ0 = y(0). Puisque la méthode est stable :

N
X −1
max |yi − y(ti )| ≤ S1 |y0,h − y(0)| + S2 |εi | .
i
i=0

Hors, les deux quantités du membre de droite tendent vers 0 quand h tend
vers 0 par consistance, donc le résultat est acquis. 

5.1.2 Condition nécessaire et suffisante de consistance

Théorème 7. La méthode à un pas est consistante si et seulement si :

Φ(t, z, 0) = f (t, z) ,

pour tous t ∈ [0, T ], z ∈ Rn .

Preuve : On ne va vraiment détailler que la condition suffisante.

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

Z ti+1
= [f (s, y(s)) − Φ(ti , y(ti ), h)]ds .
ti

Mais, pour s ∈ [ti , ti+1 ] :

|f (s, y(s)) − Φ(ti , y(ti ), h)| ≤ |f (s, y(s)) − f (ti , y(ti ))|+

|f (ti , y(ti )) − Φ(ti , y(ti ), 0)| + |Φ(ti , y(ti ), 0) − Φ(ti , y(ti ), h)|
Le premier et le troisième terme sont des termes petits (uniformément en i)
par l’uniforme continuité de f , y et Φ ; on les estime par un δ(h) qui tend
vers 0 avec h.
Si le terme du milieu est nul (condition de consistance) alors |εi | ≤ hδ(h)
et Σh ≤ N hδ(h) = T δ(h) → 0 quand h → 0. D’où la consistance.
NB : la condition suffisante se prouve en examinant la preuve d’un peu plus
près : si le terme du milieu n’est pas nul... 

35
5.1.3 Condition suffisante de stabilité
Théorème 8. La méthode à un pas est stable si Φ(t, z, h) est lipschitzien en
z pour tous t ∈ [0, T ], h ∈ [0, H] avec une constante de lipschitz indépendante
de t et h.

Preuve : On note θi = |yi − ỹi |. On a :

θi+1 = |yi + hΦ(ti , yi , h) − ỹi − hΦ(ti , ỹi , h) − εi | .

D’où :
θi+1 ≤ θi + h|Φ(ti , yi , h) − Φ(ti , ỹi , h)| + |εi | .
Si Φ est lipschitzienne en sa deuxième variable de constante de lipschitz L̃ :

θi+1 ≤ (1 + L̃h)θi + |εi | .

Une petite amélioration du Lemme de Gronwall discret nous donne :

i
X
θi+1 ≤ exp(L̃(i + 1)h)θ0 + exp(L̃(i − k)h)|εk | .
k=0

En majorant exp(L̃(i + 1)h) et exp(L̃(i − k)h) par exp(L̃T ), on a le résultat

avec S1 = S2 = exp(L̃T ). 

5.1.4 Ordre d’un schéma

La question que l’on se pose ici est la suivante : peut-on avoir une
meilleure précision que dans le cas de la méthode d’Euler en choisissant
bien la méthode à un pas et comment faut-il la choisir ?

Définition 3. On dit qu’une méthode à un pas est d’ordre p ≥ 1 si, pour

toute solution de l’EDO, il existe une constante C > 0 telle que :

|εi | = |y(ti+1 ) − y(ti ) − hΦ(ti , y(ti ), h)| ≤ Chp+1 .

Pourquoi le “p + 1” dans l’ordre p ? Deux raisons concourantes :

y(ti+1 ) − y(ti )
1. | − Φ(ti , y(ti ), h)| ≤ Chp et la quantité considérée ap-
h
proche l’équation à l’ordre p.
N
X −1
2. Σh = |y(ti+1 ) − y(ti ) − hΦ(ti , y(ti ), h)| ≤ Chp , donc ordre p =
i=0
erreur en hp .
Cette dernière idée est justifiée par le résultat suivant dont la preuve est
immédiate à partir du résultat de convergence :

36
Corollaire 1. Si la méthode à un pas est d’ordre p ≥ 1 et si |y(0) − y0,h | ≤
C̃hp alors maxi |yi − y(ti )| ≤ Ĉhp .

Nous donnons maintenant une condition nécessaire et suffisante pour

qu’une méthode à un pas soit d’ordre p dans R. Pour cela, nous commençons
par un résultat de régularité sur la solution y de l’EDO : c’est une règle
générale, on ne sait bien approcher que des fonctions régulières et donc nous
aurons besoin d’une certaine régularité de la solution. Nous formulons ce
résultat en dimension 1, c’est-à-dire quand y est à valeurs dans R.

Théorème 9. Si la fonction f est de classe C p sur [0, T ] × R alors y est de

classe C p+1 et :

y (k+1) (t) = f [k] (t, y(t)) dans ]0, T [ (k = 0, 1, · · · , p),

où les fonctions f [k] sont définies sur [0, T ] par :

f [0] (t, z) = f (t, z) ,

∂ [k] ∂
f [k+1] (t, z) = f (t, z) + f (t, z) f [k] (t, z) ,
∂t ∂y
pour k = 0, 1, · · · , p − 1.

La preuve se fait facilement et elle est laissée en exercice (où l’on pourra
aussi réfléchir à la formulation de ce résultat si y est à valeurs dans Rn ).
La condition nécessaire et suffisante pour qu’une méthode à un pas soit
d’ordre p dans R est donnée par le :

Théorème 10. On suppose que f est de classe C p et que, pour k ≤ p, les

∂kΦ
dérivées partielles existent et sont continues sur [0, T ] × R × [0, H].
∂hk
Alors la méthode à un pas est d’ordre p si, pour tous t ∈ [0, T ], z ∈ R :

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

Preuve : Comme f est de classe C p , y est de classe C p+1 et y (k) (t) =

f [k−1] (t, y(t)) par le résultat sur la régularité des solutions.

37
 Par la formule de Taylor, on a :
p
X hk
y(ti+1 ) − y(ti ) = y (k) (ti ) + O(hp+1 )
k!
k=1
p
X hk [k−1]
= f (ti , y(ti )) + O(hp+1 )
k!
k=1

et :
p
X hk ∂ k Φ
Φ(ti , y(ti ), h) = (ti , y(ti ), 0) + o(hp ) .
k! ∂hk
k=0
Il en résulte que :
p
X hk+1 ∂ k Φ
hΦ(ti , y(ti ), h) = (ti , y(ti ), 0) + o(hp+1 )
k! ∂hk
k=0
p+1
X hk ∂ k−1 Φ
= (ti , y(ti ), 0) + o(hp+1 )
(k − 1)! ∂hk−1
k=1

et :
p 
∂ k−1 Φ hk

X 1 [k−1]
εi = f (ti , y(ti )) − (ti , y(ti ), 0) + O(hp+1 ) .
k ∂hk−1 (k − 1)!
k=1

Sous les hypothèses du théorème, tous les termes entre crochets sont nuls et
le résultat est acquis.
On peut le préciser avec la formule de Taylor avec reste intégral, ce qui
donne la deuxième partie du théorème. 

5.1.5 Exemples
On va chercher la meilleure fonction Φ, c’est-à-dire celle qui donne le
meilleur ordre de convergence, parmi celles qui sont de la forme :

Φ(t, z, h) = a1 f (t, z) + a2 f (t + p1 h, z + p2 hf (t, z)) ,

où a1 , a2 , p1 , p2 sont des paramètres à fixer “au mieux”.

On a une méthode consistante d’ordre 1 si :

Φ(t, z, 0) = f (t, z) ,

donc si :
a1 f (t, z) + a2 f (t, z) = f (t, z) ,
d’où la première condition a1 + a2 = 1.

38
 Pour avoir de l’ordre 2, il faut que :
 
∂Φ 1 1 ∂f ∂f
(t, z, 0) = f [1] (t, z)) = (t, z) + f (t, z) (t, z) .
∂h 2 2 ∂t ∂y

Or :
∂Φ ∂f ∂f
(t, z, h) = a2 p1 (t+p1 h, z+p2 hf (t, z))+a2 p2 f (t, z) (t+p1 h, z+p2 hf (t, z)) .
∂h ∂t ∂y

La deuxième condition est donc a2 p1 = a2 p2 = 12 .

Par contre, on vérifie facilement (le faire !) que l’on ne peut pas aller plus
loin. En prenant, a2 = α comme paramètre, on a une famille de méthodes à
un pas d’ordre 2 :
h h
Φ(t, z, h) = (1 − α)f (t, z) + αf (t + ,z + f (t, z)) .
2α 2α
Ces méthodes sont connues sous le nom de :
— α = 1, méthode de la tangente améliorée,
— α = 1/2, méthode d’Euler modifiée,
— α = 1, méthode de Heun.

39