Linear inequality solving
The LinearSolve command in the RegularChains package, whose purpose is to solve systems of linear equations and inequalities, has been enhanced in several ways:
A more efficient algorithm based on Fourier-Motzkin elimination, has been implemented.
Two new options, canonical and projection, for controlling the form of the output have been added.
The examples below illustrate the use of this command and the two new options for computations with polyhedra.
Change of representation for polyhedra
Dual polyhedra
See Also
There are at least two ways to mathematically describe a polyhedron, either by specifying its corners, or by specifying its bounding hyperplanes. This is illustrated using a cube of side length 1 with one corner at the origin and its three adjacent edges parallel to the three coordinate axes.
Plot the cube:
withplots:withplottools:
displaycuboid0,0,0,1,1,1,axes=boxed,labels=x,y,z
Define this cube by bounding planes:
c1≔x≥0,x≤1,y≥0,y≤1,z≥0,z≤1
0≤x,x≤1,0≤y,y≤1,0≤z,z≤1
Alternatively, define the vertices of the same cube:
V≔Matrix0,0,0,0,0,1,0,1,0,0,1,1,1,0,0,1,0,1,1,1,0,1,1,1
Suppose now that you are given the vertex representation and want to automatically derive the bounding planes representation. The cube is defined as the convex hull of all its vertices:
opconvertVectorrowx,y,z=~addλi⋅Vi,i=1..8,list
x=λ5+λ6+λ7+λ8,y=λ3+λ4+λ7+λ8,z=λ2+λ4+λ6+λ8
addλi,i=1..8=1
λ1+λ2+λ3+λ4+λ5+λ6+λ7+λ8=1
seqλi≥0,i=1..8
0≤λ1,0≤λ2,0≤λ3,0≤λ4,0≤λ5,0≤λ6,0≤λ7,0≤λ8
seqλi≤1,i=1..8
λ1≤1,λ2≤1,λ3≤1,λ4≤1,λ5≤1,λ6≤1,λ7≤1,λ8≤1
sys≔,,,:
You can now use the LinearSolve command to eliminate all λi from this system of equations and inequalities. The projection option is used to specify that elimination is required, and the canonical option to eliminate redundant equations and inequalities.
withRegularChains:withSemiAlgebraicSetTools:
vars≔seqλi,i=1..8,z,y,x
λ1,λ2,λ3,λ4,λ5,λ6,λ7,λ8,z,y,x
LinearSolvesys,PolynomialRingvars,canonical,projection=3
This agrees with the representation c1 above.
For every convex 3-D polyhedron, there is a dual polyhedron with the roles of corners and faces swapped. Let us compute the dual polyhedron of the cube. Start with a representation of the cube with the origin as the center.
c2≔x≤1,−x≤1,y≤1,−y≤1,z≤1,−z≤1
x≤1,−x≤1,y≤1,−y≤1,z≤1,−z≤1
Now rewrite this in matrix form.
B≔Matrix1,0,0,−1,0,0,0,1,0,0,−1,0,0,0,1,0,0,−1
v≔Vectorx,y,z
c≔Vector1,1,1,1,1,1
B.v≤~c
In general, if a bounded convex polyhedron is given by a set of inequalities of the form B.v≤1, then the vertices of its dual are exactly the columns of the matrix B. The convex hull of these six vertices defines an octahedron. Using the same method as before, compute the bounding planes for this octahedron.
opconvertVectorrowx,y,z=~addλi⋅Bi,i=1..6,list
x=λ1−λ2,y=λ3−λ4,z=λ5−λ6
addλi,i=1..6=1
λ1+λ2+λ3+λ4+λ5+λ6=1
seqλi≥0,i=1..6
0≤λ1,0≤λ2,0≤λ3,0≤λ4,0≤λ5,0≤λ6
seqλi≤1,i=1..6
λ1≤1,λ2≤1,λ3≤1,λ4≤1,λ5≤1,λ6≤1
vars≔seqλi,i=1..6,z,y,x
λ1,λ2,λ3,λ4,λ5,λ6,z,y,x
S≔LinearSolvesys,PolynomialRingvars,canonical,projection=3
−1≤x,x≤1,−x−1≤y,−1+x≤y,y≤1+x,y≤1−x,−y−x−1≤z,−y+x−1≤z,−1−x+y≤z,−1+x+y≤z,z≤y+x+1,z≤y−x+1,z≤1+x−y,z≤1−x−y
The first two inequalities represent that projection of the octahedron onto the x axis:
S1..2
−1≤x,x≤1
Together with the following four inequalities, which involve only x and y but not z, this defines the projection of the octahedron onto the x,y-plane:
S1..2; S3..6
−x−1≤y,−1+x≤y,y≤1+x,y≤1−x
The following plot shows this projection (blue) together with the four boundary lines.
inequalS3..6,x=−2..2,y=−2..2
Finally, the last 8 inequalities, i.e., the ones containing z, define the bounding planes of the octahedron.
S−8..−1
−y−x−1≤z,−y+x−1≤z,−1−x+y≤z,−1+x+y≤z,z≤y+x+1,z≤y−x+1,z≤1+x−y,z≤1−x−y
The following 3-D plot displays the tetrahedron as well as the last of these 8 bounding planes. Rotate the plot with the mouse in order to see that the plane actually touches one of the faces.
plane≔Student:-LinearAlgebra:-PlanePlotlhs−rhsS−1=0,x,y,z,caption=,shownormal=false,showpoint=false,planeoptions=transparency=0.3:
displayoctahedron0,0,0,plane,axes=boxed,labels=x,y,z
RegularChains[SemiAlgebraicSetTools][LinearSolve]
Download Help Document