Filire : gnie mcanique/matriaux

intgration numrique
Compte rendu du TP
Novembre 2014

Auteurs :

Introduction
Aprs avoir introduit, en cours et TD, les direntes mthodes d'intgration numrique, nous allons
appliquer ces dernires sur machine par le biais du logiciel Matlab.

But du TP
Les prcdentes sances de TP nous ont permis de nous familiariser avec le langage Matlab et
ainsi acqurir les outils de base de son utilisation, ces derniers nous permettront d'crire les algorithmes
correspondant aux mthodes vu en cours/TD en langage Matlab.

f (x)dx
a

Ce TP nous permettra de comparer les mthodes d'intgration numrique savoir

les mthodes des trapzes, Simpson, Romberg, Gauss-Legendre et Gauss-Tchebychev.

Droulement du TP
Ce TP s'est droul de la manire suivante :
Rappels thoriques,
criture des algorithmes des mthodes cites,
Utilisation de la fonction
"Switch...Case",
Application des programmes sur les fonctions proposes.

Matlab

Droulement du TP
Algorithmes des mthodes
Les algorithmes des mthodes prsentes ci-dessous sont crit sous forme de fonction

Matlab.

Mthode du trapze :

L'algorithme se base sur cette formule :

n1

X
h
fi + f (b))
I = (f (a) + 2
2
1

Une simple boucle sut donc pour calculer l'approximation.

Code
:

Matlab

1
2
3
4
5
6
7
8
9
10
11
12
13
14

function [E] = trpz( a,b,n )

g=input('Donnez f : ')
f=inline (g);
h=(ba)/n;
S=(h/2)*(f(a)+f(b));
for(i=1:n1)
S=S+h*f(a+i*h);
end
I=S;
disp('I vaut : ')
I
h=int(sym(g),a,b);
E=abs(hI)
end

Mthode de Simpson :

L'algorithme se base sur cette formule :

I=

n1
X
h
i=0

(f2i + 4 f2i+1 + f2i+2 )

Une simple boucle sut.

Code
:

Matlab

1
2
3
4

function [ E ] = simpson( a,b,m )

g=input('Donnez f : ')
f=inline (g);
n=m/2;

5
6
7
8
9
10
11
12
13
14
15
16
17

h=(ba)/(2*n);
S=0;
for(i=0:n1)
S=S+(h/3)*(f(a+2*i*h)+4*f(a+(2*i+1)*h)+f(a+(2*i+2)*h))
end
I=S;
disp('I vaut : ')
I
h=int(sym(g),a,b);
E=abs(hI);
end

Mthode de Romberg :

Pour cette mthode on doit calculer les lments de la matrice triangulaire de Romberg,
Pour la premire ligne on utilise la formule suivante :
j1

T1,j

2X
1
ba
=
[f (a) + 2
f (xij ) + f (b)]
2j
i=1

O xij = a + i 2ba
j1
Pour le reste des lments :
Tk,j =

4k1 Tk1,j+1 Tk1,j

4k1 1

On construit partir de ces formules le code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

Matlab suivant :

function [ E ] = rmbrg( a,b,n )

g=input('Donnez f : ')
f=inline (g);
T=zeros(n);
T(1,1)=((ba)/2)*(f(a)+f(b));
for(j=2:n)
S=0;
for(i=1:(2^(j1)1))
x=a+i*(ba)/(2^(j1));
S=S+2*f(x);
end
T(1,j)=(ba)/(2^j)*(f(a)+f(b)+S);
end
for(k=2:n)
for(l=1:nk+1)
T(k,l)=(4^(k1)*T(k1,l+1)T(k1,l))/(4^(k1)1);
end
end
T
disp('I vaut : ')
I=T(n,1)
h=int(sym(g),a,b);
E=abs(hI)
end

Mthode de Gauss-Legendre :

La construction de l'algorithme se fait en quatre grandes tapes :

Gnration du polynme de Legendre par la formule ci-aprs et calcule de ses
racines xi :
1
dn
Pn (x) =

[(x2 1)n ]
n! 2n dxn
Gnration des coecients poids par rsolution du systme suivant :
n
X

,kpair

k+1
Ai xki = {0,kimpair

i=1

Changement de variable :

1
x = (b + a + (b a)t)
2
Calcul de l'approximation :
n

Code
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37

Matlab :

baX
b+a ba
I=
Ai f (
+
ti )
2 i=1
2
2

function [ E ] = gaussl( a,b,g )

n=input('donnez le nombre de points pour la methode : ');
x=[1 0 1];% le polynome intervenant dans la formule ( x^21)
if (n==0)
P =[0 0 1];
else
E=conv(x,x) ;%eleve le polynome au carre.
if (n>1)
q=E;
for (i=2 :n1)
q=conv(q,x) ;%eleve le polynome jusqu au degre n.
end
end
if (n==1)
D=x;
else
if (n>1)
D=polyder(q);% derive une premiere fois.
for j=1:n1
D=polyder(D);% derive (n1) autre fois.
end
end
end
P= (1/(factorial(n)*(2^n)))*D;% polynome de degre n.
end
y=roots(P);%racines du polynome de degre n.
% calcul des coefficients
C=zeros(n,n);
for(r=1:n)
for(s=1:n)

38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74

C(r,s)=y(s)^(r1) ;%calcule les valeurs de f(xi) et les mets dans une matrice
end
end
Z=zeros(n,1);
for( k =1 : (n))
if (mod((k1),2)==0)
Z(k,1)=(2/(k));
else
Z(k,1)=0;
end
end
A=zeros(n,1);
A= inv(C)*Z

;%solutions de notre systeme d equation (Ai).

f=inline (g);
I=0;
for(v=1:n)
I=(I+
( (ba)/2)*(f( ((b+a)/2)+((ba)/2)*y(v)) )*A(v));
end

I
format long
h=int(sym(g),a,b);
E=abs(hI);

end

Mthode de Gauss-Tchebychev :

Pour cette mthode on utilise la formule de rcurrence suivante :

Z

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

2i 1
X
f (x)

))
f (cos(

n i=1
2n
1 x2

function [ E ] = gausst(a,b,g)
n=input('Donnez le nombre de points : ')

f=inline (g);
if(a==0)
S=0;
for(i=1:n)
S=(S+(pi/n)*f(cos((2*i1)*pi /(2*n) ) ));
end
I=S/2;
else
S=0;
for(i=1:n)
S=S+(pi/n) *f(cos((2*i1)*pi/(2*n) ));
end
end
disp('I vaut : ')

21
22
23
24

I
h=int(sym(g),a,b);
E=abs(hI);
end

Intgration des mthodes dans un programme "Switch...Case" :

Ce programme nous ore la possibilit de choisir la mthode voulu en lui donnant le
nom de cette dernire en argument, puis le programme se charge d'appeler la fonction
associes la mthode et demande l'utilisateur les donnes ncessaires pour le lancement de celle-ci.
Ci-dessous le code du programme "Switch..Case" :

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34

function [ ] = integration(methode)
a=input('Entrez la borne inferieure : ');
b=input('Entrez la borne superieur : ');
switch (methode)
case 'trapeze'
n=input('Entrez le nombre de sub : ');
trpz(a,b,n);
case'simpson'
n=input('Entrez le nombre de sub : ');
simpson(a,b,n);
case 'romberg'
n=input('Entrez le nombre de sub : ');
rmbrg(a,b,n)
case'gausslegendre'
g=input('donnez f : ');
k=zeros(5,1);
for(i=1:5)
k(i)=gaussl(a,b,g);
end
v=[2;4;8;12;16];
plot(v,k,'b')
title('Erreur en fonction du nombre de points')
case'gausstchebychev'
g=input('Donnez la fonction f(x) qui est au nominateur ponderee par racine de (1x^2) ') ;
k=zeros(3,1);
for(i=1:3)
k(i)=gaussl(a,b,g);
end
v=[2;3;4];
plot(v,k,'b')
title('Erreur en fonction du nombre de points')
end
end

Application des programmes

Gauss Legendre

On introduit la commande suivante : integration('gausslegendre').

f (x) = sin(x) avec x [0; /2] :

Voici le rsultat obtenu aprs excution de la commande :

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33

>> integration('gausslegendre')
Entrez la borne inferieure : 0
Entrez la borne superieur : pi/2
donnez f : 'sin(x)'
donnez le nombre de points pour la methode : 2
I =
0.998472613404115
donnez le nombre de points pour la methode : 4
I =
0.999999977197114
donnez le nombre de points pour la methode : 8
I =
1
donnez le nombre de points pour la methode : 12
I =
1.000000000000398
donnez le nombre de points pour la methode : 16
I =
0.999999999996684

Figure 1  Reprsentation de l'erreur en fonction du nombre de points

h(x) = x3 avec x [0; 2] :

Voici le rsultat obtenu aprs excution de la commande :

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33

>> integration('gausslegendre')
Entrez la borne inferieure : 0
Entrez la borne superieur : 2
donnez f : 'x^3'
donnez le nombre de points pour la methode : 2
I =
3.999999999999999
donnez le nombre de points pour la methode : 4
I =
3.999999999999996
donnez le nombre de points pour la methode : 8
I =
4.000000000000002
donnez le nombre de points pour la methode : 12
I =
4.000000000000729
donnez le nombre de points pour la methode : 16
I =
3.999999999968925

g(x) = 1/x2 avec x [1; 2] :

Voici le rsultat obtenu aprs excution de la commande :

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

>> integration('gausslegendre')
Entrez la borne inferieure : 1
Entrez la borne superieur : 2
donnez f : '1/(x^2)'
donnez le nombre de points pour la methode : 2
I =
0.497041420118343
donnez le nombre de points pour la methode : 4
I =
0.499995147562620
donnez le nombre de points pour la methode : 8
I =
0.499999999992918

22
23
24
25
26
27
28
29
30
31
32
33

donnez le nombre de points pour la methode : 12

I =
0.500000000000201
donnez le nombre de points pour la methode : 16
I =
0.499999999997659

Figure 2  Reprsentation de l'erreur en fonction du nombre de points

S(x) = 1/(1 + x2 ) avec x [0; 1] :

Voici le rsultat obtenu aprs excution de la commande :

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

>> integration('gausslegendre')
Entrez la borne inferieure : 0
Entrez la borne superieur : 1
donnez f : '1/(1+x^2)'
donnez le nombre de points pour la methode : 2
I =
0.786885245901639
donnez le nombre de points pour la methode : 4
I =
0.785402976311451
donnez le nombre de points pour la methode : 8
I =
0.785398163379778

23
24
25
26
27
28
29
30
31
32
33

donnez le nombre de points pour la methode : 12

I =
0.785398163397791
donnez le nombre de points pour la methode : 16
I =
0.785398163394261

Gauss Tchebychev

On introduit la commande suivante pour les fonctions suivantes : integration('gausstchebychev').

1
f (x) = 1x
2 :
Voici le rsultat obtenu aprs excution de la commande :

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

>> integration('gausstchebychev')
Entrez la borne inferieure : 0
Entrez la borne superieur : 1
Donnez la fonction f(x) qui est au nominateur ponderee par racine de (1x^2) '1'
donnez le nombre de points pour la methode : 2
I =
1
donnez le nombre de points pour la methode : 3
I =
1.000000000000000
donnez le nombre de points pour la methode : 4
I =
0.999999999999999

f (x) =

x
1x2

Voici le rsultat obtenu aprs excution de la commande :

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

>> integration('gausstchebychev')
Entrez la borne inferieure : 0
Entrez la borne superieur : 1
Donnez la fonction f(x) qui est au nominateur ponderee par racine de (1x^2) 'x^2'
donnez le nombre de points pour la methode : 2
I =
0.333333333333333
donnez le nombre de points pour la methode : 3
I =
0.333333333333333

10

16
17
18
19
20
21

donnez le nombre de points pour la methode : 4

I =
0.333333333333333

Figure 3  Reprsentation de l'erreur en fonction du nombre de points

f (x) =

x
1x2

Voici le rsultat obtenu aprs excution de la commande :

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

>> integration('gausstchebychev')
Entrez la borne inferieure : 0
Entrez la borne superieur : 1
Donnez la fonction f(x) qui est au nominateur ponderee par racine de (1x^2) 'x^4'
donnez le nombre de points pour la methode : 2
I =
0.194444444444444
donnez le nombre de points pour la methode : 3
I =
0.200000000000000
donnez le nombre de points pour la methode : 4
I =
0.200000000000000

11

Figure 4  Reprsentation de l'erreur en fonction du nombre de points

Conclusion
Dans ce TP, en appliquant les mthodes numriques pour rsoudre les intgrales
dnies, on a ,encore une fois, prouv l'utilit et l'ecacit de l'utilisation de
,
de part sa rapidit excuter les instructions donnes et l'exactitude des rsultats
fournis.
On a aussi conrm qu'il tait prfrable de faire appel aux mthodes de quadratures
car elles convergent plus rapidement vers la valeur exacte de l'intgrale compare aux
autres mthodes, donc plus conomique pour les machines.

Matlab

12