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/