Verify the Divergence theorem in MATLAB

16 次查看（过去 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 个评论
显示 3更早的评论
Camstien
Camstien 2022-11-25
Ideally symbolically (i.e. with syms x y z)
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 个评论
显示 4更早的评论
Torsten
Torsten 2022-11-27
编辑：Torsten 2022-11-27
Ran in:
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
Paul 2022-11-27
Thanks, I knew I was missing something.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Number Theory 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by