Pràctica: INTEGRACIÓ SOBRE CORBES

Aquesta pràctica ens permetrà calcular numèricament dos tipus d'integrals sobre corbes (també anomenades integrals de línia):
i anàlogament per a corbes i camps vectorials a
En qualsevol d'aquests casos, hem d'acabar calculant una integral simple amb una funció que construïm a partir de les dades del problema. Per calcular aproximadament aquesta integral, veurem en primer lloc la regla dels trapezis com a mètode numèric molt simple, i veurem també com usar algun dels mètodes de quadratura adaptativa incorporats al Matlab.

Regla dels trapezis

Escollim un nombre i dividim l'interval en N subintervals de longitud que anomenem el pas d'integració,
i calculem els valors Aproximem la integral de g sobre cadascun dels subintervals per l'àrea del trapezi limitat pel segment que uneix els punts
Sumant aquests resultats per a obtenim una aproximació per a la integral que és la regla dels trapezis (o mètode dels trapezis):
Aquest mètode ens ve donat per la funció de Matlab trapz, que té com a arguments el vector o rang de valors format per les abscisses i un altre vector format per les ordenades
Com a exemple, calculem usant la regla dels trapezis amb N=50 subintervals, i comparem amb el valor exacte
g=@(t) t.^3.*exp(-t);
% Recordeu que es necessari utilitzar els operadors ".*", "./", ".^"
% per avaluar sobre rangs de valors les funcions que definim.
N=50;
h=(4-0)/N; %pas d'integracio
tv=0:h:4;
%vector d'abscisses, que alternativament podiem haver generat
%amb l'ordre "tv=linspace(0,4,N+1);"
gv=g(tv); %vector d'ordenades
t1=trapz(tv,gv) %aproximacio de la integral
t1 =
3.3990e+00
I=6-142*exp(-4); %valor exacte (conegut)
error1=abs(I-t1)
error1 =
1.5597e-04
L'error en la regla 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). Podem verificar aquest fet a l'exemple anterior, repetint els càlculs amb N=500 subintervals.
N=500;
h=(4-0)/N;
tv=0:h:4;
gv=g(tv);
t2=trapz(tv,gv) %nova aproximacio de la integral
t2 =
3.3992e+00
error2=abs(I-t2)
error2 =
1.5629e-06
quocient=error1/error2
quocient =
9.9793e+01
La regla dels trapezis és útil sobretot quan no tenim una fórmula explícita per a la funció, sinó que ens ve donada com una taula de valors.

Mètodes de quadratura adaptativa

Diverses funcions de Matlab incorporen mètodes d'integració més acurats que la regla dels trapezis, que podem aplicar si coneixem una fórmula explícita de la funció a integrar. En particular, podem usar quadl (quadratura adaptativa de Lobatto), en què no s'ha d'especificar un pas d'integració (que serà variable) sinó una tolerància per a l'error absolut (si no s'especifica, per defecte és 1e-6). Així, escriurem quadl(g,a,b,tol) per a calcular una integral amb una tolerància d'error tol.
Il·lustrem-ho amb la integral escollint tolerància 1e-8, i comparem amb el valor exacte
g=@(t) t.*sqrt(abs(1-t));
ql=quadl(g,0,2,1e-8)
ql =
1.3333e+00
I=4/3;
error3=abs(I-ql)
error3 =
6.2034e-10
Nota: Una opció alternativa per definir la funció a integrar és, en comptes d'usar una ordre del tipus g=@(t)..., fer-ho a través d'un fitxer, anomenat per exemple gfun.m, que contingui les línies següents,
function y=gfun(t)
y=t.*sqrt(abs(1-t));
En aquest cas, per calcular la integral escriuríem quadl(@gfun,0,2,1e-8) (posant "@" abans del nom del fitxer).
Altres funcions de Matlab que permeten calcular integrals són integral (quadratura adaptativa global), o quadgk (quadratura adaptativa de Gauss-Kronrod). Podeu consultar el "help" per a cadascuna.

Camp vectorial i parametrització de la corba

Per il·lustrar el càlcul d'integrals sobre corbes, considerem el càlcul de l'àrea d'un domini a partir del teorema de Green, integrant el camp al llarg de la frontera del domini:
Com a exemple, calcularem l'àrea del domini definit en coordenades polars per és a dir, limitat per l' espiral logarítmica i el segment entre els seus punts inicial i final (el valor exacte és ). Només caldrà integrar al llarg de l'espiral, ja que el segment és ortogonal al camp
Comencem definint el camp vectorial i la parametrització de la corba que recorre la frontera del domini en sentit antihorari.
P=@(x,y) -y;
Q=@(x,y) x;
% Nota: Si volem definir un camp que tingui alguna component constant, per
% exemple Q=1, podem escriure "Q=@(x,y) 1+0*x;" o alguna cosa semblant, per
% tal d'obtenir vectors de la mateixa dimensio quan ho apliquem sobre rangs
% de valors.
a=0.2; R=1; %per exemple
r=@(t) R*exp(-a*t);
x=@(t) r(t).*cos(t); %el parametre t coincideix amb l'angle de les polars
y=@(t) r(t).*sin(t);
Dr=@(t) -R*a*exp(-a*t); %tambe necessitarem les derivades
Dx=@(t) Dr(t).*cos(t)-r(t).*sin(t);
Dy=@(t) Dr(t).*sin(t)+r(t).*cos(t);

Dibuix del camp vectorial i la corba

Construïm una xarxa de punts a i dibuixem el camp vectorial a tots aquests punts, amb la funció de Matlab quiver, i tot seguit dibuixem la corba usant plot (a usaríem quiver3 i plot3).
close all %tanquem figures anteriors (optatiu)
figure %nova figura
hold on
[xg,yg] = meshgrid( -1.1*R:R/5:1.1*R, -1.1*R:R/5:1.1*R );
%obtenim dues matrius xg i yg, que contenen la coordenada x
%i la coordenada y de tots els punts, respectivament
Pg=P(xg,yg);
Qg=Q(xg,yg);
%hem avaluat les dues components del camp vectorial a tots els punts
quiver(xg,yg,Pg,Qg)
%opcionalment podem afegir-hi un factor per dibuixar els vectors mes
%llargs o mes curts, per exemple "quiver(xg,yg,Pg,Qg,1.25)"
axis equal %igualem les proporcions de les coordenades
N=100;
h=2*pi/N; %increment del parametre de la corba
tv=0:h:2*pi; %rang de N+1 valors del parametre t
plot(x(tv),y(tv),'r') %dibuix de la corba
plot([x(0),x(2*pi)],[y(0),y(2*pi)],'r') %tambe dibuixem el segment

Càlcul de la integral usant la regla dels trapezis

Avaluem la funció a integrar,
sobre el rang de valors del paràmetre, obtenint un altre vector amb els valors Tot seguit apliquem la regla dels trapezis per tal d'obtenir una aproximació de la integral.
N=100;
h=2*pi/N; %pas d'integracio
tv=0:h:2*pi;
xv=x(tv); %punts de la corba
yv=y(tv);
Dxv=Dx(tv); %derivades de x i y als punts de la corba
Dyv=Dy(tv);
Pv=P(xv,yv); %avaluem el camp als punts de la corba
Qv=Q(xv,yv);
gv=Pv.*Dxv+Qv.*Dyv; %valors de la funcio a integrar
intaprox1=trapz(tv,gv);
areaaprox1=intaprox1/2
areaaprox1 =
1.1488e+00
area=R^2*(1-exp(-4*pi*a))/(4*a);
errorarea1=abs(area-areaaprox1) %comparem amb el valor exacte
errorarea1 =
6.0467e-05

Càlcul de la integral usant quadratura adaptativa

En aquest cas hem de definir com una funció, a partir de les funcions P, Q, x, y, Dx, Dy ja definides.
g=@(t) P(x(t),y(t)).*Dx(t)+Q(x(t),y(t)).*Dy(t);
intaprox2=quadl(g,0,2*pi,1e-12);
areaaprox2=intaprox2/2
areaaprox2 =
1.1487e+00
errorarea2=abs(area-areaaprox2)
errorarea2 =
2.2204e-16
________
(c) Numerical Factory (by Pere Gutiérrez)