MACS 205: Analyse numérique

notes de cours (succintes)

A. Sabourin, Olivier Fercoq, François Portier

3 février 2021
Table des matières

1 Fondements du calcul scientifique 4

1.1 Exemples, sources d’erreur, problèmes bien posés, conditionnement . . . . . . 4
1.1.1 Sources d’erreur . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.1.2 Problème bien posé . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.1.3 Stabilité . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.1.4 Calcul de conditionnement lorsque G est différentiable . . . . . . . . . 8
1.2 Convergence des méthodes numériques . . . . . . . . . . . . . . . . . . . . . . 8
1.2.1 Définitions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.3 Erreurs d’arrondis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.3.1 Représentation en base β . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.3.2 Nombres à virgule flottante ("float") . . . . . . . . . . . . . . . . . . . 9
1.3.3 Opérations sur les flottants . . . . . . . . . . . . . . . . . . . . . . . . 10

2 Interpolation Polynomiale 11
2.1 Interpolation de Lagrange . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.1.1 Existence, unicité . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.1.2 Erreur d’interpolation de Lagrange . . . . . . . . . . . . . . . . . . . . 12
2.1.3 Contre exemple de Runge . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.1.4 Stabilité . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.1.5 Lien entre stabilité et convergence . . . . . . . . . . . . . . . . . . . . 15
2.1.6 Ordre de grandeur de la constantes de contrôle de l’erreur pour n grand 15
2.2 Interpolation aux noeuds de Tchebychev . . . . . . . . . . . . . . . . . . . . . 16
2.2.1 Définition et premières propriétés . . . . . . . . . . . . . . . . . . . . . 16
2.2.2 Interpolation de Tchebychev sur [a, b] . . . . . . . . . . . . . . . . . . 17
2.3 Différences divisées, Base de Newton . . . . . . . . . . . . . . . . . . . . . . . 18
2.3.1 Base de Newton . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.3.2 Différences divisées . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.3.3 Méthode de Horner pour l’évaluation . . . . . . . . . . . . . . . . . . 20

3 Intégration numérique : méthodes de quadrature 21

3.1 Introduction : méthodes de quadrature, méthodes de Newton-Cotes . . . . . . 21
3.2 Exemples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
3.2.1 Rectangles à gauche / à droite . . . . . . . . . . . . . . . . . . . . . . 23
3.2.2 Méthode du point Milieu . . . . . . . . . . . . . . . . . . . . . . . . . 23
3.2.3 Méthode des trapèzes . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
3.2.4 Méthode de Cavalieri-Simpson . . . . . . . . . . . . . . . . . . . . . . 24

1
 3.3 Ordre des méthodes de Newton-Cotes . . . . . . . . . . . . . . . . . . . . . . 24
3.4 Stabilité . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.4.1 Stabilité d’une MQS . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.4.2 Stabilité d’une MQC . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.5 Majorations de l’erreur . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.5.1 Erreur d’une MQS d’ordre N . . . . . . . . . . . . . . . . . . . . . . . 26
3.5.2 Erreur d’une MQC d’ordre N sur des intervalles de même longueur . . 27
3.5.3 Erreur des méthodes de Newton-Cotes . . . . . . . . . . . . . . . . . . 28
3.6 Évaluations de l’erreur a posteriori . . . . . . . . . . . . . . . . . . . . . . . . 28
3.7 Extrapolation de Richardson, application aux trapèzes (méthode de Romberg) 29
3.7.1 Principe de la méthode de Richardson . . . . . . . . . . . . . . . . . . 29
3.7.2 Méthode de Romberg . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

4 Équations différentielles 32
4.1 Existence et unicité des solutions . . . . . . . . . . . . . . . . . . . . . . . . . 32
4.1.1 Problème de Cauchy . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
4.1.2 Exemple : les équations du mouvement . . . . . . . . . . . . . . . . . . 32
4.1.3 Théorèmes d’existence et unicité des solution . . . . . . . . . . . . . . 33
4.2 Équations différentielles linéaires . . . . . . . . . . . . . . . . . . . . . . . . . 35
4.2.1 Géométrie de l’ensemble des solutions . . . . . . . . . . . . . . . . . . 35
4.2.2 Résolution de y 0 = Ay : équation différentielle linéaire à coefficients
constants et sans second membre . . . . . . . . . . . . . . . . . . . . . 35
4.2.3 Résolution de y 0 = Ay + b(t) : équa. diff. linéaire à coefficients constants 36
4.3 Méthodes numériques à un pas . . . . . . . . . . . . . . . . . . . . . . . . . . 37
4.3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
4.3.2 Consistance, stabilité, convergence . . . . . . . . . . . . . . . . . . . . 38
4.3.3 *Condition nécessaire et suffisante de consistance . . . . . . . . . . . . 40
4.3.4 Condition suffisante de stabilité . . . . . . . . . . . . . . . . . . . . . . 41
4.3.5 Majoration de l’erreur globale . . . . . . . . . . . . . . . . . . . . . . . 41
4.4 Méthodes de Runge-Kutta . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
4.4.1 Principe . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
4.4.2 Exemples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
4.4.3 Stabilité des méthodes de Runge-Kutta . . . . . . . . . . . . . . . . . 44
4.4.4 Ordre des méthodes de Runge-Kutta . . . . . . . . . . . . . . . . . . . 45

Ce cours d’analyse numérique (calcul scientifique) s’appuie principalement sur les ouvrages
de Quarteroni et al. [2008] et Demailly [2012]. Les deux premiers chapitres de ce document
sont issus des notes de cours rédigées par Baptiste Portenard, que nous remercions chaleu-
reusement.

2
Fonctionnement du cours

Les cours de Macs 205 (a) et (b), sont en pratique regroupés en un seul bloc Macs 205.
Celui-ci est divisé en trois partie, correspondant chacune à un intervenant :
(1) Fondements du calcul scientifique, interpolation polynomial et intégration numérique →
Anne Sabourin
(2) Méthodes numériques pour équations différentielles → Olivier Fercoq
(3) Méthodes de Monte-Carlo → François Portier.

Déroulé du cours : alternance de cours, TD et TP. La partie ’implémentation’ des mé-

thodes prend un place relativement importante (≈ la moitié du temps). Le langage utilisé
est R.

Evaluation Elle comprend une partie théorique (examens écrits) et une partie pratique
(mini-projets). Pour des raisons administratives, deux notes indépendantes doivent être at-
tribuées. Chaque note est calculée sur le barême suivants :
• MACS 205-a :
— Partie (1) : coeff 1
— Partie (2) : coeff 1/2
• MACS 205-b
— Partie (2) : coeff 1/2
— Partie (3) : coeff 1

3
Chapitre 1

Fondements du calcul scientifique

1.1 Exemples, sources d’erreur, problèmes bien posés, condi-

tionnement
But de l’analyse :
Résoudre un problème mathématique de type

F (x, d) = 0,

où F est une fonction à valeurs réelles qui représente le « modèle mathématique », x est l’in-
connue recherchée et d représente l’ensemble des données du problèmes.
On note x0 une solution en x de ce problème. On a recourt à des méthodes numériques lorsque
x0 n’a pas d’expression explicite.

Méthode générale :
On construit un algorithme produisant une suite de solutions approchées (xn )n>1 du problème
approché
Fn (x, dn ) = 0
dans l’espoir que la suite (xn ) tende vers une solution x0 du problème initial,

xn → x0 .

Pour l’instant on n’a pas précisé dans quel sens cette convergence doit avoir lieu. La notion
de convergence dépendra généralement du choix d’une norme dans l’espace des solutions.

1.1.1 Sources d’erreur

La résolution du problème initial concret (que l’on appellera problème physique passe par la
formulation successive d’un problème mathématique (censé approcher le problème physique),
puis d’un problème numérique approché, dont la résolution entraîne à son tour des erreurs de
mesure et d’arrondi. Ces différentes étapes sont illustrées dans l’exemple suivant

Exemple 1.1 (Problème du robinet):

On considère une baignoire de volume V et un robinet de débit q(t), t > 0. On désire connaitre x
le temps nécessaire pour remplir une baignoire.

4
Problème physique On peut représenter mathématiquement le problème physique sous la
forme
Fphys (x, dphys ) = 0
avec les paramètres suivants

dphys = (q : R+ → R+ , V > 0)
Z x
Fphys (x, V, q) = V − q(t)dt
0
Etape de modélisation : Modification des données du problèmes. On suppose la baignoire
est un parallélépipède rectangle (i.e., un pavé), de sorte que V ' l.L.h.

Problème mathématique Le problème devient :

F (x, d) = 0, d = {l, L, h, q}

où F s’écrit Z x
F (x, l, L, h, q) = l.L.h − q(t)dt.
0
Du fait de l’étape de modélisation, on a x0 6= xphys .
Etape de discrétisation : Modification des données du problèmes
Z x [Link]
X 1 k
q(t)dt ' q( )
0 k=0
n n

Problème numérique l’équation définissant le problème numérique est à présent

k
Fn (x, dn ) = 0, dn = {q( ), L, h, q},
n
où :
[Link]
X 1 k
Fn (x, dn ) = l.L.h − q( ).
k=0
n n

Erreur de mesure, arrondi La mesure des données est soumise a des approximations, par
exemple pour calculer un débit on moyenne sur une petite durée d’écoulement. On note avec un
chapeau les variables mesurées et on obtient la solution xcn du problème approché
[Link]
b−
X 1 k
F
cn (x, d l.L.
cn ) = b bh qb( )
k=0
n n

Dans ce cours, on ne s’intéressera pas à l’erreur venant de l’étape de modélisation. Autre-

ment dit, on considérerar que x0 est la « vraie » solution et on analysera les erreurs venant
des étapes suivantes de l’exemple ci-dessus.

Définition 1.1.1 (Erreurs). On définit les différentes erreurs suivantes

• Erreur totale = xcn − xphys
• Erreur d’arrondi = xcn − xn

5
 • Erreur de troncature = xn − x0
• Erreur de modèle = x0 − xphys .
On appelle erreur de calcul la somme des erreurs d’arrondi et de troncature (i.e. xcn − x0 )
Comme xcn − xphys = (xcn − xn ) + (xn − x0 ) + (x0 − xphys ), on en déduit que l’erreur totale
est la somme des différentes autres erreurs de la définition précédente. En particulier,
xcn − xphys = xc − x + x0 − xphys
| n {z 0} | {z }
erreur de calcul erreur de modèle

Et on cherchera à contrôler le terme d’erreur de calcul dans la décomposition ci-dessus, en

fonction de
• δd = dcn − dn : erreur d’arrondi.
[Link]
• l’erreur de troncature des données dn −d (dans l’exemple, dn −d = k=0 n1 q( nk )− 0x q).
R

La sensibilité d’une méthode numérique par rapport à δd est résumée par la notion de stabilité
et de conditionnement (voir le paragraphe 1.1.3). La sensibilité par rapport à dn −d est résumée
par la notion de consistance (voir le paragraphe 1.2).

1.1.2 Problème bien posé

Définition 1.1.2 (Problème bien posé). Un problème mathématique
F (x, d) = 0
est bien posé si
• Pour toutes données d il existe une unique solution x0 = G(d). On appelle G la fonction
résolvante.
• G est continue
Exemple 1.2:
On considère le problème suivant
F (x, d) = x2 − 2xd + 1
dont on cherche les solutions (x1 , x2 ) dans R.
Le discriminant réduit est ∆ = d2 − 1
Si |d| < 1 alors le problème est mal posé.
sinon on a √ √
G(d) = (d − ∆, d + ∆)
p p
G(d) = (d − d2 − 1, d + d2 − 1)
et G est continue sur R\ ] − 1, 1[

1.1.3 Stabilité
Définition 1.1.3 (Stabilité). Un problème F (x, d) = 0 est stable si la résolvante G(d) est
localement lipschitzienne.
∀d, ∃η > 0, ∃K0 (d) tel que khk < η ⇒ kx(d + h) − x(d)k 6 K0 (d)khk (1.1)
où x(d) = G(d)

6
Proposition 1.1.4
La condition (1) de stabilité est équivalente à
x(d + h) − x(d)
lim supkhk< < +∞ (1.2)
→0 h
Définition 1.1.5 (Conditionnement). On définit le conditionnement d’un problème stable
comme la meilleure constante de Lipschitz locale.
x(d + h) − x(d)
Kabs = lim sup (1.3)
→0 khk< h
On définit également le conditionnement relatif
||d||
Krel = Kabs (1.4)
||x(d)||
Exemple 1.3:
Soit A ∈ Mn,m (R), y ∈ Rn et y ∈ Rm . On considère le problème
F (x, d) = Ax − y avec d = (A, y)
Le problème mathématique (PM) a une unique solution si A est inversible, alors x0 = A−1 (y) =
x(A, y). Avec la norme k k2
1
x(d) = A−1 = (ComA)T y
det A
q
kAk2 = Tr(AAT )
et donc G est continue en (A, y). On en déduit que le problème est bien posé.
Maintenant on suppose que A est fixé. La donnée est donc d = y. On a
||x(d + h) − x(d)|| ||A−1 (y + h) − A−1 (y)||
=
khk khk
−1
||A (h)||
=
khk
h
 
−1
= A
khk
d’où
Kabs = sup A−1 (u) = |||A−1 |||
||u||=1

Si A est symétrique, on note ses valeurs propres (λ1 , ..., λn ). Sans perte de généralité on peut les
ordonner en valeur absolue |λ1 | > ... > |λn |. On obtient alors
1
|||A−1 ||| =
|λn |
Mais également
||d|| ||y|| 1 ||Ax||
Krel = Kabs = Kabs −1 =
||x|| ||A (y)|| |λn | ||x||
1 |λ1 |
Krel 6 sup kAxk =
|λn | ||x||=1 |λn |

7
1.1.4 Calcul de conditionnement lorsque G est différentiable
Définition 1.1.6 (Différentielle). G est différentiable en y si il existe A linéaire tel que

∀d ∈ Rn , G(y + h) = G(y) + A(h) + o(khk)

On note A = dG(y) la différentielle de G en y.

Proposition 1.1.7
On a le conditionnement
Kabs = |||dG||| (1.5)
et
||d||
Krel = |||dG|||
||G(d)||
Démonstration. On calcule le conditionnement
||x(d + h) − x(d)|| ||x(d) + dG(h) + o(khk) − x(d)||
=
khk khk
||dG(h)||
= + o(1)
khk
d’où
||x(d + h) − x(d)||
 
lim sup = |||dG|||
khk→0 khk
et donc également
||d|| ||d||
Krel = Kabs = |||dG|||
||x|| ||G(d)||


1.2 Convergence des méthodes numériques

1.2.1 Définitions
On suppose que le (PM) : F (x, d) = 0 est stable et bien posé.
Une méthode numérique est une suite de problème approchés, i.e. Fn (x, dn ) = 0 de solution
xn .

Définition 1.2.1. Une méthode numérique est convergente si

∀d ∈ N, n→∞
lim xn (dn ) = x(d) (1.6)
dn →d

En pratique la convergence est difficile à vérifier directement ; Il est parfois plus facile de
vérifier les propriétés de consistance et de stabilités définie ci-dessous.

Définition 1.2.2 (Consistance). Une méthode numérique est consistante si

∀d ∈ N, Fn (x(d), d) → 0 lorsque n → ∞. (1.7)

où x(d) est solution du problème F (x, d) = 0.

8
 La stabilité d’une méthode numérique est définie comme celle d’un problème mathéma-
tique à la section 1.1.3. Attention cependant au fait que la condition suivante est uniforme
en n

Définition 1.2.3 (stabilité d’une méthode numérique). Une méthode numérique Fn (x, dn )
est stable si pour tout d il existe K0 (d) et δ tels que :

∀n, sup kxn (d + h) − xn (d)k ≤ K0 (d)khk.

khk<δ

L’intérêt de telles définitions est le suivant :

Sous des hypothèses raisonnables, si une méthode est consistante et stable,
alors elle est également convergente.

Exemple 1.4 (Convergence dans le cas où Fn est « suffisamment régulière »):

On suppose que pour tout x, n, d,
∂Fn
1. ∂x (x, d) existe et est inversible.
2. absence de ‘plateau’ :
h ∂F i−1
n
M = sup (x, d) < +∞
x,d,n ∂x
(intuitivement : Fn n’est pas trop ’plate’).
3. La méthode numérique est stable et consistante.
On montre que sous ces hypothèses, la méthode est convergente.
Preuve : exercice.

1.3 Erreurs d’arrondis

(N.B. : arrondis = termes « δd » dans l’analyse de la stabilité).

1.3.1 Représentation en base β

Soit β ≥ 2 ∈ N. Tout nombre réel x s’écrit
n
X
x = (−1)s xk β k ,
k=−∞

pour un certain n ∈ N, avec s ∈ {−1, 1} (le signe) et xk ∈ [0, . . . , β − 1]. On note [x]β =
xn xn−1 . . . x1 x0 · x−1 x−2 · · · la représentation en base β.

1.3.2 Nombres à virgule flottante ("float")

Souvent β = 2 ou β = 16. On veut pouvoir représenter un panel le plus grand possible de
nombres en base β, avec N cases mémoires seulement.

Méthode naïve : utiliser directement la représentation en base β ci-dessus, tronquée à

k chigffres après la virgules : → réel max = β N −k − 1 : pas grand !

9
 Principe des flottants : faire ’flotter’ la virgule et réserver une case pour l’exposant e.

x = (−1)s × 0. a1 a2 . . . at ×β e
| {z }
taille t

•t : nombre de chiffres significatifs

•ai ∈ {0, . . . , β − 1}
•‘a1 a2 . . . a0t : mantisse.
•e : un entier : l’exposant. Contrainte sur e : L ≤ e ≤ U . (L et U sont des entiers
dépendent de la taille réservée pour e.
• On impose : a1 6= 0 → unicité de la représentation.
exemples :
• systeme 32 bits (’float’) : s7→ 1 bit ; e 7→ 8 bits, m 7→ 23 bits.
• systeme 64 bits (’double’) : s7→ 1 bit ; e 7→ 11 bits, m 7→ 52 bits.
on appelle F(β, t, L, U ) l’ensemble des nombres flottants représentables de cette manière.

! ! ! c’est un ensemble fini (6= R) ! ! !

réels min et max : tout flottant x vérifie

0.1β L ≤ |x| ≤ 0.(β − 1)(β − 1) . . . (β − 1) × β U = β U (1 − β −t )

epsilon machine
 = min{x > 0 : 1 + x ∈ F(β, t, L, U )}.
. . 0} ×β, le nombre suivant est 1 +  = 0. 1| . {z
Puisque 1 = 0. |1 .{z . . 01} ×β.
taille t taille t
D’où  = (1 + ) − 1 = 0.0 . . . 01 × β = β 1−t .

Dans R : Toutes ces constantes sont acessibles dans la liste .Machine

1.3.3 Opérations sur les flottants

On note f l(x) la représentation machine d’un réel x. Pour toute opération arithmétique
· ∈ {+, −, ×, ÷}, on note l’opération machine correspondante, i.e., son analogue dans le
systeme des nombres flottants. Plus précisément

x y = f l[f l(x) f l(y)].

Sur les machines suivant la dernière norme ISO en vigueur, IEE559 (i.e., toutes les machines),
on est assuré que
x y − x·y 
≤ ,
x·y 2
où  = epsilon machine.

10
Chapitre 2

Interpolation Polynomiale

Enjeux de ce chapitre
• Représenter à faible coût des fonctions définies sur un intervalle borné.
• Avoir des garanties sur la qualité de la représentation en terme de norme de l’erreur
• mathématiquement : étudier la stabilité, la convergence des méthodes d’interpolation
• Objectif pour les étudiants : savoir implémenter une méthode d’interpolation et
savoir évaluer l’erreur commise.

2.1 Interpolation de Lagrange

2.1.1 Existence, unicité
Soit a < b deux réels. On considère f une fonction de [a, b] dans R.
On cherche a construire un polynôme p de degré inférieur à n ∈ N tel que pour n + 1 points
d’évaluations a 6 x0 < x1 < ... < xn 6 b on ait
∀i ∈ [[0; n]], p(xi ) = f (xi )
Définition 2.1.1 (Base de Lagrange). Pour les points d’évaluations précédents, on définit les
n + 1 polynômes suivants
n
Y x − xj
∀i ∈ [[0; n]], li (x) =
x − xj
j=0,j6=i i

Soit Pn l’espace vectoriel (de dimension n + 1) des polynômes réels de degré n.

Proposition 2.1.2
Les polynômes de la base de Lagrange vérifient les points suivants :
1. ∀i ∈ [[0; n]], li ∈ Pn
2. ∀i, j ∈ [[0; n]], li (xj ) = δij
3. Ils constituent une base de Pn
Définition 2.1.3 (Interpolée de Lagrange). On définit l’interpolation de f pour les n + 1
points d’évaluations précédents comme le polynôme vérifiant
n
X
Ln f = f (xi )l(xi )
i=0

11
Ln est un opérateur linéaire, i.e. Ln (f + g) = Ln f + Ln g

Proposition 2.1.4
De la propriété précédente il vient que Ln f ∈ Pn et

∀i ∈ [[0; n]], f (xj ) = Ln f (xi )

De plus, comme les polynômes (li ) forment une base de Pn , on en déduit que Ln f existe et
est l’unique polynôme de degré au plus n a satisfaire cette propriété.

2.1.2 Erreur d’interpolation de Lagrange

Définition 2.1.5 (Polynôme nodal). Pour n + 1 points d’évaluations (x0 , ..., xn ), on définit
le polynôme nodal de degré n + 1 comme
n
Y
ωn+1 (x) = (x − xi )
i=0

Lemme 2.1.6
Si g : [a, b] → R est p fois dérivable sur [a, b] et si on dispose de x0 < ... < xn dans [a, b] tels
que ∀i ∈ [[0; n]]g(xi ) = 0, alors

∃y ∈]a, b[ tel que g (p) (y) = 0

Démonstration. Par récurrence sur n : Soit Hp la propriété de l’énoncé au rang p.

Pour H1 est vraie d’après le théorème de Rolle.
Supposons Hp−1 vraie. Soient x0 , . . . , xp comme dans l’énoncé. D’après le théorème de Rolle
appliqué à g entre xi et xi+1 , ∃(c0 , ..., cp−1 ) avec ci ∈]xi , xi+1 [ tel que g 0 (ci ) = 0. Par hypothèse
de récurrence appliquée à g 0 , il existe y ∈]c0 , cp−1 [ tel que (g 0 )(p−1) (y) = 0. Ainsi g (p) (y) = 0
et y ∈]a, b[, ce qui démontre Hp et achève la récurrence. 

Théorème 2.1.7
Si f est (n + 1) fois dérivable sur [a, b] alors

∀x ∈ [a, b], ∃yx ∈ Ix = [min(x, x0 ), max(x, xn )] ⊂ [a, b]

1
tel que En (x) = f (x) − Ln f (x) = ωn+1 (x)f (n+1) (yx )
(n + 1)!
Démonstration. On pose xn+1 = x et on note qn+1 = Ln+1 f − Ln f ∈ Pn+1 . On remarque que
qn+1 (x0 ) = ... = qn+1 (xn ) = 0 et donc qn+1 possède n + 1 racines. Comme qn+1 est dans Pn+1
on en déduit que
n
Y
qn+1 (X) = C (X − xi ) = Cωn+1 (X)
i=0

Or par définition de Ln+1 f , en X = xn+1 = x on a Ln+1 f (x) = f (x) et donc il vient que

En (x) = f (x) − Ln f (x) = Cωn+1 (x)

12
On pose désormais g = f − Ln+1 f qui s’annule en x0 , ..., xn+1 . Comme de plus g est n + 1 fois
dérivable on peut appliquer le lemme 2.1.6, d’où l’existence de yxn+1 = yx tel que g n+1 (yx ) = 0.
Or,
g (n+1) (y) = f (n+1) (y) − (Ln f − Cωn+1 )(n+1) (y)
= f (n+1) (y) − C(n + 1)!.
Pour conclure, remarquons que

f (n+1) (yx )
g (n+1) (yx ) = 0 ⇒ C = .
(n + 1)!


Corollaire 2.1.8
On déduit du théorème précédent
1
kf − Ln f k[a,b],∞ 6 kωn+1 k[a,b],∞ f (n+1)
(n + 1)! [a,b],∞

2.1.3 Contre exemple de Runge

Le membre de droite dans le corollaire 2.1.8 fait intervenir la norme du polynôme nodal
ωn+1 . Malheureusement, cette dernière « explose » pour les grandes valeurs de n, contraire-
ment à l’intuition qui suggère qu’une interpolation avec des polynômes de degré élevé fourni-
rait une meilleure approximation.

Exemple 2.1:
On pose sur [−1, 1]
1
f (x) =
1 + 5x2
Et on considère les noeuds equi - répartis xi = a + i b−a
n La figure 1. représente f et son interpolée
de Lagrange.

2.1.4 Stabilité
En pratique on évalue les f (xi ) par des valeurs approchées fb(xi ). On cherche à contrôler
||Ln fb − Ln f || par rapport à ||δf || = ||f − fb||. Ici et dans la suite, on choisit comme norme
sur les fonctions bornées [a, b] → R :

khk = khk∞,[a,b] = sup |h(x)|.

x∈[a,b]

En reprenant la définition du conditionnement on obtient

Ln (f + δf ) − Ln (f )
Kabs (f ) = sup
||δf ||→0 δf
Ln (δf )
= sup
||δf ||→0 δf
= sup kLn (δf )k
||δf ||=1

13
 Figure 2.1 – f et son interpolée pour n = 10

Et on en déduit le résultat suivant

Kabs (f ) = |||Ln |||.

(rappel : |||Ln ||| est la « norme de opérateur » de Ln , en prenant la norme k · k∞,[a,b] sur les
fonctions continues sur [a, b]. Autrement dit |||Ln ||| = supkf k∞,[a,b] =1 kLn f k.

Calcul de |||Ln ||| :

Dépend de X n = (x0 , ..., xn )
Soit g une fonction de [a, b] dans R
n
X
kLn gk = sup g(xi )li (x)
x∈[a,b] i=0
Xn
6 sup |g(xi )| |li (x)|
x∈[a,b] i=0
n
X
6 kgk[a,b] sup |li (x)|
x∈[a,b] i=0

Définition 2.1.9 (Constante de Lebesgue). Pour X n = (x0 , ..., xn ) donné on définit la

constante de Lebesgue
n
X
Λn = sup |li (x)|
x∈[a,b] i=0

Proposition 2.1.10 (Norme de l’opérateur d’interpolation de Lagrange et constante de

Lebesgue)

|||Ln ||| = Λn

14
Démonstration. On vient de montrer que |||Ln ||| 6 Λn , montrons que l’inégalité est atteinte
pour une certaine fonction g telle que kgk = 1.
Les li sont continues sur [a, b], donc ni=0 |li (.)| est également continue. On en déduit qu’il
P

existe y dans [a, b] tel que Λn = ni=0 |li (y)|.

P

On peut toujours construire une fonction g : [a, b] → R telle que

• Pour i ∈ {0, n}, g(xi ) = signe(li (y)) ∈ {−1, 1}
• kgk = 1.
Pour ce choix de g, on remarque enfin que
n
X n
X n
X
kLn gk = sup g(xi )li (x) > g(xi )li (y) = |li (y)| = Λn ,
x∈[a,b] i=0 i=0 i=0
d’où
Λn ≤ kLn gk ≤ |||Ln |||.


2.1.5 Lien entre stabilité et convergence

Théorème 2.1.11 (Erreur et constante de Lebesque)
Soit f : [a, b] → R et X n = (x0 , ..., xn ). On a
kf − Ln f k 6 En∗ (f )(1 + Λn )
où
En∗ (f ) = inf kf − P k
P ∈Pn

Démonstration. Soit f une fonction continue sur [a, b]. On pose

p∗n = arg min kf − P k
P ∈Pn

(L’existence de p∗npeut se montrer par le fait que K = {p ∈ Pn : kf − pk∞,[a,b] ≤ M } est

fermé, borné, donc compact dans Pn ).
Alors,
kf − Ln f k 6 kf − p∗n k + kp∗n − Ln f k
6 En∗ (f ) + kLn p∗n − Ln f k
6 En∗ (f ) + Λn kp∗n − f k
6 En∗ (f )(1 + Λn )


2.1.6 Ordre de grandeur de la constantes de contrôle de l’erreur pour n

grand
Proposition 2.1.12
Dans le cas de n + 1 points d’évaluations équidistants sur [a, b], i.e. xi = a + i b−a
n pour
i ∈ [[0; n]], on a les équivalents suivants
2n+1
Λequi
n ∼
e.n. ln(n)
b − a n+1
 
equi
ωn+1 ∼2
e

15
2.2 Interpolation aux noeuds de Tchebychev
Dans cette section une alternative aux noeuds équidistants, pour laquelle les constantes
de contrôle de l’erreur Λn et kωkn sont meilleures. Plus précisément, on va choisir les noeuds
d’interpolation (sur l’intervalle [−1, 1] comme étant les racines du n + 1eme polynôme de
Tchebychev. Les propriétés de bases de ces polynômes sont rappelées ci-dessous.

2.2.1 Définition et premières propriétés

Définition 2.2.1 (Polynômes de Tchebychev). On définit les polynômes de Tchebychev pour
x ∈ [−1, 1]
tn : x 7→ cos(n arccos x)

Proposition 2.2.2
Les (tn ) sont des polynômes définis par la relation de récurrence


 t0 (x) = 1
t1 (x) = x

t
n+1 (x) + tn−1 (x) = 2xtn (x), n > 1

Démonstration. Si le système précédent est vrai, la relation de récurrence montre que les (tn )
sont des polynômes. Montrons que le système est vrai :
On pose θ = arccos x pour x dans [−1, 1]. On a par trigonométrie

tn+1 (x) + tn−1 (x) = cos((n + 1)θ) + cos((n − 1)θ)

= 2 cos(nθ) cos(θ)
= 2xtn (x)

Corollaire 2.2.3
tn est de degré n et de coefficient dominant 2n−1 .

Proposition 2.2.4
tn a n racines distinctes dans [−1, 1]

(k + 21 )π
!
∀k ∈ [[0; n − 1]], xk = cos
n

Démonstration. Soit x dans [−1, 1]. On dispose d’un unique θ dans [0, π] tel que x = cos(θ).

tn (x) = 0 ⇔ cos(nθ) = 0
π
⇔ nθ ≡ [π]
2 
π π
⇔θ≡
2n n
d’où le résultat. 

16
 Figure 2.2 – Noeuds d’évaluation de Tchebychev pour n = 8

Définition 2.2.5 (Noeuds de Tchebychev). Les points d’interpolation de Tchebychev d’ordre

n − 1 sont les n racines de tn : {x0 , ..., xn−1 }. Il s’agit de l’ordre n − 1 car le polynôme de
Lagrange est de degré n − 1 pour n points d’interpolations.

Remarque 2.2.6. Étant donné que les (xi ) sont les racines de tn ,et que l’on connaît le
coefficient dominant de tn , on peut écrire

tn = 2n−1 ωn

2.2.2 Interpolation de Tchebychev sur [a, b]

Pour se ramener à [-1,1], on considère le changement de variable affine

ϕ[a,b] : [−1, 1] → [a, b]

a+b

b−a

(2.1)
u 7→ x = + u
2 2
On vérifie que ϕ[a,b] est bien un transformation affine entre [−1, 1] et [a, b].
On note (uk ) les racines de tn sur [−1, 1] et on pose pour k dans [[0; n]]

xk = ϕ[a,b] (uk )

[a,b]
On note par la suite ωn,T le polynôme nodal de Tchebychev pour n points d’interpolation
sur [a, b]. Rappelons que la norme du polynôme nodal est cruciale pour contrôler l’erreur
d’interpolation (c.f. Corollaire 2.1.8). Contrairement au cas des noeuds équidistants, on peut
calculer de manière simple et explicite la valeur de kωkn aux noeuds de Tchebychev :

Proposition 2.2.7
La norme du polynôme nodal aux n noeuds de Tchebychev est donnée par
n
b−a

[a,b]
ωn,T =2
4

17
Démonstration. Soit u dans [−1, 1]
n−1
Y b−a

[a,b]
ωn,T ◦ ϕ[a,b] (u) = (u − uk )
k=0
2
n n−1
b−a
 Y
= (u − uk )
2 k=0
b − a n [−1,1]
 
= ωn,T (u)
2
b−a n
 
=2 tn (u)
4

or, ktn k = kcos(n arccos)k = 1, d’où le résultat. 

Corollaire 2.2.8
On déduit du théorème précédent et de la proposition 2.1.12 :
[a,b]
ωn,T  n
e
[a,b]
∼ ( 1 lorsque n → ∞).
ωn,equi 4

La deuxième constante intervenant dans le contrôle de la stabilité et de l’erreur est la

constante de Lebesgue Λn (c.f. Proposition 2.1.10 et Théorème 2.1.11). Le résultat suivant
(admis) est à comparer avec l’ordre de grandeur de Λn pour des noeuds équidistants donnée
par la proposition 2.1.12.

Proposition 2.2.9 (Constante de Lebesgue aux noeuds de Tchebychev)

Un équivalent pour n grand de la constante de Lebesgue aux noeuds e Tchebychev sur [a, b]
est donné par
2
ΛTcheby
n ∼n→∞ ln(n)
π
2n+1
Ainsi ΛTcheby
n  Λequi
n ∼ e.n. ln(n)

2.3 Différences divisées, Base de Newton

On cherche dans cette section un moyen efficace en temps de calcul pour calculer nu-
mériquement l’interpolateur de Lagrange de f . En effet, l’évaluation ’naïve’ à partir de la
définition O(n2 ) opérations élémentaires (multiplications par (x − xi ) ou (xi − xj )) (exercice).
La procédé d’évaluation ci-dessous est en O(n) opérations élémentaires.

2.3.1 Base de Newton

Définition 2.3.1 (Base de Newton). Soit X n = (x1 , ..., xn ) dans [a, b]. On définit les poly-
nômes suivants
k−1
Y
ωk = (X − xi ) ∈ Pk
i=0

Remarque 2.3.2. Il s’agit des polynômes nodaux d’ordre k.

18
Proposition 2.3.3
La famille B = {ω0 , ..., ωn } est une base de Pn .
But :exprimer Ln f dans cette base, autrement dit trouver des coefficients a0 , . . . , an tels
que Ln f = ni=0 ai ωi .
P

2.3.2 Différences divisées

Soit k ∈ [[1; n]]. On considère le polynôme de Pk

∆k = Lk f − Lk−1 f

On remarque que l’on en connaît k racines puisque

∆k (x0 ) = ... = ∆k (xk−1 ) = 0.

On en déduit que, pour tout k ∈ [[1; n]], il existe ak ∈ R tel que

∆k (X) = ak ωk (X).

Ainsi
∀k ∈ [[1; n]], Lk f = ak ωk + Lk−1 f,
Ainsi, il existe bien (a0 , . . . , an ) tels que
n
X k
X
Ln f = aj ωj , et ∀k ≤ n, Lk f = aj ωj (2.2)
j=0 j=0

Intérêt : si on est capable de calculer les aj pour j ≤ n − 1 (pour n noeuds), et si on rajoute

un n + 1eme noeud, on n’a plus qu’à calculer le dernier coefficient an+1 (on ne doit pas tout
recalculer).
Définition 2.3.4 (Différences divisées). Les (ak ) sont appelés différences divisées de f . On
note
ak = f[x0 ,...,xk ]
Remarque 2.3.5. On a l’égalité suivante

a0 = f[x0 ] = f (x0 )

Proposition 2.3.6
On peut calculer les coefficients récursivement
f[x0 ,...,xk−1 ] − f[x1 ,...,xk ]
f[x0 ,...,xk ] =
x0 − xk
et plus généralement
f[xi ,...,xj−1 ] − f[xi+1 ,...,xj ]
∀i < j, f[xi ,...,xj ] = (2.3)
xi − xj
Démonstration. exercice de TD 

Intérêt de cette proposition : suggère une méthode itérative pour calculer les f[x0 ,...,xk ] ,
à partir d’une table triangulaire.

19
2.3.3 Méthode de Horner pour l’évaluation
c.f. TD : La méthode de Horner permet d’évaluer Ln f (x) en O(n) opérations, étant don-
nées les différences divisées f[x0 ,...,xk ] , k = 0, . . . , n.

20
Chapitre 3

Intégration numérique : méthodes

de quadrature

3.1 Introduction : méthodes de quadrature, méthodes de Newton-

Cotes
But : Approcher des intégrales de type
Z b
I(f ) = f (x) dx (a < b ∈ R)
a

Pratique courante : changement de variable pour se ramener à une intégration sur [−1, 1] :
m = (a + b)/2, h = (b − a)/2, g(u) = f (m + hu), u ∈ [−1, 1].
Z 1
h
I(f ) = g(u) du
2 −1

Les méthodes de quadratures définies ci-dessous sont des opérateurs Ib : f 7→ I(f

b ) destinés à
approcher I(f ).
Définition 3.1.1 (Méthodes de quadratures (MQ)).
1. Une MQ simple(MQS) sur [a, b] est un opérateur
n
X
f 7→ I(f
b ) = (b − a) λi f (xi ), (3.1)
i=0

avec λi > 0, ni=0 λi = 1, xi ∈ [a, b], i = 0, . . . , n.

P

2. Soit α0 = a < α1 < . . . αM = b une subdivision de [a, b]. Une MQ Composite (MQC)
est un opérateur de type
M
X −1 nm
X
f 7→ I(f
b )= (αm+1 − αm ) λm,i f (xm,i ), (3.2)
m=0 i=0

avec pour tout m ∈ {0, . . . , M − 1} :

λm,i > 0, ni=0 λm,i = 1, xm,i ∈ [αm , αm+1 ],
P m
i = 0, . . . , nm .
De manière générale on étudie les propriétés des MQS et l’on en déduit celles de MQC.

21
 Lien avec l’interpolation de Lagrange, rang d’une MQ Les MQS que nous verrons
seront toujours de type
I(f
b ) = I(Ln f )

où Ln f est le polynôme d’interpolation de Lagrange aux noeuds x0 , . . . , xn . Une telle méthode

est dite de rang n. Le membre de droite est bien de type (3.1). En effet
n
Z bX
I(Ln f ) = f (xi )li (x) dx (li : ieme polynôme de Lagrange)
a i=0
n
X Z b
= f (xi ) li (x) dx
i=0 |a {z }
λ̃i
n
X
= (b − a) λi f (xi ) ,
i=0

avec λi = λ̃i /(b − a) > 0. Reste à vérifier que = (b − a). Or

P
i λ̃i

X n
Z bX
λ̃i = li (x) dx,
i a i=0

Le polynôme ni=0 li est de degré au plus n et vaut 1 en n + 1 points. C’est donc le polynôme
P

constant égal à 1, d’où le résultat.

Cas particulier : méthodes de Newton-Cotes C’est le cas particulier où les points

d’interpolation sont équidistants. L’avantage principal est celui de la simplicité d’implémen-
tation

Définition 3.1.2 (Méthode de Newton-Cotes de rang n (N Cn )). LA M QS de Newton-

Cotes de rang n sur [a, b] est la méthode données par I(f
b ) = I(Ln f ) où Ln f est le polynôme
d’interpolation de Lagrange de f aux noeuds équidistants
b−a
(x0 = a, . . . , xi = a + i , . . . , xn = b).
n

Contrôle de l’erreur : ordre d’une MQ

Définition 3.1.3 (Ordre d’une MQ). Une MQS ( resp. une MQC) est dite d’ordre N si elle
est
1. Exacte pour tout les polynômes de degré ≤ N , i.e., si I(p)
b = I(p) pour tout p ∈ PN .
( resp. pour toute fonction polynomiale par morceaux sur la subdivision (αm )m=0:M )), et
2. Inexacte pour au moins un polynôme p ∈ PN +1 ( resp. une fonction polynomiale par
morceaux de degré n + 1).

Remarque 3.1.4 (Rang et ordre). Attention à ne pas confondre le rang d’une méthode de
quadrature(= le degré n du polynôme de Lagrange sur lequel est basée l’approximation) et son
ordre (le plus grand entier N pour lequel la méthode est exacte pour les polynômes de degré
N . On verra plus loin que de manière générale N ≥ n mais l’inégalité peut être stricte.

22
 Par linéarité de Iˆ et de I, et par changement de variable, on a le résultat suivant : notons
X0 = 1, X, X 2 , . . . , X N , . . . les polynômes de la base canonique de PN . Une MQS Ib est d’ordre
N sur [a, b] si et seulement si
b i ) = I(X i ) sur (−1, 1] pour i = 1, . . . , N et
1. I(X
b N +1 ) 6= I(X N +1 ) sur [−1, 1].
2. I(X

3.2 Exemples
3.2.1 Rectangles à gauche / à droite
1. Rectangles à Gauche :
• MQS : I(f
b ) = (b − a)f (a)
b ) = PM (αm+1 − αm )f (αm ).
• MQC : I(f m=0
2. Rectangles à droite : Idem en remplaçant f (a) par f (b) et f (αm ) par f (αm+1 ).

Ordre de la MQS des rectangles (à gauche ou à droite)

On vérifie pour [a, b] = [−1, 1], que I(1)
b = I(1) = b − a, mais que I(X)
b 6= I(X). La méthode
est donc d’ordre 0.

3.2.2 Méthode du point Milieu

a+b
On pose x0 = 2 et on prend x0 comme unique point d’interpolation. On obtient

b ) = (b − a)f a+b
I(f ).
2
On montre comme au paragraphe précédent que la méthode est d’ordre 1.

3.2.3 Méthode des trapèzes

C’est le nom donné à la méthode de Newton-Cotes de rang 1 (N C1 ). Autrement dit on
prend deux points d’interpolation x0 = a, x1 = b, de sorte que

x−b x−a (x − a)f (b) + (b − x)f (a)

L1 f (x) = f (a) + f (b) =
a−b b−a b−a
(L1 f est la la droite passant par (a, f (a)) et (b, f (b))). Ainsi
Z b
f (a) + f (b)
I(f
b )= L1 f (x) dx = . . .(calcul). . . = (b − a)
a 2
(c’est l’aire du trapèze sous le segment de droite entre les deux points d’interpolation).

Une formule simple pour la méthode composite associée

c.f. TP4.

Ordre de la méthode des trapèzes On vérifie que la méthode est d’ordre 1.

23
3.2.4 Méthode de Cavalieri-Simpson
C’est le nom donné à la méthode de Newton-Cotes de rang 2 (N C2 ) : on intègre le
polynôme de Lagrange aux points x0 = a, x1 = (a + b)/2, x2 = b. On montre que les poids
de la méthode sont
1 2 1
λ0 = , λ 1 = , λ 2 = .
6 3 6
(exercice : vérifiez-le par intégration des polynômes de Lagrange li , i = 0, 1, 2).

3.3 Ordre des méthodes de Newton-Cotes

On généralise les calculs d’ordre des méthodes données en exemple à la section précédente.
Soit Ib une méthode N Cn . Par définition I(f
b ) = I(Ln f ) avec Ln f le polynôme de Lagrange
interpolé en n + 1 points équidistants. Lorsque f est un polynôme de degré ≤ n, on a Ln f = f
ˆ ) = I(f ), la méthode N Cn est exacte. l’ordre d’une
(l(interpolation est exacte), d’où I(f
méthode N Cn est donc toujours supérieur ou égal à n. Lorsque n est pair, on gagne un ordre
de précision, comme résumé dans la proposition ci-dessous.

Proposition 3.3.1
Pour tout n ∈ N, l’ordre N de la méthode N Cn est
• N = n si n est impair
• N = n + 1 si n est pair.

Démonstration. on a montré que l’ordre est ≥ n.

• Pour n impair, on admet qu’il y a égalité.
• Pour n pair, on va montrer que l’ordre est supérieur ou égal à n + 1, et on admettra
l’égalité.
Soit n ∈ N, pair. Par linéarité / changement de variable, il suffit de montrer que sur l’intervalle
b ) = I(f ) avec f (x) = xn+1 . Par imparité, I(f ) = 0. D’autre part, par définition
[−1, 1], I(f
de I(f ), on a
b
n Z 1
xn+1
X
I(f
b )=
i li (x) dx.
i=0 −1

De plus, on a, par symétrie des xi et des li ,

n
xn/2 = 0 , et ∀i ∈ 1, . . . − 1, xn−i = −xi .
2
enfin pour tout i ∈ {0, . . . , n}, ln−i (x) = li (−x) (∀x ∈ [−1, 1]), d’où
Z Z
li = ln−i .

Ainsi
n
2
−1 Z Z
xn+1
X
I(f
b )=
i li + (−xi )n+1 li
i=0
=0

Ainsi I(f ) − I(f

b ) = 0 : la méthode est donc bien d’ordre ≥ n + 1. 

24
3.4 Stabilité
A cause des erreurs de mesure/ d’arrondi, on n’a jamais exactement accès aux valeurs
f (xi ) mais plutôt à des valeurs perturbées f˜(xi ) = f (xi ) + δf (xi ), où δf est la ’perturbation’
de la fonction f . Par linéarité de l’opérateur I,
b ceci induit une perturbation du résultat de la
méthode :
b f˜) − I(f
I( b f˜ − f ) = I(δf
b ) = I( b ).

On pose g = δf . On considère des fonctions (et des perturbations) bornées sur [a, b] et on
utilise la norme infinie sur [a, b], kf k = kf k∞,[a,b] . L’analyse de la stabilité (cf. Chapitre 1)
consiste alors à trouver une constante K (si elle existe) telle que pour toute fonction g bornée,

I(g)
b ≤ Kkgk.

Plus précisément on cherche à donner la ‘meilleure’ constante K telle que la condition ci-dessus
soit vérifiée, c’est-à-dire le conditionnement absolu Kabs = supkgk=1 |I(g)|
b

3.4.1 Stabilité d’une MQS

Par définition des méthodes de quadratures simples sur [a, b], on peut écrire
n
X Z b
I(g)
b = g(xi ) li (x) dx
i=0 a
Z b n
X
≤ kgk∞,[a,b] |li (x)| dx
a 0
| {z P
}
≤Λn =supu∈[a,b] |li (u)|

≤ kgk∞,[a,b] Λn (b − a),

où Λn est la constante de Lebesgue (norme de l’opérateur d’interpolation de lagrange, c.f.

Définition 2.1.9). La méthode est donc stable et le conditionnement absolu vérifie Kabs ≤
Λn (b − a). On peut montrer (comme pour le calcul du conditionnement de l’interpolateur de
Lagrange, Proposition 2.1.10) qu’il y a égalité. On retiendra

Proposition 3.4.1
Le conditionnement absolu d’une MQS sur [a; b] basée sur l’interpolation aux noeuds (x0 , . . . , xn )
est
Kabs = Λn (b − a)
où Λn est la constante de Lebesgue pour les noeuds (x0 , . . . , xn ).

La constante Λn croissant rapidement avec n, on voit que le conditionnement se détériore

pour des MQS de rang n élevé.

3.4.2 Stabilité d’une MQC

Soit IbM Une méthode composite avec M sous intervalles αm , αm+1 , m = 0, . . . , M − 1, sur
chacun desquels on applique une méthode simple de rang n notée Iˆ[α,m,αm+1 ] . Par définition

25
 −1 b
on a IbM = M
P
i=0 I[α,m,αm+1 . On en déduit la majoration suivante, pour toute fonction g
bornée sur [a, b] :
M
X −1
|IbM (g)| = Ib[α,m,αm+1 ] (g)
i=0
M
X −1
≤ Ib[α,m,αm+1 ] (g)
i=0
M
X −1
≤ Λn (αm − αm+1 )kgk∞,[αm ,αm+1 ]
i=0
≤ λn kgk∞ (b − a).
Comme au paragraphe précédent, il est facile de construire une fonction g telle qu’il y ait
égalité. On obtient ainsi le même résultat que pour une MQS.

Proposition 3.4.2
Le conditionnement absolu d’une MQC sur [a; b] construite avec M sous-intervalles et une
méthode simple de rang n sur chaque sous-intervalle est le même que le conditionnement
d’une MQS de rang n, c’est-à-dire

Kabs = Λn (b − a)

où Λn est la constante de Lebesgue de la MQS de rang n choisie.

Conclusion : il est beaucoup plus intéressant d’utiliser des MQC avec un grand nombre
M de sous-intervalles et un rang n modéré, plutôt qu’une méthode simple de rang O(nM ),
car la stabilité ne dépend que du rang de la méthode simple et est meilleure pour des rangs
faibles.

3.5 Majorations de l’erreur

3.5.1 Erreur d’une MQS d’ordre N
On définit l’erreur d’une méthode de quadrature simple d’ordre N (donc de rang n ≤ N )
appliquée à une fonction f sur l’intervalle [a, b], par
(a,b) b ) − I(f ).
EIN (f ) := I(f

Soit Ib une méthode d’ordre N . On suppose que f ∈ C (N +1) [a, b]. La formule de Taylor avec
reste exacte donne, pour x ∈ [a, b],

f (N +1) (θx )
f (x) = pN (x) + (x − a)N +1 , (3.3)
(N + 1)!
| {z }
gN (x)

où θx ∈ [a, x] et pN ∈ PN . Comme I et Ib sont des opérateurs linéaires, l’erreur EI l’est aussi.

Ainsi
(a,b)
EIN (f ) = EI(pN ) + EI(gN ) = EI(gN ) ,

26
puisque Ib est exacte sur PN . On obtient la majoration suivante de l’erreur
(a,b)
EIN (f ) = I(gN ) − I(g
b N ) ≤ |I(gN )| + IbN (gN ) (3.4)

D’après (3.3),

kf (N +1) k (b − a)N +2 kf (N +1) k

|I(gN )| ≤ = (b − a)N +2 .
(N + 1)! N + 2 (N + 2)!

D’autre part
n
X
b N ) = (b − a)
I(g gN (xi )λi
0
≤ (b − a)kgN k
kf (N +1) k(b − a)N +1
≤ (b − a)
(N + 1)!
kf (N +1) k
≤ (b − a)N +2
(N + 1)!

Au vu de (3.4) et des deux dernières majorations, on a finalement

(a,b)
[EIN (f )| ≤ CN kf (N +1) k(b − a)N +1 , (3.5)
avec CN = 1/(N + 1)! + 1/(N + 2)!.
Cette borne ne permet pas d’assurer que l’erreur tende vers 0 si la fonction n’est pas très
‘régulière’, au sens où la norme des dérivées d’ordre élevé ‘explose’, ou si (b − a) > 1. Ceci
suggère encore une fois d’utiliser des méthodes composites avec de petits sous-intervalles et
un ordre modéré.

3.5.2 Erreur d’une MQC d’ordre N sur des intervalles de même longueur
(a,b)
Soit EIN,M (f ) l’erreur d’une méthode composite sur M sous-intervalles [αm , αm+1 ] sur
lesquels on applique une méthode d’ordre N , Ib[αm ,αm+1 ] . On considère des αi équidistants
αm = a + m b−a b−a
M et on note h = M .
N.B. On ne suppose pas que la MQS utilisée sur chaque sous-intervalle utilise des noeuds
équirépartis.
L’erreur s’écrit alors
M −1 Z b
(a,b) X
EIN,M (f ) = Ib[αm ,αm+1 ] − f
m=0 a
M
X −1 Z αm+1
= Ib[αm ,αm+1 ] − f
m=0 αm
M
X −1
= EIn(αm ,αm+1 (f )
m=0

27
D’où
M −1
(a,b) X
EIN,M (f ) ≤ EIn(αm ,αm+1 (f )
m=0
M
X −1
≤ CN kf (N +1) k∞,[αm ,αm+1 ] (αm+1 − αm )N +2
m=0
 b − a N +2
≤ CN kf (N +1) k∞,[a,b] M
|M
{z }
h
=≤ CN kf (N +1)
k∞,[a,b] (b − a) hN +1 (3.6)
h iN +1
On gagne donc un facteur h/(b − a) = M −(N +1) par rapport à la borne d’erreur (3.5)
pour une MQS. En utilisant une MQS d’ordre N modéré (donc de faible rang et ainsi avec
une bonne stabilité), il est toujours possible d’augmenter M (diminuer h), et d’avoir un bon
contrôle de l’erreur.

3.5.3 Erreur des méthodes de Newton-Cotes

On conclut cette partie sur un résultat (admis) sur l’erreur des méthodes de Newton-Cotes
([Link] MQS avec noeuds équidistants et les MQC associées), pour lesquelles les inégalités (3.5)
et (3.6) sont des égalités.
Proposition 3.5.1
Il existe une suite de constantes (KN )n∈N > 0 indépendantes de l’intervalle d’intégration (a, b)
et de la fonction f à intégrer telles que pour tout n ∈ N, pour toute fonction f ∈ C (N +1) [a, b] :
1. L’erreur de la MQS Newton Cotes d’ordre N , est
(a,b)
EIN (f ) = KN f (N +1) (θN ) (b − a)N +2 ,
où θN ∈ [a, b],
b−a
2. La MQC de Newton Cotes d’ordre N et de pas h = M , il existe θN,M ∈ [a, b], tel que
(a,b)
EIN,M (f ) = KN f (N +1) (θN,M ) (b − a) hN +1

3.6 Évaluations de l’erreur a posteriori

L’inconvénient de l’analyse de l’erreur dans la section précédente est d’impliquer des quan-
tités le plus souvent inconnues (les dérivées de f , évaluées en des points θN inconnus) Une
autre manière dite ’a posteriori’ d’évaluation de l’erreur des méthodes composites repose
sur la comparaison entre les résultats de la méthode de pas h = (b − a)/M et celle de pas
h/2 = (b − a)/(2M ). En effet la proposition 3.5.1 2. implique
EIN,M (f ) f (N +1) (θN,M ) hN +1
= (N +1)
EIN,2M (f ) f (θN,2M ) (h/2)N +1
f (N +1) (θN,M ) N +1
= 2
f (N +1) (θN,2M )
≈ 2N +1

28
En supposant que f (n+1) (θN,2M ) ≈ f (n+1) (θN,M ). En notant IM (resp. I2M ) la MQC avec
M (resp. 2M ) sous-intervalles, la dernière expression donne une relation simple entre I, IM
et I2M . Une inversion simple permet d’exprimer l’erreur I2M − I en fonction de la différence
I2M − IM .
Exemple : Analyse a posteriori pour la méthode de Cavalieri Simpson (N C2 , ordre N = 3) :
c.f. TP.

3.7 Extrapolation de Richardson, application aux trapèzes (mé-

thode de Romberg)
3.7.1 Principe de la méthode de Richardson
L’extrapolation de Richardson est une méthode générale pour estimer efficacement la
valeur en 0 d’une fonction A : R+ → R, à partir d’évaluations en certains points bien choisis
ti > 0, i = 0, . . . , n. Dans le cas de l’intégration numérique, on prend A(h) = I(h)
b où Ib est
une MQC de pas h = (b − a)/M . On cherche alors à évaluer I(0) := limh→0 I(h) = I (la
b b
vraie intégrale), où la dernière égalité vient du fait que l’erreur tend vers 0 avec h, d’après la
majoration (3.6) de l’erreur.
Les points ti sont choisis de la manière suivante (pour des raisons qui apparaîtront plus
loin) : on fixe t > 0 et on choisit δ ∈]0, 1[, puis on pose
ti = δ i t, i = 0, . . . , n

Un estimateur ‘naïf’ est Abnaif (0) = A(δ n t). Un développement de Taylor d’ordre 1 donne
A(0) = A(δ n t) + O(δ n t). L’erreur est en O(δ n t).

L’approximation de Richardson permet de diminuer l’ordre de grandeur de l’erreur en

interpolant un polynôme de Lagrange aux noeuds xi = ti = δ i t. Soit Ltn A le polynôme de
Lagrange obtenu. On prend alors comme estimateur, la valeur de l’interpolateur en 0, soit
Abrich (0) = Ltn A(0).

Analyse de l’erreur de la méthode de Richardson

Dans la suite on appelle An (t) = Ltn A(0) l’estimateur de Richardson décrit ci-dessus. La
formule de l’erreur de l’interpolation de Lagrange (Théorème 2.1.7) sur l’intervalle [a, b] = [0, t]
donne
A(n+1) (θn )
An (t) = Ltn A(0) = A(0) + ωn+1 (0),
(n + 1)!
où θn ∈ [0, t] et ωn+1 est le polynôme nodal aux noeuds δ i t, i = 0, . . . , n. Enfin,
n Pn n(n+1)
i n+1
Y
ωn+1 (0) = (0 − δ i t) = (−1)n+1 δ 0 t = (−1)n+1 δ 2 tn+1
i=0

L’erreur de la méthode de Richardson est donc d’ordre de grandeur

n(n+1)
A(0) − An (t) = O(δ 2 tn+1 )
n(n−1)
le rapport entre l’erreur de Richardson et l’erreur de l’estimateur naïf est en O(δ 2 tn ), ce
qui rend très attractive la méthode de Richardson pour δ < 1.

29
Interprétation de la méthode comme une ‘élimination’ des termes dans le déve-
loppement de A en 0
Supposons que A admette un développement limité en 0,

A(t) = a0 + a1 t + a2 t2 + . . . + O(tk+1 ).

avec a0 = A(0) et k ≥ n.
Puisque l’estimateur de Richardson s’écrit comme une combinaison linéaire des A(δ i t), on
a aussi un développement, à n fixé, de type
(n) (n) (n)
An (t) = b0 + b1 t + b2 t2 + . . . + O(tk+1 )

Cependant l’argument du paragraphe précédent montre que

An (t) = A(0) + O(tn+1 ) = a0 + O(tn+1 ).

t→0

Par unicité on a
(n) (n) (n)
b0 = a0 = A(0) et b1 = b2 = · · · = bn(n) = 0.

Il est donc préférable d’utiliser An (t) plutôt que A(t) (où même A(δ n t)) pour approcher a0 ,
puisque le premier terme non nul dans le développement de An (t) − A(0) est d’ordre tn+1
(plutôt que t).

Implémentation : calcul de An (t)

On peut utiliser une méthode voisine de celle des différences divisées pour calculer rapidement
An (0), étant donnés les A(δ i t), c.f. TP.

3.7.2 Méthode de Romberg

C’est l’application de√la méthode de Richardson lorsque A(t) est le résultat de la MQC
des trapèzes, de pas h = t, en prenant δ = 1/4. (la raison de ce ‘changement de variable’ et
du choix de δ est expliquée plus bas).
Rappelons que la méthode des trapèzes est une Newton-Cotes de rang 1, donc d’ordre 1.
L’erreur de la MQC de pas h est en O(h2 ) d’après la Proposition 3.5.1. On note I(h)
b cette
méthode. On utilise le résultat plus fort suivant (admis)
Proposition 3.7.1 (Formule d’Euler Maclaurin)
La MQC des trapèzes de pas h sur a, b, appliquée à une fonction 2k + 2 fois dérivable, admet
un développement en 0 :
Z b
I(h)
b = f + a1 h2 + a2 h4 + . . . + ak h2k + ak+1 f 2k+2 (θ)
a

où θ ∈ [a, b] et où les ai dépendent des dérivées d’ordre 2i − 1 de f en a et b.

Ainsi, seuls les termes
√ pairs interviennent dans le développement de I(h)
b en 0. En posant
2
t = h et A(t) = I( t) = I(h), on a donc
b b

Z b k
X
A(t) = f+ ai ti + O(tk+1 ).
a i=1
| {z }
A(0)

30
Appliquer la méthode de Richardson à A(t) plutôt qu’à I(h) permet donc d’éliminer les
n premiers termes non nuls du développement, alors qu’on éliminerait seulement les n/2
premiers en appliquant Richardson à la fonction I(h). b
Concernant le choix de δ : en prenant δ = 1/4, on a δ i t = (h/2i )2 et on doit évaluer
successivement I(h),
b I(h/2),
b . . . , I(h/2
b i ), . . .. Or, il est peu coûteux d’évaluer I(h/2)
b lorsque
l’on connaît la valeur de I(h)
b (voir TP 2-3).

31
Chapitre 4

Équations différentielles

4.1 Existence et unicité des solutions

4.1.1 Problème de Cauchy
Le problème de Cauchy, que l’on cherche à résoudre, est une équation sur la fonction
y : I → Rm :

y(t0 ) = y0
y 0 (t) = f (t, y(t)), ∀t ∈ I.

Les données du problème sont :

• l’intervalle de définition I ⊂ R,
• la fonction f : R × Rm → Rm qui est continue en y et intégrable en t.
• et la condition initiale (t0 , y0 ) où t0 ∈ ˚
I.

4.1.2 Exemple : les équations du mouvement

q est la position
p = mq 0 est m fois la vitesse
2
H(q, p) = V (q) + kpk
2m est l’énergie totale
Si il y a conservation de l’énergie :
1 ∂H
q0 = p= (q, p)
m ∂p
∂H
p0 = −∇V (q) = − (q, p)
∂q
∂H
" # " #
q (y)
Avec les notations précédentes, on a y = ∈ R et f (t, y) = f (y) = ∂p
6 .
p − ∂H
∂q (y)
Résoudre le problème de Cauchy revient à simuler le mouvement du système en partant
des conditions initiales y0 en t0 .

Remarque : On aurait pu écrire q 00 = − ∇Vm(q) . La méthode pour passer d’une équation

différentielle d’ordre quelconque à une équation différentielle d’ordre 1 est générale.

32
4.1.3 Théorèmes d’existence et unicité des solution
Régularité des solutions. ∀y solution, si f est de classe C k , alors y est de classe C k+1 .
En particulier, y est dérivable.

Lemme 4.1.1
La fonction y : I → Rm est solution du problème de Cauchy si et seulement si
Z t
∀t ∈ I, y(t) = y0 + f (u, y(u))du (équation intégrale)
t0

Démonstration. ⇒ : On intègre. ⇐ : On dérive. 

Outil de preuve. Soit C(I, Rm ) l’ensemble des fonctions continues de l’intervalle I (qui
contient t0 ) vers l’ensemble Rm . Soit φ la fonction de C(I, Rm ) → C(I, Rm ) définie par
Z x
m
∀z ∈ C(I, R ) : ∀x ∈ I, φ[z](x) = y0 + f (τ, z(τ ))dτ (4.1)
t0

Lemme 4.1.2
Soit r > 0 et soit φ définie dans (4.1). Il existe T > 0 tel que
si kz − y0 kT = supt∈[t0 −T,t0 +T ] kz(t) − y0 k ≤ r, alors kφ[z] − y0 kT = supt∈[t0 −T,t0 +T ] kφ[z](t) −
y0 k ≤ r.

Démonstration.
Z t Z t
kφ[z](t) − y0 k = k f (τ, z(τ ))dτ k ≤ | kf (τ, z(τ ))kdτ |
t0 t0

Fixons T0 > 0. Si τ ∈ [t0 − T0 ; t0 + T0 ] et kz(τ ) − y0 k ≤ r, alors (τ, z(τ )) ∈ [t0 − T0 ; t0 +

T0 ] × B(y0 , r).
[t0 − T0 ; t0 + T0 ] × B(y0 , r) est un ensemble compact et f est continue donc ∃M > 0 tel que
kf (τ, z(τ ))k ≤ M, ∀τ ∈ [t0 −R T0 ; t0 + T0 ].
Ainsi, kφ[z](t) − y0 k ≤ | tt0 kf (τ, z(τ ))kdτ | ≤ |t − t0 |M
r
et il suffit de prendre T = min(T0 , M ) pour avoir le résultat. 

Soit F = C([t0 − T ; t0 + T ] × B(y0 , r)). La restriction de φ à F est telle que φ|F : F → F.

On va maintenant écrire φ = φ|F .

Théorème 4.1.3 (Cauchy-Lipschitz : existence et unicité des solutions locales)

Si f : I × Rm → Rm est localement lipschitzienne en y alors ∃T > 0 tel que le problème de
Cauchy avec condition initiale (t0 , y0 ) admet une unique solution y : [t0 − T ; t0 + T ] → Rm .

Démonstration. On va démontrer que φ est une contraction et utiliser le théorème du point

fixe de Picard. Rappelons que l’espace des fonctions continues muni de la norme sup est un
espace complet et donc que le théorème du point fixe s’applique.
Soit (z1 , z2 ) ∈ F × F et soit k > 0 tel que f est k-lipschitzienne en y sur F. Comme F
est compact et f est localement lipschitzienne, ce k existe forcément.

33
kφ[z1 ] − φ[z2 ]k = sup kφ[z1 ](t) − φ[z2 ](t)k (déf. kφ[z]k)
t∈[t0 −T ;t0 +T ]
Z t Z t
= sup ky0 + f (τ, z1 (τ ))dτ − y0 − f (τ, z2 (τ ))dτ k (déf. φ)
t∈[t0 −T ;t0 +T ] t0 t0
Z t Z Z
≤ sup | kf (τ, z1 (τ )) − f (τ, z2 (τ ))kdτ | (k fk ≤ kf k)
t∈[t0 −T ;t0 +T ] t0
Z t
≤ sup | kkz1 (τ ) − z2 (τ )kdτ | (f est k-lipschitzienne)
t∈[t0 −T ;t0 +T ] t0
Z t
≤ sup | kkz1 − z2 kdτ |
t∈[t0 −T ;t0 +T ] t0

= sup |t0 − t|kkz1 − z2 k = T kkz1 − z2 k

t∈[t0 −T ;t0 +T ]

Si T < k1 , on peut s’arrêter là. En fait, on peut raffiner un peu le résultat. En effet, par un
calcul très similaire, on trouve pour φ2 = φ ◦ φ :
Z t
kφ2 [z1 ] − φ2 [z2 ]k ≤ sup | kkφ[z1 ](τ ) − φ[z2 ](τ )kdτ |
t∈[t0 −T ;t0 +T ] t0
Z t
≤ sup | k sup kφ[z1 ] − φ[z2 ]kdτ |
t∈[t0 −T ;t0 +T ] t0 t∈[t0 ;τ ]
Z t
≤ sup | k × |τ − t0 |kkz1 − z2 kdτ |
t∈[t0 −T ;t0 +T ] t0

(t − t0 )2 2 T 2 k2
= sup k kz1 − z2 k = kz1 − z2 k
t∈[t0 −T ;t0 +T ] 2 2
p p p p
Par récurrence, on prouve que ∀p ∈ N, kφp [z1 ]−φp [z2 ]k ≤ T p!k kz1 −z2 k. Or limp→+∞ T p!k = 0,
donc ∃p, ∃α < 1 tels que ∀(z1 , z2 ) ∈ F × F, kφp [z1 ] − φp [z2 ]k ≤ αkz1 − z2 k. Par le théorème
du point fixe de Picard, φ a un point fixe, il est unique et la suite (φn [z])n∈N converge vers ce
point fixe. Notons ce point fixe y. D’après le lemme 4.1.1, y = φ(y) signifie que y est solution
du problème de Cauchy. 

Théorème 4.1.4 (Unicité globale)

Soient y1 et y2 deux solutions de l’équation différentielle sur I (leur condition initiale peut
être différente), avec f localement lipschitzienne en y. Si y1 et y2 coïncident en un point t1
de I, alors y1 (t) = y2 (t), ∀t ∈ I.

Démonstration. On considère le premier instant auquel y1 et y2 ne coïncident plus :

t̃0 = inf{t ∈ I, t ≥ t1 , y1 (t) 6= y2 (t)}.

D’après le théorème d’unicité locale, il existe T̃ tel que le problème de Cauchy avec conditions
initiales (t̃0 , y1 (t̃0 )) a une unique solution sur [t̃0 − T̃ ; t̃0 + T̃ ]. Ainsi y1 (t) = y2 (t) sur [t̃0 −
T̃ ; t̃0 + T̃ ]. Cela contredit la définition de t̃0 donc y1 et y2 coïncident sur I tout entier. 

Remarque : ce théorème ne garantit pas l’existence de solutions sur I, juste l’unicité.

34
4.2 Équations différentielles linéaires
Définition 4.2.1. Une équation différentielle y 0 (t) = f (t, y(t)) est linéaire si f est affine en
y, c’est-à-dire si il existe pour chaque t ∈ R une matrice carrée A(t) ∈ Mm (C) et un vecteur
b(t) ∈ C tels que f (t, y) = A(t)y + b(t).
f est continue dès que A et b sont des fonctions continues.
Pour tout t ∈ R, f est lipschitzienne en y de rapport k(t) = kA(t)k (norme d’opérateur).
Conséquence : Le théorème de Cauchy-Lipschitz s’applique et donc le problème de Cauchy
a une unique solution.

4.2.1 Géométrie de l’ensemble des solutions

Proposition 4.2.2
Soit S0 l’ensemble des solutions sur R de l’équation différentielle linéaire sans second membre
(E0 ) : y 0 (t) = A(t)y(t).
S0 est un espace vectoriel de dimension m.
Démonstration. Soient x et y deux solutions de (E0 ) et soient λ, µ ∈ C.

(λx + µy)0 (t) = λx0 (t) + µy 0 (t) = λA(t)x(t) + µA(t)y(t) = A(t)(λx(t) + µy(t))

donc λx + µy est une solution de (E0 ) et S0 est bien un espace vectoriel.

S0 → Cm
De plus, l’application φt0 d’évaluation au temps t0 définie par φt0 : est un
y 7→ y(t0 )
isomorphisme linéaire d’après le théorème d’existence (surjectivité) et d’unicité (injectivité).
Cela implique que dim(S0 ) = dim(Cm ) = m. 

Proposition 4.2.3
Soit S l’ensemble des solutions sur R de l’équation différentielle linéaire (E) : y 0 (t) = A(t)y(t)+
b(t) et soit y(1) une solution quelconque de (E).
Alors S = {y ∈ C(R, Cm ) : ∃z ∈ S0 , y = y(1) + z} = y(1) + S0 (espace affine).
Démonstration. Soit y ∈ S et soit z = y − y(1) . On a que z 0 (t) = y 0 (t) − y(1)
0 (t) = A(t)z(t)

donc z ∈ S0 . On en déduit que S ⊂ y(1) + S0 .

Soit y ∈ y(1) + S0 : y = y(1) + z où z ∈ S0 .

y 0 (t) = y(1)
0
(t) + z 0 (t) = A(t)y(1) (t) + b(t) + A(t)z(t) = A(t)y(t) + b(t)

et donc y(1) + S0 ⊂ S. 

4.2.2 Résolution de y 0 = Ay : équation différentielle linéaire à coefficients

constants et sans second membre
Rappel : Si m = 1, la solution du problème de Cauchy y 0 = ay, y(t0 ) = y0 est donnée par
y(t) = y0 ea(t−t0 ) , ∀t ∈ R.
Définition 4.2.4. L’exponentielle de la matrice A est définie par la somme de la série abso-
lument convergente
+∞
X 1
exp(A) = An .
n=0
n!

35
Théorème 4.2.5
La solution y du problème de Cauchy y 0 = Ay, y(t0 ) = y0 ∈ Cm est donnée par y(t) =
exp((t − t0 )A) × y0 , ∀t ∈ R.

Démonstration. y(t0 ) = exp(0 ∗ A) × y0 = I ∗ y0 = y0 (OK pour les conditions initiales).

+∞
X 1 n n
exp(tA) = t A
n=0
n!
+∞ +∞
d X 1 X 1
(exp tA) = ntn−1 An = tp Ap+1 = Aexp(tA)
dt n=0
n! p=0
p!

Par conséquent, y 0 (t) = A exp((t − t0 )A) × y0 = Ay(t) (OK pour l’équation différentielle).
D’après le théorème de Cauchy Lipschitz, on a trouvé l’unique solution. 

4.2.3 Résolution de y 0 = Ay+b(t) : équa. diff. linéaire à coefficients constants

Comme l’ensemble de solutions est un espace affine et que l’on sait résoudre y 0 = Ay, il
nous suffit de trouver une solution particulière yp .

La méthode de la variation de la constante. On cherche yp sous la forme yp (t) =

exp(tA) × v(t) où v est supposée dérivable.

yp0 (t) = A exp(tA)v(t) + exp(tA)v 0 (t) = Ayp (t) + exp(tA)v 0 (t).

Pour que yp0 (t) = Ayp (t) + b(t), il suffit de choisir v telle que pour tout t, exp(tA)v 0 (t) = b(t).
On prend donc, pour un certain t0 ∈ R, v(t) = tt0 exp(−uA)b(u)du et on obtient
R

Z t
yp (t) = exp((t − u)A)b(u)du.
t0

La solution du problème de Cauchy avec conditions initiales (t0 , y0 ) est alors

Z t
y(t) = exp((t − t0 )A)y0 + yp (t) = exp((t − t0 )A)y0 + exp((t − u)A)b(u)du.
t0

Conclusion Pour résoudre les équations différentielles linéaires à coefficients constants, il

suffit de savoir calculer :
1. l’exponentielle d’une matrice,
2. l’intégrale d’une fonction de R dans Cm .

36
4.3 Méthodes numériques à un pas
4.3.1 Introduction
On veut résoudre le problème de Cauchy
(
y 0 (t) = f (t, y(t)), ∀t ∈ I
y(t0 ) = y0

Mais très souvent, on n’arrive pas à calculer la solution exacte. De la même manière que pour
les intégrales, on va chercher des solutions approchées.

Table 4.1 – Idée générale des méthodes à 1 pas

équation différentielle méthode numérique

I ⊂ R, I = [t0 , t0 + T ] discrétisation : t0 , t1 , . . . , tN
tN = t0 + T
hn = tn+1 − tn

yn+1 − yn
y 0 (tn )
hn

f (tn , y(tn )) φ(tn , yn , hn )

( (
y 0 (t) = f (t, y(t)), ∀t ∈ I yn+1 = yn + hn φ(tn , yn , hn )
y(t0 ) = y0 y0 donné
yn est calculable par récurrence

La méthode d’Euler Elle consiste à choisir φ(tn , yn , hn ) = f (tn , yn ) et donc

yn+1 = yn + hn f (tn , yn ).

37
4.3.2 Consistance, stabilité, convergence
Définition 4.3.1 (Erreur de consistance locale). Soit z la solution du problème de Cauchy
avec condition initiale z(tn ) = yn
(
z 0 (t) = f (t, z(t))
z(tn ) = yn
L’erreur de consistance locale est définie par en = z(tn+1 ) − yn+1 = z(tn+1 ) − yn −
φ(tn , yn , hn ).

Exemple : La méthode d’Euler, yn+1 = yn + hn f (tn , yn )

Si f est C 1 ,
(tn+1 − tn )2 00
z(tn+1 ) = z(tn ) + (tn+1 − tn )z 0 (tn ) + z (tn ) + o((tn+1 − tn )2 ).
2

On a z(tn ) = yn , z 0 (tn ) = f (tn , z(tn )) = f (tn , yn ) et

∂f ∂f déf
z 00 (tn ) = (tn , z(tn )) + (tn , z(tn )) ×z 0 (tn ) = f [1] (tn , yn )
∂t ∂z
| {z }
Jacobienne
Ainsi,
h2n [1]
ken k = kz(tn+1 ) − yn+1 k = kyn + hn f (tn , yn ) + f (tn , yn ) + o(h2n ) − yn − hn f (tn , yn )k
2
h2n [1]
= kf (tn , yn )k + o(h2n ).
2
Proposition 4.3.2
L’erreur de consistance locale vérifie ken k ≤ Chp+1 n avec C > 0 indépendant de hn si et
seulement si :
• f est de classe C p
∂lφ 1 [l]
• ∀l, 0 ≤ l ≤ p − 1, l (t, y, 0) = f (t, y)
∂h l+1
 dl 
où f [l] (t, y) =


f (τ, z(τ )) z(t) = y
(t)
dτ l 0
z (τ ) = f (τ, z(τ ))
Démonstration.
yn+1 = yn + hn φ(tn , yn , hn )
h ∂φ hp−1
n ∂ p−1 φ hpn ∂ p φ i
= yn + hn φ(tn , yn , 0) + hn (tn , yn , 0) + . . . + p−1
(t ,
n n y , 0) + p
(tn , yn , 0) + o(hpn )
∂h (p − 1)! ∂h p! ∂h
h 2 p+1
hn
z(tn+1 ) = z(tn ) + hn f (tn , z(tn )) + n f [1] (tn , z(tn )) + . . . + f [p] (tn , z(tn )) + o(hp+1
n )
2 (p + 1)!
p " l l #
X hn ∂ φ hln
ken k = kyn+1 − z(tn+1 k = khn (tn , yn , 0) − f [l] (tn , z(tn )) + o(hp+1
n )k
l=0
l! ∂hl (l + 1)!
1 ∂pφ 1
= hp+1
n k p
(tn , yn , 0) − f [p] (tn , z(tn ))k + o(hp+1
n ) 
p! ∂h (p + 1)!

38
Définition 4.3.3 (Consistance). La méthode φ est consistante si pour toute solution exacte
z, la somme des erreurs de consistance relatives à z tend vers 0 quand hmax = maxn hn tend
vers 0.
Soit en = z(tn+1 ) − z(tn ) − hn φ(tn , z(tn ), hn ). La définition signifie que
N
X −1
lim ken k = 0.
hmax →0
n=0

Définition 4.3.4. La méthode φ est stable si il existe une constante S ≥ 0 indépendante de

N , appelée constante de stabilité, telle que pour toutes suites (yn ) et (ỹn ) définies par

yn+1 = yn + hn φ(tn , yn , hn ), 0≤n≤N −1

ỹn+1 = ỹn + hn φ(tn , ỹn , hn ) + n , 0≤n≤N −1
PN −1

on ait max0≤n≤N kỹn − yn k ≤ S kỹ0 − y0 k + n=0 kn k .
Remarque : Stabilité ⇒ propagation des erreurs maîtrisée
Définition 4.3.5 (Convergence). La méthode φ est convergente si pour toute solution exacte
z, la suite (yn ) telle que yn+1 = yn + hn φ(tn , yn , hn ) vérifie

max kyn − z(tn )k → 0 quand y0 → z(t0 ) et hmax → 0.

0≤n≤N

La quantité max0≤n≤N kyn − z(tn )k s’appelle l’erreur globale

Théorème 4.3.6
Si la méthode est stable et consistante, alors elle est convergente.
Démonstration. Soit z la solution exacte du problème de Cauchy.
Posons ỹn = z(tn ) : par définition de l’erreur de consistance locale en , on a

ỹn+1 = ỹn + hn φ(tn , ỹn , hn ) + en .

Comme la méthode est stable de constante de stabilité S, on a

N
X −1


max kyn − z(tn )k = max kyn − ỹn k ≤ S ky0 − z(t0 )k + ken k .
0≤n≤N 0≤n≤N
n=0

39
 PN −1
Comme la méthode est consistante, limhmax →0 n=0 ken k = 0. On obtient que

lim max kyn − z(tn )k = 0.

hmax → 0 0≤n≤N
y0 → z(t0 )

4.3.3 *Condition nécessaire et suffisante de consistance

On sait déjà que si f est C 1 et φ(t, y, 0) = f (t, y), alors ∀n, ken k ≤ Ch2n . Pour le résultat
suivant, on va supposer uniquement que f est C 0 et on va majorer N n=0 ken k.
P

Théorème 4.3.7
Supposons que f : R×Rm → Rm est continue. La méthode à 1 pas définie par φ est consistante
si et seulement si
∀t ∈ [t0 , t0 + T ], ∀y ∈ Rm , φ(t, y, 0) = f (t, y).

Démonstration. Soit z une solution exacte du problème de Cauchy et soit en l’erreur de

consistance locale en = z(tn+1 ) − z(tn ) − hn φ(tn , z(tn ), hn ).
Z tn+1
z(tn+1 ) = z(tn ) + f (t, z(t))dt
tn
ken k = kz(tn+1 ) − z(tn ) − hn φ(tn , z(tn ), hn )k
Z tn+1 Z tn+1
=k f (t, z(t))dt − φ(tn , z(tn ), hn )dtk
tn tn
Z tn+1 Z tn+1 Z tn+1
=k f (t, z(t))dt − φ(t, z(t), 0)dt + φ(t, z(t), 0) − φ(tn , z(tn ), hn )dtk
tn tn tn
Z tn+1 Z tn+1
≤ kf (t, z(t)) − φ(t, z(t), 0)kdt + kφ(t, z(t), 0) − φ(tn , z(tn ), hn )kdt
tn tn
| {z } | {z }
An Bn

PN
ke k = N A + N
P P
B
PN n
n=0 PN n
n=0 PN n R tn+1
n=0
kf (t, z(t))−φ(t, z(t), 0)kdt = tT0 kf (t, z(t))−φ(t, z(t), 0)kdt
R
An : n=0 An = n=0 tn
Pn=0
N
n=0 Bn : φ est continue donc ((t, h) 7→ φ(t, z(t), h)) l’est aussi.
Pour tout δ, [t0 , t0 + T ] × [0, δ] est compact donc ((t, h) 7→ φ(t, z(t), h)) est uniformément
continue sur cet ensemble. Cela veut dire que ∀ > 0, ∃η > 0 (η ≤ δ) tel que ∀hmax ≤ η, on a
0) − φ(tn , z(tn ), hn )k ≤ . Ainsi, N
PN
∀n, ∀t ∈ [t0 , t0 + T ], kφ(t, z(t), n=0 Bn ≤
P
PN RT n=0 hn  = T .
Au final, n=0 ken k = t0 kf (t, z(t)) − φ(t, z(t), 0)kdt + T .
Si ∃J ⊂ [0; T ] tel que f (t, z(t)) 6= φ(t, z(t), 0) sur J, alors N n=0 ken k 6→ 0.
P
PN
Si ∀t, f (t, z(t)) = φ(t, z(t), 0), alors n=0 ken k → 0. 

40
4.3.4 Condition suffisante de stabilité
Théorème 4.3.8
Si φ est lipschitzienne en y de constante Λ, alors la méthode est stable avec constante de
stabilité eΛT .
Dit autrement, si ∀t ∈ [t0 , t0 + T ], ∀y1 , y2 ∈ Rm , ∀h ≥ 0, on a kφ(t, y1 , h) − φ(t, y2 , h)k ≤
Λky1 − y2 k, alors pour yn+1 = yn + hn φ(tn , yn , hn ) et ỹn+1 = ỹn + hn φ(tn , ỹn , hn ) + n , on a

 N
X 
ΛT
max kỹn − yn k ≤ e kỹ0 − y0 k + ken k .
0≤n≤N
n=0

Lemme 4.3.9 (Gronwall, cas discret)

Soient (hn ), (θn ), (n ) des suites de R+ telles que pour tout n,

θn+1 ≤ (1 + Λhn )θn + n .

Alors pour tout n,

n−1
θn ≤ eΛ(tn −t0 ) θ0 + eΛ(tn −ti ) i .
X

i=1

Démonstration du lemme. Preuve par récurrence.

Pour n = 0, il suffit de vérifier que θ0 ≤ θ0 .
Supposons que θn ≤ eΛ(tn −t0 ) θ0 + n−1 Λ(tn −ti )  .
P
i=1 e i

 n−1 
θn+1 ≤ (1 + Λhn )θn + n ≤ (1 + Λhn ) eΛ(tn −t0 ) θ0 + eΛ(tn −ti ) i + n
X

i=1

Mais comme 1 + Λhn ≤ eΛhn = eΛ(tn+1 −tn ) , et n ≤ eΛ(tn+1 −tn ) n ,

n−1 n
θn+1 ≤ eΛ(tn+1 −tn +tn −t0 ) θ0 + eΛ(tn+1 −tn +tn −ti ) i +eΛ(tn+1 −tn ) n = eΛ(tn+1 −t0 ) θ0 + eΛ(tn+1 −ti ) i .
X X

i=1 i=1


Démonstration du théorème. Considérons deux suites (yn ) et (ỹn ) telles que

yn+1 = yn + hn φ(tn , yn , hn )
ỹn+1 = ỹn + hn φ(tn , ỹn , hn ) + n

Par différence, on obtient kyn+1 − ỹn+1 k ≤ kyn − ỹn k + hn Λkyn − ỹn k + kn k.
On peut utiliser le lemme de Gronwall avec θn = kyn − ỹn k, hn = hn et n = kn k.
Pn−1 Λ(tn −ti )
Cela donne ∀n, kyn − ỹn k ≤ eΛ(tn −t0 ) ky0 − ỹ0 k + i=0 e ki k
Pn
et donc ∀n, kyn − ỹn k ≤ eΛT ky0 − ỹ0 k + i=0 kn k , vu que eΛ(tn −ti ) ≤ eΛT . 

4.3.5 Majoration de l’erreur globale

• Si la méthode est d’ordre p et f est de classe C p , ken k ≤ Chp+1
n où C est une constante
indépendante de n et de hn ≤ δ.

41
 • Si la méthode est stable de constante de stabilité S,
 N
X 
max kyn − z(tn )k ≤ S ky0 − z(t0 )k + ken k .
0≤n≤N
n=0
PN PN p+1 PN
En combinant les deux, on obtient n=0 ken k ≤ ≤ Chpmax
n=0 Chn n=0 hn = CT hpmax
et
max kyn − z(tn )k ≤ Sky0 − z(t0 )k + SCT hpmax .
0≤n≤N
Si φ est lipschitzienne en y de constante Λ,
max kyn − z(tn )k ≤ eΛT ky0 − z(t0 )k + CT eΛT hpmax .
0≤n≤N

4.4 Méthodes de Runge-Kutta

4.4.1 Principe
Le problème de Cauchy à résoudre est
(
y 0 (t) = f (t, y(t)), t ∈ [t0 , t0 + T ]
y(t0 ) = y0
[t0 , t0 + T ] est discrétisé en t0 , t1 , . . . , tN . Rt
Soit z la solution exacte du problème de Cauchy : z(tn+1 ) = z(tn ) + tnn+1 f (t, z(t))dt
Rt
La méthode d’Euler yn+1 = yn + hn f (tn , yn ) correspond à approcher tnn+1 f (t, z(t))dt par
les rectangles à gauche, de manière successive.
Les méthodes de Runge-Kutta correspondent à schéma d’intégration quelconque.

Il faut des points intermédiaires entre tn et tn+1 :

tn,i = tn + ci hn , 1≤i≤q ci ∈ [0, 1]

Rappel : Soit g : R → Rm . 01 g(t)dt ≈ qj=1 bj g(cj ) où qj=1 bj = 1.

R P P

Les bj sont les coefficients donnés par la technique d’approximation de l’intégrale.

Difficulté : On aimerait prendre g(t) = f (t, z(t)) mais z n’est pas connu. Ainsi :
• On approche z(tn ) par yn .
Rt
• On approche z(tn,i ) = z(tn ) + tnn,i f (t, z(t))dt = z(tn ) + 0ci f (tn + uhn , z(tn + uhn ))du
R
Pi−1 Pi−1
par yn,i = yn + hn j=1 ai,j f (tn,j , yn,j ), 1 ≤ i ≤ q où j=1 ai,j = ci .
• On approche z(tn+1 ) par yn+1 = yn + hn qj=1 bj f (tn,j , bn,j ) où qj=1 bj = 1.
P P

Notez que les techniques d’intégration peuvent être toutes différentes.

On obtient l’algorithme suivant :

 
tn,i = tn + ci hn

 



∀i ∈ {1, . . . q} y = y + h Pi−1 a p

tn+1 = tn + hn
n,i n n j=1 i,j n,j


p = f (t , y )
n,i n,i n,i



 Pq
yn+1 = yn + hn j=1 bj pn,j



42
 L’algorithme est bien défini de manière explicite car la somme qui définit yn,i n’utilise que
les yn,j pour j < i. On peut le représenter par le tableau

c1 0 0 ... 0 0
c2 a2,1 0 ... 0 0
.. .. .. .. ..
. . . . .
.. .. .. .. ..
. . . . .
cq aq,1 aq,2 ... aq,q−1 0
b1 b2 ... bq−1 bq
Chaque ligne correspond à une méthode d’intégration approchée.
Pq
Il y a des contraintes sur les coefficients : ci = i−1
P
j=1 ai,j , 1 = j=1 bj , c1 = 0, tn,1 = tn ,
yn,1 = yn , pn,1 = f (tn , yn ).

4.4.2 Exemples
q=1
0 0
Le seul choix possible est .
1
c1 = 0, a1,1 = 0, b1 = 1.
L’algorithme est donné par pn,1 = f (tn , yn ), tn+1 = tn + hn , yn+1 = yn + hn pn . Il s’agit
de la méthode d’Euler.
q=2
0 0 0
Pour α ∈]0; 1], et β ∈ R, on a la méthode définie par le tableau : α α 0
1−β β
Exercice : Pour quelles valeurs de α et β la méthode est-elle d’ordre 2 ?
• Pour α = 12 et β = 1, on a la méthode du point milieu :
pn,1 = f (tn , yn )
tn,2 = tn + 12 hn
yn,2 = yn + 12 hn pn,1 
ie : yn+1 = yn + hn f (tn + 21 hn , yn + 12 hn f (tn , yn )
pn,2 = f (tn,2 , yn,2 )
tn+1 = tn + hn
yn+1 = yn + hn pn,2
• Pour α = 1 et β = 12 , on a la méthode de Heun :
pn,1 = f (tn , yn )
tn,2 = tn + hn
yn,2 = yn + hn pn,1
pn,2 = f (tn,2 , yn,2 )
tn+1 = tn + hn
yn+1 = yn + 12 hn pn,1 + 12 hn pn,1
La méthode de Heun peut aussi s’écrire 
yn+1 = yn + hn 21 f (tn , yn ) + 12 f (tn+1 , yn + hn f (tn , yn ) .
Cependant cette formule fait intervenir 3 évaluations de f alors que le fait de stocker
f (tn , yn ) dans pn,1 permet de n’en faire que 2.

43
 RK4 (Runge-Kutta classique avec q = 4)
L’algorithme s’écrit
pn,1 = f (tn , yn )
tn,2 = tn + 12 hn
yn,2 = yn + 12 hn pn,1
pn,2 = f (tn,2 , yn,2 )
0 0 0 0 0
tn,3 = tn,2 1 1
2 2 0 0 0 → rectangles à gauche
yn,3 = yn + 21 hn pn,2 1
2 0 1
2 0 0 → rectangles à droite
pn,3 = f (tn,3 , yn,3 ) 1 0 0 1 0 → point milieu
1 2 2 1
6 6 6 6 → Simpson
tn,4 = tn + hn
yn,4 = yn + hn pn,3
pn,4 = f (tn,4 , yn,4 )

tn+1 = tn + hn
1
+ 26 pn,2 + 62 pn,3 + 16 pn,4


yn+1 = yn + hn 6 pn,1

4.4.3 Stabilité des méthodes de Runge-Kutta

Ce sont des méthodes à un pas : yn+1 = yn + hn φ(tn , yn , hn )
Pi−1
où φ(t, y, h) = aj=1 bj f (t + cj h, yj ) et yi = y + h j=1
P
ai,j f (t + cj h, yj )

Pour l’étude de la stabilité, on définit pour tout i la fonction ȳi telle que
ȳi (t, y, h) = y + h i−1
P
j=1 ai,j f (t + cj h, ȳj (t, y, h)).

Lemme 4.4.1
Pi
Rm ,-


Si f est k-lipschitzienne en y alors avec α = max1≤i≤q j=1 |ai,j | , on a ∀y, z ∈
∀i ∈ {1, . . . , q}, kȳi (t, y, h) − ȳi (t, z, h)k ≤ 1 + (αkh) + (αkh)2 + . . . (αkh)i−1 ky − zk

Démonstration. On fait une récurrence sur i. Pour i = 1, on a ȳ1 (t, y, h) = y et ȳ1 (t, z, h) = z
donc c’est bon.
Supposons l’inégalité vraie pour tout j < i.

i−1
X i−1
X
kȳi (t, y, h) − ȳi (t, z, h)k ≤ y + h ai,j f (t + cj h, ȳj (t, y, h)) − z − h ai,j f (t + cj h, ȳj (t, z, h))
j=1 j=1
i−1
X
≤ ky − zk + h |ai,j | kf (t + cj h, ȳj (t, y, h)) − f (t + cj h, ȳj (t, z, h))k
j=1
i−1
X
≤ ky − zk + h |ai,j |kkȳj (t, y, h) − ȳj (t, z, h)k
j=1

44
On utilise l’hypothèse de récurrence pour obtenir
i−1
X
|ai,j |k 1 + (αkh) + (αkh)2 + . . . (αkh)j−1 ky − zk


kȳi (t, y, h) − ȳi (t, z, h)k ≤ ky − zk + h
j=1
i−1
X
|ai,j |k 1 + (αkh) + (αkh)2 + . . . (αkh)i−2 ky − zk


≤ ky − zk + h
j=1

≤ 1 + (αkh) + (αkh)2 + . . . (αkh)i−1 ky − zk



Théorème 4.4.2
Les méthodes de Runge Kutta sont stables, avec constante de stabilité S = eΛT où
Λ = k qj=1 |bj |(1 + (αkhmax ) + . . . + (αkhmax )j−1 ).
P

Démonstration. Il suffit de montrer que φ est Λ-lipschitzienne en y.

Soient y, z ∈ Rm , t ∈ R et h > 0.
q
X

kφ(t, y, h) − φ(t, z, h)k = k bj f (t + cj h, ȳj (t, y, h) − f (t + cj h, ȳj (t, z, h)) k
j=1
q
X
≤ |bj |kkȳj (t, y, h) − ȳj (t, z, h)k
j=1
Xq
≤k |bj |(1 + (αkhmax ) + . . . + (αkhmax )j−1 ) 
j=1

1
Remarquons que si hmax est petit devant αk , alors S est de l’ordre de ekT , ce qui est
comparable à la méthode d’Euler.

4.4.4 Ordre des méthodes de Runge-Kutta

On rappelle que
• φ(t, y, h) = qj=1 bj f (t + cj h, ȳj (t, y, h))
P

où ȳj (t, y, h) = y + h i−1

P
j=1 ai,j f (t + cj h, ȳj (t, y, h)).
• L’ordre de la méthode est p au moins si et seulement si ∀l ≤ p − 1, ∀t ∈ R, ∀y ∈ Rm ,

∂lφ 1 [l]
l
(t, y, 0) = f (t, y)
∂h l+1
∂0φ
= φ(t, y, 0) = f (t, y) car ȳi (t, y, 0) = y et qj=1 bj = 1. Ainsi les méthode de
P
1. ∂h0
(t, y, 0)
Runge-Kutta sont toujours d’ordre 1 et donc consistantes.
∂1φ
 
Pq ∂f ∂f ∂ ȳj
2. ∂h1
(t, y, h) = j=1 bj ∂t (t + cj h, ȳj (t, y, h)) × cj + ∂y (t + cj h, ȳj (t, y, h)) × ∂h (t, y, h)
 
∂ ȳi P P ∂f ∂f ∂ ȳj
∂h (t, y, h) = j<i ai,j f (t + cj h, ȳj (t, y, h)) + h j<i ai,j ∂t (−)cj + ∂y (−) ∂h (−)

Ainsi ∂∂h ȳi P

(t, y, 0) = j<i ai,j f (t, y) = ci f (t, y).
∂1φ
(t, y, 0) = qj=1 bj ∂f ∂f Pq
[1]
∂t (t, y) × cj + ∂y (t, y) × f (t, y) × cj =
P
∂h1 j=1 bj cj f (t, y)
On en déduit que la méthode est d’ordre 2 si et seulement si qj=1 bj cj = 12 .
P

45
3. Il faut dériver encore une fois.
q
∂2φ X ∂f 2 ∂2f ∂ ȳj
(t, y, h) = bj (t + cj h, ȳ j (t, y, h))cj + 2cj (t + cj h, ȳj (t, y, h)) (t, y, h)
∂h2 j=1
∂t2 ∂t∂y ∂y
!
∂f ∂ 2 ȳj D ∂2f ∂ ȳj ∂ ȳj E
+ (t + cj h, ȳj (t, y, h)) 2 (t, y, h) + (t + cj h, ȳ j (t, y, h)) (t, y, h), (t, y, h)
∂y ∂h ∂y 2 ∂h ∂h

∂ 2 ȳi X  ∂f ∂f ∂ ȳj  X  ∂2f 

(t, y, h) = 2 a i,j (−)cj + (−) (−) + h a i,j . . .
∂h2 j<i
∂t ∂y ∂h j<i
∂y 2
∂ 2 ȳi X
(t, y, 0) = 2 ai,j cj f [1] (t, y)
∂h2 j<i

q
∂2φ X ∂f 2
2
2 ∂ f
(t, y, 0) = bj (t, y)cj + 2cj (t, y)f (t, y)
∂h2 j=1
∂t2 ∂t∂y
!
∂f X D ∂2f E
+ (t, y)2 aj,l f [1] (t, y) + (t, y)cj f (t, y), cj f (t, y)
∂y l<j
∂y 2
q
X  ∂f ∂2f D ∂2f E
= bj c2j (t, y) + 2 (t, y)f (t, y) + (t, y)f (t, y), f (t, y)
j=1
∂t2 ∂t∂y ∂y 2
q X
q
X ∂f
+2 bj aj,l cl (t, y)f [1] (t, y)
j=1 l=1
∂y

D’un autre côté,

∂f ∂2f D ∂2f E ∂f
f [2] (t, y) = (t, y) + 2 (t, y)f (t, y) + (t, y)f (t, y), f (t, y) + (t, y)f [1] (t, y)
∂t2 ∂t∂y ∂y 2 ∂y
∂2φ Pq Pq
Pour avoir ∂h2
(t, y, 0) = 13 f [2] (t, y), il faut que 2 1
j=1 bj cj = 3 et
1
i,j=1 bi ai,j cj = 6 .
Pq 1 q 1 Pq
Ainsi la méthode est d’ordre 3 ⇔ j=1 bj cj = 2 et i,j=1 bi ai,j cj = 16 .
P
2, j=1 bj cj = 3
4. On peut obtenir des résultats similaires pour les ordres suivants. Pour l’ordre 4 :

q q q q
X 1 X 2 1 X 1 X 3 1
ordre 4 ⇔ bj cj = ; bj cj = ; bi ai,j cj = ; bj c = ;
j=1
2 j=1 3 i,j=1 6 j=1 j 4
q q q
X 1 X 1 X 1
bi ai,j c2j = ; bi ci ai,j cj = et bi ai,j aj,k ck =
i,j=1
12 i,j=1 8 i,j,k=1
12

Exemple : La méthode de Runge-Kutta classique est d’ordre 4.

0 0 0 0 0
1 1
2 2 0 0 0
1 1
2 0 2 0 0
1 0 0 1 0
1 2 2 1
6 6 6 6

46
 P
1. On a bien j bj =1
1 2 1 2 1 1 1
×0+ × × ×1=
P
2. j bj cj = 6 6 2 + 6 2 + 6 2
2 1 2 1 2 1 1 1
× 02 + × × ×1=
P
3. j bj cj = 6 6 4 + 6 4 + 6 3
2 1 2 1 1 1 1 1
× ×0+ × × ×1×
P
i,j bi ai,j cj = b2 a2,1 c1 + b3 a3,2 c2 + b4 a4,3 c3 = 6 2 6 2 2 + 6 2 = 6
P 3
4. j bj cj = . . . vérification par ordinateur.

47
Bibliographie

Jean-Pierre Demailly. Analyse numérique et équations différentielles. EDP sciences, 2012.

Alfio Maria Quarteroni, Riccardo Sacco, and Fausto Saleri. Méthodes numériques : algo-
rithmes, analyse et applications. Springer Science & Business Media, 2008.

48