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.