Hello everyone.I'm new to the community, please help me to use matlab to calculate the volume of this 3D Shape

조회 수: 9(최근 30일)
Phi Dao Anh
Phi Dao Anh 2022년 4월 24일
편집: David Goodmanson 2022년 5월 4일
clc
u=linspace(0,1,30);
v=linspace(0,2*pi,30);
[u,v]=meshgrid(u,v);
x=(1-u).*(3+cos(v)).*cos(4.*pi.*u);
y=(1-u).*(3+cos(v)).*sin(4.*pi.*u);
z=3.*u+(1-u).*sin(v);
surf(x,y,z);
xlabel('x')
ylabel('y')
zlabel('z')

채택된 답변

Bruno Luong
Bruno Luong 2022년 4월 25일
편집: Bruno Luong 2022년 4월 25일
u=linspace(0,1,31);
v=linspace(0,2*pi,37);
[u,v]=meshgrid(u,v);
x=(1-u).*(3+cos(v)).*cos(4.*pi.*u);
y=(1-u).*(3+cos(v)).*sin(4.*pi.*u);
z=3.*u+(1-u).*sin(v);
V = [x(:) y(:) z(:)];
[m,n]=size(u);
K = reshape(1:m*n,[m,n]);
V1 = K(1:m-1,1:n-1);
V2 = K(2:m, 1:n-1);
V3 = K(2:m, 2:n);
V4 = K(1:m-1,2:n);
F0 = [V3(:) V2(:) V1(:);
V1(:) V4(:) V3(:)];
V1 = K(2:m, 1);
V2 = K(1:m-1,1);
V3 = K(ones(m-1,1),1);
Fl = [V3(:) V2(:) V1(:)];
V1 = K(1:m-1, n);
V2 = K(2:m,n);
V3 = K(ones(m-1,1),n);
Fr = [V3(:) V2(:) V1(:)];
F = [F0; Fl; Fr];
% Volume estimated
VF = permute(reshape(V(F,:),[size(F) 3]),[3 1 2]);
Vol = 1/6*sum(dot(cross(VF(:,:,1),VF(:,:,2),1),VF(:,:,3),1))
Vol = 28.5915
  댓글 수: 8
Bruno Luong
Bruno Luong 2022년 4월 27일
I did also my homework. Here is the continuous surface integral to find the volume. I do numerically rather than symbolic (a little messy since I don't see clearly a simplification), but it seems working just fine:
M=3;
N=4*pi;
P=3;
x=@(u,v) (1-u).*(M+cos(v)).*cos(N*u);
y=@(u,v) (1-u).*(M+cos(v)).*sin(N*u);
z=@(u,v) P.*u+(1-u).*sin(v);
dxdu=@(u,v) (M+cos(v)).*(-cos(N*u)-N*(1-u).*sin(N*u));
dydu=@(u,v) (M+cos(v)).*(-sin(N*u)+N*(1-u).*cos(N*u));
dzdu=@(u,v) P-sin(v);
dxdv=@(u,v) (1-u).*(-sin(v)).*cos(N*u);
dydv=@(u,v) (1-u).*(-sin(v)).*sin(N*u);
dzdv=@(u,v) (1-u).*cos(v);
crossgrad=@(u,v) cross([dxdu(u,v);dydu(u,v);dzdu(u,v)],...
[dxdv(u,v);dydv(u,v);dzdv(u,v)],1);
dxyz=@(u,v) [x(u,v);y(u,v);z(u,v)]-[x(0,0);y(0,0);z(0,0)];
dVhelper=@(u,v) dot(crossgrad(u,v),dxyz(u,v),1);
dV=@(u,v) reshape(dVhelper(u(:).',v(:).'),size(u));
Volume = integral2(dV,0,1,0,2*pi)/3
Volume = 29.6088

댓글을 달려면 로그인하십시오.

추가 답변(2개)

Matt J
Matt J 2022년 4월 26일
편집: Matt J 2022년 4월 26일
The volume can be parametrized by introducing a third variable 0<=w<=1,
x=(1-u).*(3+(w*cos(v))).*cos(4.*pi.*u);
y=(1-u).*(3+(w*cos(v))).*sin(4.*pi.*u);
z=3.*u+(1-u).*(w*sin(v));
Sweeping over w produces concentric horns which partition the volume, as seen by the following,
close all
for w=0.1:0.1:1
plotIt(w);
hold on
end
hold off
view(50,15)
function plotIt(w)
u=linspace(0,1,30);
v=linspace(0,2*pi,30);
[u,v]=meshgrid(u,v);
x=(1-u).*(3+(w*cos(v))).*cos(4.*pi.*u);
y=(1-u).*(3+(w*cos(v))).*sin(4.*pi.*u);
z=3.*u+(1-u).*(w*sin(v));
surf(x,y,z,'FaceAlpha',0.3, 'EdgeColor','none','FaceColor',rand(1,3));
xlabel('x')
ylabel('y')
zlabel('z')
end
To compute the volume, we can therefore integrate the Jacobian determinant of the mapping (x,y,z)-->(u,v,w), which leads to a volume of,
Volume = 0.3333*pi*(sin(2*pi) + 9*pi)
or . This is consistent with findings mentioned by @David Goodmanson, though of course I can't know if he used the same method.
  댓글 수: 2
David Goodmanson
David Goodmanson 2022년 4월 27일
I did use the same method as Matt, where a big key to appreciating the solution is working out the Jacobian, rather than, say, letting some syms engine do it. Since u goes from 0 to 1, it's a bit cleaner to let u --> 1-u and use the general expressions
x = u(N + rcosv) cos(Ku) y = u(N + rcosv) sin(Ku) z = M(1-u) + u rsin(v)
where in this particular case K = 4pi, N=M=3. The third variable r has been inserted, necessary to define a 3d volume. I'll write out J, same as Matt probably did. It is(x,y,x across, d/du, d/dr, d/dv down)
(N+rcosv) cos(Ku) -u(N+rcosv)Ksin(Ku) (N+rcosv)sin(Ku) +u(N+rcosv)Kcos(Ku) -M+rsinv
ucosv cos(Ku) ucos(v) sin(Ku) usinv
-ursin(v) cos(Ku) -ursin(v) sin(Ku) urcosv
Multiply the first col by sin(Ku), the second col by cos(Ku), subtract them, get two zero elements, etc. The result is
J = K(N+rcosv) u^3 r
The v integral from 0 to 2pi integrates out the cosv term, and the result from the three integrations is
2pi 1/4 1/2 KN = (pi/4)KN.
which is 3pi^2 in this case.
oo From the Jacobian, the result does not depend on M. That's consistent because a translation in the z direction is in the plane of the disc generated by r and v and does not sweep out any area.
oo K is the sweep angle around the z axis and does not have to be an integral number of turns.
oo The expression N+rcosv is the distance of a given point from the z axis. The fact that the rcosv term integrates out means that the centroid of the disk is locally at x = 0, so the mean distance from the z axis is N. Then the v integration gives 2 pi N which is consistent with the theorem of Pappus. It seems that Pappus was a pretty smart guy.
oo How does one know to insert r here (same as w for Matt), instead of some other function? It doesn't matter. You could use r^n or any function of some parameter zeta that runs from a to b, such that h(zeta) is continuous, monotonic, differentiable and is 0 at a and 1 at b. In that case the Jacobian and the integration give
Int{a,b} h dh/dzeta dzeta = h^2(zeta) /2 | b,a = 1/2 same as if you use r

댓글을 달려면 로그인하십시오.


Matt J
Matt J 2022년 4월 26일
편집: Matt J 2022년 4월 26일
Just as a hint, it's pretty easy to show that the cross sections of the horn are circle with radius (1-u) and whose center is parametrized in u according to,
xc=3*(1-u)*cos(4.*pi.*u); %"center" of the horn
yc=3*(1-u)*sin(4.*pi.*u);
zc=3*u;
The distance from (xc,yc,yc) to the surface can then be calculated for each (u,v) according to,
x-xc=(1-u)*cos(v)*cos(4.*pi.*u);
y-yc=(1-u)*cos(v)*sin(4.*pi.*u);
z-zc=(1-u)*sin(v);
leading to,
which is a u-dependent constant independent of v. It is also stragithforward to show that the points corresponding to a fixed u but arbitrary v are all perpendicular to . They are therefore planar, equidistant from (xc,yc,yc), and therefore form a circle.
Therefore, the cross-sectional area is likewise a u-dependent quantity,
To get the volume, the task is then to integrate this area function along u with some kind of integration measure m(u).
The question though, is what is the functional form of m(u) that is appropriate here?
  댓글 수: 15
Matt J
Matt J 2022년 5월 2일
What we need is an even-more generalization of Pappus-Goddman-Godmann where the cross section could evolves along the curve
Seems pretty easy to get such a result for the case (like this one) where the region can be written in cylindrical coordinates. For the volume, we have
where for each fixed θ and z, the centroid with respect to r is while is the width of the region. In the special case (which holds here), where is independent of z, then, denoting the cross-sectional area as , the integration simplifies to
which is essentially the result that I was looking for in my original answer. So, this does allow for θ -varying cross sections, but admittedly requires that they have certain symmetry in r.
With this result though, it turns out to be very easy to calculate the volume of the horn. The horn has circular cross-sections of radius which are at a distance from the z-axis. We therefore have immediately that and leading to,

댓글을 달려면 로그인하십시오.

태그

Community Treasure Hunt

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

Start Hunting!

Translated by