Verify the Divergence theorem in MATLAB
25 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
I’m trying to verify the divergence theorem using MATLAB if F(x, y, z) =<x^2y^2, y^2z^2 , z^2x^2> when the solid E is bounded by x^2 + y^2 + z^2 = 2 on top and z = (x^2+y^2)^1/2 on the bottom. I know the function for divergence in MATLAB but not sure how to calculate the bounds. Thank you in advance.
5 Kommentare
William Rose
am 25 Nov. 2022
I forgot to actually show the image of the surface, which is generated by the code I posted above. Here it is.
r=0:.1:1; t=0:pi/20:2*pi; [R,T]=meshgrid(r,t); X=R.*cos(T); Y=R.*sin(T);
Ztop=(2-X.^2-Y.^2).^(.5);
Zbot=(X.^2+Y.^2).^(.5);
surf(X,Y,Ztop)
hold on
surf(X,Y,Zbot)
Antworten (2)
William Rose
am 26 Nov. 2022
I now think there was a mistake in my earlier answer. Previously I said the integral for the top surface can be written
but that is wrong, because it is not the case that dS=dy*dx. Rather, it is true that , where ϕ is the elevation angle of the point above the sphere's equatorial plane, and in this case, .
I also believe we can simplify the expression for :
Then the integral for the flux through the top surface is
The Matlab code for this integral is
syms x y z
z=sqrt(2-x^2-y^2);
S_top=int(int((x^3*y^2+y^3*z^2+z^3*x^2)/z,y,-sqrt(1-x^2),sqrt(1-x^2)),x,-1,1)
This is not simple. Matlab evaluated the inner integral symbolically, then gave up. This seems unreasonably hard for a homework problem, which I assume this is. We haven;t even looked at the surface integral over the bottom surface, or the volume integral of div(F). I probably made a mistake somehwere.
Perhaps the goal is to verify the divergence theorem for this vector field and region by numerical approximation. The numerical approach will also be non-trivial.
1 Kommentar
Torsten
am 26 Nov. 2022
Bearbeitet: Torsten
am 26 Nov. 2022
The surfaces/volumes smell as if a coordinate transformation to cylindrical coordinates would be helpful.
For the upper part of the solid:
%Volume part
x = a*sqrt(2-z^2)*cos(phi)
y = a*sqrt(2-z^2)*sin(phi)
z = z
for
0 <= a <= 1
0 <= phi <= 2*pi
1 <= z <= sqrt(2)
% Surface part
x = sqrt(2-z^2)*cos(phi)
y = sqrt(2-z^2)*sin(phi)
z = z
for
0 <= phi <= 2*pi
1 <= z <= sqrt(2)
For the lower part of the solid:
% Volume part
x = a*z*cos(phi)
y = a*z*sin(phi)
z = z
for
0 <= a <= 1
0 <= phi <= 2*pi
0 <= z <= 1
% Surface part
x = z*cos(phi)
y = z*sin(phi)
z = z
for
0 <= phi <= 2*pi
0 <= z <= 1
Torsten
am 26 Nov. 2022
Bearbeitet: Torsten
am 26 Nov. 2022
Nice exercise.
syms X Y Z x y z
syms a phi positive
F = [X^2*Y^2,Y^2*Z^2,Z^2*X^2];
divF = diff(F(1),X)+diff(F(2),Y)+diff(F(3),Z);
% Volume integrals
% Volume integral upper part
x = a*sqrt(2-z^2)*cos(phi);
y = a*sqrt(2-z^2)*sin(phi);
z = z;
assume (z>=1 & z <= sqrt(2));
J = simplify(abs(det([diff(x,a) diff(x,phi) diff(x,z);diff(y,a) diff(y,phi) diff(y,z);diff(z,a) ,diff(z,phi), diff(z,z)])));
f = J*subs(divF,[X Y Z],[x y z]);
I1 = int(f,a,0,1);
I2 = int(I1,phi,0,2*pi);
IV_upper = int(I2,z,1,sqrt(2))
%Volume integral lower part
x = a*z*cos(phi);
y = a*z*sin(phi);
z = z;
assume (z>=0 & z <= 1);
J = simplify(abs(det([diff(x,a) diff(x,phi) diff(x,z);diff(y,a) diff(y,phi) diff(y,z);diff(z,a) ,diff(z,phi), diff(z,z)])));
f = J*subs(divF,[X Y Z],[x y z]);
I1 = int(f,a,0,1);
I2 = int(I1,phi,0,2*pi);
IV_lower = int(I2,z,0,1)
% Sum of volume integrals upper and lower part
intV = IV_upper + IV_lower
% Surface integrals
% Surface integral upper part
x = sqrt(2-z^2)*cos(phi);
y = sqrt(2-z^2)*sin(phi);
z = z;
assume (z>=1 & z <= sqrt(2));
J = simplify(cross([diff(x,phi);diff(y,phi);diff(z,phi)],[diff(x,z);diff(y,z);diff(z,z)]));
J = simplify(sqrt((J(1)^2+J(2)^2+J(3)^2)));
f = J*subs(F,[X Y Z],[x y z])*[x;y;z]/simplify(sqrt(x^2+y^2+z^2));
I1 = int(f,phi,0,2*pi);
IS_upper = int(I1,z,1,sqrt(2))
% Surface integral lower part
x = z*cos(phi);
y = z*sin(phi);
z = z;
assume (z>=1 & z <= sqrt(2));
J = simplify(cross([diff(x,phi);diff(y,phi);diff(z,phi)],[diff(x,z);diff(y,z);diff(z,z)]));
J = simplify(sqrt((J(1)^2+J(2)^2+J(3)^2)));
f = J*subs(F,[X Y Z],[x y z])*[x;y;-sqrt(x^2+y^2)]/simplify(sqrt(2*(x^2+y^2)));
I1 = int(f,phi,0,2*pi);
IS_lower = int(I1,z,0,1)
% Sum of surface integrals upper and lower part
intS = IS_upper + IS_lower
6 Kommentare
Torsten
am 27 Nov. 2022
Bearbeitet: Torsten
am 27 Nov. 2022
divF has no symmetries that you could exploit. You will need to compute over 0:2*pi for theta.
syms x y z r theta real
F(x,y,z) = [x^2*y^2 , y^2*z^2 , z^2*x^2];
divF(x,y,z) = divergence(F,[x y z]);
IV_upper = int(int(int(divF(r*cos(theta),r*sin(theta),z)*r,theta,0,sym(2*pi)),r,0,sqrt(2-z^2)),z,1,sqrt(sym(2)))
Siehe auch
Kategorien
Mehr zu Number Theory finden Sie in Help Center und File Exchange
Produkte
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!