Chapitre 5

Mthodes dintgration numrique

Le but
Le but de ce chapitre est daborder le calcul gnral de lintgrale dune fonction f (x) sur un domaine
fini dlimit par des bornes finies a et b (les cas des bornes infinies nest donc pas couvert ici) :
Z b
I=
f (x)dx .
(5.1)
a

Les motivations
Dans certains cas trs limits, une telle intgrale peut tre calcule analytiquement ( la main). Cependant, ce nest que trs rarement possible, et le plus souvent un des cas suivants se prsente :
Le calcul analytique est long, compliqu et rbarbatif
Le rsultat de lintgrale est une fonction complique qui fait appel dautres fonctions elles-mme
longues valuer
Rx
02
Cette intgrale na pas dexpression analytique (par exemple la fonction erreur : Erf (x) = p2 0 e x dx0 ).
Dans tous ces cas, on prfrera calculer numriquement la valeur de lintgrale I.

Le principe
Lide principale est de trouver des mthodes qui permettent de calculer rapidement une valeur approche
Ie de lintgrale calculer :
Ie I
(5.2)
Comme toujours, un programme numrique ninvente rien, et ne fait que procder trs rapidement un calcul
que lon pourrait en principe faire la main. Une mthode bien connue consiste par exemple diviserP
laire
sous la courbe en un grand nombre de petits rectangles daire Iek et de les sommer. Le rsultat Ie = k Iek
est alors une approximation de lintgrale I. Cette approximation est dautant meilleure que la largeur h des
rectangles tend vers 0, cest dire : limh!0 Ie = I. Cette mthode dite des rectangles est un exemple parmi
dautres. Nous le reverrons, mais nous verrons aussi dautres mthodes, plus gnrales et plus performantes.
Pour presque toutes les mthodes (sauf la mthode de Monte-Carlo), lintgrale numrique est calcule
partir de lvaluation de la fonction f (x) en un nombre de point n + 1 distincts : fk = f (xk ), k 2 [0, n]. Elle
scrit alors :
n
X
Ie = (b a)
w k fk
(5.3)
k=0

33

Universit Paul Sabatier 2014-2015

CHAPITRE 5. MTHODES DINTGRATION NUMRIQUE

Dans ce cas, on parle de mthodes de quadrature.

Nous allons voir 4 types de mthodes direntes :
1.1- Les mthodes de Newton-Cotes simples
1.2- Les mthodes de Newton-Cotes composites
1.3- Les mthodes de Gauss-Legendre
1.4- Les mthodes de Monte-Carlo

Performances
La performance dune mthode se juge en comparant
la prcision du rsultat : Celle-ci se caractrise en estimant lerreur entre lapproximation et la valeur
relle de lintgrale :
= I Ie
(5.4)

La valeur de lerreur ne peut pas tre calcule exactement puisquen gnral, on ne connat pas lintgrale I que lon cherche calculer. Cependant, une majoration peut souvent tre estime en tudiant
le dveloppement en srie de Taylor de la fonction f (x).
La rapidit dexcution ncessaire pour atteindre ce rsultat. De manire gnrale, toutes les mthodes
peuvent atteindre de trs grandes prcisions. Cependant, le temps de calcul augmente avec la prcision. Ce temps naugmente pas de la mme manire pour toutes les mthodes si bien que certaines
savrent plus ecaces que dautres. En particulier, le temps de calcul des mthodes de quadrature est
proportionnel au nombre de points o la fonction f (x) est value.

5.1

Lois de Newton-Cotes simples

Comme nous allons le voir, les mthodes de Newton-Cotes simples ne permettent pas, elles-seules,
datteindre des prcisions susantes sur des intervalles [a, b] finis et ne sont donc jamais utilises dans ce
cas. En revanche, elles deviennent prcises lorsque |b a| ! 0, et elles constituent alors la base lmentaire
des mthodes composites prsentes dans la section suivante.

5.1.1

Principe

Le principe gnrale des mthodes de Newton-Cotes simples est dapproximer la fonction f (x) intgrer par un polynme P (x) f (x). Si cette approximation est susamment bonne alors, lintgrale de ce
polynme
Z b
Ie =
P (x)dx
(5.5)
Rb

sera une bonne approximation de I = a f (x)dx. Lavantage est que lon sait calculer analytiquement la
e Dans ces mthodes, on choisit des polynmes de degr p qui concident avec f (x) en p + 1
valeur exacte de I.
points distincts, espacs rgulirement entre les bornes a et b. Ces points sont situs aux positions :
{xk = a + kh, k 2 [0, p]} avec

h=

(5.6)

On a alors 8k 2 [0, p] P (xk ) = fk = f (xk ).

Des polynmes de degrs dirents dfinissent des mthodes direntes aux performances direntes.
Nous allons voir les plus courantes, cest dire les mthodes dordres les plus bas.
34

CHAPITRE 5. MTHODES DINTGRATION NUMRIQUE

5.1.2

Universit Paul Sabatier 2014-2015

Mthode du rectangle (p = 0)

Cette mthode utilise le polynme de degr le plus bas, savoir

le polynme constant :
(5.7)

P0 (x) = f (a) = f0 .

Rb
Lintgrale approche Ie0 = a P0 (x)dx se calcule alors trivialement et donne :
Ie0 = (b a)f0
(5.8)
Il sagit de laire du rectangle.

Cette intgrale numrique ncessite une unique valuation de la fonction f (en x0 = a) et reprsente donc
ce quon peut faire de plus rapide.
Lerreur peut tre estime en utilisant les dveloppements en srie de Taylor ou le thorme des accroissements finis on trouve alors pour h = b a :

9 2 [a, b]

0 =

h2 0
f ()
2

|0 |

c.a.d.

h2
Sup[a,b] (|f 0 |)
2

(5.9)

Dmonstration : Pour calculer lerreur, on peut utiliser le thorme des accroissements finis : 8x 2 [a, b], 9 2 [a, b]
tel que :
a)f 0 ()

f (x) = f (a) + (x

En remplaant dans lexpression de lintgrale et de lerreur, on trouve :

=
=
=
=

I
Z
Z

Ie

P0 (x)) dx =

(x

a)f 0 ()dx = f 0 ()

(f (x)
a
b

a)2

(b
2

f 0 () =

h2 0
f ()
2

(f (x)

f (a)) dx

b a

xdx
0

Lerreur nest pas connue car la valeur de 2 [a, b] reste indtermine. Cependant, on peut la majorer
par la plus grande valeur de la drive sur le domaine considr. Quelques remarques sur cette erreur :
Cette mthode dintgration est exacte pour toutes les fonctions f constantes (dans ce cas 0 = 0
puisque quelles vrifient f 0 = 0). Dans le cas plus gnral cette mthode est dautant plus prcise que
les variations de f sont faibles (f 0 petit).
Plus le domaine [a, b] est petit, plus lerreur est faible. Cette erreur dcroit en h2 .
35

Universit Paul Sabatier 2014-2015

5.1.3

CHAPITRE 5. MTHODES DINTGRATION NUMRIQUE

Mthode du point milieu (p = 0)

Cette mthode utilise galement le polynme constant pour

approximer la fonction f . Cependant, elle exploite mieux les
symtries du problme en choisissant la valeur milieu :

a+b
P00 (x) = f
= f0 .
(5.10)
2
Rb
Lintgrale approche Ie00 = a P0 (x)dx se calcule alors trivialement et donne :
Ie00 = (b a)f0
(5.11)

Il sagit de laire du rectangle. Cette mthode ncessite une unique valuation de la fonction f (en x0 =
(a + b)/2) et correspond donc aussi ce quon peut faire de plus rapide.
Lerreur peut tre estime en utilisant les dveloppements en srie de Taylor, ou le thorme des accroissements finis. On trouve alors pour h = b a :
9 2 [a, b]

00 =

h3 00
f ()
24

c.a.d.

|00 |

h3
Sup[a,b] (|f 00 |)
24

(5.12)

Dmonstration : Pour calculer lerreur, on peut utiliser le thorme des accroissements finis au deuxime ordre :
8x 2 [a, b], 9 2 [a, b] tel que :

2 00
f ()
a+b
a+b
a+b
a+b
f (x) = f
+ x
f0
+ x
2
2
2
2
2
En remplaant dans lexpression de lintgrale et de lerreur, on trouve :

=
=
=
=

I
Z

b
a

Ie

(f (x)
"

P0 (x)) dx =

b
a

f (x)

a+b
a+b
f0
+ x
2
2
a

Z b a
Z b a
2
2
f 00 ()
a+b
f0
xdx +
b a
b a
2
2
Z

a+b
dx
2
2 00 #
f ()
a+b
dx
2
2

3
f 00 () b a
h3 00
x2 dx = 0 +
=
f ()
3
2
3

Lerreur nest pas connue car la valeur de 2 [a, b] reste indtermine. Cependant, on peut la majorer
par la plus grande valeur de la drive seconde sur le domaine considr. Quelques remarques sur cette
erreur :
Du fait des symtries, cette mthode dintgration est exacte pour les fonctions f constante, mais aussi
pour les fonctions anes (dans ce cas 00 = 0 puisquelles vrifient f 00 = 0).
Dans le cas plus gnral, cette mthode est dautant plus prcise que les variations de f sont faibles
(f 00 petit).
Plus le domaine [a, b] est petit, plus lerreur est faible. Cette erreur dcroit en h3 , cest dire plus vite
que lerreur de la mthode prcdente : 00 /0 / h ! 0. Ainsi, pour des domaines [a, b] susamment
petits, la mthode du point milieu est toujours plus prcise que la mthode prcdente.
36

CHAPITRE 5. MTHODES DINTGRATION NUMRIQUE

5.1.4

Universit Paul Sabatier 2014-2015

Mthode du trapze (p = 1)

Pour approximer la fonction f , cette mthode utilise le polynme

dordre 1 (la droite) qui passe par f0 = f (a) et f1 = f (b) :

f0 + f1
f1 f0
a+b
P1 (x) =
+
x
(5.13)
2
b a
2
Rb
Lintgrale approche Ie1 = a P1 (x)dx se calcule alors mathmatiquement ou gomtriquement et donne :
Ie1 = (b

a)

f0 + f1
2

(5.14)

Il sagit de laire du trapze. Cette mthode ncessite deux valuations de la fonction f (en a et en b). Elle
est donc en gros deux fois plus lente que les mthodes prcdentes.
Lerreur peut tre estime en utilisant les dveloppements en srie de Taylor, ou le thorme des accroissements finis. On trouve alors 1 pour h = b a :
h3 00
h3
f ()
c.a.d.
|1 |
Sup[a,b] (|f 00 |)
(5.15)
12
12
Lerreur nest pas connue car la valeur de 2 [a, b] reste indtermine. Cependant, on peut la majorer
par la plus grande valeur de la drive seconde sur le domaine considr. Les remarques sur lerreur sont les
mmes que pour la mthode du point milieu. En prcision, cette mthode est donc quivalente celle du
point milieu (1 00 ), mais elle est deux fois plus lente.
9 2 [a, b]

5.1.5

1 =

Mthode de Simpson simple (p = 2)

Pour approximer la fonction f , cette mthode utilise le polynme de degr 2 (la parabole) qui passe par
les trois points f0 = f (a), f1 = f a+b
et f2 = f (b) :
2
P2 (x) = 2
Lintgrale approche Ie2 =
ment et donne :
Ie2 = (b

Rb
a

a)

f2 2f1 + f0
(x
(x2 x0 )2

x1 ) +

f2
x2

f0
(x
x0

x1 ) + f1

(5.16)

h5
Sup[a,b] |f (4) |
90

(5.18)

P2 (x)dx se calcule alors simplef0 + 4f1 + f2

6

(5.17)

Cette mthode ncessite trois valuations de la fonction f (en

x0 = a, x1 = (a + b)/2 et x2 = b). Elle est donc en gros 3 fois
plus lente que les mthodes 1 point.
Lerreur peut tre estime en utilisant les dveloppements en srie
de Taylor, ou le thorme des accroissements finis.
On trouve alors 2 pour h = (b
9 2 [a, b]

a)/2 :
2 =

h5 (4)
f ()
90

c.a.d.

|2 |

1. En fait, le calcul commence tre un poil compliqu et est laiss la curiosit du lecteur.
2. Cette dmonstration est aussi laisse au lecteur...

37

Universit Paul Sabatier 2014-2015

CHAPITRE 5. MTHODES DINTGRATION NUMRIQUE

Lerreur nest pas connue car la valeur de 2 [a, b] reste indtermine. Cependant, on peut la majorer par
la plus grande valeur de la drive quatrime sur lintervalle considr. Quelques remarques sur cette erreur :
Cette mthode dintgration est exacte pour les fonctions f polynomiales dordre 3 (car elles vrifient
f (4) = 0), ce qui inclut en particulier les fonctions constantes, les fonctions anes, et les paraboles
par exemple. Plus gnralement elle est dautant plus prcise que les variations de f sont faibles (f (4)
petit).
Plus lintervalle [a, b] est petit, plus lerreur est faible. Cette erreur dcroit en h5 lorsque h diminue,
cest dire beaucoup plus rapidement que les mthodes prcdentes : 2 /0 / h3 ! 0, 2 /00 / 2 /1 /
h2 ! 0. Ainsi, pour des intervalles [a, b] susamment petits, la mthode de Simpson est toujours plus
prcise que les mthodes prcdentes.

5.1.6

Mthodes dordres plus levs

Plus gnralement, on peut construire des approximations en utilisant des polynmes dordre quelconque.
Le polynme dordre p passant par p + 1 points rgulirement espacs entre a et b sexprime en fonction des
polynmes de Lagrange Lk (x), k 2 [0, p] :
Pp (x) =

p
X

fk Lk (x)

avec

Lk (x) =

k=0

Lintgrale approche Ie =

Rb
a

p
Y

j=0,j6=k

x
xk

xj
xj

(5.19)

Pp (x)dx peut alors se calculer et donne :

Ie = (b

a)

p
X

w k fk

avec wk =

k=0

Rb
a

Lk (x)dx
b a

(5.20)

Un telle mthode ncessite p + 1 valuations de la fonction f (en xk , k 2 [0, p]). Ainsi, plus le degr est lev,
plus la mthode est lente.
Lerreur peut aussi se calculer (mais a devient franchement compliqu), et on peut montrer que :
si p est impair :

9 2 [a, b]

p =

hp+2 (p+1)
f
()
Cp

c.a.d.

si p est pair :

9 2 [a, b]

p =

hp+3 (p+2)
f
()
Cp

c.a.d.

hp+2
Sup[a,b] f (p+1)
Cp

hp+3
|p |
Sup[a,b] f (p+2)
Cp
|p |

o Cp est un coecient qui dpend de lordre du polynme.

Un telle mthode est exacte pour les polynmes de degr p, car ceux-ci vrifient f (p+1) = 0 (et mme
en fait pour des polynmes de degr p + 1 si p est pair). Sinon, elle est dautant plus prcise que la
fonction varie peu sur le domaine [a, b].
Plus le domaine [a, b] est petit, plus lerreur est faible. Cette erreur dcroit en hp+2 lorsque h diminue
(et hp+3 si p est pair), cest dire beaucoup plus rapidement que les mthodes utilisant des polynmes
de degr plus faible.
On trouve ainsi parfois rfrence aux mthodes suivantes :
38

CHAPITRE 5. MTHODES DINTGRATION NUMRIQUE

Universit Paul Sabatier 2014-2015

Mthode de Simpson 3/8 (degr p = 3) :

Ie =

(b

a)

f0 + 3f1 + 3f2 + f3
8

3 5 (4)
h f ()
80

- Mthode de Boole (degr p = 4) :

Ie =

(b

a)

7f0 + 12f1 + 32f2 + 12f3 + 7f4

90

8 7 (6)
h f ()
945

En pratique cependant, les coecients deviennent de plus en plus grands lorsque lordre du polynme
augmente. De plus, partir du degr p = 7, certains coecients deviennent ngatifs, ce qui ce traduit par
des polynmes fortement oscillants. Les dirences entre des grands nombres posent des problmes darrondi
machine qui rendent trs vite ces mthodes inutilisables. Elles ne sont donc jamais utilises pour p > 7.

5.1.7

Bilan

Plus les mthodes de Newton-Cotes simples sont bases sur des polynmes de degr lev, plus elles
sont lentes et plus elles sont diciles coder, mais plus elles sont prcises. Le plus souvent en pratique, le
domaine total dintgration [a, b] est beaucoup trop grand et la fonction varie trop sur ce domaine pour que
ces mthodes donnent des rsultats satisfaisants. Elles ne sont donc quasiment jamais utilises en tant que
telles.
On subdivise donc le domaine total [a, b] en un grand nombre de petits intervalles sur chacun desquels on
peut appliquer avec succs les mthodes de Newton-Cotes simples. On parle alors de mthode de NewtonCotes composites.

5.2
5.2.1

Lois de Newton-Cotes composites

Principe

Lide est donc de dcouper le domaine total dintgration [a, b] en m intervalles. On approxime alors
laire Iek , k 2 [0, m 1] de chaque intervalle par des mthodes de Newton-Cotes simples, et on en dduit une
39

Universit Paul Sabatier 2014-2015

CHAPITRE 5. MTHODES DINTGRATION NUMRIQUE

approximation de laire totale par une simple somme :

Ie =

m
X1
k=0

Iek

(5.21)

Lorsque m est susamment grand, la largeur (b a)/m des intervalles devient aussi petite que lon veut,
si bien que ces mthodes peuvent atteindre des prcisions aussi grandes que ncessaire, sans pour autant se
heurter au problme des polynmes de grand degr.
En utilisant sur chaque intervalle des mthodes de Newton-Cotes simples de degrs dirents, on obtient
des mthodes composites aux proprits et aux performances direntes. Sur chaque intervalle, une mthode
de degr p + 1 value la fonction intgrer en p + 1 points, ce qui revient le subdiviser en p sous-intervalles.
Si on applique ce principe m intervalles contigus, la mthode de Newton-Cotes composite dfinit n = m p
sous-intervalles au total 3 et ncessite lvaluation de n + 1 points. Plus ce nombre est lev, plus la mthode
est lente mais, en gnral, plus elle est prcise. Les comparaisons de mthodes de dirents degrs se feront
donc nombre total de points n commun, cest dire rapidit quivalente. Les mthodes les plus prcises
seront les plus performantes.
Dans certains cas, ces m intervalles peuvent tre espacs de manire non rgulire pour mieux reprsenter
une zone ou la fonction f (x) varie beaucoup, mais dans cette partie, nous nous limiterons au cas dintervalles
rguliers.

5.2.2

Mthode des rectangles (p = 0, n = m)

La mthode des rectangles composite applique la mthode des rectangles simple (p = 0) sur chacun des
m intervalles. Le nombre total de sous-intervalles est donc n = m. Laire de chaque intervalle vaut :
Ie0,k = (xk+1

Si bien que lintgrale totale vaut :

Ie0

=
=

xk )fk = hfk

h (f0 + f1 + ... + fn
n
X1
h
fk

+ 0 fn )

(5.22)

(5.23)

k=0

Dans cette formule, tous les points ont le mme coecient (1), sauf le dernier point fn qui nest pas utilis
(coecient 0).
Lerreur est simplement la somme de toutes les erreurs :
0

n
X1
k=0
2

=
|0 |

0,k =

n 1
h2 X 0
f (k )
2
k=0

h
nf 0 ()
2
(b a)2
Sup[a,b] (|f 0 |)
2n

o lon a utilis le fait que h = (b a)/n pour la dernire ingalit. A nouveau, cette mthode est exacte
pour les fonctions constantes. Plus gnralement, elle est dautant plus prcise que le nombre de points est
3. pour le cas particulier des mthode dordre 0 (p = 0), on a n = m.

40

CHAPITRE 5. MTHODES DINTGRATION NUMRIQUE

Universit Paul Sabatier 2014-2015

Figure 5.1 Mthode composite des rectangles (p = 0) pour m = 6 intervalles (cest dire n = 6 sousintervalles et n + 1 = 7 points au total).

grand :

1
0 = O
= O(h)
n

(5.24)

Lerreur dcroit comme 1/n. La mthode des rectangles est une mthode dordre 1.

5.2.3

Mthode des trapzes (p = 1, n = m)

La mthode des trapzes composite applique la mthode des trapzes simple (p = 1) sur chacun des m
intervalles. Le nombre total de sous-intervalles est donc nouveau n = m. Chaque intgrale vaut :
Ie1,k = (xk+1

Si bien que lintgrale totale vaut :

Ie1

=
=

xk )

fk + fk+1
2

f0 + 2f1 + ... + 2fn 1 + fn

2
!
n
1
f0 X
fn
h
+
fk +
2
2

(5.25)

(5.26)

k=1

Dans cette formule, les points du bord du domaine ont des coecients dirents (1/2) de tous le points
intrieurs (1).
41

Universit Paul Sabatier 2014-2015

CHAPITRE 5. MTHODES DINTGRATION NUMRIQUE

Figure 5.2 Mthode composite des trapzes (p = 1), pour m = 6 intervalles (cest dire n = 6 sousintervalles et n + 1 = 7 points au total).

Lerreur est simplement la somme de toutes les erreurs :

1

1,0 /2 +

n
X1

1,k + 1,n /2 =

k=1

=
|1 |

h3
12

h3 00
nf ()
12
(b a)3
Sup[a,b] (|f 00 |)
12n2

00

f (0 )/2 +

n
X1

00

00

f (k ) + f (n )/2

k=1

o lo a utilis le fait que h = (b a)/n. A nouveau, cette mthode est exacte pour les fonctions constantes
et anes (et mme les paraboles en fait). Plus gnralement, elle est dautant plus prcise que le nombre de
points est grand :

1
1 = O
= O(h2 )
(5.27)
n2
Lerreur dcroit comme 1/n2 . La mthode des trapzes est une mthode dordre 2.

5.2.4

Mthode de Simpson (p = 2, n = 2m)

La mthode de Simpson composite applique la mthode de Simpson simple (p = 2) sur chacun des m
intervalles. Le nombre total de sous-intervalles est donc cette fois-ci n = 2m (il est forcment pair et le
nombre de points n + 1 est forcment impair !). Chaque intgrale vaut :

42

Ie2,k = (xk+2

xk )

fk + 4fk+1 + fk+2
6

(5.28)

CHAPITRE 5. MTHODES DINTGRATION NUMRIQUE

Universit Paul Sabatier 2014-2015

Si bien que (pour un nombre dintervalles n pair) lintgrale totale vaut :

Ie2

=
=

f0 + 4f1 + 2f2 + 4f3 + ... + 2fn 2 + 4fn 1 + fn

3
0
1
n/2 1
n/2 1
X
X
h@
f0 + 4
f2k+1 + 2
f2k + fn A
3

k=0

(5.29)

k=1

Dans ces formules, il y a 3 coecients dirents : 1/3 pour les points du bord, 4/3 pour le points internes
impairs, et 2/3 pour les points internes pairs.

Figure 5.3 Mthode composite de Simpson (p = 2) pour m = 3 intervalles (cest dire n = 6 sousintervalles et n + 1 = 7 points au total).
nouveau, lerreur est simplement la somme de toutes les erreurs :
2

=
=

|2 |
o lon a utilis le fait que h = (b
est prcise :

...
h5
mf (4) ()
90

(b a)5
(4)
Sup
|f
|
[a,b]
180n4

a)/n et n = 2m. Plus le nombre de points est grand, plus la mthode

2 = O

1
n4

= O h4

(5.30)

Lerreur dcroit comme 1/n4 . La mthode de Simpson est une mthode dordre 4.

5.2.5

Gnralites

Bien videmment, on peut construire des mthodes composites dordres plus levs.
43

Universit Paul Sabatier 2014-2015

CHAPITRE 5. MTHODES DINTGRATION NUMRIQUE

Les mthodes utilisant des polynmes de degr plus lev (tout en restant <8) sont plus prcises (
nombre de points gal). Plus prcisment, lerreur associe lutilisation dun polynme de degr p
dcrot en n (p+1) si p est impair et n (p+2) si p est pair. Autrement dit, lordre des mthodes
associes est p + 1 si p est impair et p + 2 si p est pair.
La prcision de toutes les mthodes de Newton-Cote augmente avec le nombre de points utiliss.
Numriquement, chaque addition gnre une petite erreur darrondi machine (lerreur relative est de
lordre de 10 17 pour des rels doubles prcision). Lorsque lon somme beaucoup de nombres, les erreurs de chaque addition sajoutent et lerreur relative totale augmente. Dans le cas des mthodes de
Newton-Cotes composites, une erreur due aux arrondis machine sajoute donc aux erreurs de troncations tudies prcdemment. Et cette erreur augmente avec le nombre de points. En pratique, on ne
peut donc pas augmenter infiniment le nombre dintervalles des mthodes de Newton-Cotes, et il existe
un nombre de points particulier qui permet une prcision optimale.

Figure 5.4 Erreur des direntes mthodes de Newton-Cotes en fonction du nombre dintervalles utiliss
pour intgrer la fonction f (x) = sin x sur lintervalle [0, /2]. Les points correspondent aux erreurs mesures.
Les droites dcroissantes sont les erreurs maximales thoriques associes. La droite en pointills correspond
aux erreurs darrondi machine.

44

CHAPITRE 5. MTHODES DINTGRATION NUMRIQUE

5.3

Universit Paul Sabatier 2014-2015

Exercices

Dfinir en C une fonction simple intgrer sur un intervalle [a, b], par exemple f (x) = sin x sur [0, /2],
f (x) = e x sur [0, 1] , ou autre (choisir une fonction dont on connait la valeur exacte de lintgrale).
En utilisant la mthode des rectangles et un nombre dintervalles n donns, intgrer numriquement
cette fonction par la mthode des rectangles. Essayer plusieurs valeurs de n et observer lvolution de
lerreur mesure.
Intgrer cette mme fonction avec les mthodes des trapzes et de Simpson (on pourra par exemple
dfinir une fonction pour chaque mthode). Observer qualitativement lvolution de lerreur mesure
pour les direntes mthodes.
Raliser une boucle et enregistrer dans un fichier lerreur des trois mthodes pour toutes les valeurs du
nombre dintervalles n entre 1 et 108 (on pourra multiplier n par 2 chaque itration pour limiter le
temps de calcul et la taille du fichier). Visualiser le rsultat en chelle log-log avec gnuplot.
Aprs avoir calcul analytiquement le maximum thorique des drives successives de f sur lintervalle
considr, tracer galement la majoration thorique de lerreur pour ces trois mthodes. Commenter
Recommencer avecpf (x) = sin x sur lintervalle [0, ], ou encore f (x) = sin2 x sur lintervalle [ , ].
Idem pour f (x) = x sur [0, 1]. Interprter...
Intgrer la fonction f avec la mthode de Gauss (voir section suivante) et comparer lerreur celle des
mthodes de Newton-Cotes composites.

45

Universit Paul Sabatier 2014-2015

5.4
5.4.1

CHAPITRE 5. MTHODES DINTGRATION NUMRIQUE

Mthode de Gauss-Legendre
Principe

Pn
Nous avons vu que les mthodes de Newton-Cotes simples pouvaient scrire Ie = (b a) k=0 wk fk o la
fonction f (x) est value en n + 1 points xk , k 2 [0, n]. Avec ces mthodes, les points taient rgulirement
espacs sur le domaine et on approximait la fonction f (x) pour un polynme de degr n qui concidait avec
f en ces n + 1 points.
Lide de la quadrature de Gauss-Legendre est de gnraliser cette mthode pour des points espacs de
manire non-rgulire sur lintervalle dintgration. Si on ne fixe pas a priori les positions des n + 1 points,
cela laisse n + 1 degrs de liberts supplmentaires et on peut peut choisir les positions de ces points de
manire obtenir une mthode optimale.
Quelle que soit la position des n + 1 points, la mthode qui consiste calculer lintgrale du polynme
de degr n passant par ces points est exacte pour tous les polynmes de degr d n. On peut donc utiliser
les n + 1 degrs de libert qui correspondent la position des points et calculer leur positions de manire
ce que la valeur de lintgrale soit exacte galement pour tous les polynmes de n + 1 degrs de plus, cest
dire pour tous les polynme de degr d 2n + 1 4 . Lerreur de cette mthode est alors beaucoup plus petite
que celle des mthodes de Newton-Cotes simples utilisant le mme nombre de points.
Afin de trouver lexpression des mthodes de Gauss pour dirents nombres de points, on cherche
dterminer les positions des xk et les coecients wk associs de manire ce que cette mthode soit exacte
pour les polynmes 1, x, x2 , ... et x2n+1 .
Le plus souvent, la mthode de Gauss-Legendre est prsente sur un intervalle [ 1, 1]. Une intgration
b a
sur un intervalle plus gnral [a, b] peut nanmoins sobtenir par le changement de variable : x ! a+b
2 + 2 x.

5.4.2

Mthode 1 point (n = 0)

lordre le plus bas (n + 1 = 1 point), on cherche la position de lunique point x0 ainsi que lunique
coecient associ w0 de manire ce que la mthode soit exacte pour tous les polynmes de degr 2n+1 = 1.
En loccurrence, il faut ici quelle soit exacte pour les polynmes P0 = 1 etR P1 = x.
1
Pour le polynme constant : P0 = 1, lintgrale exacte est I0 =
1 dx = 2. Et lintgrale par
1
quadrature scrit :
0
0
X
X
Ie = 2
wk P0 (xk ) = 2
wk 1 = 2w0
(5.31)
k=0

k=0

si bien que Ie = I donne directement : w0 = 1.

R1
Pour le polynme de degr 1 : P1 = x, lintgrale exacte est I1 =
x dx = 0. Et lintgrale par
1
quadrature scrit :
0
0
X
X
Ie1 = 2
wk P1 (xk ) = 2
wk xk = 2w0 x0
(5.32)
k=0

k=0

si bien que Ie = I donne directement : w0 x0 = 0, soit x0 = 0.

Lunique point de quadrature correspond au point milieu et
Ie0 = 2f (0)

(5.33)

La mthode de Gauss-Legendre un point correspond en fait la mthode du point milieu (Newton-Cotes).

4. Cela revient dire que lintgrale de tous les polynmes de degr n 2n + 1 qui concident avec f (x) aux n + 1 points de
cette quadrature possdent la mme intgrale. Lintgrale du polynme de degr n tant exacte, elle lest aussi pour tous ces
autres polynmes de degr plus lev.

46

CHAPITRE 5. MTHODES DINTGRATION NUMRIQUE

5.4.3

Universit Paul Sabatier 2014-2015

Mthode 2 points (n = 1)

lordre suivant (n + 1 = 2 points), on cherche la position des deux points x0 et x1 ainsi que les deux
coecients associs w0 et w1 de manire ce que la mthode soit exacte pour tous les polynmes de degr
2n + 1 = 3. En loccurrence, il faut ici quelle soit exacte pour les polynmes P0 = 1, P1 = x, P1 = x2 et
P1 = x 3 .
R1
Pour le polynme constant : P0 = 1, lintgrale exacte est I0 =
1 dx = 2. Et lintgrale par
1
quadrature scrit :
1
1
X
X
Ie = 2
wk P0 (xk ) = 2
wk 1 = 2(w0 + w1 )
(5.34)
k=0

k=0

si bien que Ie = I donne : w0 + w1 = 1.

R1
Pour le polynme de degr 1 : P1 = x, lintgrale exacte est I1 =
x dx = 0. Et lintgrale par
1
quadrature scrit :
1
1
X
X
Ie = 2
wk P1 (xk ) = 2
wk xk = 2(w0 x0 + w1 x1 )
(5.35)
k=0

k=0

si bien que Ie = I donne : w0 x0 + w1 x1 = 0.

R1
Pour le polynme de degr 2 : P2 = x2 , lintgrale exacte est I2 = 1 x2 dx = 2/3 Et lintgrale par
quadrature scrit :
1
1
X
X
e
I=2
wk P2 (xk ) = 2
wk x2k = 2(w0 x20 + w1 x21 )
(5.36)
k=0

k=0

si bien que Ie = I donne : w0 x20 + w1 x21 = 1/3.

R1
Pour le polynme de degr 3 : P3 = x3 , lintgrale exacte est I3 = 1 x3 dx = 0. Et lintgrale par
quadrature scrit :
1
1
X
X
Ie = 2
wk P3 (xk ) = 2
wk x3k = 2(w0 x30 + w1 x31 )
(5.37)
k=0

k=0

si bien que Ie = I donne : w0 x30 + w1 x31 = 0.

Un peu de travail pour rsoudre ce systme donne finalement :
p
p
w0 = w1 = 1/2, x0 = 1/ 3, x1 = +1/ 3

(5.38)

La mthode de Gauss deux points scrit donc :

5.4.4

Gnralisation

p
p
Ie = f ( 1/ 3) + f (1/ 3)

(5.39)

Principe
Cette mthode peut tre applique un nombre de points n + 1 quelconque. On peut alors montrer
quune bonne approximation de lintgrale calculer est :

o :

Ie = 2

n
X

w k fk

(5.40)

k=0

47

Universit Paul Sabatier 2014-2015

CHAPITRE 5. MTHODES DINTGRATION NUMRIQUE

Figure 5.5 Polynmes de Legendre Pn de degrs n = 1 6 et leurs racines.

les positions xk des points de quadrature sont les racines du polynme de Legendre Pn+1 de degr
n + 1. Les polynmes de Legendre peuvent se calculer par rcurrence avec la relation (voir Fig. 5.5) :
P0

P1

Pn

x
(2n

(n 1)Pn 2 (x)
n
Leurs racines nont pas toujours dexpression analytique, et le plus souvent, elles doivent tre calcules
numriquement.
Les n + 1 coecients associs aux points de quadrature valent respectivement :
8k 2 [0, n]

wk =

1)xPn

1 (x)

(1

0
x2k )Pn+1
(xk )2

(1 x2k )
(n + 1)2 Pn2 (xk )

(5.41)

Pn
On peut montrer que 8n > 0, k=0 wk = 1.
On trouve facilement ces valeurs dans des livres ou sur internet. Quelques exemples sont donns dans la table
5.1.
Integrale sur un intervalle quelconque
Dans le cas o lintgrale nest pas calcule sur [ 1, 1], mais sur un domaine quelconque [a, b], alors la
mthode est la mme avec :
n
X
Ie = (b a)
w k fk
(5.42)
k=0

o les poids wk sont les mmes que prcdemment et la fonction f est value aux points de quadrature
b a
suivants : fk = f a+b
2 + 2 xk (o les xk sont les points de quadrature racines des polynmes de Legendre
calculs sur [ 1, 1]).
48

CHAPITRE 5. MTHODES DINTGRATION NUMRIQUE

n+1
1
2
3
4
5
...

2wk
2
1, 1
5/9, 8/9, 5/9

p
18+ 30
,
36
36
p
p
18+ 30 18
, 36 30
36
p
p
322 13 70 322+13 70 128
,
, 255 ,
900
900
p
p
322+13 70 322 13 70
,
900
900
18

30

....

Universit Paul Sabatier 2014-2015

xk
p0
p
1/
3,
1/
p
p3
3/5, 0,q 3/5

p
p
(3 + 2 6/5)/7,
(3 2 6/5)/7,
q
q
p
p
(3 2 6/5)/7, (3 + 2 6/5)/7
q
q
p
p
1
1
5
+
2
10/7,
2 10/7, 0,
3 q
3q 5
p
p
1
5 2 10/7, 13 5 + 2 10/7
3
...

Table 5.1 Position

xk et coecients wk est points associs pour la mthode de Gauss-Legendre n + 1
Pn
points : Ie = 2 k=0 wk f (xk ).
Erreur
Contrairement aux mthodes de Newton-Cotes qui dcroissent en 1/nd+1 ou 1/nd+2 (o n + 1 est le
nombre de points o f (x) est value et d le degr du polynme utilis), les mthodes de Gauss dcroissent
exponentiellement avec d, ce qui savre extrmement rapide. De plus, on peut montrer que les coecients
wk sont toujours positifs, si bien quil ny a aucune soustraction de grands nombres comme il peut y avoir
pour les mthodes de Newton-Cotes de degr d > 8. Les mthodes de Gauss peuvent donc tre utilises avec
un nombre de points arbitrairement grand. Nanmoins, la mise en place de ces mthodes est beaucoup plus
complexe et il nest pas facile de changer le nombre de points une fois le programme crit. En gnral, un
nombre de points de 4 6 est largement susant pour la plupart des applications.

49

Universit Paul Sabatier 2014-2015

CHAPITRE 5. MTHODES DINTGRATION NUMRIQUE

Figure 5.6 Erreur de la mthode de Gauss en fonction du nombre de utiliss pour intgrer fonction
f (x) = sin x sur lintervalle [0, /2]. Les erreurs de quelques mthodes de Newton-Cotes composites sont
rappeles pour comparaison. La droite en pointills tout en bas correspond aux erreurs darrondi machine.

50