Re: Drawing polytopes

*To*: mathgroup at smc.vnet.net*Subject*: [mg24267] Re: [mg24236] Drawing polytopes*From*: Andrzej Kozlowski <andrzej at tuins.ac.jp>*Date*: Wed, 5 Jul 2000 23:10:48 -0400 (EDT)*Sender*: owner-wri-mathgroup at wolfram.com

on 7/4/00 9:39 AM, Matthijs Sypkens Smit at matthijs at helena.tux.nu wrote: > I'm wondering if it's possible to draw (convex) polytopes defined by > systems of linear inequalities. For example the system: > > x + 2y + z <= 3 > 2x - y >= 0 > -x + 2y + z >= 0 > y + z >= 0 > 2x + 3y + 3z >= 1 > > Can I get Mathematica to draw the polytope defined on R^3 by this system? > > I'm able to solve this system with > InequalitySolve[{x+2y+2z<=3,2x-y>=0,-x+2y+z>=0,y+z>=0,2x+3y+3z>=1},{x,y,z}] > but the data is hard to interpret this way. > > btw; the same problem, but a little more intuitive, would be to draw a > cube with the system defined below: > > {x >= -1, x <= 1, y >= -1, y <= 1, z >= -1, z <= 1} > > Any hints are appreciated, > There is no built in way to do this. I can think of two reasonable approaches, in one the code is simple but it is computation intensive and unlikely to produce an enlightening result, and the second in which the code is tedious to write out but which gives accurate results and does it fast. Let m estart with the second, (which is actually pretty obvious). Let me start with the second. First we load In[1]:= << Algebra`InequalitySolve` We define the inequalitis In[2]:= ineqs = {x + 2y + z <= 3, 2x - y >= 0, -x + 2y + z >= 0, y + z >= 0, 2x + 3y + 3z >= 1}; and solve them In[3]:= ineqs2 = LogicalExpand[InequalitySolve[ineqs, {x, y, z}]] Out[3]= 3 3 3 x == - && z == - - 2 y && y <= - || 2 2 2 3 y == 3 - x && z == 3 - x - 2 y && x < - && 1 <= x || 2 x <= -1 && x - 2 y <= z && y <= 2 x && z <= 3 - x - 2 y || 1 1 -1 < x && x < - && x - 2 y <= z && y <= - (-1 + 5 x) && 2 3 1 z <= 3 - x - 2 y || x < 1 && y < x && - <= x && 2 x - 2 y <= z && z <= 3 - x - 2 y || 3 x < - && y < x && 1 <= x && x - 2 y <= z && 2 1 z <= 3 - x - 2 y || -1 < x && x < - && 2 1 1 - (-1 + 5 x) < y && - (1 - 2 x - 3 y) <= z && y <= 2 x && 3 3 1 z <= 3 - x - 2 y || x < 1 && - <= x && x <= y && 2 -y <= z && y <= 2 x && z <= 3 - x - 2 y || 3 x < - && y < 3 - x && 1 <= x && x <= y && -y <= z && 2 z <= 3 - x - 2 y (Applying LogicalExpand makes the output clearer). Next the tedious part. First we shall define some global variables, that we shall use in the plot: In[4]:= Minx = -10; Miny = -10; Minz = -10; rx = 0.1; ry = 0.1; Next, we take each inequality and simply create the corresponding plot (using our parameters Minx,Miny and Minz whenever the lower bounds are not determined by the inequalities). We notice that the first two inequalities determine straight lines which can be ploted using ParametricPlot. The other inequalities determine 2 dimensional linear subspaces. In[7]:= plot[1] = ParametricPlot3D[{3/2, y, 3/2 - 2y}, {y, Miny, 3/2}, DisplayFunction -> Identity]; In[8]:= plot[2] = ParametricPlot3D[{x, 3 - x, 3 - x - 2}, {x, 1, 3/2}, DisplayFunction -> Identity]; In[9]:= plot[3] = Table[Line[{{x, y, x - 2y}, {x, y, 3 - x - 2y}}], {x, Minx, -1, rx}, {y, Miny, 2x, ry}]; In[10]:= plot[4] = Table[Line[{{x, y, x - 2y}, {x, y, 3 - x - 2y}}], {x, -1, 0.5, rx}, {y, Miny, (5x - 1)/3, ry}]; In[11]:= plot[5] = Table[Line[{{x, y, x - 2y}, {x, y, 3 - x - 2y}}], {x, 0.5, 1, rx}, {y, Miny, x, ry}]; In[12]:= plot[6] = Table[Line[{{x, y, x - 2y}, {x, y, 3 - x - 2y}}], {x, 1, 3/2, rx}, {y, Miny, x, ry}]; In[13]:= plot[7] = Table[Line[{{x, y, (1 - 2 x - 3 y)/3}, {x, y, 3 - x - 2 y}}], {x, -1, 0.5, rx}, {y, (5 x - 1)/3, 2 x, ry}]; In[14]:= plot[8] = Table[Line[{{x, y, -y}, {x, y, 3 - x - 2 y}}], {x, 0.5, 1, rx}, {y, x, 2x, ry}]; In[15]:= plot[9] = Table[Line[{{x, y, -y}, {x, y, 3 - x - 2 y}}], {x, 1, 3/2, rx}, {y, x, 3 - x, ry}]; The polytope (which is really of course unbounded) can now be seen with: In[16]:= polytope = Show[Graphics3D[Table[{RGBColor[i/10, 0, 0], plot[i]}, {i, 3, 9}]]] (It is easier to see it properly if one loads in real-time 3d grpahics using <<RealTime3D`). The one dimensional pieces can be obtained in the same way: In[18]:= lines = Show[plot /@ {1, 2}, DisplayFunction -> $DisplayFunction] Of course one can display them together, though the result is not very impressive: In[19]:= Show[polytope, lines] Next, briefly, about the second approach. This method works very nicely with the cube example, but less well with the more complex one. I 'll describe it only for the cube: We only the package: In[1]:= << Graphics`ContourPlot3D` define the test function: In[2]:= test[x_, y_, z_] = x >= -1 && x <= 1 && y >= -1 && y <= 1 && z >= -1 && z <= 1; and the corresponding numerical function: In[3]:= f[x_, y_, z_] /; test[x, y, z] := 1; f[x_, y_, z_] := 0 Unfortunately ContourPlot3D is useless here so we need to use ListContourPLot3D. We make a table of data In[5]:= data = Table[f[x, y, z], {z, -2, 2, .25}, {y, -2, 2, .25}, {x, -2, 2, .25}]; now: In[6]:= ListContourPlot3D[data, Contours -> {1}, Lighting -> False, Axes -> True, ContourStyle -> {{RGBColor[0, 1, 0]}}] will produce a nice cube. Unfortunately when we apply the same method to th eother problem the computation takes far longer and the resulting polytope is very inaccurate. -- Andrzej Kozlowski Toyama International University JAPAN http://platon.c.u-tokyo.ac.jp/andrzej/ http://sigma.tuins.ac.jp/