Verify the Divergence theorem in MATLAB
17 views (last 30 days)
Show older comments
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 Comments
William Rose
on 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)
Answers (2)
William Rose
on 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 Comment
Torsten
on 26 Nov 2022
Edited: Torsten
on 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
on 26 Nov 2022
Edited: Torsten
on 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 Comments
Torsten
on 27 Nov 2022
Edited: Torsten
on 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)))
See Also
Categories
Find more on Number Theory in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!