Integració Numèrica de funcions

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 precissió
N=10; %número de divisions de l'interval
a=0;
b=pi/2;
h=(b-a)/N; %pas per cada punt
x=a:h:b; %tenim N+1 punts
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 f(x)
per punts equidistants on per
Rectangles superiors i inferiors:
Obs.- Com que la funcio f és decreixent el màxim en l'interval s'assoleix en l'extrem inferior de l'interval, mentre que el mínim es fa en l'extrem superior.
AreaSup=0;
AreaInf=0;
for i=1:N
AreaSup=AreaSup+(x(i+1)-x(i))*y(i); %base per altura
AreaInf=AreaInf+(x(i+1)-x(i))*y(i+1); %base per altura
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)); %area 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 precissió 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(2019)