x=@(th,s) (R+s.*cos(th/2)).*cos(th);
y=@(th,s) (R+s.*cos(th/2)).*sin(th);
Dxth=@(th,s) -s/2.*sin(th/2).*cos(th)-(R+s.*cos(th/2)).*sin(th);
Dyth=@(th,s) -s/2.*sin(th/2).*sin(th)+(R+s.*cos(th/2)).*cos(th);
Dzth=@(th,s) s/2.*cos(th/2);
Dxs=@(th,s) cos(th/2).*cos(th);
Dys=@(th,s) cos(th/2).*sin(th);
nx=@(th,s) Dyth(th,s).*Dzs(th,s)-Dzth(th,s).*Dys(th,s);
ny=@(th,s) Dzth(th,s).*Dxs(th,s)-Dxth(th,s).*Dzs(th,s);
nz=@(th,s) Dxth(th,s).*Dys(th,s)-Dyth(th,s).*Dxs(th,s);
nn=@(th,s) sqrt(nx(th,s).^2+ny(th,s).^2+nz(th,s).^2);
area=integral2(nn,0,2*pi,-A/2,A/2);
fprintf('area = %.6f\n',area)
% rangs de valors dels parametres i xarxa de punts
N1=80; thv=linspace(0,2*pi,N1+1);
N2=10; sv=linspace(-A/2,A/2,N2+1);
[thm,sm]=meshgrid(thv,sv);
xm=x(thm,sm); ym=y(thm,sm); zm=z(thm,sm);
% representacio grafica de la superficie
figure, hold on, axis equal
xc1=x(thv,-A/2); yc1=y(thv,-A/2); zc1=z(thv,-A/2);
plot3(xc1,yc1,zc1,'r','LineWidth',3)
xc2=x(thv,A/2); yc2=y(thv,A/2); zc2=z(thv,A/2);
plot3(xc2,yc2,zc2,'r','LineWidth',3)
% Nota: podem dibuixar la vora com una unica corba, considerant
% per a theta un rang de valors de 0 a 4*pi
%vectors normals sobre la circumferencia central (un de cada dos punts)
thv0=linspace(0,2*pi,N1/2+1);
xv0=x(thv0,0); yv0=y(thv0,0); zv0=z(thv0,0);
nxv0=nx(thv0,0); nyv0=ny(thv0,0); nzv0=nz(thv0,0);
quiver3(xv0,yv0,zv0,nxv0,nyv0,nzv0,0.5,'Color','k')
% funcio que defineix la grafica i element de superficie
h=@(x,y) (x.^2-y.^2)/(2*a);
Dhx=@(x,y) x/a; Dhy=@(x,y) -y/a;
nn=@(x,y) sqrt(1+Dhx(x,y).^2+Dhy(x,y).^2);
% (per simetria, ens restringim al primer quadrant i multipliquem per 4)
g=@(x,y) rho(x,y,h(x,y)).*nn(x,y); %funcio a integrar
gpolars=@(r,th) g(r.*cos(th),r.*sin(th)).*r;
%hem passat a polars i multiplicat pel jacobia
massa=4*integral2(gpolars,0,1,0,pi/2);
fprintf('massa = %.6f\n',massa)
% xarxa de punts sobre el rectangle
N=50; xv=linspace(-1.5*a,1.5*a,N+1); yv=xv;
% funcio logica que ens diu si un punt es troba dins el cilindre
chi=@(x,y) ( x.^2+y.^2<=a^2 );
% fem una copia zm0 de la matriu zm, a la qual canviem a NaN
% els termes que corresponen a punts del parabolide hiperbolic
% que no pertanyen al cilindre
% (guardem zm per dibuixar el paraboloide sencer)
% representem el parabolide hiperbolic dos cops,
% primer usant la matriu zm0, i despres usant la matriu original zm
% pero aquesta darrera amb linies de punts i un grau de transparencia
% definit per 'FaceAlpha'
figure, hold on, axis equal
surf(xm,ym,zm,'FaceAlpha',0.2,'LineStyle',':')
N2=100; thv2=linspace(0,2*pi,N2+1);
xc=a*cos(thv2); yc=a*sin(thv2);
plot3(xc,yc,zc,'r','LineWidth',3)
% dibuixem el cilindre, tambe amb linies de punts i transparent
xc=@(th,z) a*cos(th); yc=@(th,z) a*sin(th);
N3=20; thcv=linspace(0,2*pi,N3+1);
N4=20; zcv=linspace(-1.5*a,1.5*a,N4+1);
[thcm,zcm]=meshgrid(thcv,zcv);
xcm=xc(thcm,zcm); ycm=yc(thcm,zcm);
surf(xcm,ycm,zcm,'FaceAlpha',0.2,'LineStyle',':')
% parametritzacio de la superficie i les seves derivades
x=@(r,th) r.*cos(th); y=@(r,th) r.*sin(th); z=@(r,th) 0.75*th;
Dxr=@(r,th) cos(th); Dyr=@(r,th) sin(th); Dzr=@(r,th) 0;
Dxth=@(r,th) -r.*sin(th); Dyth=@(r,th) r.*cos(th); Dzth=@(r,th) 0.75;
nx=@(r,th) Dyr(r,th).*Dzth(r,th)-Dzr(r,th).*Dyth(r,th);
ny=@(r,th) Dzr(r,th).*Dxth(r,th)-Dxr(r,th).*Dzth(r,th);
nz=@(r,th) Dxr(r,th).*Dyth(r,th)-Dyr(r,th).*Dxth(r,th);
P=@(x,y,z) x.^2+y+z; Q=@(x,y,z) cos(x+z); R=@(x,y,z) y.^2+2*x.*z;
% definim la funcio g a integrar i calculem el flux aproximat
P(x(r,th),y(r,th),z(r,th)).*nx(r,th)...
+Q(x(r,th),y(r,th),z(r,th)).*ny(r,th)...
+R(x(r,th),y(r,th),z(r,th)).*nz(r,th);
flux=integral2(g,0,5,0,2*pi);
% calculem el vector normal en el punt (r,theta)=(2.5,pi) per tal
% de verificar-ne l'orientacio i, en cas que no sigui la que es demana,
% canviarem el signe del flux
x0=x(r0,th0); y0=y(r0,th0); z0=z(r0,th0);
nx0=nx(r0,th0); ny0=ny(r0,th0); nz0=nz(r0,th0);
fprintf('flux = %.6f\n',flux)
N1=20; rv=linspace(0,5,N1+1);
N2=40; thv=linspace(0,2*pi,N2+1);
[rm,thm]=meshgrid(rv,thv);
xm=x(rm,thm); ym=y(rm,thm); zm=z(rm,thm);
figure, hold on, axis equal
% representem el vector normal nomes en el punt (r,theta)=(2.5,pi)
% per tal de verificar-ne l'orientacio
% (com podem veure, ja tenia l'orientacio demanada)
quiver3(x0,y0,z0,nx0,ny0,nz0,'Color','r')
% la superficie es una grafica z=z(x,y), definim aquesta funcio
% i les seves derivades parcials
Dzx=@(x,y) 4*x; Dzy=@(x,y) 2*y;
P=@(x,y) x.*y.^2; Q=@(x,y) -2*x.^2.*y; R=@(x,y) x.^2+y.^2;
% definim la funcio a integrar com a producte escalar del vector F=(P,Q,R)
% pel vector normal (-Dzx,-Dzy,1), el qual te l'orientacio demanada
g=@(x,y) -P(x,y).*Dzx(x,y)-Q(x,y).*Dzy(x,y)+R(x,y);
flux=integral2(g,0,1,0,ymax);
fprintf('flux = %.6f\n',flux) %(el valor exacte es 1/6)
% discretitzacio del domini dels parametres
N=50; xv=linspace(0,1,N+1); yv=xv;
chi=@(x,y) ( x>=0 & y>=0 & x+y<=1 );
% trobem la coordenada z dels punts de la superficie discretitzada
% i la redefinim com a NaN fora del domini
% dibuix de la superfície
figure, hold on, axis equal
% dibuix de la part x=0 del contorn
plot3(xc1,yc1,zc1,'r','LineWidth',3)
% dibuix de la part y=0 del contorn
plot3(xc2,yc2,zc2,'r','LineWidth',3)
% dibuix de la part x+y=1 del contorn
plot3(xc3,yc3,zc3,'r','LineWidth',3)
x=@(th,z1) A*cos(th); y=@(th,z1) A*sin(th); z=@(th,z1) z1;
Dxth=@(th,z1) -A*sin(th); Dyth=@(th,z1) A*cos(th); Dzth=@(th,z1) 0;
Dxz=@(th,z1) 0; Dyz=@(th,z1) 0; Dzz=@(th,z1) 1;
% tot i que no es estrictament necessari, hem diferenciat entre
% la component z de la parametritzacio (que hem anomenat 'z')
% i el parametre z (que hem anomenat 'z1')
nx=@(th,z1) Dyth(th,z1).*Dzz(th,z1)-Dzth(th,z1).*Dyz(th,z1);
ny=@(th,z1) Dzth(th,z1).*Dxz(th,z1)-Dxth(th,z1).*Dzz(th,z1);
nz=@(th,z1) Dxth(th,z1).*Dyz(th,z1)-Dyth(th,z1).*Dxz(th,z1);
nn=@(th,z1) sqrt(nx(th,z1).^2+ny(th,z1).^2+nz(th,z1).^2);
% de fet, podiem haver escrit directament "nn=@(th,z1) A;",
% que es l'element de superficie del cilindre, prou conegut
(x(th,z1).^2+y(th,z1).^2).*rho(x(th,z1),y(th,z1),z(th,z1)).*nn(th,z1);
Iz=2*integral2(g,0,pi,zmin,zmax);
fprintf('Iz = %.6f\n',Iz)
% rangs de valors dels parametres i xarxa de punts
N1=60; thv=linspace(0,2*pi,N1+1);
N2=50; zv=linspace(-2*A,2*A,N2+1);
[thm,zm]=meshgrid(thv,zv);
xm=x(thm,zm); ym=y(thm,zm);
% funcio logica corresponent a S1
chi=@(x,y,z) ( x.^2+z.^2<=A^2 & y>=0 );
% redefinim una copia de la matriu zm, per a S1
zm0( ~chi(xm,ym,zm) )=NaN;
% dibuix de S1, S2 i els dos cilindres
% (notem que podem dibuixar el cilindre horitzontal a partir del vertical,
% simplement intercanviant les coordenades y,z)
figure, hold on, axis equal
surf(xm,ym,zm0) %dibuix de S1
surf(xm,-ym,zm0) %dibuix de S2
surf(xm,ym,zm,'FaceAlpha',0.05,'LineStyle',':') %cilindre vertical
surf(xm,zm,ym,'FaceAlpha',0.05,'LineStyle',':') %cilindre horitzontal
% dibuix de les vores de S1 i S2
zc=A*sin(thv); xc=x(thv,zc); yc=y(thv,zc);
plot3(xc,yc,-zc,'r','LineWidth',3) %vora inferior de S1
plot3(xc,yc,zc,'r','LineWidth',3) %vora superior de S1
plot3(xc,-yc,-zc,'r','LineWidth',3) %vora inferior de S2
plot3(xc,-yc,zc,'r','LineWidth',3) %vora superior de S2
% camp vectorial F=(P,Q,R)
P=@(x,y,z) x.*z; Q=@(x,y,z) y.*z; R=@(x,y,z) z.^2;
% funcio a integrar, com a producte escalar del camp vectorial,
% avaluat sobre cada punt (x,y,h(x,y)) de la grafica,
% pel vector normal (-Dzx,-Dzy,1), el qual te l'orientacio demanada
g=@(x,y) -P(x,y,h(x,y)).*Dhx(x,y)-Q(x,y,h(x,y)).*Dhy(x,y)+R(x,y,h(x,y));
% calculem el flux com la integral de sobre el domini D=[0,1]x[0,1]
flux1=integral2(g,0,1,0,1);
fprintf('flux1 = %.6f\n',flux1) %(el valor exacte es 1/450)
% primera corba (la corba generatriu)
x1=@(t) 3+cos(t); y1=@(t) 0; z1=@(t) sin(t); %orientacio no compatible
Dx1=@(t) -sin(t); Dy1=@(t) 0; Dz1=@(t) cos(t);
x2=@(t) 0; y2=@(t) 3+cos(t); z2=@(t) sin(t); %orientacio compatible
Dx2=@(t) 0; Dy2=@(t) -sin(t); Dz2=@(t) cos(t);
P=@(x,y,z) x.^2.*z; Q=@(x,y,z) -2*y.^2.*z; R=@(x,y,z) x.*y;
% funcions a integrar (entre 0 i 2*pi)
h1=@(t) P(x1(t),y1(t),z1(t)).*Dx1(t)+Q(x1(t),y1(t),z1(t)).*Dy1(t)...
+R(x1(t),y1(t),z1(t)).*Dz1(t);
h2=@(t) P(x2(t),y2(t),z2(t)).*Dx2(t)+Q(x2(t),y2(t),z2(t)).*Dy2(t)...
+R(x2(t),y2(t),z2(t)).*Dz2(t);
% calcul del flux com a suma de les dues circulacions
% (canviem el signe de la primera, perque l'orientacio no es compatible
% amb la de la superficie)
flux1=-integral(h1,0,2*pi)+integral(h2,0,2*pi);
fprintf('flux1 = %.6f\n',flux1) %(el valor exacte es 111*pi/4)
% parametritzacio de la superficie
x=@(th,t) x1(t).*cos(th);
y=@(th,t) x1(t).*sin(th);
% derivades parcials de la parametritzacio
Dxth=@(th,t) -x1(t).*sin(th);
Dyth=@(th,t) x1(t).*cos(th);
Dxt=@(th,t) Dx1(t).*cos(th);
Dyt=@(th,t) Dx1(t).*sin(th);
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);
Prot=@(x,y,z) x+2*y.^2; Qrot=@(x,y,z) x.^2-y; Rrot=@(x,y,z) 0;
Prot(x(th,t),y(th,t),z(th,t)).*nx(th,t)...
+Qrot(x(th,t),y(th,t),z(th,t)).*ny(th,t)...
+Rrot(x(th,t),y(th,t),z(th,t)).*nz(th,t);
flux2=integral2(g,0,pi/2,0,2*pi);
fprintf('flux2 = %.6f\n',flux2)
N1=15; thv=linspace(0,pi/2,N1+1);
N2=40; tv=linspace(0,2*pi,N2+1);
[thm,tm]=meshgrid(thv,tv);
xm=x(thm,tm); ym=y(thm,tm); zm=z(thm,tm);
figure, hold on, axis equal
surf(xm,ym,zm,'LineStyle',':')
xc1=x1(tv); yc1=0*tv; zc1=z1(tv);
plot3(xc1,yc1,zc1,'r','LineWidth',2)
xc2=0*tv; yc2=y2(tv); zc2=z2(tv);
plot3(xc2,yc2,zc2,'r','LineWidth',2)
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.5,...
'Color','k','LineWidth',2,'MaxHeadSize',0.8)
x01=x1(pi/2); y01=y1(pi/2); z01=z1(pi/2);
tx01=Dx1(pi/2); ty01=Dy1(pi/2); tz01=Dz1(pi/2);
quiver3(x01,y01,z01,-tx01,-ty01,-tz01,1.25,...
'Color','b','LineWidth',2,'MaxHeadSize',1)
x02=x2(pi/2); y02=y2(pi/2); z02=z2(pi/2);
tx02=Dx2(pi/2); ty02=Dy2(pi/2); tz02=Dz2(pi/2);
quiver3(x02,y02,z02,tx02,ty02,tz02,1.25,...
'Color','b','LineWidth',2,'MaxHeadSize',1)