126
4 Intégration et différentiation numérique
Faisons maintenant l’hypothèse que f
(4) (x) est approximativement
constante sur l’intervalle [α, β]. Dans ce cas, f
(4) (ξ) f
(4) (η). On peut
calculer f
(4) (η) à partir de (4.31) puis, injectant cette valeur dans l’équation (4.30), on obtient cette estimation de l’erreur
β
α
f(x) dx − I
c
s (f)
1
15
ΔI.
Le pas (β−α)/2 (qui est le pas utilisé pour calculer I
c
s (f)) sera accepté
si |ΔI|/15 < ε(β − α)/[2(b − a)]. La formule de quadrature qui utilise ce
critère dans le procédé d’adaptation décrit ci-dessus est appelée formule
de Simpson adaptative. Elle est implémentée dans le Programme 4.3.
Parmi les paramètres d’entrée, f est la chaîne de caractère qui définit la
fonction f, a et b sont les extrémités de l’intervalle d’intégration, tol est
la tolérance fixée sur l’erreur et hmin est la longueur minimale admise
pour le pas d’intégration (afin d’assurer que le procédé d’adaptation ne
boucle pas indéfiniment).
Programme 4.3. simpadpt : formule de Simpson adaptative
function [ JSf , nodes ]= simpadpt ( fun ,a ,b , tol , hmin , varargin )
% SIMPADPT calcul n u m ériq ue de l ’ i n t egra le avec la
% méthode de Simpson a d a pt ative.
% JSF = SIMPADPT (FUN ,A ,B , TOL , HMIN ) tente d ’ a p p r oche r
% l ’ i n t égr ale de la fonction FUN de A à B avec une
% erreur i n f é rie ure à TOL en u t i l isant par r é c u rr ence
% la méthode a d a pta tive de Simpson avec H >= HMIN .
% La fonction Y = FUN( X) doit accepter en
% entrée un vecteur X et r e t ou rner dans un vecteur Y ,
% les valeurs de l ’ i n t ég rande en chaque c o m p osante de V.
% FUN peut être une fonction inline , une fonction
% anonyme ou définie par un m - file .
% JSF = SIMPADPT (FUN ,A ,B , TOL , HMIN ,P1 , P2 ,...) appelle la
% fonction FUN en passant les p a r amè tres o p t i onn els
% P1 , P2 ,... de la maniere suivante : FUN (X ,P1 , P2 ,...).
% [ JSF , NODES ] = SIMPADPT (...) renvoie la
% d i s t ribu tion des noeuds .
A =[a , b ]; N =[]; S =[]; JSf = 0; ba = 2*(b - a ); nodes =[];
while ~ isempty ( A ) ,
[ deltaI , ISc ]= c a l delt ai(A , fun , varargin {:});
if abs( deltaI ) < 15* tol *(A (2) -A (1))/ ba;
JSf = JSf + ISc ;
S = union (S , A );
nodes = [ nodes , A (1) (A (1)+ A (2))*0.5 A (2)];
S = [S(1), S(end)]; A = N; N = [];
elseif A (2) -A (1) < hmin
JSf= JSf+ ISc ;
S = union (S , A );
S = [S (1) , S( end )]; A =N ; N =[];
warning ( ’Pas d ’’ i n t e grati on trop petit ’);
else
Am = ( A (1)+ A ( 2 ) )*0. 5;
A = [A (1) Am ];
N = [Am , b ];
end
end
4 Intégration et différentiation numérique
Faisons maintenant l’hypothèse que f
(4) (x) est approximativement
constante sur l’intervalle [α, β]. Dans ce cas, f
(4) (ξ) f
(4) (η). On peut
calculer f
(4) (η) à partir de (4.31) puis, injectant cette valeur dans l’équation (4.30), on obtient cette estimation de l’erreur
β
α
f(x) dx − I
c
s (f)
1
15
ΔI.
Le pas (β−α)/2 (qui est le pas utilisé pour calculer I
c
s (f)) sera accepté
si |ΔI|/15 < ε(β − α)/[2(b − a)]. La formule de quadrature qui utilise ce
critère dans le procédé d’adaptation décrit ci-dessus est appelée formule
de Simpson adaptative. Elle est implémentée dans le Programme 4.3.
Parmi les paramètres d’entrée, f est la chaîne de caractère qui définit la
fonction f, a et b sont les extrémités de l’intervalle d’intégration, tol est
la tolérance fixée sur l’erreur et hmin est la longueur minimale admise
pour le pas d’intégration (afin d’assurer que le procédé d’adaptation ne
boucle pas indéfiniment).
Programme 4.3. simpadpt : formule de Simpson adaptative
function [ JSf , nodes ]= simpadpt ( fun ,a ,b , tol , hmin , varargin )
% SIMPADPT calcul n u m ériq ue de l ’ i n t egra le avec la
% méthode de Simpson a d a pt ative.
% JSF = SIMPADPT (FUN ,A ,B , TOL , HMIN ) tente d ’ a p p r oche r
% l ’ i n t égr ale de la fonction FUN de A à B avec une
% erreur i n f é rie ure à TOL en u t i l isant par r é c u rr ence
% la méthode a d a pta tive de Simpson avec H >= HMIN .
% La fonction Y = FUN( X) doit accepter en
% entrée un vecteur X et r e t ou rner dans un vecteur Y ,
% les valeurs de l ’ i n t ég rande en chaque c o m p osante de V.
% FUN peut être une fonction inline , une fonction
% anonyme ou définie par un m - file .
% JSF = SIMPADPT (FUN ,A ,B , TOL , HMIN ,P1 , P2 ,...) appelle la
% fonction FUN en passant les p a r amè tres o p t i onn els
% P1 , P2 ,... de la maniere suivante : FUN (X ,P1 , P2 ,...).
% [ JSF , NODES ] = SIMPADPT (...) renvoie la
% d i s t ribu tion des noeuds .
A =[a , b ]; N =[]; S =[]; JSf = 0; ba = 2*(b - a ); nodes =[];
while ~ isempty ( A ) ,
[ deltaI , ISc ]= c a l delt ai(A , fun , varargin {:});
if abs( deltaI ) < 15* tol *(A (2) -A (1))/ ba;
JSf = JSf + ISc ;
S = union (S , A );
nodes = [ nodes , A (1) (A (1)+ A (2))*0.5 A (2)];
S = [S(1), S(end)]; A = N; N = [];
elseif A (2) -A (1) < hmin
JSf= JSf+ ISc ;
S = union (S , A );
S = [S (1) , S( end )]; A =N ; N =[];
warning ( ’Pas d ’’ i n t e grati on trop petit ’);
else
Am = ( A (1)+ A ( 2 ) )*0. 5;
A = [A (1) Am ];
N = [Am , b ];
end
end