Verify the Divergence theorem in MATLAB

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.

[Torsten]

You know the formulae on how to compute volume and surface integrals in three dimensions ?

[Camstien]

Yes,

I have for this problem, the surface is the solid E, which is bounded by x^2 + y^2 + z^2 = 2 on top and z = (x^2+y^2)^1/2 on the bottom.

The unit normal to the surface is:

n=<x,y,z>/sqrt(x^2+y^2+z^2)

The dot product of the vectors F and n is:

F dot n=x^2y^2x/sqrt(x^2+y^2+z^2)+y^2z^2y/sqrt(x^2+y^2+z^2)+z^2x^2z/sqrt(x^2+y^2+z^2)

The divergence of the vector F is:

div F=2xy^2+2yz^2+2xz^2

The surface integral is:

int int (F dot n) dS= int int int (x^2y^2x/sqrt(x^2+y^2+z^2)+y^2z^2y/sqrt(x^2+y^2+z^2)+z^2x^2z/sqrt(x^2+y^2+z^2)) dS

The volume integral is:

int int int div F dV= int int int (2xy^2+2yz^2+2xz^2) dV

The limits of integration for the surface integral are:

0<=z<=(x^2+y^2)^1/2

-sqrt(2)<=x<=sqrt(2)

-sqrt(2)<=y<=sqrt(2)

The limits of integration for the volume integral are:

0<=z<=2

-sqrt(2)<=x<=sqrt(2)

-sqrt(2)<=y<=sqrt(2)

After calculating the integrals, we find that:

int int (F dot n) dS=int int int div F dV

I need the functions to test this in MATLAB

[William Rose]

@Camstien,

The top and bottom meet where : where , i.e. where , which is the unit circle, which is at the level z=1. Here is an image of the volume:

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)

Are you supposed to do this symbolically (i.e. with syms x y z) or numerically?

You gave the unit normal for the top but not for the bottom. YOu will need to do the bottom as well as the top.

The integral for the top surface can be written

Symbolic, top surface:

syms x y z
Fdotn=x^3*y^2/sqrt(x^2+y^2+z^2)+y^3*z^2/sqrt(x^2+y^2+z^2)+z^3*x^2/sqrt(x^2+y^2+z^2);
z=sqrt(2-x^2-y^2);
S=int(int(Fdotn,y,-sqrt(1-x^2),sqrt(1-x^2)),x,-1,1)

I tried the code above. Matlab has been working on it for many minutes and is still busy. I am not sure that my code above, for S, is correct. Maybe we need to try numerical integration.

[Camstien]

Ideally symbolically (i.e. with syms x y z)

[William Rose]

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)

Answer 1

[William Rose]

@Camstien

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)

S_top = 

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.

[Torsten]

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

Answer 2

[Torsten]

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))

IV_upper = 

%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)

IV_lower = 

% Sum of volume integrals upper and lower part
intV = IV_upper + IV_lower

intV = 

% 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))  

IS_upper = 

% 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)

IS_lower = 

% Sum of surface integrals upper and lower part
intS = IS_upper + IS_lower

intS = 

[William Rose]

@Torsten, wow.

[Paul]

@Torsten,

Could save a couple of keystrokes using divergence and jacobian.

More importantly to me, can you explain the coordinates you're using for the top half of the volume integral? If I use cylindrical coordinates x = r*cos(theta), y = r*sin(theta), z = z, I get

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 = 4*int(int(int(divF(r*cos(theta),r*sin(theta),z)*r,theta,0,sym(pi)/2),r,0,sqrt(2-z^2)),z,1,sqrt(sym(2)))

IV_upper = 

[Torsten]

More importantly to me, can you explain the coordinates you're using for the top half of the volume integral?

I wrote them explicitly in my comment above. The coordinate transformations are only similar to cylindrical coordinates.

[Paul]

Yes, I saw the definitions.

When I use cylindrical coordinates, I get a different result for IV_upper and I'd like to reconcile the results

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 = 4*int(int(int(divF(r*cos(theta),r*sin(theta),z)*r,theta,0,sym(pi)/2),r,0,sqrt(2-z^2)),z,1,sqrt(sym(2)))

IV_upper = 

[Torsten]

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)))

IV_upper = 

[Paul]

Thanks, I knew I was missing something.