Integració Numèrica de funcions

Table of Contents
Exemple del càlcul de la integral d'una funció Regla dels rectangles (suma de Riemann) Regla dels trapezis Mètode de Simpson Funcions del Matlab Exercici 1 Exercici 2 (c) NumFactory 2022
En aquesta pràctica veurem com es pot aproximar la integral definida d'una funció a partir del càlcul de l'area entre la funció i l'eix x.

Exemple del càlcul de la integral d'una funció

format long e; %per veure els valors numèrics amb el màxim de precisió
N=10; %número de divisions de l'interval
a=0;
b=pi/2;
h=(b-a)/N; % pas entre cada punt
x=a:h:b; % tenim N+1 punts (inclou els extrems)
y=cos(x);
%-------- fem el dibuix -------
xplot=a:0.1*h:b;
yplot=cos(xplot);
figure()
plot(xplot,yplot);
hold on;
plot(x,y,'ro')
hold off;

Regla dels rectangles (suma de Riemann)

Consisteix en aproximar l'area de la funció pel rectangle definit per les imatges de
on . La subdivisió de l'interval es fa amb pas i els punts es defineixen com , .
Així també satisfà que .
Rectangles superiors i inferiors:
Si prenem com un dels extrems de l'interval, o respectivament, com que la funció de l'exemple és decreixent el màxim en l'interval s'assoleix en l'extrem inferior de l'interval, mentre que el mínim s'assoleix en l'extrem superior.
AreaSup=0;
AreaInf=0;
for i=1:N
AreaSup=AreaSup+(x(i+1)-x(i))*y(i); %base per altura: h*f(x(i))
AreaInf=AreaInf+(x(i+1)-x(i))*y(i+1); %base per altura: h*f(x(i+1))
end
fprintf('Arees calculades: Asup= %24.16e, Ainf= %24.16e \n',AreaSup,AreaInf)
Arees calculades: Asup= 1.0764828026941020e+00, Ainf= 9.1940317001461258e-01
 
%-------- fem el dibuix dels rectangles ------
figure()
set(gcf, 'Position', [200,200, 900, 400])
subplot(1,2,1) %superior en vermell
plot(xplot,yplot);
hold on;
for i=1:N
plot([x(i),x(i)],[0,y(i)],'r')
plot([x(i),x(i+1)],[y(i),y(i)],'r')
plot([x(i+1),x(i+1)],[0,y(i)],'r')
end
hold off;
subplot(1,2,2) %inferior en verd
plot(xplot,yplot);
hold on;
for i=1:N
plot([x(i),x(i)],[0,y(i+1)],'g')
plot([x(i),x(i+1)],[y(i+1),y(i+1)],'g')
plot([x(i+1),x(i+1)],[0,y(i+1)],'g')
end
hold off

Regla dels trapezis

Consisteix en aproximar l'area de la funció pel trapezi definit per les imatges de f(x)
La fórmula general (composta) és:
on per i=0...N, on N és el número de divisions (npunts-1) i
AreaTrap=0; %calculem fent un bucle sobre cada trapezi
for i=1:N
AreaTrap=AreaTrap+0.5*(x(i+1)-x(i))*(y(i)+y(i+1)); %àrea del trapezi
end
fprintf('Arees calculades: Asup= %24.16e, Ainf= %24.16e, ATrap= %24.16e \n',AreaSup,AreaInf,AreaTrap)
Arees calculades: Asup= 1.0764828026941020e+00, Ainf= 9.1940317001461258e-01, ATrap= 9.9794298635435719e-01
% Formula composta: (equivalent a l'anterior, però més eficient)
AreaTrapComp= ((b-a)/(2*N))*(y(1)+2*sum(y(2:end-1))+y(end));
%veiem que dona el mateix:
[AreaTrap, AreaTrapComp, AreaTrap-AreaTrapComp]
ans = 1×3
9.979429863543572e-01 9.979429863543573e-01 -1.110223024625157e-16
% error comès sobre el valor de la integral (recordem que dona 1)
er=abs(AreaTrapComp-1);
fprintf('Trapezis: Error comès amb N= %d divisions er=%24.16e \n',N,er);
Trapezis: Error comès amb N= 10 divisions er= 2.0570136456427024e-03
 
%-------- fem el dibuix dels trapezis ------
figure()
plot(xplot,yplot);
hold on;
for i=1:N
plot([x(i),x(i)],[0,y(i)],'r')
plot([x(i),x(i+1)],[y(i),y(i+1)],'r')
plot([x(i+1),x(i+1)],[0,y(i+1)],'r')
end
hold off;

Mètode de Simpson

Una millora notable en la precisió consisteix en aproximar l'area de la funció, per l'area calculada en aproximar cada 3 punts per un segment de paràbola (noteu que el cas dels trapezis és el mateix però aproximant per una recta cada 2 punts)
La fórmula composta és:
Obs.- Pel mètode de Simpson el múmero de divisions, N, ha de ser parell (npunts imparell !!)
AreaSimp=((b-a)/(3*N))*(y(1)+4*sum(y(2:2:(end-1)))+2*sum(y(3:2:(end-1)))+y(end));
er=abs(AreaSimp-1);
fprintf('Arees calculades: Asup= %16.10e, Ainf= %16.10e, ATrap= %16.10e, ASimp= %16.10e \n',AreaSup,AreaInf,AreaTrap,AreaSimp)
Arees calculades: Asup= 1.0764828027e+00, Ainf= 9.1940317001e-01, ATrap= 9.9794298635e-01, ASimp= 1.0000033922e+00
fprintf('Simpson: Error comès amb N= %d divisions er=%24.16e \n',N,er);
Simpson: Error comès amb N= 10 divisions er= 3.3922209004000337e-06

Funcions del Matlab

El Matlab ja te implementada una funció trapezis: trapz
AreaTrapzMatlab = trapz(x,y); %noteu que el valor coincideix amb el càlcul anterior
er=abs(AreaTrapzMatlab-1);
fprintf('Trapz Matlab: Error comès amb N= %d, valor=%24.16e divisions er=%24.16e \n',length(x)-1,AreaTrapzMatlab,er);
Trapz Matlab: Error comès amb N= 10, valor= 9.9794298635435719e-01 divisions er= 2.0570136456428134e-03
 
% Obs: Si es vol millorar la precisió s'han de tenir més punts per la funció (només en
% tenim 11), per exemple si prenem els punts que fem servir per dibuixar (en tindrem 101)
 
AreaTrapzMatlab2 = trapz(xplot,yplot);
er=abs(AreaTrapzMatlab2-1);
fprintf('Trapz Matlab: Error comès amb N= %d, valor=%24.16e divisions er=%24.16e \n',length(x)-1,AreaTrapzMatlab2,er);
Trapz Matlab: Error comès amb N= 10, valor= 9.9997943823960755e-01 divisions er= 2.0561760392445727e-05
 
% També el Matlab té la funció *integral* que és un mètode de Simpson adaptatiu.
% Requereix definir la funció de forma inline o en un fitxer de funció i dona el resultat amb un error relatiu per defecte de 10^-6
 
f=@(x) cos(x);
AreaInt = integral(f,a,b);
er=abs(AreaInt-1);
fprintf('Trapz Matlab: valor= %24.16e, error comès error=%24.16e \n',AreaInt,er);
Trapz Matlab: valor= 9.9999999999999989e-01, error comès error= 1.1102230246251565e-16

Exercici 1

(by P. Gutiérrez)
Ens proposem calcular numèricament la integral
el valor exacte de la qual és , aplicant els mètodes de trapezis i Simpson amb diferents passos d'integració i veure l'error comès en cada cas.
(1.1) Calculeu els valors T1, T2, T3 obtinguts aplicant el mètode dels trapezis amb N=50,500,5000 subintervals respectivament (és a dir, 51,501,5001 punts), usant la funció trapz.
f=@(x) (5-x.^2).*sin(x);
a=0;
b=7;
 
N=50;
h1=(b-a)/N;
x1=a:h1:b;
y1=f(x1);
T1=trapz(x1,y1);
fprintf('T1=%.12f\n',T1) %amb 12 decimals i saltem linia
T1=29.388676198730
 
N=500;
h2=(b-a)/N;
x2=a:h2:b;
y2=f(x2);
T2=trapz(x2,y2);
fprintf('T2=%.12f\n',T2)
T2=29.465308594746
 
N=5000;
h3=(b-a)/N;
x3=a:h3:b;
y3=f(x3);
T3=trapz(x3,y3);
fprintf('T3=%.12f\n',T3)
T3=29.466074563335
(1.2) L'error en el mètode dels trapezis és . Això vol dir que si dividim el pas d'integració per 10, l'error es dividirà aproximadament per 100 (és a dir, guanyem 2 xifres decimals correctes).
Per verificar aquest fet, calculeu el valor exacte I i calculeu els errors err1=|I-T1|, err2=|I-T2|, err3=|I-T3| obtinguts amb cada aproximació. A continuació, calculeu els quocients err1/err2, err2/err3 (han de ser valors propers a 100).
I=42*cos(7)-14*sin(7)+7;
fprintf('I=%.12f\n\n',I)
I=29.466082300356
 
err1=abs(I-T1);
fprintf('err1=%.6e\n',err1) %format exponencial amb 6 decimals
err1=7.740610e-02
 
err2=abs(I-T2);
fprintf('err2=%.6e\n',err2)
err2=7.737056e-04
 
err3=abs(I-T3);
fprintf('err3=%.6e\n\n',err3)
err3=7.737021e-06
 
q1=err1/err2;
fprintf('q1=%.6f\n', q1)
q1=100.045936
 
q2=err2/err3;
fprintf('q2=%.6f\n', q2)
q2=100.000459
(1.3) Calculeu els valors S1, S2, S3 obtinguts aplicant el mètode de Simpson amb N=50,500,5000 subintervals respectivament, aplicant el sumatori de coeficients 1,4,2,4,...,2,4,1.
N=50;
S1=h1/3*(y1(1)+4*sum(y1(2:2:N))+2*sum(y1(3:2:N-1))+y1(N+1));
%no cal repetir el calcul del pas h1 i dels vectors x1,y1
fprintf('S1=%.12f\n',S1)
S1=29.466226253913
 
N=500;
S2=h2/3*(y2(1)+4*sum(y2(2:2:N))+2*sum(y2(3:2:N-1))+y2(N+1));
fprintf('S2=%.12f\n',S2)
S2=29.466082314707
 
N=5000;
S3=h3/3*(y3(1)+4*sum(y3(2:2:N))+2*sum(y3(3:2:N-1))+y3(N+1));
fprintf('S3=%.12f\n',S3)
S3=29.466082300357
(1.4) L'error en el mètode de Simpson és . Això vol dir que si dividim el pas d'integració per 10, l'error es dividirà aproximadament per 10000 (és a dir, guanyem 4 xifres decimals correctes).
Per verificar aquest fet, calculeu el valor exacte I i calculeu els errors err1=|I-S1|, err2=|I-S2|, err3=|I-S3| obtinguts amb cada aproximació. A continuació, calculeu els quocients err1/err2, err2/err3 (han de ser valors propers a 10000).
err1=abs(I-S1);
fprintf('err1=%.6e\n',err1)
err1=1.439536e-04
 
err2=abs(I-S2);
fprintf('err2=%.6e\n',err2)
err2=1.435131e-08
 
err3=abs(I-S3);
fprintf('err3=%.6e\n\n',err3)
err3=1.421085e-12
 
q1=err1/err2;
fprintf('q1=%.6f\n', q1)
q1=10030.687754
 
q2=err2/err3;
fprintf('q2=%.6f\n', q2)
q2=10098.840000
(1.5) Per calcular l'aproximació amb el mètode de Simpson, una alternativa al sumatori és aplicar la fórmula d'extrapolació que la relaciona amb el mètode dels trapezis aplicat amb dos passos d'integració diferents: Per calcular el valor de Simpson amb pas d'integració h, primer calculem els valors i de trapezis amb passos h i respectivament, i apliquem la fórmula
Useu aquesta fórmula d'extrapolació per calcular Simpson amb N=50 subintervals (ha de coincidir amb el valor S1 calculat anteriorment).
N=50;
h=(b-a)/N;
 
x50=a:h:b;
y50=f(x50);
T50=trapz(x50,y50); %trapezis amb pas h
 
x25=a:2*h:b;
y25=f(x25);
T25=trapz(x25,y25); %trapezis amb pas 2*h
 
S50=(4*T50-T25)/3; %formula d'extrapolacio
fprintf('S50=%.12f\n',S50)
S50=29.466226253913
 
dif=S50-S1 %comparem amb el valor S1 ja calculat amb el sumatori
dif =
3.552713678800501e-15

Exercici 2

Calculeu la integral definida
el valor exacte de la qual és . Feu el càlcul amb:
Sol: T(32)= 2.500970952556458e-01, error= 9.709525564577381e-05
Sol: T(64)= 2.500242744345849e-01, error= 2.427443458485889e-05
Sol: S(64)= 2.500000008275646e-01, error= 8.275645724253877e-10

(c) NumFactory 2022