       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[] - pts[]];
>   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:= fourtet = Nest[tetrahedrize, {{0}, {1}}, 2]
>
> Out= {{0, 0, 0}, {1, 0, 0}, {1/2, Sqrt/2, 0}, {1/2, 1/(
>   2 Sqrt), Sqrt[2/3]}}
>
> In:= base = Map[Append[#, 0] &, fourtet]
> len = Length[base];
> vars = Array[x, len]
> heights = Array[h, len]
>
> Out= {{0, 0, 0, 0}, {1, 0, 0, 0}, {1/2, Sqrt/2, 0, 0}, {1/2,
>   1/(2 Sqrt), Sqrt[2/3], 0}}
>
> Out= {x, x, x, x}
>
> Out= {h, h, h, h}
>
> In:= upperpts = Apply[vertical, Transpose[{base, heights}], {1}]
>
> Out= {{0, 0, 0, h}, {1, 0, 0, h}, {1/2, Sqrt/2, 0,
>   h}, {1/2, 1/(2 Sqrt), Sqrt[2/3], h}}
>
> In:= hyperplanesets =
>  Join[{base},
>   Table[Append[Delete[base, j], upperpts[[Mod[j + 1, len, 1]]]], {j,
>     len}], {upperpts}]
>
> Out= {{{0, 0, 0, 0}, {1, 0, 0, 0}, {1/2, Sqrt/2, 0, 0}, {1/2,
>    1/(2 Sqrt), Sqrt[2/3], 0}}, {{1, 0, 0, 0}, {1/2, Sqrt/2, 0,
>    0}, {1/2, 1/(2 Sqrt), Sqrt[2/3], 0}, {1, 0, 0, h}}, {{0, 0,
>    0, 0}, {1/2, Sqrt/2, 0, 0}, {1/2, 1/(2 Sqrt), Sqrt[2/3],
>    0}, {1/2, Sqrt/2, 0, h}}, {{0, 0, 0, 0}, {1, 0, 0, 0}, {1/2,
>    1/(2 Sqrt), Sqrt[2/3], 0}, {1/2, 1/(2 Sqrt), Sqrt[2/3],
>    h}}, {{0, 0, 0, 0}, {1, 0, 0, 0}, {1/2, Sqrt/2, 0, 0}, {0, 0,
>     0, h}}, {{0, 0, 0, h}, {1, 0, 0, h}, {1/2, Sqrt/2, 0,
>    h}, {1/2, 1/(2 Sqrt), Sqrt[2/3], h}}}
>
> In:= midpt = Mean[Join[base, upperpts]]
>
> Out= {1/2, 1/8 (1/Sqrt + Sqrt), 1/(2 Sqrt),
>  1/8 (h + h + h + h)}
>
> In:= assumed = LessEqual @@ Prepend[heights, 0]
>
> Out= 0 <= h <= h <= h <= h
>
> In:= hspaces =
>  Map[halfSpace[#, vars, midpt, assumed] &, hyperplanesets]
>
> Out= {x >= 0,
>  Sqrt (-1 + x) + Sqrt x + x <=
>   0, -Sqrt x + Sqrt x + x <=
>   0, -2 Sqrt x + x <= 0,
>  x >= 0, -h + (h - h) x + ((h + h - 2 h) x[
>      2])/Sqrt[
>    3] + ((Sqrt h + Sqrt h + Sqrt h -
>       3 Sqrt h) x)/(3 Sqrt) + x <= 0}
>
> To find the mass, and center of mass coordinates, we'll need to integrate=
over this region.
>
> In:= maxes =
>  Append[Table[base[[j + 1, j]], {j, Length[base] - 1}], Last[heights]]
>
> Out= {1, Sqrt/2, Sqrt[2/3], h}
>
> In:= ranges = Thread[{vars, 0, maxes}]
>
> Out= {{x, 0, 1}, {x, 0, Sqrt/2}, {x, 0, Sqrt[2/
>   3]}, {x, 0, h}}
>
> Here we need to make explicit values for the heights. I simply set h[j]=
=j for j from 1 to 4.
>
> In:= reps = Thread[heights -> Range[len]]
>
> Out= {h -> 1, h -> 2, h -> 3, h -> 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:= Timing[
>  mass = Integrate[Boole[And @@ hspaces] /. reps,
>    Sequence @@ (ranges /. reps)]]
>
> Out= {7.114, 5/(12 Sqrt)}
>
> 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), (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.

```

• Prev by Date: Re: command to save as .m file
• Next by Date: Re: issue with LinearSolve[] when using SparseArray when
• Previous by thread: Re: Balance point of a solid
• Next by thread: Re: Balance point of a solid