Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
calculating the volume of a given shape

Subject: calculating the volume of a given shape

From: David Edson

Date: 15 Apr, 2009 14:17:02

Message: 1 of 20


Hello Dears,

I have a shape defined by the following formula?
-pi/2<theta<pi/2; 0<beta<2pi

x=rx.*cos(theta).*cos(beta);
y=ry.*sin(theta)+(2.^((-2*theta))+2.^((0.5*theta)));
z=rz.*cos(theta).*sin(beta);

how could I calculate the volume?

thanks

Subject: calculating the volume of a given shape

From: John D'Errico

Date: 15 Apr, 2009 14:49:01

Message: 2 of 20

"David Edson" <tojadeb@example.com> wrote in message <gs4q8u$kg1$1@fred.mathworks.com>...
>
> Hello Dears,
>
> I have a shape defined by the following formula?
> -pi/2<theta<pi/2; 0<beta<2pi
>
> x=rx.*cos(theta).*cos(beta);
> y=ry.*sin(theta)+(2.^((-2*theta))+2.^((0.5*theta)));
> z=rz.*cos(theta).*sin(beta);
>
> how could I calculate the volume?
>
> thanks

It depends. Is y correct as written? Or are you
missing a pair of parens? Which one is it?

 y=ry.*sin(theta)+(2.^((-2*theta))+2.^((0.5*theta)));
 y=ry.*(sin(theta)+(2.^((-2*theta))+2.^((0.5*theta))));

If it is the first case (as you have written it), then
how do you define the volume of a manifold that
intersects itself for some values of {rx,ry,rz}?

If it is the second case, a good guess is pi*rx*ry*rz.

John

Subject: calculating the volume of a given shape

From: Roger Stafford

Date: 16 Apr, 2009 06:04:09

Message: 3 of 20

"David Edson" <tojadeb@example.com> wrote in message <gs4q8u$kg1$1@fred.mathworks.com>...
> I have a shape defined by the following formula?
> -pi/2<theta<pi/2; 0<beta<2pi
>
> x=rx.*cos(theta).*cos(beta);
> y=ry.*sin(theta)+(2.^((-2*theta))+2.^((0.5*theta)));
> z=rz.*cos(theta).*sin(beta);
>
> how could I calculate the volume?

  There is a way to obtain a precise answer for your volume using just a single integration, David. Its validity depends, as John has pointed out, on your surface not intersecting itself and that depends on the value of ry. You can tell if that happens by plotting cos(theta) against y(theta) as theta varies from -pi/2 to +pi/2 to see if the curve crosses itself.

  Cross sections orthogonal to the y-axis yield ellipses of area pi*rx*rz*cos(theta)^2, so you just need to integrate this area along the y-axis. In some places the total cross section actually consists of an ellipse with a similar ellipse of smaller size cut out of its center. However, if you express your integral in terms of the parameter theta, this will automatically handle that situation (again provided the surface doesn't intersect itself.)

  If we abbreviate theta with simply t, you can express the volume as the definite integral of the function

 f(t) = pi*rx*rz*cos(t)^2 * dy/dt

taken with respect to t, from t equals -pi/2 to +pi/2.

  As it turns out, this integral can be evaluated in terms of elementary functions, so numerical methods are not really necessary. This integrand has three parts: the first part is a constant times cos(t)^3, and the second and third parts are both constants times cos(t)^2*exp(k*t) for two appropriate values of k. Each of these is of a form that can be evaluated using either integral tables or matlab's 'int' function in the Symbolic Toolbox.

  The above is the way I would tackle the problem, but I won't spoil things for you by carrying this any further (and anyway it's my bedtime.)

Roger Stafford

Subject: calculating the volume of a given shape

From: Roger Stafford

Date: 16 Apr, 2009 16:59:01

Message: 4 of 20

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <gs6hop$mfh$1@fred.mathworks.com>...
> There is a way to obtain a precise answer for your volume using just a single integration, David. Its validity depends, as John has pointed out, on your surface not intersecting itself and that depends on the value of ry. You can tell if that happens by plotting cos(theta) against y(theta) as theta varies from -pi/2 to +pi/2 to see if the curve crosses itself.
>
> Cross sections orthogonal to the y-axis yield ellipses of area pi*rx*rz*cos(theta)^2, so you just need to integrate this area along the y-axis. In some places the total cross section actually consists of an ellipse with a similar ellipse of smaller size cut out of its center. However, if you express your integral in terms of the parameter theta, this will automatically handle that situation (again provided the surface doesn't intersect itself.)
>
> If we abbreviate theta with simply t, you can express the volume as the definite integral of the function
>
> f(t) = pi*rx*rz*cos(t)^2 * dy/dt
>
> taken with respect to t, from t equals -pi/2 to +pi/2.
>
> As it turns out, this integral can be evaluated in terms of elementary functions, so numerical methods are not really necessary. This integrand has three parts: the first part is a constant times cos(t)^3, and the second and third parts are both constants times cos(t)^2*exp(k*t) for two appropriate values of k. Each of these is of a form that can be evaluated using either integral tables or matlab's 'int' function in the Symbolic Toolbox.
>
> The above is the way I would tackle the problem, but I won't spoil things for you by carrying this any further (and anyway it's my bedtime.)

  After a good night's sleep I decided to give you a little more help, David. The following identities are needed in solving your integration problem along the lines I outlined last night.

 1. 2^u = exp(log(2)*u)
 2. int('cos(t)^3','t',-pi/2,pi/2) = 4/3
 3. int('cos(t)^2*exp(k*t)','t',-pi/2,pi/2) =
      2/k/(k^2+4)*(exp(k*pi/2)-exp(-k*pi/2))

As I asserted earlier, these can be obtained either from a table of integrals or from matlab's symbolic toolbox.

  From these, you should be able to find a precise expression for the volume of your "shape" using the integral of f(t) defined earlier. Actually you may prefer to use the symbolic toolbox to evaluate the integral of f(t) directly, though the results may look somewhat messier.

Roger Stafford

Subject: calculating the volume of a given shape

From: David Edson

Date: 20 Apr, 2009 07:56:02

Message: 5 of 20


Thank you Roger,

I can calculate the integral using summation by taking small t intervals. I have one question though, What would be therCross sections orthogonal to the y-axis. The aim is still to calculate the volume.

If the shape is defined as follows. I introduced a sin(beta) in the y direction.


 x=rx.*cos(theta).*cos(beta);
 y=ry.*(sin(theta)+ (abs(sin(k1*beta))+ 2.^(k2.*theta)- 2.^(k3*theta)));
 z=rz.*cos(theta).*sin(beta);


Thanks in advance



"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <gs6hop$mfh$1@fred.mathworks.com>...
> "David Edson" <tojadeb@example.com> wrote in message <gs4q8u$kg1$1@fred.mathworks.com>...
> > I have a shape defined by the following formula?
> > -pi/2<theta<pi/2; 0<beta<2pi
> >
> > x=rx.*cos(theta).*cos(beta);
> > y=ry.*sin(theta)+(2.^((-2*theta))+2.^((0.5*theta)));
> > z=rz.*cos(theta).*sin(beta);
> >
> > how could I calculate the volume?
>
> There is a way to obtain a precise answer for your volume using just a single integration, David. Its validity depends, as John has pointed out, on your surface not intersecting itself and that depends on the value of ry. You can tell if that happens by plotting cos(theta) against y(theta) as theta varies from -pi/2 to +pi/2 to see if the curve crosses itself.
>
> Cross sections orthogonal to the y-axis yield ellipses of area pi*rx*rz*cos(theta)^2, so you just need to integrate this area along the y-axis. In some places the total cross section actually consists of an ellipse with a similar ellipse of smaller size cut out of its center. However, if you express your integral in terms of the parameter theta, this will automatically handle that situation (again provided the surface doesn't intersect itself.)
>
> If we abbreviate theta with simply t, you can express the volume as the definite integral of the function
>
> f(t) = pi*rx*rz*cos(t)^2 * dy/dt
>
> taken with respect to t, from t equals -pi/2 to +pi/2.
>
> As it turns out, this integral can be evaluated in terms of elementary functions, so numerical methods are not really necessary. This integrand has three parts: the first part is a constant times cos(t)^3, and the second and third parts are both constants times cos(t)^2*exp(k*t) for two appropriate values of k. Each of these is of a form that can be evaluated using either integral tables or matlab's 'int' function in the Symbolic Toolbox.
>
> The above is the way I would tackle the problem, but I won't spoil things for you by carrying this any further (and anyway it's my bedtime.)
>
> Roger Stafford

Subject: calculating the volume of a given shape

From: Roger Stafford

Date: 20 Apr, 2009 12:41:01

Message: 6 of 20

"David Edson" <tojadeb@example.com> wrote in message <gsh9qi$lmh$1@fred.mathworks.com>...
> I can calculate the integral using summation by taking small t intervals. I have one question though, What would be therCross sections orthogonal to the y-axis. The aim is still to calculate the volume.
>
> If the shape is defined as follows. I introduced a sin(beta) in the y direction.
>
> x=rx.*cos(theta).*cos(beta);
> y=ry.*(sin(theta)+ (abs(sin(k1*beta))+ 2.^(k2.*theta)- 2.^(k3*theta)));
> z=rz.*cos(theta).*sin(beta);

  The change you have made creates a much more difficult problem, David. The y cross section would no longer be an ellipse, since holding y fixed is now not a matter of merely holding theta fixed, but of varying beta and theta subject to an implicit relationship. This will cause x and z to vary in other than a simple ellipse. Taking y cross sections is definitely not the way to solve this altered problem!

  I'll outline a way you could solve it in terms of the parameters beta and theta. Call these parameters b and t respectively for simplification. Define the vectors P, Q, and R as:

 P(b,t) = [dx/db,dy/db,dz/db]
 Q(b,t) = [dx/dt,dy/dt,dz/dt]
 R(b,t) = [x,y,z]

where each of these quantities is expressed as functions of b and t. Then the volume enclosed by this surface can be expressed as the double integral with respect to b and t of the function

 f(b,t) = 1/3*dot(cross(P,Q),R)

as b varies from 0 to 2*pi and t from -pi/2 to +pi/2. As you can see, this integrand is a rather complicated function of b and t (beta and theta) and it is doubtful if an integral solution can be found in terms of elementary functions, so it would probably have to be evaluated numerically. Matlab's 'dblquad' function should be quite suitable for that task.

  You can test this formula using just the first term in y which defines an ellipsoidal surface. Its volume should turn out to be 4/3*pi*rx*ry*rz (within round-off tolerance of course.)

  Note that the above depends on the fact that x and z become zero when theta is at its two extremes -pi/2 and +pi/2, so the surface is thereby a closed surface. You will of course still have the same concerns as before about the surface intersecting itself and this time it may be more difficult to determine if that happens. You want to avoid situations where two different pairs of the parameters yield the same values for x, y, and z, (except of course at the endpoints of the parametric ranges where the surface closes on itself.)

Roger Stafford

Subject: calculating the volume of a given shape

From: Roger Stafford

Date: 20 Apr, 2009 17:55:03

Message: 7 of 20

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <gshqgt$clk$1@fred.mathworks.com>...
> ......
> I'll outline a way you could solve it in terms of the parameters beta and theta. Call these parameters b and t respectively for simplification. Define the vectors P, Q, and R as:
>
> P(b,t) = [dx/db,dy/db,dz/db]
> Q(b,t) = [dx/dt,dy/dt,dz/dt]
> R(b,t) = [x,y,z]
>
> where each of these quantities is expressed as functions of b and t. Then the volume enclosed by this surface can be expressed as the double integral with respect to b and t of the function
>
> f(b,t) = 1/3*dot(cross(P,Q),R)
>
> as b varies from 0 to 2*pi and t from -pi/2 to +pi/2.
> ......

  I have a couple of additional observations to make about that formula I gave you.

  First, the addition of the extra term introduces a strong discontinuity to your defined surface where the parameter theta approaches pi/2 and -pi/2. As beta sweeps through the range 0 to 2*pi, this term, abs(sin(k1*beta)), continues its full variation in spite of the fact that the variation is now right along the y-axis. It is not a natural-looking surface. It is a fortunate thing to be evaluating the volume in terms of the beta and theta parameters. Otherwise the numerical errors would probably be excessively large.

  The second is that the integral formula I gave you is a slightly disguised form of Gauss's Theorem, which should increase your confidence in it. The theorem states that the volume integral of the divergence of a field vector is equal to the surface integral of the dot product between that field vector at the surface and the normal vector to the surface. In this case the field vector would be 1/3*R whose divergence is a constant 1, and the surface normal elements are cross(P,Q)*db*dt, thus giving just the volume contained within the surface.

Roger Stafford

Subject: calculating the volume of a given shape

From: David Edson

Date: 21 Apr, 2009 07:29:02

Message: 8 of 20

Hello Dear Roger,

Thank you for your help. I am an applied scientist and some of the formulas you gave me were difficult to program in matlab. For example the cross products hold for only three element vectors.

C = CROSS(A,B) returns the cross product of the vectors
    A and B. That is, C = A x B. A and B must be 3 element
    vectors.
 
To follow the rule, I did the following:


[theta, beta]=meshgrid((-pi/2):(pi/100):(pi/2),0:(pi/50):2*pi);

t1=[theta(:,1);theta(:,2)]; b1=[beta(:,1);beta(:,2)];
  
  dxb=[ ]; dyb=[ ]; dzb=[ ];
 dxt=[ ]; dyt=[ ]; dzt=[ ];
 
 xx=[ ]; yy=[]; zz=[ ];
 
for k=1:length(t1)

t=t1(k); b=b1(k);

dxb1=-rx*cos(t)*sin(b);
dyb1=ry*sin(t)+(ry/10)*(k1*abs(cos(k1*b))+ 2.^(k2*t)- 2.^(k3*t));
dzb1=rz*cos(t)*cos(b);

dxt1=rz*-sin(t)*cos(b);
dyt1=ry*cos(t)+(ry/10)*(abs(sin(k1.*b)) +log(2)*k2*(2.^(k2*t))- log(2)*k3*(2^(k3.*t)));
dzt1=rz*-sin(t)*sin(b);

dxb=[dxb;dxb1];
dyb=[dyb;dyb1];
dzb=[dzb;dzb1];

dxt=[dxt;dxt1];
dyt=[dyt;dyt1];
dzt=[dzt;dzt1];

x1=rx*cos(t).*cos(b);
y1=ry*sin(t)+ (ry/10)*(abs(sin(k1*b))+ 2.^(k2*t)- 2.^(k3*t)); % including the concave shape of a tomato
z1=rz*cos(t).*sin(b);

xx=[xx;x1];
yy=[yy;y1];
zz=[zz;z1];

end


P=[dxb,dyb,dzb];
Q=[dxt,dyt,dzt];
R=[xx,yy,zz];
size(R)


cross(P,Q);

f=(1/3)*dot(cross(P,Q),R)


But, from f I can see that this is not the way.


I appreciate your smart mind; I am trying to model different shapes and it would be more interesting if you can tell me , or send me, articles or books that can help me define shapes, and calculate volumes and surface areas. Thank you for your help and time.

DE


"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <gsictn$flu$1@fred.mathworks.com>...
> "Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <gshqgt$clk$1@fred.mathworks.com>...
> > ......
> > I'll outline a way you could solve it in terms of the parameters beta and theta. Call these parameters b and t respectively for simplification. Define the vectors P, Q, and R as:
> >
> > P(b,t) = [dx/db,dy/db,dz/db]
> > Q(b,t) = [dx/dt,dy/dt,dz/dt]
> > R(b,t) = [x,y,z]
> >
> > where each of these quantities is expressed as functions of b and t. Then the volume enclosed by this surface can be expressed as the double integral with respect to b and t of the function
> >
> > f(b,t) = 1/3*dot(cross(P,Q),R)
> >
> > as b varies from 0 to 2*pi and t from -pi/2 to +pi/2.
> > ......
>
> I have a couple of additional observations to make about that formula I gave you.
>
> First, the addition of the extra term introduces a strong discontinuity to your defined surface where the parameter theta approaches pi/2 and -pi/2. As beta sweeps through the range 0 to 2*pi, this term, abs(sin(k1*beta)), continues its full variation in spite of the fact that the variation is now right along the y-axis. It is not a natural-looking surface. It is a fortunate thing to be evaluating the volume in terms of the beta and theta parameters. Otherwise the numerical errors would probably be excessively large.
>
> The second is that the integral formula I gave you is a slightly disguised form of Gauss's Theorem, which should increase your confidence in it. The theorem states that the volume integral of the divergence of a field vector is equal to the surface integral of the dot product between that field vector at the surface and the normal vector to the surface. In this case the field vector would be 1/3*R whose divergence is a constant 1, and the surface normal elements are cross(P,Q)*db*dt, thus giving just the volume contained within the surface.
>
> Roger Stafford

Subject: calculating the volume of a given shape

From: Roger Stafford

Date: 21 Apr, 2009 19:32:00

Message: 9 of 20

"David Edson" <tojadeb@example.com> wrote in message <gsjsju$fcg$1@fred.mathworks.com>...
> .......
> But, from f I can see that this is not the way.

  David, most of what you wrote looks valid, but I had in mind your using 'dblquad' which requires a function handle, rather than any use of a discrete mesh from 'meshgrid'. You could write it as a subfunction so as to be able to pass through the parameters, k1, k2, and k3. Based on your most recent shape, it might look something like this:

function f = tomato(b,t) % Named after one of your shape descriptions
% Must be able to handle vectors for b and scalars for t
% Uses parameters k1, k2, and k3
x = cos(t)*cos(b);
y = sin(t)+abs(sin(k1*b))+2^(k2*t)-2^(k3*t);
z = cos(t)*sin(b);
xb = -cos(t)*sin(b); % <-- dx/db
yb = sign(sin(k1*b)).*k1*cos(k1*b); % <-- dy/db
zb = cos(t)*cos(b); % <-- dz/db
xt = -sin(t)*cos(b); % <-- dx/dt
yt = cos(t)+log(2)*(k2*2^(k2*t)-k3*2^(k3*t)); % <-- dy/dt
zt = -sin(t)*sin(b); % <-- dz/dt
f = (yb.*zt-zb.*yt).*x+(zb.*xt-xb.*zt).*y+(xb.*yt-yb.*xt).*z;
return

Here xb, yb, and zb are the partial derivatives with respect to b of x, y, and z, respectively, and similarly with xt, yt, and zt as partial derivatives with respect to t. Note that the dot and cross products are written out explicitly in the formula for f here rather than calling on the 'dot' and 'cross' functions.

  Then outside this function you could call on it with:

V = rx*ry*rz/3*dblquad(@tomato,0,2*pi,-pi/2,pi/2,tol); % Volume

to obtain the shape's volume. Note that 'dblquad' requires that 'tomato' be able to handle vectors for its first argument, b, and scalars for the second one, t.

  Notice that I have factored out rx, ry, rz, and 1/3 here since it is senseless to force numerous computations with them inside 'tomato' when a single operation will suffice after performing 'dblquad'. This double integration operation does require a comparatively long time, especially if the tolerance is set for very accurate computation.

Roger Stafford

Subject: calculating the volume of a given shape

From: Roger Stafford

Date: 21 Apr, 2009 22:14:02

Message: 10 of 20

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <gsl6vg$hjo$1@fred.mathworks.com>...
> .....
> function f = tomato(b,t) % Named after one of your shape descriptions
> % Must be able to handle vectors for b and scalars for t
> % Uses parameters k1, k2, and k3
> x = cos(t)*cos(b);
> y = sin(t)+abs(sin(k1*b))+2^(k2*t)-2^(k3*t);
> z = cos(t)*sin(b);
> xb = -cos(t)*sin(b); % <-- dx/db
> yb = sign(sin(k1*b)).*k1*cos(k1*b); % <-- dy/db
> zb = cos(t)*cos(b); % <-- dz/db
> xt = -sin(t)*cos(b); % <-- dx/dt
> yt = cos(t)+log(2)*(k2*2^(k2*t)-k3*2^(k3*t)); % <-- dy/dt
> zt = -sin(t)*sin(b); % <-- dz/dt
> f = (yb.*zt-zb.*yt).*x+(zb.*xt-xb.*zt).*y+(xb.*yt-yb.*xt).*z;
> return
> .....

  As you may recall, I mentioned earlier that the newly added term,

 abs(sin(k1*b))

in your surface creates a discontinuity in that surface. One consequence is a jump discontinuity in the partial derivative, dy/db, (called yb in the code.) It occurs to me that I may have written the 'tomato' function with its arguments in the less favorable order for double integration. I have reason to believe that 'dblquad' treats the first argument it sees (the one it sends out as a vector) as the first integration variable in an iterated integration process. This means that for every scalar value of t used, the integration with respect to b will (probably) encounter a discontinuity and thereby be slowed down in adjusting to it because of the need for a dense concentration of b values in the vicinity of the discontinuity.

  If you encounter such a difficulty, you should rewrite 'tomato' with its b and t arguments in the opposite order. The only consequence that I can see is that the two expressions, 2^(k2*t) and 2^(k3*t) would then need a dot before the '^' operator, giving 2.^(k2*t) and 2.^(k3*t), since t would then be a vector. The function 'dblquad' will still encounter a discontinuity when it crosses the critical value of b, now in the second integration variable, but this would presumably occur only once (or at worst, a few times, depending on k1) thereby (hopefully) causing less delay all told.

  In your endeavor to "model different shapes" the changes you would need for differing shapes would only be to the x, y, and z functions and their six derivatives. The quantity, f, would presumably remain the same since it is directly based on Gauss's divergence integral theorem. (I'm ignoring the ability to factor out rx, ry, and rz in making that statement. With some changes you might have to put them back in the x, y, z expressions.)

  As for a recommendation of books or articles, this use of Gauss's divergence integral theorem can be found in a great many textbooks in advanced calculus as well as in many engineering and physics texts on electric and magnetic field theory. It is usually accompanied by a description of Stokes Theorem. You will not find it in tables of volumes for common objects, since your objects seem to be distinctly uncommon. I think you should resign yourself to working out such problems the hard way.

  (At this point I am reminded of a story told about Thomas Edison who needed to find the volume of the light bulbs he was getting ready to sell. When told by his advisors that it was a difficult problem mathematically, he promptly removed one of the bulbs' base and filled it with water whose volume he could then measure easily.)

  You mentioned calculating surface areas, too. You could obtain those using 'dblquad' parametrically by using as an integrand the function

 f=sqrt((yb.*zt-zb.*yt).^2+(zb.*xt-xb.*zt).^2+(xb.*yt-yb.*xt).^2)

which is just the length of the cross product vector, P x Q. It does not involve the R field. You would still be integrating with respect to the b and t parameters here.

Roger Stafford

Subject: calculating the volume of a given shape

From: David Edson

Date: 22 Apr, 2009 16:03:03

Message: 11 of 20


Dear Roger,

Thank you indeed.

The Volume I calculated, though, is most of the time a negative value( I changed k1, k2, k3, to see if it works). It looks that the Volume is not affected by the change in the value of k1, which I could not accept. I am trying to work it out.

Thanks you Roger... You are really a smart guy
















"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <gslgfa$hrf$1@fred.mathworks.com>...
> "Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <gsl6vg$hjo$1@fred.mathworks.com>...
> > .....
> > function f = tomato(b,t) % Named after one of your shape descriptions
> > % Must be able to handle vectors for b and scalars for t
> > % Uses parameters k1, k2, and k3
> > x = cos(t)*cos(b);
> > y = sin(t)+abs(sin(k1*b))+2^(k2*t)-2^(k3*t);
> > z = cos(t)*sin(b);
> > xb = -cos(t)*sin(b); % <-- dx/db
> > yb = sign(sin(k1*b)).*k1*cos(k1*b); % <-- dy/db
> > zb = cos(t)*cos(b); % <-- dz/db
> > xt = -sin(t)*cos(b); % <-- dx/dt
> > yt = cos(t)+log(2)*(k2*2^(k2*t)-k3*2^(k3*t)); % <-- dy/dt
> > zt = -sin(t)*sin(b); % <-- dz/dt
> > f = (yb.*zt-zb.*yt).*x+(zb.*xt-xb.*zt).*y+(xb.*yt-yb.*xt).*z;
> > return
> > .....
>
> As you may recall, I mentioned earlier that the newly added term,
>
> abs(sin(k1*b))
>
> in your surface creates a discontinuity in that surface. One consequence is a jump discontinuity in the partial derivative, dy/db, (called yb in the code.) It occurs to me that I may have written the 'tomato' function with its arguments in the less favorable order for double integration. I have reason to believe that 'dblquad' treats the first argument it sees (the one it sends out as a vector) as the first integration variable in an iterated integration process. This means that for every scalar value of t used, the integration with respect to b will (probably) encounter a discontinuity and thereby be slowed down in adjusting to it because of the need for a dense concentration of b values in the vicinity of the discontinuity.
>
> If you encounter such a difficulty, you should rewrite 'tomato' with its b and t arguments in the opposite order. The only consequence that I can see is that the two expressions, 2^(k2*t) and 2^(k3*t) would then need a dot before the '^' operator, giving 2.^(k2*t) and 2.^(k3*t), since t would then be a vector. The function 'dblquad' will still encounter a discontinuity when it crosses the critical value of b, now in the second integration variable, but this would presumably occur only once (or at worst, a few times, depending on k1) thereby (hopefully) causing less delay all told.
>
> In your endeavor to "model different shapes" the changes you would need for differing shapes would only be to the x, y, and z functions and their six derivatives. The quantity, f, would presumably remain the same since it is directly based on Gauss's divergence integral theorem. (I'm ignoring the ability to factor out rx, ry, and rz in making that statement. With some changes you might have to put them back in the x, y, z expressions.)
>
> As for a recommendation of books or articles, this use of Gauss's divergence integral theorem can be found in a great many textbooks in advanced calculus as well as in many engineering and physics texts on electric and magnetic field theory. It is usually accompanied by a description of Stokes Theorem. You will not find it in tables of volumes for common objects, since your objects seem to be distinctly uncommon. I think you should resign yourself to working out such problems the hard way.
>
> (At this point I am reminded of a story told about Thomas Edison who needed to find the volume of the light bulbs he was getting ready to sell. When told by his advisors that it was a difficult problem mathematically, he promptly removed one of the bulbs' base and filled it with water whose volume he could then measure easily.)
>
> You mentioned calculating surface areas, too. You could obtain those using 'dblquad' parametrically by using as an integrand the function
>
> f=sqrt((yb.*zt-zb.*yt).^2+(zb.*xt-xb.*zt).^2+(xb.*yt-yb.*xt).^2)
>
> which is just the length of the cross product vector, P x Q. It does not involve the R field. You would still be integrating with respect to the b and t parameters here.
>
> Roger Stafford

Subject: calculating the volume of a given shape

From: Roger Stafford

Date: 22 Apr, 2009 20:40:18

Message: 12 of 20

"David Edson" <tojadeb@example.com> wrote in message <gsnf3n$fmv$1@fred.mathworks.com>...
> The Volume I calculated, though, is most of the time a negative value( I changed k1, k2, k3, to see if it works). It looks that the Volume is not affected by the change in the value of k1, which I could not accept. I am trying to work it out.

  I have oriented the cross product P x Q in the formula I gave you so that ordinarily it points outward away from the origin as in the case of a centered ellipsoid. Some of your "tomato" shapes can legitimately cause it to point with an inward component relative to the origin but that is needed to subtract out regions that may be unoccupied between the origin and the shape. However, to get on overall negative value to the volume indicates to me that your shape has seriously intersected itself or is entirely turned inside out with respect to the beta and theta parameters. The part of the volume region in which inside/outside become reversed will count as negative volumes. You don't want that!

  Have you ever seen pictures of what mathematicians call a "klein bottle" where the neck of a vase is turned down where it reenters the side of the bottle and connects inside with its base. There is ambiguity in such a surface as to which side is outside and which is inside. With the brief studies I made of your earlier shapes I (and John D'Errico too) could see that for some values of the parameters you were using, your surfaces intersected themselves in a manner analogous to a klein bottle.

  As your shape was before the abs(sin(k1*beta)) term was added, I made plots of cos(t) versus y when ry was not included in the second two terms of y. When the ry value was sufficiently small this curve did indeed intersect itself and showed that you would then have an anomalous surface in which some of its volume receives a negative value. I suspect the same thing is happening with your later version. To find it with y now depending on beta will be a more difficult task, but this is certainly what you need to do. Remember I warned you to avoid any case where you could get the same x, y, z value for more than one pair of beta and theta values. This is how it happens. It is important to avoid such intersections if you are to obtain valid values for volumes. It is NOT enough to check that your overall volume is positive.

  As a test for the general validity of the method, you should remove the extra terms on y and use only

 y = ry*sin(t)

which makes the shape a pure centered ellipsoid, to see if the resulting volume is the positive value 4/3*pi*rx*ry*rz. If it comes out negative, then there would be something wrong in the formula.

Roger Stafford

Subject: calculating the volume of a given shape

From: David Edson

Date: 23 Apr, 2009 09:58:02

Message: 13 of 20


Dear Roger,

with only y=sin(t), I get a negative value. The magnetude of the area is ok. I think the problem is with the cross product. Cross(P,Q) is different from cross(Q,P).

thanks,

DE










"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <gsnvbi$jl$1@fred.mathworks.com>...
> "David Edson" <tojadeb@example.com> wrote in message <gsnf3n$fmv$1@fred.mathworks.com>...
> > The Volume I calculated, though, is most of the time a negative value( I changed k1, k2, k3, to see if it works). It looks that the Volume is not affected by the change in the value of k1, which I could not accept. I am trying to work it out.
>
> I have oriented the cross product P x Q in the formula I gave you so that ordinarily it points outward away from the origin as in the case of a centered ellipsoid. Some of your "tomato" shapes can legitimately cause it to point with an inward component relative to the origin but that is needed to subtract out regions that may be unoccupied between the origin and the shape. However, to get on overall negative value to the volume indicates to me that your shape has seriously intersected itself or is entirely turned inside out with respect to the beta and theta parameters. The part of the volume region in which inside/outside become reversed will count as negative volumes. You don't want that!
>
> Have you ever seen pictures of what mathematicians call a "klein bottle" where the neck of a vase is turned down where it reenters the side of the bottle and connects inside with its base. There is ambiguity in such a surface as to which side is outside and which is inside. With the brief studies I made of your earlier shapes I (and John D'Errico too) could see that for some values of the parameters you were using, your surfaces intersected themselves in a manner analogous to a klein bottle.
>
> As your shape was before the abs(sin(k1*beta)) term was added, I made plots of cos(t) versus y when ry was not included in the second two terms of y. When the ry value was sufficiently small this curve did indeed intersect itself and showed that you would then have an anomalous surface in which some of its volume receives a negative value. I suspect the same thing is happening with your later version. To find it with y now depending on beta will be a more difficult task, but this is certainly what you need to do. Remember I warned you to avoid any case where you could get the same x, y, z value for more than one pair of beta and theta values. This is how it happens. It is important to avoid such intersections if you are to obtain valid values for volumes. It is NOT enough to check that your overall volume is positive.
>
> As a test for the general validity of the method, you should remove the extra terms on y and use only
>
> y = ry*sin(t)
>
> which makes the shape a pure centered ellipsoid, to see if the resulting volume is the positive value 4/3*pi*rx*ry*rz. If it comes out negative, then there would be something wrong in the formula.
>
> Roger Stafford

Subject: calculating the volume of a given shape

From: Roger Stafford

Date: 24 Apr, 2009 15:09:01

Message: 14 of 20

"David Edson" <tojadeb@example.com> wrote in message <gspe3a$3dg$1@fred.mathworks.com>...
> ........
> with only y=sin(t), I get a negative value. The magnetude of the area is ok. I think the problem is with the cross product. Cross(P,Q) is different from cross(Q,P).
> ........

  Yes, you are right, David. What threw me off is that with your definition of the angle beta, it represents a left-hand rather than right-hand rotation about the y-axis. That is, as viewed from the positive y-axis, beta advances in a clockwise sense in the x-z plane, which is opposite to the way azimuthal-type angles are usually defined in spherical coordinate systems. So, yes, for your coordinate definitions you should be using cross(Q,P) rather than cross(P,Q) in computing volumes and surface areas, (where P = dR/db and Q = dR/dt.)

  It still remains true that you should guard against having your parametrically defined surfaces intersect themselves. If that were to happen, even though you may obtain positive values for the total computed volume this would be in error because some volume regions are given a negative value and reduce rather than increase the total volume.

Roger Stafford

Subject: calculating the volume of a given shape

From: David Edson

Date: 27 Apr, 2009 14:06:02

Message: 15 of 20


Roger,

I want your idea again. I want to have my geometry zero y values at t values where abs(sin(k1t)) equals one. how could I approach it?

thanks,

DE








"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <gsskmd$ghg$1@fred.mathworks.com>...
> "David Edson" <tojadeb@example.com> wrote in message <gspe3a$3dg$1@fred.mathworks.com>...
> > ........
> > with only y=sin(t), I get a negative value. The magnetude of the area is ok. I think the problem is with the cross product. Cross(P,Q) is different from cross(Q,P).
> > ........
>
> Yes, you are right, David. What threw me off is that with your definition of the angle beta, it represents a left-hand rather than right-hand rotation about the y-axis. That is, as viewed from the positive y-axis, beta advances in a clockwise sense in the x-z plane, which is opposite to the way azimuthal-type angles are usually defined in spherical coordinate systems. So, yes, for your coordinate definitions you should be using cross(Q,P) rather than cross(P,Q) in computing volumes and surface areas, (where P = dR/db and Q = dR/dt.)
>
> It still remains true that you should guard against having your parametrically defined surfaces intersect themselves. If that were to happen, even though you may obtain positive values for the total computed volume this would be in error because some volume regions are given a negative value and reduce rather than increase the total volume.
>
> Roger Stafford

Subject: calculating the volume of a given shape

From: Roger Stafford

Date: 27 Apr, 2009 19:40:19

Message: 16 of 20

"David Edson" <tojadeb@example.com> wrote in message <gt4e4a$3vo$1@fred.mathworks.com>...
> I want your idea again. I want to have my geometry zero y values at t values where abs(sin(k1t)) equals one. how could I approach it?

  I'm afraid I don't understand your question, David. What exactly do you mean by "I want to have my geometry zero y values at t values where abs(sin(k1t)) equals one"? I believe your most recent equation for y was

 y = sin(t)+abs(sin(k1*b))+2^(k2*t)-2^(k3*t)

where b = beta varies from 0 to 2*pi. How did abs(sin(k1t)) get into the act? What significance is there to having abs(sin(k1t)) equal to one?

  If this is somehow related to the question of whether your surface intersects itself, for that you need to make studies of radial cross sections of planes containing the y-axis for various values of beta. In each section, as theta varies from -pi/2 to +pi/2, you would be plotting the radial distance cos(theta) (ignoring rx and rz) against y. Earlier on before y depended on beta there were values of ry where this did occur, as I recall.

Roger Stafford

Subject: calculating the volume of a given shape

From: Roger Stafford

Date: 28 Apr, 2009 00:42:01

Message: 17 of 20

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <gt51n3$3cc$1@fred.mathworks.com>...
> "David Edson" <tojadeb@example.com> wrote in message <gt4e4a$3vo$1@fred.mathworks.com>...
> > I want your idea again. I want to have my geometry zero y values at t values where abs(sin(k1t)) equals one. how could I approach it?
>
> I'm afraid I don't understand your question, David. ......

  Your mention of the term abs(sin(k1t)) has triggered an important thought, David. I should have thought of it before this. Presumably this is a term you wish to add to the y coordinate of your surface. I make the claim that this addition will have no effect whatsoever on the total volume enclosed by the surface, nor will it have a bearing on the question of whether the surface will intersect itself or not. Also the same holds for the additive term abs(sin(k1*b)). (Of course either of these would have a profound effect on the shape itself.)

  See if you agree with my reasoning. In the case of abs(sin(k1t)), for any given beta its value will be the same at one value of t = theta as at its negative. Moreover the values of x and z are the same for the two opposite-signed values of theta since cos(theta) is the same in both cases. What this means is that y is being moved upwards by the same amount at the two ends of a rod extending through the volume parallel to the y-axis. Since the volume can be considered as the double integral over the x-z plane of the length of that rod throughout the volume region, simply moving the whole rod farther up by variable amounts without changing its length will not affect the total volume occupied. The same would be all the more true for the additive term abs(sin(k1*b)) since for any given b = beta, it is a fixed value unchanging with theta.

  In looking for surface intersections, as I have stated, you consider the y value plotting against the value cos(theta) for any fixed beta and look for cases of the curve crossing itself. That would occur with two different values of theta with the same y and the same cos(theta) These thetas can only be negatives of one another. With each of the two above additions, abs(sin(k1t)) and abs(sin(k1*b)), the curve will therefore have such intersections at points where both ends of the "rod" above (zero length in this case) are moved upwards by the same amount. Hence they will still intersect in the same way.

  I recall you mentioning that the volume calculation was independent of k1. This is the explanation for that phenomenon! It is also independent of the two entire additive terms. Put them both in, multiply them by any desired constant factor, and they will make no difference to either the volume or the property of self intersection.

  This fact should greatly simplify any study of such intersections for the added term abs(sin(k1*b)). It should also make it again possible to find volumes using the analytic solution for the volume with the methods I mentioned early on in this thread and not require the use of numerical integration with either 'dblquad' or 'quad2d'. You are back to the equivalent (in the volume sense) equations:

 x=rx.*cos(theta).*cos(beta);
 y=ry.*sin(theta)+2.^(k2*t)-2.^(k3*t));
 z=rz.*cos(theta).*sin(beta);

which have an exact analytic expression for their volume.

Roger Stafford

Subject: calculating the volume of a given shape

From: David Edson

Date: 28 Apr, 2009 06:42:05

Message: 18 of 20

Hello Roge,

Let me rephrase my question. I want to have a geometry defined as follows:


x=rx.*cos(theta).*cos(beta);
y=ry.*(sin(theta)+ (abs(sin(k1.*beta))+ 2.^(k2.*theta)- 2.^(k3.*theta)));
z=rz.*cos(theta).*sin(beta);
 
But, I want also to have x and z to be zero and the geometry to have a y value different from zero when abs(sin(k1t)) is equal to 1.

for example when k1=2, I want to have a cross with curved tips when seen from the top.

I do not know if I make it clearer?






> > I'm afraid I don't understand your question, David. ......
>
> Your mention of the term abs(sin(k1t)) has triggered an important thought, David. I should have thought of it before this. Presumably this is a term you wish to add to the y coordinate of your surface. I make the claim that this addition will have no effect whatsoever on the total volume enclosed by the surface, nor will it have a bearing on the question of whether the surface will intersect itself or not. Also the same holds for the additive term abs(sin(k1*b)). (Of course either of these would have a profound effect on the shape itself.)
>
> See if you agree with my reasoning. In the case of abs(sin(k1t)), for any given beta its value will be the same at one value of t = theta as at its negative. Moreover the values of x and z are the same for the two opposite-signed values of theta since cos(theta) is the same in both cases. What this means is that y is being moved upwards by the same amount at the two ends of a rod extending through the volume parallel to the y-axis. Since the volume can be considered as the double integral over the x-z plane of the length of that rod throughout the volume region, simply moving the whole rod farther up by variable amounts without changing its length will not affect the total volume occupied. The same would be all the more true for the additive term abs(sin(k1*b)) since for any given b = beta, it is a fixed value unchanging with theta.
>
> In looking for surface intersections, as I have stated, you consider the y value plotting against the value cos(theta) for any fixed beta and look for cases of the curve crossing itself. That would occur with two different values of theta with the same y and the same cos(theta) These thetas can only be negatives of one another. With each of the two above additions, abs(sin(k1t)) and abs(sin(k1*b)), the curve will therefore have such intersections at points where both ends of the "rod" above (zero length in this case) are moved upwards by the same amount. Hence they will still intersect in the same way.
>
> I recall you mentioning that the volume calculation was independent of k1. This is the explanation for that phenomenon! It is also independent of the two entire additive terms. Put them both in, multiply them by any desired constant factor, and they will make no difference to either the volume or the property of self intersection.
>
> This fact should greatly simplify any study of such intersections for the added term abs(sin(k1*b)). It should also make it again possible to find volumes using the analytic solution for the volume with the methods I mentioned early on in this thread and not require the use of numerical integration with either 'dblquad' or 'quad2d'. You are back to the equivalent (in the volume sense) equations:
>
> x=rx.*cos(theta).*cos(beta);
> y=ry.*sin(theta)+2.^(k2*t)-2.^(k3*t));
> z=rz.*cos(theta).*sin(beta);
>
> which have an exact analytic expression for their volume.
>
> Roger Stafford

Subject: calculating the volume of a given shape

From: David Edson

Date: 29 Apr, 2009 13:53:01

Message: 19 of 20

Hello Roger,

The two dimensional version of what I want in three dimension is the following:

e=0:20/pi:2*pi;
r=0;
x=r*cos(e); y=r*sin(e);

x(t==3*pi/4)=0;
x(t==3*pi/4)=0;
 x(t==5*pi/4)=0;
 x(t==7*pi/4)=0;

x(t==3*pi/4)=0;
 y(t==3*pi/4)=0;
 y(t==5*pi/4)=0;
 y(t==7*pi/4)=0;

plot(x,y)

Now, I want to have such a shape in 3D. I thought that would be achieved by making x,and z values zero at abs(sin(k1.beta))=1. The shape was defined as follows:

x=rx.*cos(theta).*cos(beta);
y=ry.*(sin(theta)+ (abs(sin(k1.*beta))+ 2.^(k2.*theta)- 2.^(k3.*theta)));
z=rz.*cos(theta).*sin(beta);

is it solvable?

thanks

Subject: calculating the volume of a given shape

From: David Edson

Date: 7 Jun, 2009 14:59:02

Message: 20 of 20

Roger,

please give your hand how I could get a 2D contour of half of my shape I defined previously. I set rz=0, but the contour I got was not a simplified contour I wnted.

thanks,

TE




"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <gs6hop$mfh$1@fred.mathworks.com>...
> "David Edson" <tojadeb@example.com> wrote in message <gs4q8u$kg1$1@fred.mathworks.com>...
> > I have a shape defined by the following formula?
> > -pi/2<theta<pi/2; 0<beta<2pi
> >
> > x=rx.*cos(theta).*cos(beta);
> > y=ry.*sin(theta)+(2.^((-2*theta))+2.^((0.5*theta)));
> > z=rz.*cos(theta).*sin(beta);
> >
> > how could I calculate the volume?
>
> There is a way to obtain a precise answer for your volume using just a single integration, David. Its validity depends, as John has pointed out, on your surface not intersecting itself and that depends on the value of ry. You can tell if that happens by plotting cos(theta) against y(theta) as theta varies from -pi/2 to +pi/2 to see if the curve crosses itself.
>
> Cross sections orthogonal to the y-axis yield ellipses of area pi*rx*rz*cos(theta)^2, so you just need to integrate this area along the y-axis. In some places the total cross section actually consists of an ellipse with a similar ellipse of smaller size cut out of its center. However, if you express your integral in terms of the parameter theta, this will automatically handle that situation (again provided the surface doesn't intersect itself.)
>
> If we abbreviate theta with simply t, you can express the volume as the definite integral of the function
>
> f(t) = pi*rx*rz*cos(t)^2 * dy/dt
>
> taken with respect to t, from t equals -pi/2 to +pi/2.
>
> As it turns out, this integral can be evaluated in terms of elementary functions, so numerical methods are not really necessary. This integrand has three parts: the first part is a constant times cos(t)^3, and the second and third parts are both constants times cos(t)^2*exp(k*t) for two appropriate values of k. Each of these is of a form that can be evaluated using either integral tables or matlab's 'int' function in the Symbolic Toolbox.
>
> The above is the way I would tackle the problem, but I won't spoil things for you by carrying this any further (and anyway it's my bedtime.)
>
> Roger Stafford

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us