[Date Index] [Thread Index] [Author Index]

Re: Balance point of a solid

• To: mathgroup at smc.vnet.net
• Subject: [mg113754] Re: Balance point of a solid
• From: Ray Koopman <koopman at sfu.ca>
• Date: Wed, 10 Nov 2010 06:30:20 -0500 (EST)
• References: <ib5u0i$c11$1@smc.vnet.net>

On Nov 7, 2:11 am, Daniel Lichtblau <d... at wolfram.com> wrote:
> ----- Original Message -----
> > From: "Andreas" <aa... at ix.netcom.com>
> > To: mathgr... at smc.vnet.net
> > Sent: Saturday, November 6, 2010 5:01:14 AM
> > Subject:  Re: Balance point of a solid
> > Ray,
>
> > I think that sounds right and follows the I found in the wikipedia
> > entry for Simplexhttp://en.wikipedia.org/wiki/Simplex. The entry has
> > a discussion of
> > "Cartesian coordinates for regular n-dimensionsl simplex in Rn" that
> > aplies to this.
>
> > To relate this to what I described earlier... for n = 4 dimensions, I
> > need a tetrahedron as my "base" it needs 4 points to describe it and
> > as you stated their coordinates must be mutually equidistant. In my
> > application the distance between them would always equal 1 for all of
> > these "bases". Their must be another word to describe this besides
> > base, but base seems intuitive to me.
>
> If I follow correctly, the main issue is setting up the constraints that =
define the region. If so, one might proceed as follows.
>
> First, we want to obtain a regular tetrahedron (not sure if that's the co=
rrect term. I mean a tetrahedron with all sides of same length, in this cas=
e 1). Also we'll orient so that all points have nonnegative coordinate valu=
es, and one is the origin, another is one unit along the x axis. From this =
point onward, it is not hard to go to higher dimensions. One sets up and so=
lves equations for which we increase dimension by one, and find points of u=
nit distance from all prior points. There are two solutions; we discard the=
 one that has a negative coordinate value.
>
> Next we will need to find the "upper" points (the ones that form the cap)=
. Here we again need to add one to the dimension of the base points (the te=
trahedron we just described).
>
> After this we need to find hyperplanes that define the region boundaries.=
 A bit of thought shows that these are formed by the base, the upper cap, a=
nd sets of n points (where n is the dimension of our space) comprised of n-=
1 base points and one upper cap point that lies above one of the base point=
s in the selection. Given the right number of points on the hyperplane, a n=
ull vector to certain differences gives a perpendicular vector, and basic l=
inear algebra gives a linear polynomial that defines said hyperplane.
>
> We also need to know on which side of each such boundary hyperplane our r=
egion lies. This can be determined by deciding on which side the barycenter=
 (average of the base and upper cap coordinates) lies. We use this to give =
an appropriate linear inequality for each boundary hyperplane.
>
> Here is the code I use for all this.
>
> normSquare[vec_] := vec.vec
>
> tetrahedrize[pts_] := Module[
>   {dim = Length[pts], newp, len, x, vars, polys, sol},
>   newp = Map[Append[#, 0] &, pts];
>   len = Norm[pts[[1]] - pts[[2]]];
>   vars = Array[x, dim];
>   polys = Map[normSquare[vars - #] &, newp] - len;
>   sol = RootReduce[vars /. Solve[polys == 0, vars]];
>   Append[newp, First[Select[sol, Last[#] > 0 &]]]
>   ]
>
> vertical[basept_, h_] :=
>  basept + UnitVector[Length[basept], Length[basept]]*h
>
> hyperplane[pts_,
>   vars_] := (vars - First[pts]).First[
>    NullSpace[Map[# - First[pts] &, Rest[pts]]]]
>
> halfSpace[pts_, vars_, pt_, assumed_] := Module[
>   {hp, ineq},
>   hp = hyperplane[pts, pt];
>   ineq = If[Reduce[hp > 0 && assumed] === False, LessEqual,
>     GreaterEqual];
>   ineq[hyperplane[pts, vars], 0]]
>
> Here we do the case of four dimensions. I'll keep it symbolic until the e=
nd. We'll assume, without loss of generality, that the heights are in incre=
asing order.
>
> In[100]:= fourtet = Nest[tetrahedrize, {{0}, {1}}, 2]
>
> Out[100]= {{0, 0, 0}, {1, 0, 0}, {1/2, Sqrt[3]/2, 0}, {1/2, 1/(
>   2 Sqrt[3]), Sqrt[2/3]}}
>
> In[142]:= base = Map[Append[#, 0] &, fourtet]
> len = Length[base];
> vars = Array[x, len]
> heights = Array[h, len]
>
> Out[142]= {{0, 0, 0, 0}, {1, 0, 0, 0}, {1/2, Sqrt[3]/2, 0, 0}, {1/2,
>   1/(2 Sqrt[3]), Sqrt[2/3], 0}}
>
> Out[144]= {x[1], x[2], x[3], x[4]}
>
> Out[145]= {h[1], h[2], h[3], h[4]}
>
> In[146]:= upperpts = Apply[vertical, Transpose[{base, heights}], {1}]
>
> Out[146]= {{0, 0, 0, h[1]}, {1, 0, 0, h[2]}, {1/2, Sqrt[3]/2, 0,
>   h[3]}, {1/2, 1/(2 Sqrt[3]), Sqrt[2/3], h[4]}}
>
> In[110]:= hyperplanesets =
>  Join[{base},
>   Table[Append[Delete[base, j], upperpts[[Mod[j + 1, len, 1]]]], {j,
>     len}], {upperpts}]
>
> Out[110]= {{{0, 0, 0, 0}, {1, 0, 0, 0}, {1/2, Sqrt[3]/2, 0, 0}, {1/2,
>    1/(2 Sqrt[3]), Sqrt[2/3], 0}}, {{1, 0, 0, 0}, {1/2, Sqrt[3]/2, 0,
>    0}, {1/2, 1/(2 Sqrt[3]), Sqrt[2/3], 0}, {1, 0, 0, h[2]}}, {{0, 0,
>    0, 0}, {1/2, Sqrt[3]/2, 0, 0}, {1/2, 1/(2 Sqrt[3]), Sqrt[2/3],
>    0}, {1/2, Sqrt[3]/2, 0, h[3]}}, {{0, 0, 0, 0}, {1, 0, 0, 0}, {1/2,
>    1/(2 Sqrt[3]), Sqrt[2/3], 0}, {1/2, 1/(2 Sqrt[3]), Sqrt[2/3],
>    h[4]}}, {{0, 0, 0, 0}, {1, 0, 0, 0}, {1/2, Sqrt[3]/2, 0, 0}, {0, 0,
>     0, h[1]}}, {{0, 0, 0, h[1]}, {1, 0, 0, h[2]}, {1/2, Sqrt[3]/2, 0,
>    h[3]}, {1/2, 1/(2 Sqrt[3]), Sqrt[2/3], h[4]}}}
>
> In[107]:= midpt = Mean[Join[base, upperpts]]
>
> Out[107]= {1/2, 1/8 (1/Sqrt[3] + Sqrt[3]), 1/(2 Sqrt[6]),
>  1/8 (h[1] + h[2] + h[3] + h[4])}
>
> In[147]:= assumed = LessEqual @@ Prepend[heights, 0]
>
> Out[147]= 0 <= h[1] <= h[2] <= h[3] <= h[4]
>
> In[141]:= hspaces =
>  Map[halfSpace[#, vars, midpt, assumed] &, hyperplanesets]
>
> Out[141]= {x[4] >= 0,
>  Sqrt[6] (-1 + x[1]) + Sqrt[2] x[2] + x[3] <=
>   0, -Sqrt[6] x[1] + Sqrt[2] x[2] + x[3] <=
>   0, -2 Sqrt[2] x[2] + x[3] <= 0,
>  x[3] >= 0, -h[1] + (h[1] - h[2]) x[1] + ((h[1] + h[2] - 2 h[3]) x[
>      2])/Sqrt[
>    3] + ((Sqrt[3] h[1] + Sqrt[3] h[2] + Sqrt[3] h[3] -
>       3 Sqrt[3] h[4]) x[3])/(3 Sqrt[2]) + x[4] <= 0}
>
> To find the mass, and center of mass coordinates, we'll need to integrate=
 over this region.
>
> In[149]:= maxes =
>  Append[Table[base[[j + 1, j]], {j, Length[base] - 1}], Last[heights]]
>
> Out[149]= {1, Sqrt[3]/2, Sqrt[2/3], h[4]}
>
> In[150]:= ranges = Thread[{vars, 0, maxes}]
>
> Out[150]= {{x[1], 0, 1}, {x[2], 0, Sqrt[3]/2}, {x[3], 0, Sqrt[2/
>   3]}, {x[4], 0, h[4]}}
>
> Here we need to make explicit values for the heights. I simply set h[j]=
=j for j from 1 to 4.
>
> In[154]:= reps = Thread[heights -> Range[len]]
>
> Out[154]= {h[1] -> 1, h[2] -> 2, h[3] -> 3, h[4] -> 4}
>
> We now compute the mass. This is simply the integral over the region, und=
er the assumption that our object in question has unit density.
>
> In[158]:= Timing[
>  mass = Integrate[Boole[And @@ hspaces] /. reps,
>    Sequence @@ (ranges /. reps)]]
>
> Out[158]= {7.114, 5/(12 Sqrt[2])}
>
> The center of mass coordinates are similarly obtained. For the jth coordi=
nate, one integrates x[j] over the region and divides by the mass.

cmass = With[{h = heights/.reps},
        Simplify[#/Tr@#&[h+Tr at h].Most/@base]]

{51/100, 53/(100*Sqrt[3]), (7*Sqrt[2/3])/25}

That's the same as the result given by

Integrate[Most@vars * Boole[And @@ hspaces] /. reps,
          Sequence @@ (ranges /. reps)] / mass

I've checked some other heights, and both methods always agree.
I have not checked higher dimensional problems,
because the integration takes too long.