Chapitre 6

Différentiation et intégration
numériques

6.1 Introduction
Le contenu de ce chapitre prolonge celui du chapitre 5 sur l’interpolation.
À peu de choses près, on y manie les mêmes outils d’analyse. Dans le cas de
l’interpolation, on cherchait à évaluer une fonction f (x) connue seulement en
quelques points. Dans le présent chapitre, le problème consiste à obtenir des
approximations des différentes dérivées de cette fonction de même que de :
! xn
f (x)dx
x0

On parle alors de dérivation numérique et d’intégration numérique. On fait

face à ce type de problèmes lorsque, par exemple, on connaît la position d’une
particule à intervalles de temps réguliers et que l’on souhaite obtenir sa vi-
tesse. On doit alors effectuer la dérivée de la position connue seulement en
quelques points. De même, l’accélération de cette particule nécessite le calcul
de la dérivée seconde.
Si, à l’inverse, on connaît la vitesse d’une particule à certains intervalles de
temps, on obtient la distance parcourue en intégrant la vitesse dans l’intervalle
[x0 , xn ].
Nous avons vu au chapitre précédent que la fonction f (x) peut être conve-
nablement estimée à l’aide d’un polynôme de degré n avec une certaine erreur.
En termes concis :
f (x) = pn (x) + En (x) (6.1)
où En (x) est le terme d’erreur d’ordre (n + 1) donné par la relation 5.21.
L’expression 6.1 est à la base des développements de ce chapitre.

6.2 Différentiation numérique

On peut aborder la différentiation numérique d’au moins deux façons. La
première approche consiste à utiliser le développement de Taylor et la seconde
292 Chapitre 6

est fondée sur l’égalité 6.1. Nous utiliserons un mélange des deux approches,
ce qui nous permettra d’avoir un portrait assez complet de la situation.
Commençons d’abord par l’équation 6.1. Si l’on dérive de chaque côté de
l’égalité, on obtient successivement :

f ′ (x) = p′n (x) + En′ (x)

f ′′ (x) = p′′n (x) + En′′ (x)
f ′′′ (x) = p′′′ ′′′
n (x) + En (x)
(6.2)
.. ..
. = .

Ainsi, pour évaluer la dérivée d’une fonction connue aux points ((xi , f (xi ))
pour i = 0, 1, 2, · · · , n), il suffit de dériver le polynôme d’interpolation passant
par ces points. De plus, le terme d’erreur associé à cette approximation de la
dérivée est tout simplement la dérivée de l’erreur d’interpolation. Ce résultat
est vrai quel que soit l’ordre de la dérivée.
Remarque 6.1
Bien qu’en théorie on soit en mesure d’estimer les dérivées de tout ordre, sur
le plan pratique, on dépasse rarement l’ordre 4. Cela s’explique par le fait que
la différentiation numérique est un procédé numériquement instable. !

6.2.1 Dérivées d’ordre 1

Commençons par faire l’approximation des dérivées d’ordre 1, ce qui re-
vient à évaluer la pente de la fonction f (x). Tout comme pour l’interpolation,
nous avons le choix entre plusieurs polynômes de degré plus ou moins élevé.
De ce choix dépendent l’ordre et la précision de l’approximation. Nous avons
rencontré un problème semblable dans le cas de l’interpolation : si un poly-
nôme de degré n est utilisé, on obtient une approximation d’ordre (n + 1) de
la fonction f (x) (voir la relation 5.26).
Il est également utile de se rappeler que l’erreur d’interpolation s’écrit :

f (n+1) (ξ(x))
En (x) = [(x − x0 )(x − x1 ) · · · (x − xn )] (6.3)
(n + 1)!

pour un certain ξ compris dans l’intervalle [x0 , xn ]. En dérivant l’expression

précédente, tout en tenant compte de la dépendance de ξ envers x, on obtient
une relation pour la dérivée de l’erreur d’interpolation :

f (n+2) (ξ(x))ξ ′ (x)

En′ (x) = [(x − x0 )(x − x1 ) · · · (x − xn )]
(n + 1)!

f (n+1) (ξ(x))
+ [(x − x0 )(x − x1 ) · · · (x − xn )]′
(n + 1)!

La dérivée du produit apparaissant dans le deuxième terme de droite est plus

délicate. Cette dérivée débouche sur une somme de produits où, tour à tour,
 Différentiation et intégration numériques 293

l’un des facteurs (x − xi ) est manquant. Il est facile de se convaincre, en repre-

nant ce développement avec n = 2 par exemple, que l’on obtient :

f (n+2) (ξ(x))ξ ′ (x)

En′ (x) = [(x − x0 )(x − x1 ) · · · (x − xn )]
(n + 1)!
⎛ ⎞ (6.4)
(n+1) $ n %n
f (ξ(x)) ⎝
+ (x − xj )⎠
(n + 1)!
k=0 j=0(j̸=k)

On peut simplifier cette expression quelque peu complexe en choisissant l’un

ou l’autre des points d’interpolation. En effet, en x = xi , le premier terme de
droite s’annule, faisant disparaître la dérivée de ξ(x), qui est inconnue. De la
somme, il ne reste qu’un seul terme puisque tous les autres contiennent un
facteur (x − xi ) et s’annulent. Il reste :
⎛ ⎞
n
%
f (n+1) (ξ(xi )) ⎝
En′ (xi ) = (xi − xj )⎠
(n + 1)!
j=0(j̸=i)

Si l’on suppose de plus que les xi sont également distancés, c’est-à-dire :

xi+1 − xi = h

ce qui signifie que xi − xj = (i − j)h, on obtient :

⎛ ⎞
n
%
f (n+1) (ξi )hn ⎝
En′ (xi ) = (i − j)⎠ (6.5)
(n + 1)!
j=0(j̸=i)

où ξi est simplement une notation différente de ξ(xi ). En particulier, si i = 0,

on trouve :
⎛ ⎞ ⎛ ⎞
(n+1) n n
% (n+1) n n
%
f (ξ0 )h ⎝ f (ξ0 )h ⎝
En′ (x0 ) = (−j)⎠ = (−j)⎠
(n + 1)! (n + 1)!
j=0(j̸=0) j=1

c’est-à-dire :
(−1)n hn f (n+1) (ξ0 )
En′ (x0 ) = (6.6)
(n + 1)
pour un certain ξ0 compris dans l’intervalle [x0 , xn ].
Remarque 6.2
L’équation 6.5 montre que, si l’on utilise un polynôme d’interpolation de degré
n (c’est-à-dire d’ordre (n + 1)), la dérivée de ce polynôme évaluée en x = xi
est une approximation d’ordre n de f ′ (xi ). !
294 Chapitre 6

Définition 6.3
Aux points d’interpolation, on a :

f ′ (xi ) = p′n (xi ) + En′ (xi ) (6.7)

Le terme p′n (xi ) dans l’équation 6.7 est une formule aux différences finies
ou plus simplement une formule aux différences. Nous proposons plus loin
plusieurs formules aux différences finies pour évaluer les différentes dérivées
de f (x). Elles se distinguent principalement par le degré du polynôme et par
les points d’interpolation retenus.

Appproximations d’ordre 1
Si l’on choisit le polynôme de degré 1 passant par les points (x0 , f (x0 )) et
(x1 , f (x1 )), on a, grâce à la formule d’interpolation de Newton :

p1 (x) = f (x0 ) + f [x0 , x1 ](x − x0 )

et donc :
f ′ (x) = p′1 (x) + E1′ (x) = f [x0 , x1 ] + E1′ (x) (6.8)
En vertu de la relation 6.6 avec n = 1 et puisque (x1 − x0 ) = h, on arrive à :

f (x1 ) − f (x0 ) f (x1 ) − f (x0 ) (−1)1 h1 f (2) (ξ0 )

f ′ (x0 ) = + E1′ (x0 ) = +
x1 − x0 h 2

qui peut encore s’écrire :

′ f (x1 ) − f (x0 ) hf (2) (ξ0 )

f (x0 ) = − pour ξ0 ∈ [x0 , x1 ] (6.9)
h 2
qui est la différence avant d’ordre 1. On l’appelle différence avant car, pour
évaluer la dérivée en x = x0 , on cherche de l’information vers l’avant (en
x = x1 ).
De la même manière, si l’on évalue l’équation 6.8 en x = x1 , la relation 6.5
avec (i = 1) donne :

f (x1 ) − f (x0 )
f ′ (x1 ) = + E1′ (x1 )
x1 − x0
⎛ ⎞
1
%
f (x1 ) − f (x0 ) h1 f (2) (ξ1 ) ⎝
= + (1 − j)⎠
h 2!
j=0(j̸=1)

ou encore :

′ f (x1 ) − f (x0 ) hf (2) (ξ1 )

f (x1 ) = + pour ξ1 ∈ [x0 , x1 ] (6.10)
h 2
qui est la différence arrière d’ordre 1.
 Différentiation et intégration numériques 295

On constate ainsi que la même différence divisée est une approximation de

la dérivée à la fois en x = x0 et en x = x1 . On remarque cependant que le
terme d’erreur est différent aux deux endroits.
Appproximations d’ordre 2
Passons maintenant aux polynômes de degré 2. Soit les points (x0 , f (x0 )),
(x1 , f (x1 )) et (x2 , f (x2 )). Le polynôme de degré 2 passant par ces trois points
est :

p2 (x) = f (x0 ) + f [x0 , x1 ](x − x0 ) + f [x0 , x1 , x2 ](x − x0 )(x − x1 )

dont la dérivée est :

p′2 (x) = f [x0 , x1 ] + f [x0 , x1 , x2 ](2x − (x0 + x1 ))

Lorsque x prend successivement les valeurs x0 , x1 et x2 , il est facile de montrer

que l’on obtient des approximations d’ordre 2 de la dérivée.

Formules de différences d’ordre 2 pour f ′ (x)

−f (x2 ) + 4f (x1 ) − 3f (x0 ) h2 f ′′′ (ξ0 )

f ′ (x 0) = +
2h 3
Différence avant d’ordre 2

f (x2 ) − f (x0 ) h2 f ′′′ (ξ1 ) (6.11)

f ′ (x1 ) = −
2h 6
Différence centrée d’ordre 2

3f (x2 ) − 4f (x1 ) + f (x0 ) h2 f ′′′ (ξ2 )

f ′ (x2 ) = +
2h 3
Différence arrière d’ordre 2

Les termes d’erreur de ces formules aux différences finies découlent tous de
la relation 6.5 lorsqu’on pose successivement i = 0, 1 et 2. Pour i = 0, on peut
utiliser directement l’équation 6.6. Les points ξ0 , ξ1 et ξ2 sont situés quelque
part dans l’intervalle [x0 , x2 ] et sont inconnus (voir les exercices de fin de
chapitre).
Remarque 6.4
Toutes ces formules aux différences sont d’ordre 2. Les mentions avant, centrée
et arrière renvoient au point où l’on calcule la dérivée et aux points utilisés
pour la calculer. Ainsi, la différence avant est évaluée en x0 sur la base des
valeurs situées vers l’avant, soit en x1 et en x2 . La différence arrière fixe la
dérivée en x = x2 avec l’appui des valeurs de la fonction en x0 et en x1 . La
différence centrée, quant à elle, fait intervenir des valeurs situées de part et
d’autre de x1 .
La figure 6.1 illustre les différentes possibilités. Pour les différences
d’ordre 1, on estime la dérivée par la pente du segment de droite joignant
les points (x0 , f (x0 )) et (x1 , f (x1 )). Dans le cas des différences d’ordre 2, on
296 Chapitre 6
Différence avant d’ordre 2 (en x0 )

Différence centrée d’ordre 2 (en x1 )

p2 (x)
Différence arrière d’ordre 2 (en x2 )

Différence arrière (en x1 )

Différence avant (en x1 )

x0 x1 x2

Figure 6.1 – Interprétation géométrique des formules aux différences

détermine un polynôme de degré 2 dont la pente en x0 , en x1 et en x2 donne

respectivement les différences avant, centrée et arrière. !

On peut aussi convenir de toujours évaluer la dérivée en x. Dans ce cas,

on utilise les valeurs de f (x + h) et de f (x + 2h) pour la différence avant et
les valeurs de f (x + h) et de f (x − h) pour la différence centrée. En ce qui
concerne le terme d’erreur, on ne retient que son ordre. Les tableaux suivants
résument la situation.

Formules de différences finies d’ordre 1 pour f ′ (x)

f (x + h) − f (x)
f ′ (x) = + O(h)
h
Différence avant d’ordre 1 (6.12)

f (x) − f (x − h)
f ′ (x) = + O(h)
h
Différence arrière d’ordre 1
 Différentiation et intégration numériques 297

Formules de différences finies d’ordre 2 pour f ′ (x)

−f (x + 2h) + 4f (x + h) − 3f (x)
f ′ (x) = + O(h2 )
2h
Différence avant d’ordre 2

f (x + h) − f (x − h) (6.13)
f ′ (x) = + O(h2 )
2h
Différence centrée d’ordre 2

3f (x) − 4f (x − h) + f (x − 2h)
f ′ (x) = + O(h2 )
2h
Différence arrière d’ordre 2

Exemple 6.5
On tente d’évaluer la dérivée de f (x) = ex en x = 0. La solution exacte est
dans ce cas f ′ (0) = e0 = 1. On peut dès lors comparer ce résultat avec ceux
que l’on obtient par les différentes formules aux différences. Par exemple, la
différence avant d’ordre 1 donne pour h = 0,1 :

′ e0+h − e0 e 0, 1 − 1
f (0) ≃ = = 1,051 709 18
h 0,1

Une valeur plus petite de h conduit à un résultat plus précis. Ainsi, si h = 0,05 :

′ e0,05 − 1
f (0) ≃ = 1,025 4219
0,05

On obtient ainsi une erreur à peu près deux fois plus petite, ce qui confirme
que cette approximation est d’ordre 1. Si l’on utilise cette fois une différence
centrée d’ordre 2, on obtient avec h = 0,05 :

′ e0,05 − e−0,05
f (0) ≃ = 1,000 4167
2(0,05)

qui est un résultat beaucoup plus précis. Avec h = 0,025, on obtient :

′ e0,025 − e−0,025
f (0) ≃ = 1,000 104 18
2(0,025)

soit une erreur à peu près 4 fois plus petite qu’avec h = 0,05. On obtiendrait
des résultats similaires avec les différences avant et arrière d’ordre 2. "
298 Chapitre 6

6.2.2 Dérivées d’ordre supérieur

Avec les dérivées d’ordre supérieur, on agit à peu près de la même manière
qu’avec les dérivées d’ordre 1, c’est-à-dire que l’on dérive un polynôme d’in-
terpolation aussi souvent que nécessaire. Les dérivées d’ordre supérieur posent
toutefois une difficulté supplémentaire, qui provient principalement de l’ana-
lyse d’erreur. En effet, dériver plusieurs fois le terme d’erreur 5.21 est long et
fastidieux. Nous préférons suivre une approche légèrement différente basée sur
le développement de Taylor.
Reprenons le polynôme de degré 2 déjà utilisé pour calculer la dérivée
première. Ce polynôme s’écrit :
p2 (x) = f (x0 ) + f [x0 , x1 ](x − x0 ) + f [x0 , x1 , x2 ](x − x0 )(x − x1 )
et sa dérivée seconde est :
f (x2 ) − 2f (x1 ) + f (x0 )
p′′2 (x) = 2f [x0 , x1 , x2 ] = (6.14)
h2
qui constitue une approximation de la dérivée seconde f ′′ (x) partout dans
l’intervalle [x0 , x2 ]. Il reste à en déterminer l’ordre. Cet ordre dépend du
point retenu pour l’approximation.
– Premier cas : On fait l’approximation de la dérivée en x0 .
L’équation 6.14 peut alors s’écrire :
f (x0 + 2h) − 2f (x0 + h) + f (x0 )
f ′′ (x0 ) ≃ p′′2 (x0 ) =
h2
On remarque immédiatement qu’il s’agit d’une formule aux différences
avant. Pour déterminer l’ordre de l’erreur liée à cette approximation,
on utilise le développement de Taylor 1.20. Dans un premier temps,
on a :
f ′′ (x0 )
f (x0 + 2h) = f (x0 ) + f ′ (x0 )(2h) + (2h)2
2!

f ′′′ (x0 ) f (4) (x0 )

+ (2h)3 + (2h)4 + · · ·
3! 4!
et de même :
f ′′ (x0 ) 2
f (x0 + h) = f (x0 ) + f ′ (x0 )h + h
2!

f ′′′ (x0 ) 3 f (4) (x0 ) 4

+ h + h + ···
3! 4!
On parvient alors à :
f (x0 + 2h) − 2f (x0 + h) + f (x0 ) f ′′ (x0 )h2 + f ′′′ (x0 )h3 + O(h4 )
=
h2 h2

= f ′′ (x0 ) + f ′′′ (x0 )h + O(h2 )

= f ′′ (x0 ) + O(h)
 Différentiation et intégration numériques 299

Cette différence avant est donc une approximation d’ordre 1 de la dé-

rivée seconde. C’est cette approximation que l’on a utilisée pour évaluer
l’erreur d’interpolation 5.21 à l’aide de la formule 5.22 au chapitre pré-
cédent.
– Deuxième cas : On fait l’approximation de la dérivée en x1 .
L’équation 6.14 peut alors s’écrire :

f (x1 + h) − 2f (x1 ) + f (x1 − h)

f ′′ (x1 ) ≃ p′′2 (x1 ) =
h2
qui est une différence centrée. Pour en déterminer l’ordre, on fait appel,
comme dans le cas précédent, aux développements de Taylor, mais cette
fois autour de x1 . On a :

f ′′ (x1 ) 2
f (x1 + h) = f (x1 ) + f ′ (x1 )h + h
2!

f ′′′ (x1 ) 3 f (4) (x1 ) 4

+ h + h + ···
3! 4!
En remplaçant h par (−h), on obtient également :

f ′′ (x1 ) 2
f (x1 − h) = f (x1 ) − f ′ (x1 )h + h
2!

f ′′′ (x1 ) 3 f (4) (x1 ) 4

− h + h + ···
3! 4!
Une fois combinées, ces deux relations deviennent :

f (4) (x1 ) 4
f (x1 + h) − 2f (x1 ) + f (x1 − h) f ′′ (x1 )h2 + 12 h + O(h6 )
=
h2 h2

f (4) (x1 ) 2
= f ′′ (x1 ) + h + O(h4 )
12
= f ′′ (x1 ) + O(h2 )

c’est-à-dire une approximation d’ordre 2 de la dérivée.

– Troisième cas : On fait l’approximation de la dérivée en x2 .
En reprenant un raisonnement similaire à celui du premier cas, on pour-
rait montrer que la relation 6.14 est une approximation d’ordre 1 de la
dérivée seconde en x = x2 .
Remarque 6.6
Il peut sembler surprenant de constater que la même équation aux différences,
obtenue à partir d’un polynôme de degré 2, soit d’ordre 1 en x = x0 et en x = x2
et soit d’ordre 2 en x = x1 . Cela s’explique par la symétrie des différences
centrées, qui permet de gagner un ordre de précision. !
300 Chapitre 6

On peut obtenir toute une série de formules aux différences finies en uti-
lisant des polynômes de degré plus ou moins élevé et en choisissant les déve-
loppements de Taylor appropriés pour en obtenir l’ordre de convergence. Les
tableaux suivants présentent les principales d’entre elles.

Formules de différences finies pour f ′′ (x)

f (x − 2h) − 2f (x − h) + f (x)
f ′′ (x) = + O(h)
h2
Différence arrière d’ordre 1

f (x + 2h) − 2f (x + h) + f (x)
f ′′ (x) = + O(h)
h2
Différence avant d’ordre 1 (6.15)

f (x + h) − 2f (x) + f (x − h)
f ′′ (x) = + O(h2 )
h2
Différence centrée d’ordre 2

−f (x+2h)+16f (x+h)−30f (x)+16f (x−h)−f (x−2h)

f ′′ (x) = 12h2
+ O(h4 )

Différence centrée d’ordre 4

Formule de différences finies pour f (4) (x)

f
(4)
(x) = f (x+2h)−4f (x+h)+6f (x)−4f (x−h)+f (x−2h)
+ O(h2 ) (6.16)
h4

Différence centrée d’ordre 2

Pour terminer, nous démontrons que la différentiation est un procédé nu-

mériquement instable. Toutes les formules de différences finies dépendent d’un
paramètre h qui est la distance entre les points d’interpolation. On pourrait
croire, de façon intuitive, que la précision du résultat augmente à mesure que
diminue la valeur de h. Dans le cas de la différentiation numérique, il y a une
limite aux valeurs de h qui peuvent être utilisées. En effet, si l’on prend, par
exemple, une différence centrée pour estimer la dérivée première, c’est-à-dire :

f (x0 + h) − f (x0 − h)
f ′ (x0 ) ≃
2h
on constate que, lorsque h tend vers 0, le numérateur contient la soustraction
de deux termes très proches l’un de l’autre. Cela résulte en l’élimination par
soustraction (voir la section 1.5.2) de plusieurs chiffres significatifs lorsque h
est trop petit. À quoi s’ajoute une division par un nombre très petit. L’exemple
suivant illustre ce phénomène.
 Différentiation et intégration numériques 301

Exemple 6.7
On considère les différences centrées d’ordre 2 pour le calcul des dérivées pre-
mière et deuxième de la fonction f (x) = ex en x = 0. Ces deux calculs, qui
doivent normalement aboutir à 1, permettent d’apprécier la précision des ré-
sultats. Le tableau suivant rassemble les résultats en simple précision, ce qui
correspond à peu près à travailler avec une mantisse de 7 chiffres décimaux.

Instabilité numérique (simple précision)

f (x+h)−f (x−h) f (x+h)−2f (x)+f (x−h)

h f ′ (x) ≃ 2h f ′′ (x) ≃ h2

10−1 1,175 201 178 1,086 161 137

10−2 1,001 667 619 1,000 839 472
10−3 1,000 016 928 1,000 165 939
10−4 1,000 017 047 1,013 279 080
10−5 1,000 166 059 0,000 000 000
10−6 1,001 358 151 0,000 000 000
10−7 0,983 476 758 −59 604,6601

La valeur de h est successivement réduite d’un facteur de 10 à partir de h = 1,0.

On constate que, lorsque h diminue, la précision liée aux dérivées augmente
dans un premier temps, puis se dégrade brusquement pour les valeurs de h plus
faibles. Cela est particulièrement évident pour la dérivée seconde (troisième
colonne). Si l’on passe en double précision (l’équivalent d’environ 15 chiffres
décimaux dans la mantisse), on observe un comportement similaire, mais qui
se produit à des valeurs de h plus faibles.

Instabilité numérique (double précision)

f (x+h)−f (x−h) f (x+h)−2f (x)+f (x−h)

h f ′ (x) ≃ 2h f ′′ (x) ≃ h2

10−02 1,000 016 666 749 992 12 1,000 008 333 360 558 05
10−03 1,000 000 166 666 681 34 1,000 000 083 406 504 81
10−04 1,000 000 001 666 889 74 1,000 000 005 024 759 28
10−05 1,000 000 000 012 102 32 0,999 998 972 517 346 26
10−06 0,999 999 999 973 244 40 0,999 977 878 279 878 16
10−07 0,999 999 999 473 643 93 0,999 200 722 162 640 44
10−08 0,999 999 993 922 528 80 0,000 000 000 000 000 00
10−09 1,000 000 027 229 219 55 111,022 302 462 515 597
10−10 1,000 000 082 740 370 78 0,000 000 000 000 000 00
10−11 1,000 000 082 740 370 78 0,000 000 000 000 000 00
10−13 0,999 755 833 674 953 35 −11 102 230 246,251 562
10−15 1,054 711 873 393 898 71 111 022 302 462 515,641
10−17 0,000 000 000 000 000 00 0,000 000 000 000 000 00

Lorsque h est trop petit, l’élimination par soustraction des chiffres significatifs
302 Chapitre 6

a un impact dévastateur sur la précision des résultats. Il est donc recommandé

d’être très prudent dans le choix de h et d’éviter des valeurs trop petites. "

6.3 Extrapolation de Richardson

La méthode d’extrapolation de Richardson est valable non seulement pour
la différentiation et l’intégration numériques, mais aussi pour l’interpolation, la
résolution numérique des équations différentielles, etc. Cette technique permet
d’augmenter la précision d’une méthode d’approximation par une technique
d’extrapolation que nous décrivons dans cette section.
Prenons comme point de départ une approximation numérique, notée Qapp (h),
d’une certaine quantité exacte Qexa inconnue. L’approximation numérique dé-
pend d’un paramètre h, comme c’est souvent le cas. Généralement, plus h est
petit, plus l’approximation est précise. On suppose de plus que cette approxi-
mation est d’ordre n, c’est-à-dire :

Qexa = Qapp (h) + O(hn )

La notation O(hn ) signifie en fait que l’on a :

Qexa = Qapp (h) + cn hn + cn+1 hn+1 + cn+2 hn+2 + · · · (6.17)

où les constantes cn dépendent de la méthode numérique utilisée. La tech-

nique d’extrapolation de Richardson consiste à obtenir, à partir de l’approxi-
mation 6.17 d’ordre n, une nouvelle approximation d’ordre au moins (n + 1).
Pour ce faire, il suffit de remplacer h par h2 dans l’équation 6.17, ce qui conduit
à la relation :
( ) ( )n ( )n+1 ( )n+2
h h h h
Qexa = Qapp + cn + cn+1 + cn+2 + · · · (6.18)
2 2 2 2

L’approximation Qapp ( h2 ) est généralement plus précise que Qapp (h). On peut
cependant se servir de ces deux approximations pour en obtenir une nouvelle,
encore plus précise. L’idée consiste à combiner les relations 6.17 et 6.18 de
telle sorte que le terme d’ordre n (cn hn ) disparaisse. Cela est possible si l’on
multiplie l’équation 6.18 par 2n pour obtenir :
( ) ( n+1 ) ( n+2 )
h h h
2n Qexa = 2n Qapp + cn hn + cn+1 + cn+2 + ···
2 2 22
En soustrayant l’énoncé 6.17 de cette dernière relation, on obtient :
( )
n n h 1 3
(2 − 1)Qexa = 2 Qapp − Qapp (h) − cn+1 hn+1 − cn+2 hn+2 + · · ·
2 2 4
d’où :
2n Qapp ( h2 ) − Qapp (h) − 12 cn+1 hn+1 − 34 cn+2 hn+2 + · · ·
Qexa = + (6.19)
(2n − 1) (2n − 1)
 Différentiation et intégration numériques 303

qui s’écrit plus simplement :

2n Qapp ( h2 ) − Qapp (h)

Qexa = n
+ O(hn+1 )
(2 − 1)

L’expression de droite est donc une approximation d’ordre au moins (n+1)

de Qexa . L’extrapolation de Richardson permet donc de gagner au moins un
ordre de convergence. En fait, on peut en gagner davantage si, par exemple, on
a cn+1 = 0 dès le départ. Dans ce cas, la nouvelle approximation est d’ordre
(n+2). Cette situation se produit fréquemment, notamment avec les différences
centrées et la méthode d’intégration dite des trapèzes que nous verrons plus
loin.

Exemple 6.8
On a vu qu’en utilisant une différence avant d’ordre 1 pour calculer la dérivée
de ex en x = 0 on obtient :
e0+h − e0 e 0, 1 − 1
f ′ (0) ≃ = = 1,051 709 18 = Qapp (0,1)
h 0,1
pour h = 0,1 et :

e0,05 − 1
f ′ (0) ≃ = 1,025 4219 = Qapp (0,05)
0,05
pour h = 0,05. On peut maintenant faire le calcul à l’aide de l’équation 6.19
avec n = 1 :
21 Qapp (0,05) − Qapp (0,1)
f ′ (0) ≃
21 − 1
= (2)(1,025 421 9) − 1,051 709 18 = 0,999 134 62

qui est une approximation d’ordre 2 et donc plus précise de f ′ (0). De même,
si l’on utilise une différence centrée d’ordre 2, on obtient pour h = 0,05 :

′ e0,05 − e−0,05
f (0) ≃ = 1,000 4167
2(0,05)

et avec h = 0,025 :

e0,025 − e−0,025
f ′ (0) ≃ = 1,000 104 18
2(0,025)

Dans ce cas, l’extrapolation de Richardson permet de gagner 2 ordres de pré-

cision puisque seules les puissances paires de h apparaissent dans le terme
d’erreur (voir les exercices de fin de chapitre). Plus précisément, on a :

f (x + h) − f (x − h) f ′′′ (x)h2 f (5) (x)h4

= f ′ (x) + + + O(h6 )
2h 3! 5!
304 Chapitre 6

La différence centrée étant d’ordre 2, l’extrapolation de Richardson avec n = 2

donne :
22 Qapp (0,025) − Qapp (0,05)
f ′ (0) ≃
22 − 1
(4)(1,000 104 18) − 1,000 4167
= = 1,000 000 007
3
qui est une approximation d’ordre 4 de la solution exacte. "

Remarque 6.9
Des exemples précédents, on conclut qu’il vaut mieux éviter d’utiliser des va-
leurs de h très petites pour calculer une dérivée à l’aide d’une formule de
différences finies. Il est en effet préférable de choisir une valeur de h pas trop
petite et de faire des extrapolations de Richardson. !

6.4 Intégration numérique

L’intégration numérique est basée principalement sur la relation :
! xn ! xn ! xn
f (x)dx = pn (x)dx + En (x)dx (6.20)
x0 x0 x0

où pn (x) est un polynôme d’interpolation et En (x) est l’erreur qui y est asso-
ciée. En faisant varier la valeur de n, on obtient les formules de Newton-Cotes.
En principe, plus n est élevé, plus grande est la précision liée à la valeur de
l’intégrale recherchée. En pratique cependant, on emploie rarement des valeurs
de n supérieures à 4.
Par ailleurs, l’extrapolation de Richardson, alliée judicieusement à l’une
des formules de Newton-Cotes, conduit à la méthode de Romberg, l’une des
techniques d’intégration numérique les plus précises. Enfin, nous traitons des
quadratures de Gauss-Legendre, très fréquemment utilisées dans les méthodes
numériques plus avancées comme celle des éléments finis (voir Reddy, réf. [31]).

6.4.1 Formules de Newton-Cotes simples et composées

Méthode des trapèzes
Commençons par la méthode la plus simple. On souhaite évaluer :
! x1
f (x)dx
x0

où f (x) est une fonction connue seulement en deux points ou encore une fonc-
tion n’ayant pas de primitive. La solution qui vient tout de suite à l’esprit
consiste à remplacer f (x) par le polynôme de degré 1 passant par les points
(x0 , f (x0 )) et (x1 , f (x1 )) comme l’illustre la figure 6.2.
 Différentiation et intégration numériques 305

f (x)
f (x1 )

f (x0 )

x0 x1 x

Figure 6.2 – Méthode du trapèze

La valeur approximative de l’intégrale correspond à l’aire sous la courbe

du polynôme. Cette aire forme un trapèze qui donne son nom à la méthode du
trapèze. Évidemment, l’approximation est grossière et l’on peut d’ores et déjà
soupçonner que le résultat sera peu précis. Le polynôme de Newton 5.6 et la
relation 5.21 conduisent à :
! x1 ! x1 ! x1
f (x)dx = p1 (x)dx + E1 (x)dx
x0 x0 x0
! x1
= {f (x0 ) + f [x0 , x1 ](x − x0 )}dx
x0
! x1
f ′′ (ξ(x))
+ (x − x0 )(x − x1 )dx
x0 2!

ce qui peut également s’écrire, si l’on intègre le polynôme :

! x1
(x1 − x0 )
f (x)dx = (f (x0 ) + f (x1 ))
x0 2
! (6.21)
x1
f ′′ (ξ(x))
+ (x − x0 )(x − x1 )dx
x0 2!

Le premier terme de droite n’est rien d’autre que l’aire du trapèze de la fi-
gure 6.2, tandis que le deuxième terme est l’erreur commise. Le changement
de variable 5.25 permet d’écrire :
x − x0
s=
h
d’où l’on tire que (x−xi ) = (s−i)h et que dx = hds. Le terme d’erreur devient
306 Chapitre 6

alors :
! x1 ! 1
f ′′ (ξ(x)) f ′′ (ξ(s))
(x − x0 )(x − x1 )dx = s(s − 1)h3 ds
x0 2! 0 2!

On peut encore simplifier cette expression en faisant appel au second théo-

rème de la moyenne.

Théorème 6.10
Soit f1 (x), une fonction continue dans l’intervalle [a , b] et f2 (x), une fonction
intégrable qui ne change pas de signe dans l’intervalle [a , b]. Il existe alors
η ∈ [a , b] tel que :
! b ! b
f1 (x)f2 (x)dx = f1 (η) f2 (x)dx (6.22)
a a

⋆
Comme la fonction (s(s − 1)) ne change pas de signe dans [0 , 1], on peut
mettre à profit ce théorème, ce qui donne :
! 1 ! 1
f ′′ (ξ(s)) 3 f ′′ (η) 3 f ′′ (η) 3
s(s − 1)h ds = h s(s − 1)ds = − h
0 2! 2! 0 12

La méthode du trapèze se résume donc à l’égalité :

! x1
h f ′′ (η) 3
f (x)dx = (f (x0 ) + f (x1 )) − h pour η ∈ [x0 , x1 ] (6.23)
x0 2 12

La méthode du trapèze demeure peu précise, comme en témoigne l’exemple

suivant.

Exemple 6.11
Il s’agit d’évaluer numériquement :
! π
2
sin xdx
0

dont la valeur exacte est 1. La méthode du trapèze donne dans ce cas :

! * π+ π
π π
2
2
sin xdx ≃ sin 0 + sin = = 0,785 398 164
0 2 2 4

qui est une piètre approximation de la valeur exacte 1. Ce résultat peu im-
pressionnant vient du fait que l’on approche la fonction sin x dans l’intervalle
[0 , π2 ] au moyen d’un polynôme de degré 1. Cette approximation est assez
médiocre, comme en témoigne la figure 6.3. "
 Différentiation et intégration numériques 307
1,0
0,9
0,8
0,7
0,6
0,5
0,4

0,3
0,2

0,1
0
π
0 2

Figure 6.3 – Méthode du trapèze, f (x) = sin x

Une meilleure stratégie consiste à décomposer l’intervalle où l’on doit faire

l’intégration, soit l’intervalle [a , b], en n sous-intervalles de longueur (fig. 6.4) :

b−a
h= (6.24)
n
Les différents points engendrés sont notés xi pour i = 0, 1, 2, · · · , n. Les valeurs
aux extrémités sont a = x0 et b = xn . Dans chaque sous-intervalle [xi , xi+1 ],
on peut utiliser la méthode du trapèze. On a alors :

! b $ ! xi+1
n−1 n−1
$ h
f (x)dx = f (x)dx ≃ [f (xi ) + f (xi+1 )]
a xi 2
i=0 i=0

h
= ([f (x0 ) + f (x1 )] + [f (x1 ) + f (x2 )] + · · ·
2
+ [f (xn−2 ) + f (xn−1 )] + [f (xn−1 ) + f (xn )])

On remarque que tous les termes f (xi ) sont répétés deux fois, sauf le premier
et le dernier. On en conclut que :
! b
h
f (x)dx ≃ (f (x0 ) + 2 [f (x1 ) + f (x2 ) + · · · + f (xn−1 )] + f (xn )) (6.25)
a 2

qui est la formule des trapèzes composée. Qu’en est-il du terme d’erreur ? Dans
chacun des n sous-intervalles [xi , xi+1 ], on commet une erreur liée à la méthode
du trapèze. Puisque :

(b − a) (b − a)
h= et donc n =
n h
308 Chapitre 6

f (x)

···

a = x0 x1 x2 x3 · · · xn−2 xn−1 xn = b x

Figure 6.4 – Méthode des trapèzes composée

l’erreur totale commise est :

( ′′ )
f (η) 3 (b − a) f ′′ (η) 3 (b − a) ′′
n − h =− h =− f (η)h2
12 h 12 12

Remarque 6.12
Le raisonnement précédent n’est pas parfaitement rigoureux, même si le ré-
sultat final est juste. En effet, dans chaque sous-intervalle [xi , xi+1 ], l’erreur
liée à la méthode du trapèze simple devrait faire intervenir f ′′ (ηi ), c’est-à-dire
une valeur de η différente pour chaque sous-intervalle. Un autre théorème de
la moyenne est alors nécessaire pour conclure (voir Burden et Faires, réf. [5]).
L’erreur globale étant donnée par :

(b − a) ′′
− f (η)h2 pour η ∈ [a , b] (6.26)
12
la méthode des trapèzes composée est d’ordre 2. !

Exemple 6.13
On reprend le calcul de :
! π
2
I= sin xdx
0

mais cette fois à l’aide de la méthode des trapèzes composée. Soit d’abord
4 intervalles de longueur :

( π2 − 0) π
h= =
4 8
 Différentiation et intégration numériques 309

tels que les montre la figure 6.5a). On a alors :

π ( , - )
8 π π 3π π
I ≃ sin 0 + 2 sin + sin + sin + sin = 0,987 1158
2 8 4 8 2

soit une erreur absolue d’environ 0,012 88 par rapport à la solution exacte. On
constate une nette amélioration en comparaison du résultat obtenu avec un seul
intervalle. Il est intéressant de refaire ce calcul avec 8 intervalles (fig. 6.5b). La
valeur de h est maintenant 16 π
et l’on a :
! π π ( ,
2
16 π π 3π
sin xdx ≃ sin 0 + 2 sin
+ sin + sin
0 2 16 8 16
- )
π 5π 3π 7π π
+ sin + sin + sin + sin + sin
4 16 8 16 2

= 0,996 7852

L’erreur absolue a été réduite à 0,0032. Cette erreur absolue est environ 4 fois
plus petite que l’erreur obtenue avec 4 intervalles, ce qui confirme que cette
méthode est d’ordre 2. On peut de plus utiliser l’extrapolation de Richardson
pour améliorer la précision de ces deux résultats. En utilisant l’équation 6.19
avec n = 2, on obtient l’approximation d’ordre au moins 3 suivante :
! π
2 22 (0,996 7852) − 0,987 1158
sin xdx ≃ = 1,000 008 33
0 22 − 1

ce qui s’approche de plus en plus de la valeur exacte. Comme il sera démontré

un peu plus loin, il s’agit en fait d’une approximation d’ordre 4. "

Remarque 6.14
La méthode du trapèze avec un seul intervalle est également connue sous le
nom de méthode des trapèzes simple. !

Remarque 6.15
La méthode des trapèzes composée est d’ordre 2. La méthode des trapèzes
simple, bien que d’ordre 3, est rarement utilisée, car elle est trop imprécise. !

Remarque 6.16
La méthode des trapèzes composée donne un résultat exact si la fonction f (x)
est un polynôme de degré inférieur ou égal à 1. Cela s’explique par la présence
de la dérivée seconde de f (x) dans le terme d’erreur : celle-ci s’annule dans le
cas de polynômes de degré 1. !
310 Chapitre 6
1,0
0,9

0,8
0,7
0,6
a) 0,5
0,4
0,3
0,2
0,1
0 π π 3π π
0 8 4 8 2

1,0
0,9
0,8
0,7
0,6
b) 0,5
0,4
0,3
0,2
0,1
0 π π 3π π 5π 3π 7π π
0 16 8 16 4 16 8 16 2

Figure 6.5 – Méthode des trapèzes composée, f (x) = sin x (4 et 8 intervalles)

Définition 6.17
Les formules d’intégration numérique sont également appelées formules de
quadrature.

Définition 6.18
Le degré d’exactitude ou encore le degré de précision d’une formule de qua-
drature est le plus grand entier n pour lequel la formule de quadrature est
exacte pour tout polynôme de degré inférieur ou égal à n.

Remarque 6.19
Le degré d’exactitude de la formule des trapèzes est 1. !
 Différentiation et intégration numériques 311

Formule de Simpson 1/3

Reprenons le raisonnement utilisé avec la méthode des trapèzes, mais cette
fois en utilisant un polynôme de degré 2 dont la courbe passe par les points
(x0 ,f (x0 )), (x1 , f (x1 )) et (x2 , f (x2 )). Ce polynôme est donné par la formule de
Newton :

p2 (x) = f (x0 ) + f [x0 , x1 ](x − x0 ) + f [x0 , x1 , x2 ](x − x0 )(x − x1 )

On se sert ensuite de l’approximation :

! x2 ! x2
f (x)dx ≃ p2 (x)dx
x0 x0
! x2
= {f (x0 ) + f [x0 , x1 ](x − x0 ) + f [x0 , x1 , x2 ](x − x0 )(x − x1 )}dx
x0

On se place de nouveau dans le cas où les abscisses sont également distan-

cées. On pose encore (x−xh
0)
= s, ce qui entraîne que (x − xi ) = (s − i)h. La
dernière expression devient :
! 2. /
f (x0 ) + f [x0 , x1 ]hs + f [x0 , x1 , x2 ]h2 s(s − 1) hds
0

h
= (f (x0 ) + 4f (x1 ) + f (x2 ))
3
où l’on a remplacé les différences divisées par leur valeur respective :

f (x1 ) − f (x0 ) f (x2 ) − 2f (x1 ) + f (x0 )

f [x0 , x1 ] = et f [x0 , x1 , x2 ] =
h 2h2
En résumé, on a :
! x2
h
f (x)dx ≃ (f (x0 ) + 4f (x1 ) + f (x2 ))
x0 3

qui est la formule de Simpson 1/3 simple. Cette terminologie est due au facteur
de 13 qui multiplie h.
L’analyse de l’erreur est plus délicate dans ce cas. On s’est vite rendu
compte que la méthode de Simpson 1/3 était plus précise que ce que l’on
escomptait. Une analyse plus fine est donc nécessaire. Cette méthode est basée
sur l’utilisation d’un polynôme de degré 2 et l’on devrait s’attendre à ce que
l’erreur soit donnée par : ! x2
E2 (x)dx
x0
312 Chapitre 6

On peut pousser plus loin l’analyse de l’erreur en introduisant un quatrième

point (x3 , f (x3 )) quelconque et le polynôme de degré 3 correspondant :

(f (x3 ) − p2 (x3 ))
p3 (x) = p2 (x) + (x − x0 )(x − x1 )(x − x2 ) (6.27)
(x3 − x0 )(x3 − x1 )(x3 − x2 )

qui n’est rien d’autre que le polynôme de degré 2 déjà utilisé auquel on ajoute
une correction de degré 3 permettant au polynôme de passer également par le
point (x3 , f (x3 )). Or :
! x2 ! 2
(x − x0 )(x − x1 )(x − x2 )dx = s(s − 1)(s − 2)h4 ds = 0
x0 0

comme on peut le vérifier facilement. Il s’ensuit que :

! x2 ! x2
p2 (x)dx = p3 (x)dx
x0 x0

En utilisant un polynôme de degré 2, on obtient en fait la même précision

qu’avec un polynôme de degré 3. Le terme d’erreur est donc de ce fait :
! x2 ! x2 (4)
f (ξ)
E3 (x)dx = (x − x0 )(x − x1 )(x − x2 )(x − x3 )dx
x0 x0 4!

Il n’est pas possible à ce stade-ci d’appliquer le théorème de la moyenne,

comme nous l’avons fait pour la méthode du trapèze. En effet, la fonction
(x−x0 )(x−x1 )(x−x2 )(x−x3 ) peut changer de signe dans l’intervalle [x0 , x2 ],
à moins de choisir judicieusement x3 . Comme le choix de x3 est arbitraire, on
peut poser x3 = x1 . Le terme d’erreur devient alors :
! x2 ! x2 (4)
f (ξ)
E3 (x)dx = (x − x0 )(x − x1 )(x − x2 )(x − x1 )dx
x0 x0 4!
! 2
f (4) (ξ)
= s(s − 1)2 (s − 2)h5 ds
0 4!

On remarque que la fonction s(s − 1)2 (s − 2) ne change pas de signe dans

l’intervalle [0 , 2]. La figure 6.6 illustre cette fonction.
On peut maintenant se servir du théorème de la moyenne pour obtenir :
! x2 !
f (4) (η) 5 2 f (4) (η) 5
E3 (x)dx = h s(s − 1)2 (s − 2)ds = − h
x0 4! 0 90

La méthode de Simpson 1/3 simple se résume donc à :

! x2
h f (4) (η) 5
f (x)dx = (f (x0 ) + 4f (x1 ) + f (x2 )) − h (6.28)
x0 3 90

où η ∈ [x0 , x2 ].
 Différentiation et intégration numériques 313
0,00

-0,05

-0,10

-0,15

-0,20

-0,25
0 0,2 0,4 0,6 0,8 1,0 1,2 1,4 1,6 1,8 2,0

Figure 6.6 – Fonction s(s − 1)2 (s − 2)

Remarque 6.20
Choisir comme nous l’avons fait x3 = x1 dans l’équation 6.27 n’est pas tout
à fait immédiat. En effet, on constate facilement que le dernier terme de la
relation 6.27 devient singulier. En fait, il faut faire un passage à la limite
comme suit :
lim p3 (x) = lim p2 (x)+
x3 →x1 x3 →x1
(f (x3 ) − p2 (x3 ))(x − x0 )(x − x1 )(x − x2 )
lim
x3 →x1 (x3 − x0 )(x3 − x1 )(x3 − x2 )

= p2 (x)+
(x − x0 )(x − x1 )(x − x2 ) (f (x3 ) − p2 (x3 ))
lim lim
x3 →x1 (x3 − x0 )(x3 − x2 ) x3 →x1 (x3 − x1 )

(x − x0 )(x − x1 )(x − x2 ) (f (x3 ) − p2 (x3 ))

= p2 (x) + lim
(x1 − x0 )(x1 − x2 ) x3 →x1 (x3 − x1 )
Étudions donc maintenant la seule limite qui reste que nous noterons L. On
peut utiliser un développement de Taylor autour de x1 en posant :
x3 = x1 + (x3 − x1 ) = x1 + h′
de sorte que :
(f (x1 + (x3 − x1 )) − p2 (x1 + (x3 − x1 )))
L = lim
x3 →x1 (x3 − x1 )

(f (x1 + h′ ) − p2 (x1 + h′ ))
= lim
h →0′ h′
. / . /
f (x1 ) + f ′ (x1 )h′ + O(h′2 ) − p2 (x1 ) + p′2 (x1 )h′ + O(h′2 )
= lim
h′ →0 h′
314 Chapitre 6

Puisque f (x1 ) = p2 (x1 ), on obtient finalement que L = f ′ (x1 ) − p′2 (x1 ), ce qui
entraîne que :

(x − x0 )(x − x1 )(x − x2 )
lim p3 (x) = p2 (x) + (f ′ (x1 ) − p′2 (x1 ))
x3 →x1 (x1 − x0 )(x1 − x2 )

ce qui montre que p3 (x) est bien défini. !

Remarque 6.21
La valeur de h exprime toujours la distance entre les points xi , c’est-à-dire
qu’elle équivaut dans ce cas à la longueur de l’intervalle divisée par 2. !

La méthode de Simpson 1/3 simple est peu précise, tout comme la méthode
du trapèze, comme en témoigne l’exemple suivant.

Exemple 6.22
On reprend une fois de plus le calcul des exemples précédents. Pour la fonction
f (x) = sin x dans l’intervalle [0 , π2 ], on a :
! * π+
π π
2
4 π
sin xdx ≃ sin 0 + 4 sin + sin = 1,002 2799
0 3 4 2

Ce résultat est plus précis que l’approximation obtenue par la méthode du

trapèze simple, mais il demeure peu satisfaisant. "

On peut encore une fois améliorer la précision de la formule de Simpson

1/3 en la composant. Puisque la méthode simple requiert deux intervalles,
il semble souhaitable de diviser l’intervalle d’intégration [a , b] en 2n sous-
intervalles et d’utiliser la méthode de Simpson 1/3 simple dans chaque paire
de sous-intervalle. La figure 6.7 illustre cette approche. On a alors :

! b $ ! x2i+2
n−1 n−1
$ h
f (x)dx = f (x)dx ≃ (f (x2i ) + 4f (x2i+1 ) + f (x2i+2 ))
a x2i 3
i=0 i=0

h
= ((f (x0 ) + 4f (x1 ) + f (x2 )) + (f (x2 ) + 4f (x3 ) + f (x4 )) + · · ·
3
+ (f (x2n−4 ) + 4f (x2n−3 ) + f (x2n−2 ))

+ (f (x2n−2 ) + 4f (x2n−1 ) + f (x2n )))

h
= (f (x0 ) + 4f (x1 ) + 2f (x2 ) + 4f (x3 ) + 2f (x4 ) + · · ·
3
+ 4f (x2n−3 ) + 2f (x2n−2 ) + 4f (x2n−1 ) + f (x2n ))
 Différentiation et intégration numériques 315

f (x)

...

a = x0 x1 x2 x3 x4 x2n−4 x2n−3 x2n−2 x2n−1 x2n = b x

Figure 6.7 – Méthode de Simpson 1/3 composée

Tous les termes de rang impair sont multipliés par 4 tandis que ceux de rang
pair sont multipliés par 2, sauf le premier (f (x0 )) et le dernier (f (x2n )).
L’analyse de l’erreur liée à la méthode de Simpson 1/3 composée est simi-
laire à celle qui s’applique à la méthode des trapèzes composée. En divisant
[a , b] en 2n intervalles, on utilise n fois la méthode de Simpson 1/3 simple et
l’on commet donc n fois l’erreur liée à cette méthode. On a alors :
b−a b−a
h= et donc n =
2n 2h
et l’erreur totale est :
0 1 0 1
f (4) (η) 5 (b − a) f (4) (η) 5 (b − a) (4)
n − h =− h =− f (η)h4
90 2h 90 180

Remarque 6.23
Le terme d’erreur de la méthode de Simpson 1/3 composée est :

(b − a) (4)
− f (η)h4 pour un certain η ∈ [a , b] (6.29)
180
ce qui en fait une méthode d’ordre 4. De plus, en raison de la présence de la
dérivée quatrième de f (x), cette méthode est exacte dans le cas des polynômes
de degré 3. Le degré d’exactitude de cette méthode est donc 3. !
316 Chapitre 6

Exemple 6.24
On divise l’intervalle [0 , π2 ] en 4 sous-intervalles de longueur h = π8 . On a
alors :
! π π ( )
2 π π 3π π
sin xdx ≃ 8 sin 0 + 4 sin + 2 sin + 4 sin + sin
0 3 8 4 8 2

= 1,000 1346
Pour une quantité de travail similaire, on obtient une précision supérieure à
π
celle de la méthode des trapèzes. Avec 8 sous-intervalles de longueur 16 , on a :
! π π (
2
16 π π 3π
sin xdx ≃ sin 0 + 4 sin + 2 sin + 4 sin
0 3 16 8 16
)
π 5π 3π 7π π
+ 2 sin + 4 sin + 2 sin + 4 sin + sin
4 16 8 16 2

= 1,000 008 296

Cette plus grande précision vient du fait que cette méthode est d’ordre 4. On
constate qu’en passant de 4 à 8 intervalles (c’est-à-dire en divisant h par 2) on
divise l’erreur par un facteur d’environ 16,22, ce qui confirme l’ordre 4 de la
méthode. On peut également utiliser l’extrapolation de Richardson 6.19 avec
n = 4 à partir de ces deux valeurs. On obtient ainsi l’approximation :
24 (1,000 008 296) − 1,000 1346
= 0,999 999 876
24 − 1
qui est d’ordre au moins 5. On verra plus loin qu’elle est en fait d’ordre 6. "

Exemple 6.25
On doit calculer : ! 1
2
e−x dx
0
à l’aide de la méthode de Simpson 1/3 composée avec 8 intervalles de longueur :
1−0 1
=
8 8
2
Comme la fonction e−x n’a pas de primitive, il faut absolument utiliser une
méthode numérique. Dans ce cas :
! 1 1
−x2
(e0 + 4e−0,125 + 2e−0,25 + 4e−0,375 + 2e−0,5
2 2 2 2
8
e dx ≃
0 3

+4e−0,625 + 2e−0,75 + 4e−0,875 + e−1,0 )

2 2 2

= 0,746 826 1205

 Différentiation et intégration numériques 317

Il est intéressant de poursuivre les calculs un peu plus loin et de comparer une
fois de plus les méthodes des trapèzes et de Simpson 1/3 composées. En pre-
nant 64 intervalles et en travaillant en double précision, on obtient les valeurs
suivantes.

21 2
Calcul de 0 e−x dx
Méthode des trapèzes composée 0,746 809 163
Méthode de Simpson 1/3 composée 0,746 824 133
Solution exacte à 9 chiffres 0,746 824 133

ce qui démontre la supériorité de la méthode de Simpson. "

On peut poursuivre dans la même voie et développer des formules de

Newton-Cotes basées sur des polynômes de degré de plus en plus élevé. Nous
ne présentons ci-dessous que les formules de Simpson 3/8 et de Boole sans les
démontrer.

Formule de Simpson 3/8

Si l’on utilise un polynôme de degré 3 dans l’intervalle [x0 , x3 ] et passant
par les points ((xi , f (xi )) pour i = 0, 1, 2, 3), on obtient la formule de Simpson
3/8 simple qui s’écrit :
! x3
3h 3f (4) (η) 5
f (x)dx = (f (x0 ) + 3f (x1 ) + 3f (x2 ) + f (x3 )) − h (6.30)
x0 8 80

pour un certain η ∈ [x0 , x3 ]. On peut également composer cette méthode en

divisant l’intervalle d’intégration [a , b] en 3n sous-intervalles de longueur :

b−a
h=
3n
et en utilisant la formule de Simpson 3/8 simple dans chaque triplet de sous-
intervalle. On obtient alors :
! b $ ! x3i+3
n−1
f (x)dx = f (x)dx
a i=0 x3i

n−1
$ 3h
≃ (f (x3i ) + 3f (x3i+1 ) + 3f (x3i+2 ) + f (x3i+3 ))
8
i=0

3h
= (f (x0 ) + 3f (x1 ) + 3f (x2 ) + 2f (x3 ) + 3f (x4 ) + · · ·
8
+ 2f (x3n−3 ) + 3f (x3n−2 ) + 3f (x3n−1 ) + f (x3n ))
318 Chapitre 6

et le terme d’erreur :
0 1
3f (4) (η) 5 (b − a) 3f (4) (η) 5 (b − a)f (4) (η) 4
n − h =− h =− h
80 3h 80 80

Remarque 6.26
La méthode de Simpson 3/8 composée a le même ordre de convergence (4) et le
même degré d’exactitude (3) que la méthode de Simpson 1/3 composée. Pour
cette raison, on lui préfère souvent la méthode de Simpson 1/3. !

Formule de Boole
Si l’on a au départ un polynôme de degré 4 dans l’intervalle [x0 , x4 ] dont
la courbe passe par les points ((xi , f (xi )) pour i = 0, 1, 2, 3, 4), la formule de
Boole simple s’écrit :
! x4
2h
f (x)dx = (7f (x0 ) + 32f (x1 ) + 12f (x2 ) + 32f (x3 ) + 7f (x4 ))
x0 45

8f (6) (η) 7
− h
945
(6.31)
pour un certain η ∈ [x0 , x4 ]. On compose cette méthode en divisant cette fois
l’intervalle d’intégration [a , b] en 4n sous-intervalles de longueur :
b−a
h=
4n
et en utilisant la formule de Boole simple dans chaque quadruplet de sous-
intervalle. On obtient alors :
! b $ ! x4i+4
n−1
f (x)dx = f (x)dx
a i=0 x4i

n−1
$ 2h
≃ (7f (x4i ) + 32f (x4i+1 ) + 12f (x4i+2 )
45
i=0

+ 32f (x4i+3 ) + 7f (x4i+4 ))

2h
= (7f (x0 ) + 32f (x1 ) + 12f (x2 ) + 32f (x3 ) + 14f (x4 ) + · · ·
45
+ 32f (x4n−5 ) + 14f (x4n−4 ) + 32f (x4n−3 ) + 12f (x4n−2 )

+ 32f (x4n−1 ) + 7f (x4n ))

et le terme d’erreur :
0 1
8f (6) (η) 7 (b − a) 8f (6) (η) 7 2(b − a)f (6) (η) 6
n − h =− h =− h
945 4h 945 945
 Différentiation et intégration numériques 319

Remarque 6.27
En ce qui concerne l’erreur, il se produit un phénomène déjà observé avec la
formule de Simpson 1/3 en ce sens que la formule de Boole conduit à une
approximation d’ordre 6 au lieu de 5. La méthode de Boole a de plus un
degré d’exactitude de 5 puisqu’elle est exacte pour tous les polynômes de degré
inférieur ou égal à 5. !

6.4.2 Méthode de Romberg

La méthode de Romberg est une méthode d’intégration qui permet d’at-
teindre des résultats très précis. Elle est basée sur une utilisation très astucieuse
de la méthode des trapèzes composée (d’ordre 2) et de la technique d’extrapo-
lation de Richardson 6.19. On peut en effet démontrer, sous des hypothèses de
régularité suffisante de la fonction f (x), que le terme d’erreur de la méthode
des trapèzes composée s’écrit :

(b − a) ′′
− f (η)h2 = c2 h2 + c4 h4 + c6 h6 + · · ·
12
où les ci sont des constantes. L’information supplémentaire que l’on tire de
cette relation est que seuls les termes d’ordre pair sont présents. L’absence
des puissances impaires de h permet, du point de vue de l’extrapolation de
Richardson, de gagner deux ordres de convergence à chaque extrapolation.
De plus, les valeurs extrapolées, qui sont d’ordre 4, peuvent à leur tour être
extrapolées pour passer à l’ordre 6, et ainsi de suite. Cette utilisation systé-
matique de l’extrapolation de Richardson permet d’obtenir successivement des
approximations de :
! b
f (x)dx
a

d’ordre 2, 4, 6, 8 et plus. Sur le plan pratique, on obtient généralement des

résultats extrêmement précis.
Dans un premier temps, introduisons quelques notations. On note T1,i le
résultat obtenu à l’aide de la méthode des trapèzes composée avec 2i−1 inter-
valles. Les T1,i sont des approximations d’ordre 2. Pour passer de T1,i à T1,i+1 ,
on doit doubler le nombre de sous-intervalles, ce qui revient à diviser la valeur
de h par deux. Au moyen de l’extrapolation de Richardson 6.19 avec n = 2,
on définit alors :
22 T1,i+1 − T1,i
T2,i = (6.32)
22 − 1
et les T2,i sont des approximations d’ordre 4. On pose ensuite successivement :

24 T2,i+1 − T2,i 26 T3,i+1 − T3,i

T3,i = T4,i =
24 − 1 26 − 1
(6.33)
28 T4,i+1 − T4,i 210 T
5,i+1 − T5,i
T5,i = T6,i =
28 − 1 210 − 1
320 Chapitre 6

et ainsi de suite, ce qui définit un tableau triangulaire de la forme :

Méthode de Romberg
T1,1 T1,2 T1,3 T1,4 T1,5 T1,6 (ordre 2)

T2,1 T2,2 T2,3 T2,4 T2,5 (ordre 4)

T3,1 T3,2 T3,3 T3,4 (ordre 6)

T4,1 T4,2 T4,3 (ordre 8)

T5,1 T5,2 (ordre 10)

T6,1 (ordre 12)

Chaque ligne de ce triangle est de deux ordres de convergence plus précis
que la ligne précédente. La première ligne est tout simplement constituée des
approximations obtenues à l’aide de la méthode des trapèzes composée avec
1, 2, 4, 8, 16 · · · intervalles. Pour passer d’une ligne à l’autre, on utilise l’extra-
polation de Richardson par le biais des relations 6.32 et 6.33.

Remarque 6.28
On peut montrer (voir les exercices de fin de chapitre) que la deuxième ligne
de ce tableau n’est autre que le résultat de la méthode de Simpson 1/3 avec
respectivement 2, 4, 8 · · · intervalles. On pourrait donc éliminer la première
ligne et commencer directement avec la méthode de Simpson. !

Exemple 6.29
On a déjà obtenu, lors de calculs précédents, les valeurs T1,1 = 0,785 3982,
T1,3 = 0,987 1158 et T1,4 = 0,996 7852 correspondant à la formule des trapèzes
composée avec respectivement 1, 4 et 8 intervalles pour évaluer :
! π
2
sin xdx
0
Il est alors possible de remplir la première ligne du tableau en calculant :
π *
π π+
T1,2 = 4 sin 0 + 2 sin + sin = 0,948 0594
2 4 2
On peut ensuite effectuer les différentes extrapolations de Richardson.
! π
2
Méthode de Romberg : sin xdx
0
0,785 3982 0,948 0594 0,987 1158 0,996 7852 (ordre 2)

1,002 2799 1,000 1346 1,000 0083 (ordre 4)

0,999 9916 0,999 9999 (ordre 6)

1,000 0000 (ordre 8)

 Différentiation et intégration numériques 321

La première ligne du tableau étant d’ordre 2, la deuxième ligne est donnée par :
(22 )(0,948 0594) − 0,785 3982
= 1,002 2799
22 − 1

(22 )(0,987 1158) − 0,948 0594

= 1,000 1346
22 − 1

(22 )(0,996 7852) − 0,987 1158

= 1,000 0083
22 − 1
qui sont toutes des approximations d’ordre 4. La troisième ligne devient alors :

(24 )(1,000 1346) − 1,002 2799

= 0,999 9916
24 − 1

(24 )(1,000 0083) − 1,000 1346

= 0,999 9999
24 − 1
d’ordre 6. Puis enfin :
(26 )(0,999 9999) − 0,999 9916
= 1,000 0000
26 − 1
Il en résulte une approximation d’ordre 8 ayant plus de 7 chiffres significatifs.
On remarque que la précision augmente à mesure que l’on se déplace vers le
bas (car l’ordre de l’approximation augmente) et vers la droite sur une même
ligne (car h est divisé par 2 entre chaque valeur). "

Exemple 6.30
Soit une fonction f (x) connue seulement pour quelques valeurs de x.

Valeurs d’une fonction tabulée

x f (x) x f (x)
0,00 0,3989 0,75 0,3011
0,25 0,3867 1,00 0,2420
0,50 0,3521
On tente d’évaluer : ! 1
f (x)dx
0
selon la méthode de Romberg. Puisqu’il y a en tout 5 points, on peut utiliser
la méthode des trapèzes composée avec 1, 2 et 4 intervalles seulement. On a
respectivement :
1
T1,1 = (0,3989 + 0,2420) = 0,320 45
2
1
2
T1,2 = (0,3989 + 2(0,3521) + 0,2420) = 0,336 275
2
1
4
T1,3 = (0,3989 + 2(0,3867 + 0,3521 + 0,3011) + 0,2420) = 0,340 0875
2
322 Chapitre 6

On peut dès lors remplir la première ligne du tableau de la méthode de Rom-

berg.

Méthode de Romberg
0,320 4500 0,336 2750 0,340 0875 (ordre 2)
0,341 5500 0,341 3583 (ordre 4)
0,341 3456 (ordre 6)
Les autres lignes du tableau sont tirées elles aussi des relations 6.32 et 6.33 :
(22 )(0,336 275) − 0,320 450
T2,1 = = 0,341 5500
22 − 1

(22 )(0,340 0875) − 0,336 275

T2,2 = = 0,341 3583
22 − 1

(24 )(0,341 3583) − 0,341 550

T3,1 = = 0,341 3456
24 − 1
On obtient ainsi une approximation d’ordre 6 de l’intégrale. "

Remarque 6.31
Dans le cas d’une fonction connue seulement en certains points, comme dans
l’exemple précédent, le nombre de points doit être de la forme 2n + 1 pour
que la méthode de Romberg puisse s’appliquer. En effet, il faut que le nombre
de sous-intervalles soit une puissance de 2. Dans l’exemple précédent, on avait
22 + 1 points et 4 intervalles. !

6.4.3 Quadratures de Gauss-Legendre

Les quadratures de Gauss-Legendre reposent sur un raisonnement différent
de celui qui est à la base des méthodes de Newton-Cotes. D’une certaine façon,
on cherche à optimiser les schémas d’intégration numérique en choisissant plus
judicieusement les points où est évaluée la fonction f (x). Dans le cas où l’éva-
luation de f (x) est coûteuse en temps de calcul, ces quadratures permettent
d’atteindre une grande précision avec relativement peu d’évaluations de f (x).
Par exemple, la méthode du trapèze requiert l’évaluation de la fonction f (x)
aux deux extrémités de l’intervalle sous la forme :
! b
(b − a)
f (x)dx ≃ (f (a) + f (b))
a 2
Nous avons vu que le degré d’exactitude de cette méthode est 1, car cette
quadrature est exacte dans le cas de tout polynôme de degré inférieur ou égal
à 1. On peut se demander s’il est possible de trouver deux points situés dans
l’intervalle d’intégration ainsi que des coefficients appropriés de telle sorte que
l’expression :
! b
f (x)dx ≃ w1 f (t1 ) + w2 f (t2 )
a
 Différentiation et intégration numériques 323

ait un degré d’exactitude supérieur à celui de la méthode du trapèze. Bien sûr,

si :
(b − a)
w1 = w2 = , t1 = a et t2 = b
2
on retrouve la formule du trapèze. Mais est-ce un choix optimal ?
Pour répondre à cette question, nous allons dans un premier temps nous
restreindre à l’intervalle [−1 , 1], où nous ferons tout le développement. Pour
un intervalle quelconque, il suffira d’effectuer le changement de variable :
(b − a)t + (a + b) (b − a)
x= et dx = dt (6.34)
2 2
qui envoie l’intervalle [−1 , 1] sur un intervalle quelconque [a , b]. En effet, le
changement de variable 6.34 permet d’écrire que :
! b ! 1 ( ) !
(b − a)t + (a + b) (b − a) (b − a) 1
f (x)dx = f dt = g(t)dt
a −1 2 2 2 −1

où : ( )
(b − a)t + (a + b)
g(t) = f
2
Il est donc toujours possible de revenir à l’intervalle [−1 , 1]. De manière
générale, on cherche des expressions de la forme :
! 1 $ n
g(t)dt ≃ wi g(ti ) (6.35)
−1 i=1

dont le degré d’exactitude soit le plus élevé possible.

Définition 6.32
L’expression 6.35 est appelée quadrature de Gauss-Legendre à n points. Les ti
sont appelés points d’intégration, tandis que les coefficients wi sont les poids
d’intégration.

On choisit les points et les poids d’intégration de façon à ce que la quadra-

ture 6.35 soit exacte dans le cas des polynômes de degré le plus élevé possible.
De toute évidence, les points d’intégration ti doivent tous être distincts les
uns des autres et les poids d’intégration doivent être non nuls. Puisque tout
polynôme de degré m peut s’écrire :
m
$
pm (t) = c k tk
k=0

il suffit que la relation 6.35 soit exacte successivement pour les monômes g(t) =
tk , pour k = 0, 1, 2, · · · , m qui constituent une base de l’espace des polynômes
de degré m. On gagne alors à accroître le plus possible le degré m. Le degré
maximal atteint dépend du nombre de points d’intégration n. Puisqu’il y a 2n
coefficients à déterminer dans l’équation 6.35, il est raisonnable de penser que
l’on peut atteindre le degré m = (2n − 1). La valeur de k varie donc entre 0 et
2n − 1.
324 Chapitre 6

Quadrature de Gauss-Legendre à 1 point

Cherchons donc une expression de la forme :

! 1
g(t)dt = w1 g(t1 ) (6.36)
−1

qui soit exacte dans le cas des polynômes de degré le plus élevé possible. Com-
mençons par les polynômes de degré 0. La formule 6.36 doit être exacte pour
g(t) = 1, ce qui donne une première équation :
! 1
1dt = 2 = w1
−1

et l’unique poids d’intégration est déjà déterminé. L’équation 6.35 doit de plus
être exacte pour g(t) = t. On trouve donc :
! 1
tdt = 0 = w1 t1 = 2t1
−1

ce qui entraîne que t1 = 0. Ainsi, la quadrature de Gauss-Legendre à 1 point

s’écrit :
! 1
g(t)dt ≃ 2g(0)
−1

et est exacte pour tout polynôme de degré 1.

Remarque 6.33
La quadrature de Gauss-Legendre à 1 point a le même degré d’exactitude (1)
que la méthode du trapèze, qui est une formule à 2 points. La quadrature de
Gauss-Legendre à 1 point est également connue sous le nom de formule du
point milieu. !

Quadrature de Gauss-Legendre à 2 points

On doit maintenant déterminer les 4 coefficients inconnus de l’expression :

! 1
g(t)dt ≃ w1 g(t1 ) + w2 g(t2 ) (6.37)
−1

On remarque immédiatement que t1 doit être différent de t2 et que les deux

wi doivent être non nuls. Sinon, on se retrouve avec une formule à 1 point. Il
nous faut alors 4 équations qui proviendront de la relation 6.37, où l’on choisit
successivement g(t) = 1, g(t) = t, g(t) = t2 et g(t) = t3 . Les 4 équations
 Différentiation et intégration numériques 325

résultantes sont :
! 1
1dt = 2 = w1 + w2 (6.38)
−1
! 1
tdt = 0 = w1 t1 + w2 t2 (6.39)
−1
! 1
2
t2 dt = = w1 t21 + w2 t22 (6.40)
−1 3
! 1
t3 dt = 0 = w1 t31 + w2 t32 (6.41)
−1

et forment un système non linéaire qu’il est heureusement possible de résoudre

analytiquement. On multiplie l’équation 6.39 par t21 et l’on soustrait du résultat
l’équation 6.41 pour obtenir :
w2 t2 (t21 − t22 ) = 0
Pour que ce produit soit nul, il faut que l’un ou l’autre des facteurs s’annule,
c’est-à-dire :
– w2 = 0.
Cette possibilité doit être écartée, car dans ce cas la formule de Gauss-
Legendre à 2 points 6.37 dégénère en une formule à 1 seul point.
– t2 = 0.
De l’équation 6.39, on tire que w1 = 0 ou t1 = 0, ce qui conduit de
nouveau à une formule à 1 point.
– t21 = t22 .
On en conclut que t1 = −t2 , puisque le cas t1 = t2 conduit encore à une
formule à 1 point.
Cette conclusion permet d’obtenir les poids d’intégration. En effet, en vertu de
l’équation 6.39 :
t1 (w1 − w2 ) = 0
et puisque t1 ne peut être nul, w1 = w2 et la relation 6.38 entraîne que :
w1 = w2 = 1
Enfin, selon l’équation 6.40, on a :
2
= t21 + t22 = t21 + (−t1 )2 = 2t21
3
3 3
1 1
ce qui entraîne que t1 = − et t2 = . La formule de Gauss-Legendre à 2
3 3
points s’écrit donc :
! 1 0 3 1 03 1
1 1
g(t)dt ≃ g − +g
−1 3 3

et est exacte dans le cas des polynômes de degré inférieur ou égal à 3.