# Verify the Divergence theorem in MATLAB

14 views (last 30 days)

Show older comments

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

### See Also

### Categories

### Products

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!