MATLAB Answers

Verify the Divergence theorem in MATLAB

12 次查看（过去 30 天）
Camstien
Camstien ，2022-11-24
评论： Paul ，2022-11-27
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 个评论
显示 4更早的评论
William Rose
William Rose 2022-11-25
Ran in:
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)

请先登录，再进行评论。

请先登录，再回答此问题。

回答（2 个）

William Rose
William Rose 2022-11-26
Ran in:
@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.
  1 个评论
Torsten
Torsten 2022-11-26
编辑：Torsten 2022-11-26
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
Torsten 2022-11-26
编辑：Torsten 2022-11-26
Ran in:
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 = 
  6 个评论
显示 5更早的评论
Paul
Paul 2022-11-27
Thanks, I knew I was missing something.

请先登录，再进行评论。

类别

Find more on Calculus in Help Center and File Exchange

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by