Pràctica: INTEGRACIÓ SOBRE SUPERFÍCIES

Aquesta pràctica ens permetrà calcular numèricament dos tipus d'integrals sobre superfícies de :
En qualsevol d'aquests casos, hem d'acabar calculant una integral doble amb una funció de 2 variables que construïm a partir de les dades del problema. Per calcular aproximadament aquesta integral, usarem la funció integral2 incorporada al Matlab (vegeu-ne el funcionament a l'apèndix).

Representació gràfica d'una superfície i càlcul d'una integral

Veurem amb dos exemples com graficar una superfície amb Matlab, i calcular una integral sobre la superfície. En el primer exemple la superfície ve parametritzada amb un domini rectangular, i en el segon exemple amb un altre tipus de domini.
Exemple 1. Donada la porció de tor S parametritzada per
amb domini (recordem que el tor complet ve donat per ), representem S gràficament i, a més, calculem el flux del camp vectorial a través de que suposem orientada per la normal exterior.
Descrivim la resolució en diferents etapes. Si només ens interessa calcular la integral i no la representació gràfica, després de fer (a) podem passar directament a (f).
(a) Definicions de la parametrització i del camp vectorial. Comencem definint les components de la parametrització Φ com a funcions dels paràmetres
x=@(th,t) (1+cos(t)/2).*cos(th);
y=@(th,t) (1+cos(t)/2).*sin(th);
z=@(th,t) sin(t)/2;
Ara definim les components del vector com a derivades parcials de les funcions anteriors respecte
Dxth=@(th,t) -(1+(cos(t)/2)).*sin(th);
Dyth=@(th,t) (1+cos(t)/2).*cos(th);
Dzth=@(th,t) 0;
i anàlogament les components de com a derivades parcials respecte
Dxt=@(th,t) -sin(t).*cos(th)/2;
Dyt=@(th,t) -sin(t).*sin(th)/2;
Dzt=@(th,t) cos(t)/2;
i també, fent el producte vectorial, les components del vector normal (en general, no unitari)
nx=@(th,t) Dyth(th,t).*Dzt(th,t)-Dzth(th,t).*Dyt(th,t);
ny=@(th,t) Dzth(th,t).*Dxt(th,t)-Dxth(th,t).*Dzt(th,t);
nz=@(th,t) Dxth(th,t).*Dyt(th,t)-Dyth(th,t).*Dxt(th,t);
Finalment, definim les components del camp vectorial com a funcions de
P=@(x,y,z) x;
Q=@(x,y,z) y;
R=@(x,y,z) z;
(b) Rangs de valors dels paràmetres. Per tal de recórrer tot el domini, escollim i subintervals en els paràmetres θ i t respectivament, i definim rangs de valors thv i tv.
N1=40; thv=linspace(0,pi,N1+1);
N2=20; tv=linspace(0,2*pi/3,N2+1);
% Nota: tambe ho podem definir a partir dels increments en cada parametre,
% delta1=(pi-0)/N1; thv=0:delta1:pi;
% delta2=(2*pi/3-0)/N2; tv=0:delta2:2*pi/3;
Amb la funció meshgrid següent, obtenim dues matrius thm i tm amb la primera i segona coordenada dels punts de la xarxa determinada pels rangs de valors thv i tv.
[thm,tm]=meshgrid(thv,tv);
(c) Dibuix de la superfície. A continuació, avaluant les components de la parametrització Φ als punts d'aquesta xarxa obtenim tres noves matrius, que utilitzem per fer la representació gràfica amb la funció surf. Aquesta representació estarà formada per petits paral·lelograms.
xm=x(thm,tm);
ym=y(thm,tm);
zm=z(thm,tm);
 
close all
figure, hold on, axis equal
 
surf(xm,ym,zm)
Sovint, el punt de vista que tenim sobre la superfície no serà el més adequat. Normalment el corregirem de manera interactiva aplicant una rotació sobre la figura obtinguda. Una altra possibilitat és, un cop s'ha creat el gràfic, modificar la gca (get current axis) especificant un nou punt de vista. Ho podem fer donant dos angles AZ i EL anomenats azimut i elevació (similars a la longitud i la latitud de l'esfera). Cal expressar aquests angles en graus (no en radians). Podem usar la instrucció set(gca, 'View', [AZ,EL]) o, simplement, view(AZ,EL).
view(75,20) %per exemple
(d) Dibuix de la vora de la superfície. Podem completar el gràfic amb la vora o contorn. Per a la nostra superfície, són 4 corbes que obtenim fixant, en cada cas, un dels paràmetres com el seu valor inicial o final, i que podem dibuixar amb la funció plot3.
% fixant th=0
xc1=x(0,tv);
yc1=y(0,tv);
zc1=z(0,tv);
plot3(xc1,yc1,zc1,'r','LineWidth',3)
%amb 'r' dibuixem la corba de color vermell
%i amb 'LineWidth' indiquem el gruix (3) de la corba
 
% fixant th=pi
xc2=x(pi,tv);
yc2=y(pi,tv);
zc2=z(pi,tv);
plot3(xc2,yc2,zc2,'r','LineWidth',3)
 
% fixant t=0
xc3=x(thv,0);
yc3=y(thv,0);
zc3=z(thv,0*thv);
%com que la component z no depen de theta, fem aixo per obtenir
%un vector zc3 de la mateixa longitud que xc3 i yc3
plot3(xc3,yc3,zc3,'r','LineWidth',3)
 
% fixant t=2*pi/3
xc4=x(thv,2*pi/3);
yc4=y(thv,2*pi/3);
zc4=z(thv,2*pi/3+0*thv);
plot3(xc4,yc4,zc4,'r','LineWidth',3)
(e) Dibuix del camp vectorial i del vector normal. El pas següent és representar el camp vectorial sobre una xarxa de punts de la superfície, usant la funció quiver3 (ens trobem a l'espai 3D, al pla usaríem quiver).
% definim nous rangs de valors per tal que el grafic no sigui massa dens,
% considerant per exemple 1 de cada 5 punts considerats anteriorment
thv5=linspace(0,pi,N1/5+1);
tv5=linspace(0,2*pi/3,N2/5+1);
[thm5,tm5]=meshgrid(thv5,tv5);
 
% nova xarxa de punts, menys densa
xm5=x(thm5,tm5);
ym5=y(thm5,tm5);
zm5=z(thm5,tm5);
 
% avaluem les components del camp vectorial a la xarxa obtinguda
Pm5=P(xm5,ym5,zm5);
Qm5=Q(xm5,ym5,zm5);
Rm5=R(xm5,ym5,zm5);
 
% dibuixem els vectors, de longitud multiplicada pel factor 0.5
% (tambe en podem variar el gruix amb 'LineWidth')
quiver3(xm5,ym5,zm5,Pm5,Qm5,Rm5,0.5,'Color','g')
Recordem que el flux a través d'una superfície depèn de com estigui orientada, en aquest cas per la normal exterior al tor. Perquè la fórmula del flux sigui vàlida cal que el vector normal tingui aquesta mateixa orientació. Si no, haurem de canviar el signe al flux obtingut.
Així doncs, representem el vector només en un punt concret, per tal de verificar-ne l'orientació. Escollim, per exemple, el punt "central" de la parametrització,
th0=pi/2; t0=pi/3;
x0=x(th0,t0);
y0=y(th0,t0);
z0=z(th0,t0);
nx0=nx(th0,t0);
ny0=ny(th0,t0);
nz0=nz(th0,t0);
quiver3(x0,y0,z0,nx0,ny0,nz0,0.75,'Color','k')
(si cal, podem tornar a fer una rotació de la figura o usar view, per tal de visualitzar millor el vector normal)
Així doncs, aquesta parametrització orienta el tor per la normal exterior, com volíem.
(f) Càlcul del flux. Hem de calcular la integral doble, sobre el domini D considerat, de la funció
Recordem que a l'apartat (e) hem vist que la parametrització Φ ja té l'orientació demenada (si no fos així, hauríem de canviar el signe al flux obtingut).
Així doncs, definim g i en calculem la integral amb la funció integral2, la qual cosa ens dóna el flux aproximat.
g=@(th,t)...
P(x(th,t),y(th,t),z(th,t)).*nx(th,t)...
+Q(x(th,t),y(th,t),z(th,t)).*ny(th,t)...
+R(x(th,t),y(th,t),z(th,t)).*nz(th,t);
 
fluxaprox=integral2(g,0,pi,0,2*pi/3)
fluxaprox = 3.9978
Podem veure l'error comès comparant amb el valor exacte:
flux=(8*pi^2+9*pi*sqrt(3))/32;
errflux=abs(flux-fluxaprox)
errflux = 3.7346e-11
Exemple 2. Calculem la massa de la porció de paraboloide hiperbòlic definida per amb pertanyent al domini suposant una funció densitat
Com que es tracta d'una gràfica podem parametritzar-la per les pròpies coordenades és a dir El vector normal és que té norma
Comencem definint la funció i les seves derivades parcials, la norma del vector normal, i la densitat
z=@(x,y) x.^2-y.^2;
 
Dzx=@(x,y) 2*x;
Dzy=@(x,y) -2*y;
 
nn=@(x,y) sqrt(1+Dzx(x,y).^2+Dzy(x,y).^2);
 
rho=@(x,y,z) z;
Tant per a la representació gràfica com per a la integral, hem de tenir en compte que en aquest cas el domini D no és un rectangle, sinó un triangle.
Per a la gràfica, començarem definint una xarxa de punts sobre el rectangle que conté Per als punts d'aquesta xarxa que pertanyen a calcularem segons la fórmula. I, per als que no pertanyen a donarem a z un "valor" NaN (not a number), amb la qual cosa aquests punts no seran representats.
Per diferenciar els punts que pertanyen al domini i els que no, definim una "funció lògica" que val 1 als punts de D i val 0 als punts que no hi pertanyen (s'escriu i rep el nom de funció característica de ),
chi=@(x,y) ( x<=2 & -x<=y & y<=x );
%podem comprovar, per exemple, els valors de chi(1,0.9) i chi(1,1.1)
Ara escollim rangs de valors per a les coordenades i definim una xarxa de punts sobre el rectangle, obtenint dues matrius xm i ym,
N1=30; xv=linspace(0,2,N1+1);
N2=60; yv=linspace(-2,2,N2+1);
 
[xm,ym]=meshgrid(xv,yv);
Avaluem a tots els punts de la xarxa obtenint una matriu zm però, tot seguit, redefinim com a NaN els elements de la matriu que corresponen a punts que NO pertanyen al domini.
zm=z(xm,ym);
zm( ~chi(xm,ym) )=NaN;
Ara representem gràficament la superfície,
close all
figure, hold on, axis equal
 
surf(xm,ym,zm)
 
view(15,40)
Nota: Els talls que podem observar a les vores inferiors de la superfície són deguts al fet que el domini no és rectangular. Una manera de disminuir aquest efecte és considerar considerar una xarxa de punts més densa escollint i més grans.
I dibuixem també la vora de la superfície, formada per 3 corbes,
xc1=xv;
yc1=xc1;
zc1=z(xc1,yc1);
plot3(xc1,yc1,zc1,'r','LineWidth',3)
 
xc2=xv;
yc2=-xc2;
zc2=z(xc2,yc2);
plot3(xc2,yc2,zc2,'r','LineWidth',3)
 
yc3=yv;
xc3=2+0*yv;
%com que la component x es constant sobre el segment vertical,
%fem aixo per obtenir un vector xc3 de la mateixa longitud que yc3
zc3=z(xc3,yc3);
plot3(xc3,yc3,zc3,'r','LineWidth',3)
Finalment calculem, en aquest exemple, la massa de la superfície com la integral de la funció densitat sobre la superfície, Així, definirem la funció i en calculem la seva integral doble sobre
g=@(x,y) rho(x,y,z(x,y)).*nn(x,y);
 
ymin=@(x) -x;
ymax=@(x) x;
massaaprox=integral2(g,0,2,ymin,ymax)
massaaprox = 19.4215

Apèndix: Integració en 2 variables

Veurem com calcular una integral doble
amb dominis D de diferents tipus, usant la funció integral2 del Matlab (per a integrals triples tenim la funció integral3 que funciona de manera anàloga).
Cas 1: Domini rectangular
La funció integral2 de Matlab requereix com a inputs la funció a integrar, i els límits d'integració definits pel domini, que en tractar-se d'un domini rectangular serien respectivament. A més, opcionalment podem afegir-hi altres opcions com toleràncies d'error absolut i relatiu, 'AbsTol', 'RelTol' (si no les afegim s'empraran els valors per defecte i respectivament),
g=@(u,v) ...;
I=integral2(g,a,b,c,d, 'AbsTol',1e-12, 'RelTol',1e-12)
Nota: Una opció alternativa per definir la funció a integrar és, en comptes d'usar una ordre del tipus g=@(u,v) ..., fer-ho a través d'un fitxer de Matlab, anomenat per exemple gfun.m, que contingui un codi tipus
function z=gfun(u,v)
z=...;
En aquest cas, per calcular la integral escriuríem integral2(@gfun,a,b,c,d,opcions) (afegint "@" abans del nom del fitxer).
Exemple: Calculem el valor de la integral següent, sobre el rectangle
També compararem amb el valor exacte, que és
g1=@(u,v) 1./(u+2*v).^2;
I1aprox=integral2(g1,2,5,1,3)
I1aprox = 0.1206
%o bé, si volem especificar les tolerancies d'error,
% I1aprox=integral2(g1,2,5,1,3,'AbsTol',1e-12,'RelTol',1e-12)
 
I1=log(14/11)/2; %valor exacte
err1=abs(I1-I1aprox)
err1 = 3.4577e-13
Cas 2: Domini de la forma que correspon a una integral iterada
Per a aquest cas, haurem de començar definint els extrems d'integració que depenen de u com a funcions,
vmin=@(u) ...;
vmax=@(u) ...;
g=@(u,v) ...;
I=integral2(g,a,b,vmin,vmax)
Exemple: Calculem la integral doble sobre el domini El valor exacte és
g2=@(u,v) (u.^2)./(v.^2);
vmin=@(u) 1./u;
vmax=@(u) u;
I2aprox=integral2(g2,1,2,vmin,vmax)
I2aprox = 2.2500
 
I2=9/4;
err2=abs(I2-I2aprox)
err2 = 1.3823e-09
S'ha de tenir en compte que integral2 només permet que els extrems de la segona variable siguin funció de la primera variable, i no al revés. Aleshores, en el cas d'un domini de la forma haurem de permutar les variables a la funció inicial.
Exemple: Calculem la integral doble sobre el domini El valor exacte és
g3=@(v,u) u+2*v; %definim la funcio amb les variables permutades
umin=@(v) v.^2-4;
I3aprox=integral2(g3,-3,3,umin,5)
I3aprox = 50.4000
 
I3=252/5;
err3=abs(I3-I3aprox)
err3 = 5.8608e-09
________
(c) Numerical Factory (by Ana Domínguez & Pere Gutiérrez)