|
Visites depuis le 14 juin 2002 |
Chapitre IV : Intégration
numériques des fonctions
1. Introduction
2. Méthodes d’intégrations numériques
2.1. Méthode des trapèzes
2.2. Méthode de Simpson
3. Fonctions MATLAB utilisées pour l'intégration numérique
Dans la plupart des cas, les fonctions analytiques, du fait de leurs complexités, ne sont pas intégrables analytiquement. Dans d'autres cas, on a des fonctions qui sont évaluées numériquement en différents points de l’intervalle où ces dernières sont données, et l’intégrale de ces types de fonctions ne peut être obtenue que par des approches numériques. Dans ce chapitre, on s’intéresse aux méthodes utilisées fréquemment ; à savoir la méthode des trapèzes, la méthode de Simpson, les formules de Newton-Cotes et la méthode de Gauss. Nous étudierons également les méthodes de calcul d’intégrales doubles et les intégrales des fonctions non définies sur leurs bornes.
2. Méthodes d’intégrations numériques
Soit 
la fonction à intégrer sur 
.
L’intégrale 
de 
s’écrit
en utilisant la méthode des trapèzes :
où 
; 
; 
et 
.
Le terme représentant l’erreur est :
est la
moyenne de 
sur l’intervalle 
.
L’erreur 
est inversement proportionnelle à la valeur de 
.
Si la fonction 
est donnée sur les intervalles réguliers (
),
la méthode des trapèzes peut être écrite dans
Matlab
sous la forme :
avec :
On peut également faire un programme pour calculer l'intégrale 
.
Celui-ci appelé 'trapez_v.m' par exemple, est listé
ci-dessous :
sur
un intervalle 
où le pas est "h" égal à 1. En mode
interactif dans MATLAB (mode commande), avant de lancer le programme 'trapez_v.m',
on donne la fonction
ainsi que ses bornes :
>>x=-10:1:8;
>>f=x.^2+2*x-1;
>>h=1;
On exécute ensuite le programme 'trapez_v.m', et on obtient :
>>I=trapez_v(f,h)
I =
2265
En utilisant les programmes 'trapez_n.m' et 'trapez_g.m' listés ci-après, on peut calculer aussi l'intégrale I.
>>I=trapez_n('fonction_f',a,b,n)
ou
>>I=trapez_g('fonction_f',a,b,n)
Liste du programme 'trapez_n.m' :
Liste du programme 'trapez_g.m' :
h2=(b-a)/100;
xc=a+(0:100)*h2;
fc=feval(‘fonction_f’,xc);
plot(xc,fc,'r');
hold on;
title('Méthode
des trapèzes');
xlabel('x');
ylabel('y');
grid on;
plot(x,f,'m');
plot(x,zeros(size(x)),'c')
for
i=1:n;
plot([x(i),x(i)],[0,f(i)],'g');
end
La fonction 'fonction_f' est donnée par le programme 'fonction_f.m'.
Remarque :
La fonction y=feval('sin',x) est équivalente à y=sin(x).
Exemples :
1°) On donne la fonction 'fonction_f' par le sous-programme 'fonction_f.m' listé ci-dessous :
En exécutant le programme trapez_g('fonction_f',1,50,60),
où les bornes 
et 
sont
égales respectivement à 1 et 50 et le nombre d'intervalles 
est égal 60, on obtient l'intégrale de 
sur 
pour
60 pas :
>>trapez_g('fonction_f',1,50,60)
ans =
279.3889
2°) Un véhicule de masse m=2000 kg se déplace
à la vitesse de V=30 m/s. Soudainement, le moteur est débrayé
à 
.
L'équation du mouvement après l'instant
est donnée par :
(a)
où x est la distance linéaire mesurée à
partir de .
Dans cette équation (a), le terme de gauche représente la force d'accélération, le 1er terme de droite représente la résistance aérodynamique exercée par le vent sur le véhicule et le second terme est le coefficient de frottement. Calculer la distance au delà de laquelle la vitesse de la voiture se réduit à 15 m/s.
Solution :
L'équation (a) peut être écrite sous la forme :
L'intégration de cette équation donne :
Cette intégrale peut être évaluée par la méthode des trapèzes.
Si on se donne 60 intervalles (ou 61 points), on peut écrire :
où 
et 
.
En définissant :
et en appliquant la méthode des trapèzes pour l'intégration, on obtient :
(b)
Le programme suivant (trapeze.m)permet de calculer la valeur 
donnée par la formulation précédente (b) :
>>trapeze
x =
127.5066
En comparant cette solution à la solution exacte (=127,51 m), on remarque que l'erreur relative induite par cette méthode est de l'ordre de 2,7.10-3 %.
La liste du programme trapeze.m est la suivante :
3°) Connaissant l'intégrale exacte de la fonction 
sur l'intervalle 
,
on a 
.
Étudier l'effet du nombre d'intervalles 
sur l'erreur de calcul en adoptant la méthode trapézoïdale
:
Solution :
Le programme suivant 'trz.m' permet de calculer l'intégrale 
:
Les résultats obtenus sont les suivants :
Méthode Trapézoïdale étendue
n I Pourc. erreur relat.
---------------------------------------------
2 4.08358
0.01911429
4 4.02619
0.00478998
8 4.01180
0.00119825
16 4.00819
0.00029965
32 4.00729
0.00007496
64 4.00707
0.00001878
128 4.00701 0.00000474
256 4.00700 0.00000123
512 4.00700 0.00000035
1024 4.00699 0.00000013
D'après ces résultats, on montre que quand
est doublé, l'erreur relative décroît par un facteur
4.
Soit 
l'intégrale de 
sur l'intervalle 
.
Par la méthode de Simpson, 
s'écrit :
où 
; 
et .
Le terme 
est donné par :
et 
est la moyenne de 
sur l'intervalle ouvert 
.
Exemple :
Évaluer l'intégrale de 
sur l'intervalle 
avec la méthode de Simpson pour 
et 
.
Solution :
La liste du programme 'msimp.m' est la suivante :
Les résultats de cet algorithme sont donnés ci-dessous :
Méthode de Simpson
n I Pourc. erreur relat.
---------------------------------------------
2 4.00791
0.00022935
4 4.00705
0.00001521
8 4.00700
0.00000101
16 4.00699
0.00000012
On montre qu'à partir de 
,
l'erreur relative devient presque nulle (de l'ordre de 
).
Les programmes 'simp_v' et 'simp_n' peuvent être
utilisés pour intégrer la fonction 
.
Dans l'algorithme 'simp_v(f,h)',
est un vecteur contenant les ordonnées de l'intégrante et 
représente le pas de l'intervalle d'intégration.
D'où :
>>I=simp_n('fonction_f',a,b,n)
permet d'évaluer l'intégrale
dans laquelle 'fonction_f' est le nom du sous-programme contenant
la fonction 
, 
et 
sont
les limites de l'intervalle d'intégration et 
.
Exemple :
Soit 
où 
et 
.
è
Ecrire un algorithme utilisant le sous-programme simp_v(f,h)
(à écrire) et qui permet d'évaluer
sur 100 intervalles.
Solution :
Appelons ce programme simp.m par exemple, et donnons sa liste
:
Le sous-programme simp_v(f,h) est :
function
I=simp_v(f,h);
n=length(f)-1; % Nombre d’intervalles
if n==1
fprintf('Données
à un seul intervalle'\n');
return;
end
if n==2
I=h/3*(f(1)+4*f(2)+f(3));
return; % Revenir et continuer le calcul
end
if n==3
I=(3/8)*h*(f(1)+3*f(2)+3*f(3)+f(4));
return;
end
I=0;
if
2*floor(n/2)~=n
; %Le nombre d’intervalle doit être pair
I=3*h/8*(f(n-2)+3*f(n-1)+3*f(n)+f(n+1));
m=n-3;
else
m=n;
end
I=I+h/3*(f(1)+4*sum(f(2:2:m))+f(m+1));
if m>2
I=I+(h/3)*2*sum(f(3:2:m))
end
En exécutant le programme simp.m, on trouve l'intégrale
par la méthode de Simpson.:
>>simp
I =
1.8522e+003
3. Fonctions MATLAB utilisées pour l'intégration numérique
Dans Matlab (Toolbox), il existe 2 fonctions appelées 'quad'
et
'quad8'
pour l'intégration numérique. La fonction
'quad' utilise
la méthode de Simpson et la fonction 'quad8' utilise les
formules de Newton-Cotes à l'ordre 8. Ces deux fonctions
sont utilisées de la façon suivante :
Dans la première forme, la tolérance 'tol' qui
correspond à l'erreur relative 
,
est considérée égale à 0,001 par défaut.
Ainsi, le calcul quadratique de l'intégrale est réitéré
jusqu'à ce que la tolérance soit satisfaite. Si la 3ème
forme est utilisée avec une valeur non nulle de 'trace',
un graphique d'évolution des itérations sera affiché
sur l'écran. Quand à la fonction 'quad8', elle est
utilisée de la même manière que la fonction 'quad'
dans les 3 formes d'expressions précédentes.
Exemple :
Soit à calculer l'intégrale suivante :
avec la transformation suivante :
Dans ce cas, on a :
On obtient donc :
Or, on connaît la fonction 
(fonction Gamma) eulérienne du deuxième espèce qui
est définie par :
où 
est un réel.
Pour un entier 
strictement positif, 
a la propriété suivante :
On pourra donc s'en servir pour calculer la factorielle d'un entier
naturel 
:
è 
Dans Matlab, cette fonction est notée 'gamma'.
Exemple :
La factorielle de 10 est donnée par gamma(10+1)=gamma(11) :
>> fact=gamma(11)
fact =
3628800
Cette propriété est à l'origine du nom de la 'fonction factorielle', attribuée souvent à la fonction 'Gamma'. La fonction 'gamma' est donc prédéfinie dans Matlab. Cette fonction permet, de part sa définition pour tous les arguments réels positifs, d'introduire 'la factorielle' des arguments non entiers.
Exemple :
>> format long
>>gamma(0.5)
ans =
1.77245385090552
>>racine_de_pi=sqrt(pi)
racine_de_pi =
1.77245385090552
>>gamma(0)
Warning: Divide by zero.
ans =
Inf
>>gamma(-2)
Warning: Divide by zero.
ans =
-Inf
Tracé de la fonction Gamma
Tracer la fonction Gamma pour 
.
>>x=-3:0.1:3;
>>gx=gamma(x);
>>plot(x,gx);grid;title('Allure de la fonction Gamma');
>>xlabel('x');ylabel('Gamma(x)');
>>hold on;plot(x,zeros(size(x)))
Fig. 2 : Fonction 
Revenons maintenant à l'intégration de la fonction
:
D'après le changement de variable effectué, cette intégrale s'écrit :
Calculons cette intégrale par les fonctions Matlab 'quad8' et 'gamma'. La fonction 'quad8' nécessite l'écriture de la fonction à intégrer dans un fichier contenant l'extension '.m' (programme Matlab).
Liste du programme fct.m :
Tracé de la fonction à intégrer :
>>x=0:0.01:5;
>>ft=fct(x);
>>plot(x,ft);
>>grid; title('fonction à intégrer');
>>xlabel('x');ylabel('f(x)');
Dans le programme donné ci-dessous (quadgamm.m), on calcule l'intégrale précédente par les fonctions 'quad8' et 'gamma', pour lesquelles on compare les performances en terme de temps de calcul et du nombre d'opérations en virgule flottante.
%***************************************
% Calcul
de l'intégrale de f(x) par la *
% la
fonction 'gamma' *
%***************************************
flops(0);
tic;
Igamma=gamma(1)/2
Nbre_op2=flops
tps_g=toc
En exécutant le programme quadgamm.m, on obtient les résultats suivants :
>> format long
>> quadgamm
Iquad =
0.49999987773178
Nbre_op1 =
686
tps_q =
0.06000000000000
Igamma =
0.50000000000000
Nbre_op2 =
40
tps_g =
0.05000000000000
D'après ces résultats, on remarque bien que le nombre d'opérations en virgule flottante et le temps d'exécution dans le cas de l'utilisation de la fonction 'gamma' sont nettement inférieurs à ceux que l'on obtient en utilisant la fonction 'quad8'.
La fonction 'gamma' prend des valeurs très grandes au voisinage de zéro et des entiers négatifs. Afin d'éviter un dépassement de capacité, MATLAB dispose de la fonction 'gammaln' qui retourne le logarithme népérien de ces valeurs.
Exemple :
>>format
>>x=5e-200;
>>gamma(x)
ans =
2.0000e+199
>>gammaln(x)
ans =
458.9076