Re: Balance point of a solid
- To: mathgroup at smc.vnet.net
- Subject: [mg113673] Re: Balance point of a solid
- From: Daniel Lichtblau <danl at wolfram.com>
- Date: Sun, 7 Nov 2010 05:11:02 -0500 (EST)
----- Original Message ----- > From: "Andreas" <aagas at ix.netcom.com> > To: mathgroup at smc.vnet.net > Sent: Saturday, November 6, 2010 5:01:14 AM > Subject: [mg113660] Re: Balance point of a solid > Ray, > > I think that sounds right and follows the I found in the wikipedia > entry for Simplex http://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 correct term. I mean a tetrahedron with all sides of same length, in this case 1). Also we'll orient so that all points have nonnegative coordinate values, 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 solves equations for which we increase dimension by one, and find points of unit 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 tetrahedron 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, and 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 points in the selection. Given the right number of points on the hyperplane, a null vector to certain differences gives a perpendicular vector, and basic linear algebra gives a linear polynomial that defines said hyperplane. We also need to know on which side of each such boundary hyperplane our region 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 end. We'll assume, without loss of generality, that the heights are in increasing 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, under 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 coordinate, one integrates x[j] over the region and divides by the mass. Daniel Lichtblau Wolfram Research