Cours d’Analyse Numérique II

Exemples et applications sous Matlab

Licence3 et Mastères

Slim Chaabane : Professeur à la Faculté des Sciences de Sfax

Année universitaire 2019-2020

Table des matières

1 Résolution des équations et des systèmes d’q́uations non linéaires 4

2 Intégration numérique 5
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.2 Dérivation numérique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.2.1 Approximation de la dérivée première . . . . . . . . . . . . . . . . . 7
2.2.2 Approximation de la dérivée seconde . . . . . . . . . . . . . . . . . 11
2.3 Quelques méthodes usuelles d’ntégration numérique . . . . . . . . . . . . . 13

2
3
Chapitre 1

Résolution des équations et des

systèmes d’q́uations non linéaires

4
Chapitre 2

Intégration numérique

2.1 Introduction
Soient a < b et f : [a, b] −→ R une fonction continue. Désignons par F une primitive de f
Z b
sur [a, b] et par I(f ) = f (x)dx. On a donc :
a

I(f ) = F (b) − F (a).

Dans la plus part des cas et sauf dans des cas très particuliers où on peut calculer explici-
tement une primitive de f , il est impossible de calculer la valeur exacte de I(f ), alors qu’en
pratique le calcul de l’intégral est un outil essentiel et très utile dans plusieurs domaines
comme l’ingénierie, la physique, la mécanique, la finance, ...
L’objet de ce cours est de donner des méthodes numériques permettant de calculer une
valeure approchée de I(f ) sans passer par le calcul de la primitive de f . L’idée de base est
de subdiviser l’intervalle [a, b] en N intervalles [ak , ak+1 ], k ∈ {0, 1, ...., N −Z1} de tailles
ak+1
suffisament petites (a0 = a < a1 < a2 ... < aN −1 < aN = b) et de remplacer f (x)dx
Z ak+1 ak

par Pk (x)dx, où Pk est une fonction qui approche f sur l’intervalle [ak , ak+1 ] et dont
ak
le calcul de sa primitive est simple.

5
 6

Nous savons que par la relation de Chasles, on a :

Z b N
X −1 Z ak+1
I(f ) = f (x)dx = f (x)dx
a k=0 ak

ce qui nous donne donc une approximation IN (f ) de I(f ) par la formule suivante :

N
X −1 Z ak+1
I(f ) ' IN (f ) = Pk (x)dx
k=0 ak

Dans ce cours, nous limitons notre étude uniquement au cas où les fonction Pk sont des
fonctions polynômiales et nous étudions en particulier les 3 types de questions suivants :

• Définir la méthode d’intégration numérique et calculer pour cette méthode l’approxi-

mation IN (f ) de I(f ).

• Etudier la Convergence de cette méthode, c’est à dire, étudier si la limite de IN (f )

quand N tend vers +∞ est égale à I(f ).

• Pour chaque méthode étudiée, déterminer une estimation d’erreur entre IN (f ) et

I(f ) permettant donc de comparer la rapidité de ces différentes méthodes.

Nous commençons ce cours par rappeler quelques techniques de dérivation numérique,

puis nous présentons quelques méthodes usuelles d’integration numérique, notamment la
méthode des rectangles, la méthode des points mileux, la méthode des trapèzes et la
méthode de Simpson.

2.2 Dérivation numérique

L’objet de cette partie est d’approximer une ou plusieurs dérivées d’une fonction donnée
f , supposée suffisament régulière, dans les deux cas suivants :

• le calcul exact des dérivées de f est très difficile voir même impossible

• la fonction f n’est connue qu’en un nombre fini de points.

 7

L’outil principal utilisé est la formule de Taylor-Lagrange que nous rappelons dans le
théorème suivant :

Théorème 2.2.1. Soit f une fonction de classe Cn+1 sur un intervalle contenant a et x.
Alors , il existe ξ entre a et x telle que :

n
X f (k) (a) f (n+1) (ξ)
f (x) = (x − a)k +
k=0
k! (n + 1)!

2.2.1 Approximation de la dérivée première

Soit f : [a, b] −→ R une fonction de classe C2 et M2 = max |f 00 (t)|. Soient x ∈]a, b[ et
t∈[a,b]
h > 0 suffisament petit, alors d’après la formule de Taylor-Lagrange, il existe c ∈ [x, x+h],
telle que :
h2 00
f (x + h) = f (x) + hf 0 (x) + f (c)
2
ce qui nous donne donc la formule d’approximation suivante :

f (x + h) − f (x)
( F1 ) : f 0 (x) = + O(h)
h

et que l’on a l’estimation d’erreure suivante :

f (x + h) − f (x) M2
f 0 (x) − ≤ h
h 2

de la même manière, en utilisant la formule de Taylor-Lagrange sur l’intervalle [x − h, x],

il existe d ∈]x − h, x[ :

h2 00
f (x − h) = f (x) − hf 0 (x) + f (d)
2

ce qui nous donne donc la formule d’approximation suivante :

f (x) − f (x − h)
( F2 ) : f 0 (x) = + O(h)
h
 8

et que que l’on a l’estimation d’erreure suivante :

f (x) − f (x − h) M2
f 0 (x) − ≤ h
h 2

On a alors la proposition suivante :

Proposition 1. Soit f : [a, b] −→ R une fonction de classe C2 et M2 = max |f 00 (t)|.

t∈[a,b]
Soient x ∈]a, b[ et h > 0 suffisament petit, on a alors :

f (x + h) − f (x) M2
f 0 (x) − ≤ h.
h 2

f (x) − f (x − h) M2
f 0 (x) − ≤ h.
h 2
Supposons maintenant que la fonction f est de classe C3 , alors en écrivant la formule de
Taylor-Lagrange, il existe c1 ∈ [x, x + h] et c2 ∈ [x − h, x], telle que :

h2 00 h3
f (x + h) = f (x) + hf 0 (x) + f (x) + f 000 (c1 )
2 6

h2 00 h3
f (x − h) = f (x) − hf 0 (x) + f (x) − f 000 (c2 )
2 6
et par suite,
h3 000
f (x + h) − f (x − h) = 2hf 0 (x) + [f (c1 ) + f 000 (c2 )]
6
ce qui nous donne donc la formule d’approximation suivante :

f (x + h) − f (x − h)
( F3 ) : f 0 (x) = + O(h2 )
2h

ce qui nous donne donc la proposition suivante :

Proposition 2. Soit f : [a, b] −→ R une fonction de classe C3 et M3 = max |f 000 (t)|.

t∈[a,b]
Soient x ∈]a, b[ et h > 0 suffisament petit, on a alors :

f (x + h) − f (x − h) h2
f 0 (x) − ≤ M3
2h 6
 9

Exemple : Calculer avec les trois formules précédentes une approximation de la dérivée
première de la fonction : f (x) = x3 au point x0 = 1 en prenant h = 0.01. Conclure.
Réponse : En prenant h = 0.01, on a donc :
(1.01)3 − 1
• avec la première formule : f 0 (1) ' ' 3.0301
0.01
1 − (0.99)3
• avec la deuxième formule : f 0 (1) ' ' 2.9701
0.01
(1.01)3 − (0.99)3
• avec la troixième formule : f 0 (1) ' ' 3.0001
0.02
Comme f 0 (1) = 3, les deux première formules donnent une erreur de l’ordre de 3 × 10−2
tandis que la troisième formule présente une erreur de l’ordre de 10−4 , ce qui montre bien
l’avantage de la troisième formule d’approximation par rapport aux deux autres formules.
L’execice suivant permettra d’obtenir une autre formule d’approximation de la dérivée
première encore plus précise que ces 3 formules.
Exercice : Soit f : [a, b] −→ R une fonction de classe C5 et M5 = max |f (5) (t)|.
t∈[a,b]
1/ Montrer que l’on a :

[f (x − 2h) − f (x + 2h)] + 8 [f (x + h) − f (x − h)]

( F4 ) : f 0 (x) = + O(h4 ),
12h

puis trouver en fonction de h et M5 une estimation d’erreur.

2/ calculer avec les formules ( Fk ), k = 1, 2, 3, 4, une approximation de la dérivée première
π
de la fonction : f (x) = sin(x) au point x0 = 4
en prenant h = 0.01. Trouver les erreurs
correspondantes et interpréter les résultats obtenus.
Réponse : Pour h > 0 suffisament petit, il existe c3 , c4 , d3 , d4 , telle que :

h2 00 h3 h4 (4) h5 (5)
f (x + h) = f (x) + hf 0 (x) + f (x) + f 000 (x) + f (x) + f (c3 )
2 6 120 120

h2 00 h3 h4 (4) h5 (5)
f (x − h) = f (x) − hf 0 (x) + f (x) − f 000 (x) + f (x) − f (c4 )
2 6 120 120
et par suite,

h3 000 h5  (5)
f (x + h) − f (x − h) = 2hf 0 (x) + f (c3 ) + f (5) (c4 )

f (x) + (2.1)
3 120
 10

de même, on a :

0 4h3 000
2 00 2h4 (4) 4h5 (5)
f (x + 2h) = f (x) + 2hf (x) + 2h f (x) + f (x) + f (x) + f (d3 )
3 15 15

4h3 000 2h4 (4) 4h5 (5)

f (x − 2h) = f (x) − 2hf 0 (x) + 2h2 f 00 (x) − f (x) + f (x) − f (d4 )
3 15 15
et par suite,

8h3 000 4h5  (5)

f (x + 2h) − f (x − 2h) = 4hf 0 (x) + f (d3 ) + f (5) (d4 )

f (x) + (2.2)
3 15

Les deux equations (2.1) et (2.2) nous donnent donc :

h5
[f (x − 2h) − f (x + 2h)] + 8 [f (x + h) − f (x − h)] = 12hf 0 (x) +

15
f (5) (c3 ) + f (5) (c4 )

−4f (5) (d3 ) − 4f (5) (d4 )
(2.3)
on a donc la formule d’approximation suivante :

[f (x − 2h) − f (x + 2h)] + 8 [f (x + h) − f (x − h)]

( F4 ) : f 0 (x) = + O(h4 ),
12h

avec une estimation de l’erreur donnée par :

[f (x − 2h) − f (x + 2h)] + 8 [f (x + h) − f (x − h)] h4

f 0 (x) − ≤ M5
12h 18

π
Pour x0 = 4
et h = 0.01, on obtient avec les 4 formules précédentes :
π f (x0 + 0.01) − f (x0 )
• ( F1 ) : f 0 ( ) ' ' 0.703559491689210
4 0.01
π f (x0 ) − f (x0 − 0.01)
• ( F2 ) : f 0 ( ) ' ' 0.710630500575693
4 0.01
π f (x0 + 0.01) − f (x0 − 0.01)
• ( F3 ) : f 0 ( ) ' ' 0.707094996132451
4 0.02
π [f (x0 − 0.02) − f (x0 + 0.02)] + 8 [f (x0 + 0.01) − f (x0 − 0.01)]
• ( F4 ) : f 0 ( ) ' '
4 12 ∗ 0.01
0.707106780950846.
√
0 π 2
et comme f ( ) = ' 0.707106781186548, on a donc le tableau d’erreur suivant :
4 2
 11

Formule ( Fk ) ( F1 ) ( F2 ) ( F3 ) ( F4 )
Erreur associée 3.55 × 10−3 3.52 × 10−3 1.17 × 10−5 2.35 × 10−10

Le tableau précédent montre que la meilleure méthode d’approximation de la dérivée

première est donnée par la formule ( F4 ) où l’erreur entre la dérivée exacte et la dérivée
approchée est 2.35 × 10−10 , suivis de la formule ( F3 ) avec une erreur de 1.17 × 10−5 , tandis
que les deux premières formules d’approximation ( F1 ) et ( F2 ) sont moins performantes
puisque l’erreur est de l’ordre de 3.5 × 10−3 .

2.2.2 Approximation de la dérivée seconde

Supposons que la fonction f est de classe C4 sur [a, b] et soit M4 = | max |f (4) (t)|, alors
t∈[a,b]
en écrivant la formule de Taylor-Lagrange, il existe c3 ∈ [x, x + h] et c4 ∈ [x − h, x], telle
que :
h2 00 h3 h4
f (x + h) = f (x) + hf 0 (x) + f (x) + f 000 (x) + f (4) (c3 )
2 6 24
h2 00 h3 h4
f (x − h) = f (x) − hf 0 (x) + f (x) − f 000 (x) + f (4) (c4 )
2 6 24
ce qui nous donne donc :

h4  (4)
f (x + h) + f (x − h) = 2f (x) + h2 f 00 (x) + f (c3 ) + f (4) (c4 )

24

et donc,
f (x + h) − 2f (x) + f (x − h) 00 h2 000
= f (x) + [f (c3 ) + f 000 (c4 )]
h2 24
ce qui nous donne donc la proposition suivante :

Proposition 3. Soit f : [a, b] −→ R une fonction de classe C4 et M4 = max |f (4) (t)|.

t∈[a,b]
Soient x ∈]a, b[ et h > 0 suffisament petit, on a alors la formule d’approximation suivante :

f (x + h) − 2f (x) + f (x − h)
( F5 ) : f 00 (x) = + O(h2 )
h2
 12

avec une estimation de l’erreur donnée par :

00 f (x + h) − 2f (x) + f (x − h) h2
f (x) − ≤ M4
h2 12

Exercice : Calculer avec la formule précédente une approximation de la dérivée seconde

de la fonction : f (x) = x3 au point x0 = 1 en prenant h = 0.01. Expliquer pourquoi cette
formule est exacte dans ce cas et même dans le cas plus général où f ∈ R3 [x] ?
Réponse : En utilisant la formule précédente pour h = 0.01 :

1.013 − 2 × 1 + 0.993
f 00 (1) ' =6
10−4

Comme M4 = 0, donc d’après l’estimation d’erreur précédente, on a :

1.013 − 2 × 1 + 0.993
f 00 (1) − ≤0
10−4

donc
1.013 − 2 × 1 + 0.993
f 00 (1) = =6
10−4
et même si f ∈ R3 [x], alors sa dérivée 4eme est nulle, donc la formule précédente est aussi
exacte dans ce cas.
Exercice : Soit f : [a, b] −→ R une fonction de classe C6 et M6 = max |f (6) (t)|.
t∈[a,b]
1/ Montrer la formule d’approximation suivante :

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

( F6 ) : f 00 (x) = 2
+ O(h4 ),
12h

puis, trouver en fonction de h et M6 une estimation d’erreur.

2/ Calculer avec cette nouvelle formule une approximation de la dérivée seconde de la
π
fonction : f (x) = sin(x) au point x0 = en prenant h = 0.01, puis comparer ce résultat
4
avec celle obtenue avec la formule ( F5 ).
Réponse : 1/ Voir TD
2/ En prenant h = 0.01, on a les deux résultats suivants :

• En utilisant la formule ( F5 ) : f 00 ( π4 ) ' −0.707100888648338

 13

• En utilisant la formule ( F6 ) : f 00 ( π4 ) ' −0.707106781104490

√
π 2
comme f 00 ( ) = − ' −0.707106781186548, donc avec la formule ( F5 ), on a une erreur
4 2 √
2
e = | − 0.707100888648338 + | ' 5.89 × 10−6 et avec la formule ( F6 ), on a une erreur
√ 2
2
e0 = |0.707106781104490 + | ' 8.2 × 10−11 , ce qui montre bien l’avantage de la formule
2
d’approximation ( F6 ) par rapport à ( F5 ).

2.3 Quelques méthodes usuelles d’ntégration numérique

Nous nous sommes intéressés dans cette parrtie à donner quelques méthodes numériques
usuelles permettant de calculer :
Z b
I(f ) = f (x)dx,
a

où f : [a, b] −→ une fonction continue.

La méthode des rectangles :

La méthode des rectangles à gaucheZou tout simplement la méthode des rectangles consiste
Z b b
à approcher I(f ) = f (x)dx par P0 (x)dx, où P0 désigne le polynôme d’interpolation
a a
de Lagrange de f au point a.
Un simple calcul nous donne que P0 (x) = f (a) et donc :
Z b
I(f ) ' P0 (x)dx = (b − a)f (a).
a

Géometriquement, la méthode des rectangles à gauche consiste à approcher I(f ) par l’air
du rectangle hacheré en rouge de la figure suivante :
 14

Figure 2.1 – Méthode des rectangles

Remarques 1.

1. Si nous désignons par Q0Zle polynôme d’interpolation de Lagrange de f au point b et

b Z b
nous approchons I(f ) = f (x)dx par Q0 (x)dx = (b − a)f (b), nous obtenons ce
a a
qu’on appelle la méthode des rectangles à droite qui est semblable à la méthode des
rectangles à gauche et présente les mêmes propriétés de calcul, précision et rapidité,
c’est pour cette raison, que nous étudions uniquement la méthode des rectangles à
gauche que l’on appellera dans la suite la méthode des rectangles.

2. Sans précision, la méthode des rectangles désigne la méthode des rectangles à gauche.

3. Nous remarquons que d’autant que la longueur de l’intervalle [a, b] est ”grande” d’au-
Z b
tant que l’erreur entre la valeur exacte de f (x)dx et la valeur approchée est ”gran-
a
de”, c’est pour cette raison qu’il ne faut appliquer cette formule que sur des intervalles
de ”petites” tailles : il s’agit donc de la méthode des rectangles composée.
 15

Méthode des rectangles à gauche Méthode des rectangles à droite

Méthode des rectangles composée :

Le principe est le même que l’exemple ?? du cas d’approximation des fonctions par la
méthode d’interpolation locale. Il s’agit de subdiviser l’intervalle [a, b] en N intervalles
N[−1
Ik = [ak , ak+1 ]; k ∈ {0, 1..., N −1} de tailles suffisament petites, tel que [a, b] = Ik . Sur
k=0
Z ak+1
chaque élément Ik , nous utilisons la méthode des rectangles pour approcher f (x)dx,
ak
c’est à dire : Z ak+1
f (x)dx ' (ak+1 − ak )f (ak ).
ak

Dans ce cours, nous utilisons uniquement une subdivision régulière de l’intervalle [a, b], ce
b−a
qui nous donne donc un maillage régulier de l’intervalle [a, b] de pas h = . Les noeuds
N
ak et les éléments Ik de ce maillage sont donnés par :

ak = a + kh; k ∈ {0, . . . , N } ; Ik = [ak , ak + h]; k ∈ {0, . . . , N − 1}

Sur chaque élément Ik ; k ∈ {0, . . . , N −1}, nous obtenons donc l’approximation suivante :
Z ak+1
f (x)dx ' hf (ak ),
ak
 16

or d’après la relation de Chasles, on a :

Z b N
X −1 Z ak+1
I(f ) = f (x)dx = f (x)dx
a k=0 ak

En faisant donc la sommation de k = 0 à N − 1, nous obtenons une approximation RN (f )

de I(f ) par la formule suivante :

N
X −1
I(f ) ' RN (f ) = h f (ak )dx
k=0

Figure 2.2 – Méthode des rectangles composée

Question : A-t-on lim RN (f ) = I(f ) ?

N →+∞

Estimation d’erreur

Pour tout k ∈ {0, . . . , N − 1}, on désigne par ek l’erreur commise sur l’intervalle Ik :
Z ak+1
ek (f ) = f (x)dx − hf (ak )
ak
 17

et par EN (f ) l’erreur globale sur tout l’intervalle [a, b] :

N
X −1
EN (f ) = I(f ) − RN (f ) = ek (f ).
k=0

Pour tout k ∈ {0, . . . , N − 1}, on a :

Z ak+1 Z ak+1
|ek (f )| = f (x) − f (ak )dx ≤ |f (x) − f (ak )| dx
ak ak

Supposons que la fonction f est de classe C1 sur [a, b] et notons par : M1 = max |f 0 (x)|.
x∈[a,b]
En utilisant donc le théorème des accroissements finis, on obtient :

|f (x)| ≤ M1 |x − ak | = M1 (x − ak ) ∀ x ∈ [ak , ak+1 ],

et par suite : ak+1

(x − ak )2 h2

|ek (f )| ≤ M1 = M1 .
2 ak 2
il en résulte donc que :

N −1
X h2 (b − a)2
|E(f )| ≤ |ek (f )| ≤ N M1 = M1 .
k=0
2 2N

et par suite lim RN (f ) = I(f ) et que l’on a le théorème suivant :

N →+∞

Théorème 2.3.1. Soit f : [a, b] −→ R une fonction de classe C1 sur [a, b] et M1 =

max |f 0 (x)|. On a donc :
x∈[a,b]

(b − a)2
|EN (f )| = |I(f ) − RN (f )| ≤ M1 .
2N

1
Exemple 1. Soit f la fonction définie sur [0, 1] par : f (x) = .
1+x
En prenant ZN = 5, calculer avec la méthode des rectangles composée une valeur approchée
1
de I(f ) = f (x)dx, puis déduire E5 (f ). Indiquer ensuite, comment faut-il choisir N
0
pour calculer I(f ) à 10−4 près avec la méthode des rectangles ?
 18

Réponse :
On a : h = 0.2 et que :
 
5 5 5 5
R5 (f ) = h[f (0)+f (0.2)+f (0.4)+f (0.6)+f (0.8)] = 0.2 1 + + + + ' 0.745634920634921
6 7 8 9
Z 1
Comme f (x)dx = [Log(1 + x)]10 = Log(2), donc :
0

|E5 (f )| = |Log(2) − R5 (f ) ≤ 5.24877 × 10−2 .

1
d’autre part : M1 = max |f 0 (x)| = max = 1, donc, d’après le théorème 2.3.1, on
x∈[0,1] x∈[0,1] (1 + x)2
a:
1
|EN (f )| ≤
2N
Pour calculer I(f ) à 10−4 près, il suffit de chercher le premier entier N qui réalise :
1
≤ 10−4 . Un simple calcul, nous donne que pour N = 5000, et donc : R5000 (f ) ' I(f )
2N
à 10−4 près, ce qui est assez coûteux.

Nous avons développé un code Matlab permettant de calculer avec la méthode des rec-
tangles, le nombre d’itŕations nécessaire pour calculer l’intégrale de n’importe quel fonction
de classe C1 sur un intervalle [a, b] à une précision ε donné ainsi que la valeur de cette
intégrale. Les entrés et les sortie de ce code sont :

• Les entrés : a et b, ε, la fonction f et M1 = sup |f 0 (x)|

x∈[a,b]

• Les sorties
Z b : N : le nombre nécessaire d’itérations et I : une valeur approchée de
I(f ) = f (x)dx à ε près.
a
 19

function [N,l]=rect(a,b,eps,M1)
N=floor(M1*(b-a)^2./(2.*eps))+1;
h=(b-a)/N;
l=f(a);
for i=1:N-1
l=l+f(a+i*h);
end
l=h*l;
end
function y=f(x)
y=1./(1+x);
end

En exécutant ce code, on obtient :

>> [N,l]=rect(0,1,10^[-4],1)
N = 5001
l = 0.693197173060946

Remarques 2.
h2
1. En prenant f (x) = x, nous remarquons par un simple calcul que : ek (f ) = M1 , et
2
(b − a)2
en faisant donc la sommation de k = 0 à N − 1, on obtient : EN (f ) = M1 ,
2N
ce qui montre bien que l’estimation d’erreur obtenue dans le théorème 2.3.1 est op-
timale.

2. Si on remplace le point d’interpolation a par le milieu de l’intervalle [a, b], nous

obtenons une nouvelle méthode dite méthode des points milieux qui fera l’objet de la
section suivante :
 20

La méthode des points milieux :

Cette méthode estZ légérement différente de la méthode des recatngles, elle consiste à
b Z b
approcher I(f ) = f (x)dx par : P˜0 (x)dx, où P˜0 désigne le polynôme d’interpolation
a a
a+b
de Lagrange de f au point : α = qui représente en fait le mileu de l’intervalle [a,b].
2
Comme : P˜0 (x) = f ( a+b
2
), donc :

Z b
a+b
I(f ) ' P˜0 (x)dx = (b − a)f ( ).
a 2

Géometriquement, la méthode des points milieux consiste à approcher I(f ) par l’air du
rectangle hacheré en rouge de la figure suivante :

Figure 2.3 – Méthode des points milieux

Comme pour le cas de la méthode des rectangles, nous remarquons que d’autant que
la longueur de l’intervalle [a, b] est ”grande” d’autant que l’erreur entre la valeur exacte
Z b
de f (x)dx et la valeur approchée est ”grande”, c’est pour cette raison qu’il ne faut
a
appliquer cette formule que sur des intervalles de ”petites” tailles : il s’agit donc de la
méthode des points milieux composée.

Méthode des points milieux composée :

 21

Nous utilisons dans cette partie, une subdivision régulière de l’intervalle [a, b] de pas :
b−a
h= ; N ∈ N∗ et nous concervons les mêmes notations des noeuds ak et des éléments
N
h
Ik que la section précédente. Pour tout k ∈ {0, . . . , N − 1}, on note par : xk = ak + le
2
milieu de l’inetrvalle Ik = [ak , ak+1 ]. Z ak+1
Sur chaque élément Ik ; k ∈ {0, . . . , N − 1}, nous approchons f (x)dx par la formule
ak
suivante : Z ak+1
f (x)dx ' hf (xk),
ak

En faisant donc la sommation de k = 0 à N − 1, nous obtenons une approximation IN (f )

de I(f ) par la formule suivante :

N
X −1
I(f ) ' IN (f ) = h f (xk )dx.
k=0

Figure 2.4 – Méthode des points mileux composée

Question : A-t-on lim IN (f ) = I(f ) ?

N →+∞

Estimation d’erreur
 22

Pour tout k ∈ {0, . . . , N − 1}, on désigne par ek l’erreur commise sur l’intervalle Ik :
Z ak+1
ek (f ) = f (x)dx − hf (xk ),
ak

et par EN (f ) l’erreur globale sur tout l’intervalle [a, b] :

N
X −1
EN (f ) = I(f ) − IN (f ) = ek (f ).
k=0

Supposons que la fonction f est de classe C2 sur [a, b]. En utilisant la formule de Taylor-
Lagrange à l’ordre 2, on déduit que pour tout x ∈ [ak , ak+1 ], il existe ck,x entre xk et x,
telle que :
(x − xk )2 00
f (x) = f (xk ) + (x − xk )f 0 (xk ) + f (ck,x )
2
donc :
ak+1 ak+1
(x − xk )2 00
Z Z
ek (f ) = f (x) − f (xk )dx = (x − xk )f 0 (xk ) + f (ck )dx
ak ak 2

Comme : " #ak+1

ak+1
(x − xk ) 2 h2 h2
Z
(x − xk )dx = = − = 0.
ak 2 2 2
ak

donc :
ak+1
(x − xk )2
Z
|ek (f )| ≤ M2 dx
ak 2
où : M2 = max |f 00 (x)|. Il en résulte donc :
x∈[a,b]

h3 (−h)3 h3
 
|ek (f )| ≤ M2 − =
48 48 24

en faisant donc la sommation de k = 0 à N − 1, on obtient

N −1
X h3 (b − a)3
|EN (f )| ≤ |ek (f )| ≤ N M2 = M2 .
k=0
24 24N 2

et par suite : lim IN (f ) = I(f ) et que l’on a le théorème suivant :

N →+∞
 23

Théorème 2.3.2. Soit f : [a, b] −→ R une fonction de classe C2 sur [a, b] et M2 =

max |f 00 (x)|. On a donc :
x∈[a,b]

(b − a)3
|EN (f )| = |I(f ) − IN (f )| ≤ M2 .
24N 2

1
Exemple 2. Soit f la fonction définie sur [0, 1] par : f (x) = .
1+x
En prenant N = 5, Z calculer avec la méthode des points milieux composée une valeur
1
approchée de I(f ) = f (x)dx, puis déduire E5 (f ). Indiquer ensuite, comment faut-il
0
choisir N pour calculer I(f ) à 10−4 près avec la méthode des points milieux ?

Réponse :
On a : h = 0.2 et que :
 
10 10 10 10 10
I5 (f ) = h[f (0.1)+f (0.3)+f (0.5)+f (0.7)+f (0.9)] = 0.2 + + + + ' 0.69190788
11 13 15 17 19
Z 1
Comme f (x)dx = [Log(1 + x)]10 = Log(2), donc :
0

|E5 (f )| = |Log(2) − I5 (f ) ≤ 1.23929 × 10−3 .

2
d’autre part : M2 = max |f 00 (x)| = max = 2, donc, d’après le théorème 2.3.2,
x∈[0,1] x∈[0,1] (1 + x)3
on a :
1
|EN (f )| ≤
12N 2
Pour calculer I(f ) à 10−4 près, il suffit donc de chercher le premier entier N qui réalise :
1
≤ 10−4 , ou encore :
12N 2
 12
104

100
N≥ = √ ' 28.86
12 2 3

donc pour N = 29, on a bien : I29 (f ) ' I(f ) à 10−4 près.

 24

Comme pour la méthode des rectangles, nous avons développé un code Matlab permet-
tant de calculer avec la méthode des points milieux le nombre d’itŕations nécessaire pour
calculer l’intégrale de n’importe quel fonction de classe C2 sur un intervalle [a, b] à une
précision ε donné ainsi que la valeur de cette intégrale. Les entrés et les sortie de ce code
sont :

• Les entrés : a et b, ε, la fonction f et M2 = max |f 00 (x)|

x∈[a,b]

• Les sorties
Z b : N : le nombre nécessaire d’itérations et I : une valeur approchée de
I(f ) = f (x)dx à ε près.
a

function [N,l]=ptmil(a,b,eps,M2)
N=floor(0.5*(b-a)*sqrt(M2*(b-a)/(6.*eps)))+1;
h=(b-a)/N;
l=f(a+0.5*h);
for i=1:N-1
l=l+f(a+(0.5+i)*h);
end
l=h*l;
end
function y=f(x)
y=1./(1+x);
end

En exécutant ce code, on obtient :

>> [N,l]=ptmil(0,1,10^[-4],2)
N = 29
l = 0.693110032073620

Remarques 3. [ Comparaison des deux méthodes étudiées]

 25

h3
1. En prenant f (x) = x2 , un simple calcul nous donne : ek (f ) = M2 , et en faisant
24
(b − a)3
donc la sommation de k = 0 à N − 1, on obtient : EN (f ) = M2 , ce qui
24N 2
montre bien que l’estimation d’erreur obtenue dans le théorème 2.3.2 est optimale.

cte
2. D’après les théorèmes 2.3.1 et 2.3.2, l’estimation d’erreur EN (f ) est d’ordre N
pour
cte
la méthode des rectangles et N2
pour la méthode des points milieux, ce qui montre
bien que la méthode des points milieux est plus rapide et précise que la méthode des
rectangles.

3. Les résultats numériques obtenus dans les exemples 1 et 2 montrent clairement

l’avantage de la méthode des points milieux par rapport à la méthode des rectangles.

La méthode des trapèzes :

Z b Z b
Cette méthode consiste à approcher I(f ) = f (x)dx par P1 (x)dx, où P1 désigne
a a
le polynôme d’interpolation de Lagrange de f aux points a et b. Un simple calcul nous
donne :
f (b) − f (a)
P1 (x) = f (a) + (x − a)
b−a
et que :
b
[f (a) + f (b)](b − a)
Z
I(f ) ' P1 (x)dx = . (2.4)
a 2
Géometriquement, la méthode des trapèzes consiste à approcher I(f ) par l’air du trapèze
T hacheré en rouge de la figure suivante :
 26

Figure 2.5 – Méthode des trapèzes

Rappelons que l’équation (2.4) n’est autre que l’air A du trapèze obtenue par la formule
classique :
[petite base+grande base] × hauteur
A=
2
et que comme les deux autres méthodes, d’autant que la longueur de l’intervalle [a, b] est
Z b
”grande” d’autant que l’erreur entre la valeur exacte de f (x)dx et la valeur approchée
a
est ”grande”, c’est pour cette raison qu’il ne faut appliquer cette formule que sur des
intervalles de ”petites” tailles : il s’agit donc de la méthode des trapèzes composée.

Méthode des trapèzes composée :

Nous utilisons dans cette partie, une subdivision régulière de l’intervalle [a, b] de pas :
b−a
h= ; N ∈ N∗ et nous concervons les mêmes notations des noeuds ak et des éléments
N
Ik que les deux sections précédentes. Sur chaque élément Ik ; k ∈ {0, . . . , N − 1}, nous
Z ak+1
approchons f (x)dx par la formule (2.4), on a donc :
ak

Z ak+1
[f (ak ) + f (ak+1 )]h
f (x)dx ' ,
ak 2

En faisant donc la sommation de k = 0 à N − 1, nous obtenons une approximation TN (f )

 27

de I(f ) par la formule suivante :

"N −1 N −1
#
h X X
I(f ) ' TN (f ) = f (ak ) + f (ak+1 ) .
2 k=0 k=0

Comme :

N
X −1 N
X −1 N
X −1 N
X −1
f (ak ) = f (a) + f (ak ) et f (ak+1 ) = f (b) + f (aj ),
k=0 k=1 k=0 j=1

" N −1
#
h X
donc : I(f ) ' TN (f ) = f (a) + f (b) + 2 f (ak ) , ou encore :
2 k=1

" N −1
#
f (a) + f (b) X
I(f ) ' TN (f ) = h + f (ak ) .
2 k=1

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

Question : A-t-on lim TN (f ) = I(f ) ?

N →+∞

Estimation d’erreur
 28

Pour tout k ∈ {0, . . . , N − 1}, on désigne par P1k le polynôme d’interpolation de Lagrange
de f aux points (ak , ak+1 ) et par ek l’erreur commise sur l’intervalle Ik :
Z ak+1 Z ak+1
[f (ak ) + f (ak+1 )]h
ek (f ) = f (x)dx − = f (x) − P1k (x)dx
ak 2 ak

Supposons que la fonction f est de classe C2 sur [a, b] et soit M2 = max |f 00 (x)|. D’après
k∈[a,b]
le corollaire ??, pour tout x ∈ [ak , ak+1 ], on a :

(x − ak )(x − ak+1 ) (x − ak )(ak+1 − x)

|f (x) − P1k (x)| ≤ M2 = M2
2 2

donc :
Z ak+1 Z ak+1
M2 M2
|ek (f )| ≤ (x − ak )(ak+1 − ak + ak − x) = h(x − ak ) − (x − ak )2 dx
2 ak 2 ak

et par suite :
M2 h3 h3 h3
 
|ek (f )| ≤ − = M2
2 2 3 12
en faisant donc la sommation de k = 0 à N − 1, L’erreur globale :

N
X −1
EN (f ) = I(f ) − TN (f ) = ek (f ),
k=0

vérifie donc :
N −1
X h3 (b − a)3
|EN (f )| ≤ |ek (f )| ≤ N M2 = M2 .
k=0
12 12N 2

et par suite lim TN (f ) = I(f ) et que l’on a le théorème suivant :

N →+∞

Théorème 2.3.3. Soit f : [a, b] −→ R une fonction de classe C2 sur [a, b] et M2 =

max |f 00 (x)|. On a donc :
x∈[a,b]

(b − a)3
|EN (f )| = |I(f ) − TN (f )| ≤ M2 .
12N 2
 29

1
Exemple 3. Soit f la fonction définie sur [0, 1] par : f (x) = .
1+x
En prenant N = 5, calculer avec la méthode des trapèzes une valeur approchée de I(f ) =
Z 1
f (x)dx, puis déduire E5 (f ). Indiquer ensuite, comment faut-il choisir N pour calculer
0
I(f ) à 10−4 près avec la méthode des trapèzes ?

Réponse :
On a : h = 0.2 et que :
   
f (0) + f (1) 3 1 1 1 1
T5 (f ) = h + f (0.2) + f (0.4) + f (0.6) + f (0.8) = 0.2 + + + +
2 4 1.2 1.4 1.6 1.8
' 0.695634920634921

et par suite :
|E5 (f )| = |Log(2) − T5 (f )| ≤ 2.4877 × 10−3 .

Comme M2 = 2, donc, d’après le théorème 2.3.3, on a :

1
|EN (f )| ≤
6N 2

Pour calculer I(f ) à 10−4 près, il suffit donc de chercher le premier entier N qui réalise :
1
≤ 10−4 , ou encore :
6N 2  4  21
10 100
N≥ = √ ' 40.82
6 6
donc pour N = 41, on a bien T41 (f ) ' I(f ) à 10−4 près.

Comme pour les deux autres méthodes étudiées, nous avons développé un code Matlab
permettant de calculer avec la méthode des trapèzes le nombre d’itŕations nécessaire pour
calculer l’intégrale de n’importe quel fonction de classe C2 sur un intervalle [a, b] à une
précision ε donné ainsi que la valeur de cette intégrale. Les entrés et les sortie de ce code
sont :

• Les entrés : a et b, ε, la fonction f et M2 = max |f 00 (x)|

x∈[a,b]

• Les sorties : N : le nombre nécessaire d’itérations et I : une valeur approchée de

 30

Z b
I(f ) = f (x)dx à ε près.
a

function [N,l]=trapeze(a,b,eps,M2)
N=floor(0.5*(b-a)*sqrt(M2*(b-a)/(3.*eps)))+1;
h=(b-a)/N;
l=0.5*(f(a)+f(b));
for i=1:N-1
l=l+f(a+i*h);
end
l=h*l;
end
function y=f(x)
y=1./(1+x);
end

En exécutant ce code, on obtient :

>> [N,l]=trapeze(0,1,10^[-4],2)
N = 41
l = 0.693184358045877

Remarques 4. [ Comparaison des trois méthodes étudiées]

h3
1. En prenant f (x) = x2 , on remarque que ek (f ) = M2 pour tout k ∈ {0, . . . , N − 1}
12
(b − a)3
et donc : EN (f ) = M2 , ce qui montre bien que l’estimation d’erreur du
12N 2
théorème 2.3.3 est optimale.

2. Dans le cas où la fonction f est de classe C2 sur [a, b], l’estimation d’erreur pour la
cte
méthode des rectangles est de l’ordre de , alors que pour les deux autres méthodes
N
cte
étudiées cette estimation est de l’ordre de 2 , donc, les deux méthodes des points
N
milieux et des trapèzes sont plus rapides et précises que la méthode des rectangles.
 31

(b − a)3 (b − a)3
3. Comme M2 ≤ M 2 et que les estimations d’erreurs obtenus dans les
24N 2 12N 2
théorèmes 2.3.2 et 2.3.3 sont optimaux, donc la méthode des points milieux est plus
rapide et précise que la méthode des trapèzes.

4. Les résultats numériques obtenus dans les exemples 1, 2 et 3 montrent qu’avec un

même nombre de points N=5, on a une erreur de : 5.24877 × 10−2 pour la méthode
des rectangles, 2.4877 × 10−3 pour la méthode des trapèzes et 1.23929 × 10−3 pour la
méthode des points milieux, ce qui confirme bien que la méthode des points milieux
est plus précise que la méthode des trapèzes, laquelle est plus précise que la méthode
des rectangles.

5. Nous avons montré dans les exemples 1, 2 et 3 que pour atteindre une précision de
10−4 , il faut choisir : N = 29 pour la méthode des points milieux, N = 41 pour la
méthode des trapèzes et N = 5000 pour la méthode des rectangles, ce qui confirme
bien que la méthode des points milieux est plus rapide que la méthode des trapèzes,
laquelle est plus rapide que la méthode des rectangles.

La méthode de Simpson :

Nous commençons cette partie par énoncer la proposition suivante :

Proposition 4. Pour tout polynôme P ∈ R3 [X], on a :

b  
(b − a)
Z
a+b
P (x)dx = P (a) + 4P ( ) + P (b) .
a 6 2

Preuve : (voir correction de la série 2)

Z b Z b
La méthode de Simpson consiste à approcher I(f ) = f (x)dx par P2 (x)dx, où P2
a a
a+b
désigne le polynôme d’interpolation de Lagrange de f aux points : a, et b.
2
a+b
Comme f et P2 coïncident aux points a, et b, donc, d’après la proposition 4, on a :
2
b  
(b − a)
Z
a+b
I(f ) ' P2 (x)dx = f (a) + 4f ( ) + f (b) . (2.5)
a 6 2
 32

Géometriquement, la méthode de Simpson consiste à approcher I(f ) par l’air de la para-

a+b
bole P qui coïncide avec f aux points a, et b (voir la partie hacherée en rouge de la
2
figure suivante)

Figure 2.7 – Méthode de Simpson

Méthode de Simpson composée :

Nous utilisons dans cette partie, une subdivision régulière de l’intervalle [a, b] avec les
mêmes notations que les trois sections précédentes.
Sur chaque élément Ik ; k ∈ {0, . . . , N − 1}, nous utilisons la formule (2.5) pour approcher
Z ak+1
f (x)dx. On a donc :
ak

Z ak+1  
h h
f (x)dx ' f (ak ) + 4f (ak + ) + f (ak+1 ) .
ak 6 2

En faisant donc la sommation de k = 0 à N − 1, nous obtenons une approximation SN (f )

de I(f ) par la formule suivante :
"N −1 N −1 N −1
#
h X X X h
I(f ) ' SN (f ) = f (ak ) + f (ak+1 ) + 4 f (ak + ) .
6 k=0 k=0 k=0
2
 33

donc : " −1 −1
N N
#
h X X h
I(f ) ' SN (f ) = f (a) + f (b) + 2 f (ak ) + 4 f (ak + )
6 k=1 k=0
2
ou encore :
" N −1 N −1
#
h f (a) + f (b) X X h
I(f ) ' SN (f ) = + f (ak ) + 2 f (ak + )
3 2 k=1 k=0
2

Remarques 5. D’après l’équation précédente, il est évident que :

TN (f ) + 2IN (f )
SN (f ) = .
3

Concernant l’estimation d’erreur de la méthode de Simpson, nous admettons le théorème

suivant :

Théorème 2.3.4. Soit f : [a, b] −→ R une fonction de classe C4 sur [a, b] et M4 =

max |f (4) (x)|. On a donc :
x∈[a,b]

(b − a)5
|EN (f )| = |I(f ) − SN (f )| ≤ M4 .
2880N 4

1
Exemple 4. Soit f la fonction définie sur [0, 1] par : f (x) = .
1+x
Z 1 prenant N = 5, calculer avec la méthode de Simpson une valeur approchée de I(f ) =
En
f (x)dx, puis déduire E5 (f ). Indiquer ensuite, comment faut-il choisir N pour calculer
0
I(f ) à 10−4 près avec la méthode de Simpson ?

Réponse :
TN (f ) + 2IN (f )
d’après la remarque 5, on a : SN (f ) = , donc :
3

0.695634920634921 + 2 × 0.691907885715935
S5 (f ) ' ' 0.693150230688930
3

et donc :
|E5 (f )| = |Log(2) − S5 (f )| ' 3.05 × 10−6 .
 34

24
Comme f (4) (x) = , donc M4 = 24 et d’après le théorème 2.3.4, on obtient :
(1 + x)5

1
|EN (f )| ≤
120N 4

Pour calculer I(f ) à 10−4 près, il suffit donc de chercher le premier entier N qui réalise :
1
≤ 10−4 , ou encore :
120N 4  1
250 4
N≥ ' 3.02
3
donc pour N = 4, on a bien S4 (f ) ' I(f ) à 10−4 près.

Comme pour les trois autres méthodes étudiées, nous avons développé un code Matlab
permettant de calculer avec la méthode de Simpson le nombre d’itŕations nécessaire pour
calculer l’intégrale de n’importe quel fonction de classe C4 sur un intervalle [a, b] à une
précision ε donné ainsi que la valeur approchée de cette intégrale. Les entrés et les sortie
de ce code sont :

• Les entrés : a et b, ε, la fonction f et M4 = max |f 00 (x)|

x∈[a,b]

• Les sorties
Z b : N : le nombre nécessaire d’itérations et I : une valeur approchée de
I(f ) = f (x)dx à ε près.
a
 35

function [N,l]=simp(a,b,eps,M4)
N=floor(0.5*(b-a)*((M4*(b-a)/(180.*eps)))^(0.25))+1;
h=(b-a)/N;
h=(b-a)/N;
l=f(a)+f(b)+4*f(b-0.5*h);
for i=1:N-1
l=l+2.*f(a+i*h)+4.*f(a+(i-1/2)*h);
end
l=h*(l)/6.;
end
function y=f(x)
y=1./(1.+x);
end

En exécutant ce code, on obtient :

>> [N,l]=simp(0,1,10^[-4],24)
N = 4

Remarques 6. [ Comparaison des quatre méthodes étudiées]

h5
1. En prenant f (x) = x4 , on remarque que ek (f ) = M4 pour tout k ∈ {0, . . . , N −
2880
(b − a)5
1} et donc : EN (f ) = M4 , ce qui montre bien que l’estimation d’erreur du
2880N 4
théorème 2.3.4 est optimale.

2. Dans le cas où la fonction f est de classe C4 sur [a, b], l’estimation d’erreur pour
cte
la méthode de Simpson est de l’ordre de 4 , donc elle est la plus rapide et précise
N
parmis les trois autres méthodes étudiées.

3. Les résultats numériques obtenus dans les exemples 1, 2, 3 et 4 montrent clairement

que la méthode de Simpson est la plus rapide et précise que les trois autres méthodes
étudiées. Le tableau suivant confirme bien cette observation.
 36

Méthode Erreur E5 (f ) Nombre d’itérations nécéssaire pour

atteindre une précision de 10−4
Rectangle 5.24 × 10−2 5000
Trapèze 2.4877 × 10−3 41
Points milieux 1.239 × 10−3 29
Simpson 3.05 × 10−6 4

Remarques 7. A partir du chapitre précédent, nous pouvons construire d’autres méthodes

d’integration
Z anumérique en approchant
Z ak+1 sur chaque élément Ik = [ak , ak+1 ]; k = 0, . . . N −1,
k+1
l’integrale f (x)dx par Pk (x)dx où Pk est un polynôme d’interpolation de f
ak ak
construit à partir d’un certain nombre de points de l’élément Ik . L’augmentation du nombre
de points d’interpolation sur chaque élément Ik permet d’obtenir des formules d’intégration
numérique plus précises et plus performantes que les méthodes usuelles.