Tracer les inégalités linéaires dans Mathematica

Question

J'ai des systèmes linéaires d'inégalités dans 3 variables et j'aimerais tracer ces régions. Idéalement, j'aimerais quelque chose qui ressemble à des objets dans PolyhedronData. J'ai essayé RegionPlot3D, mais les résultats sont visuellement pauvres et trop polygone lourd pour tourner en temps réel

Voici ce que je veux dire, le code ci-dessous génère 10 ensembles de contraintes linéaires et les parcelles

 
randomCons := Module[{}, 
    hadamard = KroneckerProduct @@ Table[{{1, 1}, {1, -1}}, {3}]; 
    invHad = Inverse[hadamard]; 
    vs = Range[8]; 
    m = mm /@ vs; 
    sectionAnchors = Subsets[vs, {1, 7}]; 
    randomSection := 
    Mean[hadamard[[#]] & /@ #] & /@ 
    Prepend[RandomChoice[sectionAnchors, 3], vs]; {p0, p1, p2, p3} = 
    randomSection; 
    section = 
    Thread[m -> 
     p0 + {x, y, z}.Orthogonalize[{p1 - p0, p2 - p0, p3 - p0}]]; 
    And @@ Thread[invHad.m >= 0 /. section] 
    ]; 
Table[RegionPlot3D @@ {randomCons, {x, -3, 3}, {y, -3, 3}, {z, -3, 
    3}}, {10}] 

Toutes les suggestions ?

Mise à jour: L'intégration des suggestions ci-dessous, voici la version que je fini par utiliser pour tracer la région réalisable d'un système d'inégalités linéaires

 
(* Plots feasible region of a linear program in 3 variables, \ 
specified as cons[[1]]>=0,cons[[2]]>=0,... 
Each element of cons must \ 
be an expression of variables x,y,z only *) 

plotFeasible3D[cons_] := 
Module[{maxVerts = 20, vcons, vertCons, polyCons}, 
    (* find intersections of all triples of planes and get rid of \ 
intersections that aren't points *) 

    vcons = Thread[# == 0] & /@ Subsets[cons, {3}]; 
    vcons = Select[vcons, Length[Reduce[#]] == 3 &]; 
    (* Combine vertex constraints with inequality constraints and find \ 
up to maxVerts feasible points *) 
    vertCons = Or @@ (And @@@ vcons); 
    polyCons = And @@ Thread[cons >= 0]; 
    verts = {x, y, z} /. 
    FindInstance[polyCons && vertCons, {x, y, z}, maxVerts]; 
    ComputationalGeometry`Methods`ConvexHull3D[verts, 
    Graphics`Mesh`FlatFaces -> False] 
    ] 

code pour tester

 
randomCons := Module[{}, 
    hadamard = KroneckerProduct @@ Table[{{1, 1}, {1, -1}}, {3}]; 
    invHad = Inverse[hadamard]; 
    vs = Range[8]; 
    m = mm /@ vs; 
    sectionAnchors = Subsets[vs, {1, 7}]; 
    randomSection := 
    Mean[hadamard[[#]] & /@ #] & /@ 
    Prepend[RandomChoice[sectionAnchors, 3], vs]; {p0, p1, p2, p3} = 
    randomSection; 
    section = 
    Thread[m -> 
     p0 + {x, y, z}.Orthogonalize[{p1 - p0, p2 - p0, p3 - p0}]]; 
    And @@ Thread[invHad.m >= 0 /. section] 
    ]; 
Table[plotFeasible3D[List @@ randomCons[[All, 1]]], {50}]; 

Answer 1

Voici un petit programme qui semble faire ce que vous voulez:

rstatic = randomCons;     (* Call your function *) 
randeq = rstatic /. x_ >= y_ -> x == y; (* make a set of plane equations 
              replacing the inequalities by == *) 

eqset = Subsets[randeq, {3}];   (* Make all possible subsets of 3 planes *) 

             (* Now find the vertex candidates 
              Solving the sets of three equations *) 
vertexcandidates =  
    Flatten[Table[Solve[eqset[[i]]], {i, Length[eqset]}], 1]; 

             (* Now select those candidates 
              satisfying all the original equations *)   
vertex = Union[Select[vertexcandidates, rstatic /. # &]]; 

             (* Now use an UNDOCUMENTED Mathematica 
              function to plot the surface *) 

gr1 = ComputationalGeometry`Methods`ConvexHull3D[{x, y, z} /. vertex]; 
             (* Your plot follows *) 
gr2 = RegionPlot3D[rstatic, 
     {x, -3, 3}, {y, -3, 3}, {z, -3, 3}, 
     PerformanceGoal -> "Quality", PlotPoints -> 50] 

Show[gr1,gr2] (*Show both Graphs superposed *) 

Le résultat est:

Inconvénient: la fonction non documentée est pas parfait .Lorsque le visage est pas un triangle, il affichera une triangulation:

Modifier

Il y a une option pour se débarrasser de la triangulation faute

ComputationalGeometry`Methods`ConvexHull3D[{x, y, z} /. vertex, 
           Graphics`Mesh`FlatFaces -> False] 

-ce que le la magie. Exemple:

Edit 2

Comme je l'ai mentionné dans un commentaire, voici deux ensembles de sommet dégénéré généré par vos randomCons

{{x -> Sqrt[3/5]}, 
    {x -> -Sqrt[(5/3)] + Sqrt[2/3] y}, 
    {x -> -Sqrt[(5/3)], y -> 0}, 
    {y -> -Sqrt[(2/5)], x -> Sqrt[3/5]}, 
    {y -> 4 Sqrt[2/5], x -> Sqrt[3/5]} 
} 

et

{{x -> -Sqrt[(5/3)] + (2 z)/Sqrt[11]}, 
{x -> Sqrt[3/5], z -> 0}, 
{x -> -Sqrt[(5/3)], z -> 0}, 
{x -> -(13/Sqrt[15]), z -> -4 Sqrt[11/15]}, 
{x -> -(1/Sqrt[15]), z -> 2 Sqrt[11/15]}, 
{x -> 17/(3 Sqrt[15]), z -> -((4 Sqrt[11/15])/3)} 
} 

Toujours essayer de voir comment faire face doucement avec les ...

Edit 3

Ce code ne suffit pas en général pour le problème complet, mais élimine le problème de degenerancy cylindres pour votre générateur de données d'échantillon. J'ai utilisé le fait que les cas pathogènes sont toujours des cylindres avec leur axe parallèle à l'un des axes de coordonnées, et ensuite utilisé RegionPlot3D pour les tracer. Je ne sais pas si cela sera utile pour votre cas général :(.

For[i = 1, i <= 160, i++, 
rstatic = randomCons; 
r[i] = rstatic; 
s1 = Reduce[r[i], {x, y, z}] /. {x -> var1, y -> var2, z -> var3}; 
s2 = Union[StringCases[ToString[FullForm[s1]], "var" ~~ DigitCharacter]]; 

If [Dimensions@s2 == {3}, 

    (randeq = rstatic /. x_ >= y_ -> x == y; 
    eqset = Subsets[randeq, {3}]; 
    vertexcandidates = Flatten[Table[Solve[eqset[[i]]], {i, Length[eqset]}], 1]; 
    vertex = Union[Select[vertexcandidates, rstatic /. # &]]; 
    a[i] = ComputationalGeometry`Methods`ConvexHull3D[{x, y, z} /. vertex, 
      Graphics`Mesh`FlatFaces -> False, Axes -> False, PlotLabel -> i]) 
    , 

    a[i] = RegionPlot3D[s1, {var1, -2, 2}, {var2, -2, 2}, {var3, -2, 2}, 
      Axes -> False, PerformanceGoal -> "Quality", PlotPoints -> 50, 
      PlotLabel -> i, PlotStyle -> Directive[Yellow, Opacity[0.5]], 
      Mesh -> None] 
    ]; 
] 

GraphicsGrid[Table[{a[i], a[i + 1], a[i + 2]}, {i, 1, 160, 4}]] 

Here vous pouvez trouver une image de la sortie générée, les cas dégénérés (tous les cylindres) sont en jaune transparent

HTH

Answer 2

Triplet choisi parmi votre ensemble d'inégalités déterminera généralement un point obtenu en résolvant le triplet d'équations correspondant. Je crois que vous voulez la coque convexe de cet ensemble de points. Vous pouvez générer cela comme ça.

cons = randomCons; (* Your function *) 
eqs = Apply[Equal, List @@@ Subsets[cons, {3}], {2}]; 
sols = Flatten[{x, y, z} /. Table[Solve[eq, {x, y, z}], {eq, eqs}], 1]; 
pts = Select[sols, And @@ (NumericQ /@ #) &]; 
ComputationalGeometry`Methods`ConvexHull3D[pts] 

Bien sûr, certains triplés pourrait en fait être et conduire à sous-déterminé lignes ou evan tout un plan. Ainsi, le code émettra une plainte dans ces cas.

Cela semble fonctionner dans les quelques cas aléatoires que j'ai essayés mais, comme Yaro le fait remarquer, cela ne fonctionne pas du tout. L'image suivante illustre exactement pourquoi:

{p0, p1, p2, 
    p3} = {{1, 0, 0, 0, 0, 0, 0, 0}, {1, 1/2, -(1/2), 0, -(1/2), 0, 
    0, -(1/2)}, {1, 0, 1/2, 1/2, 0, 0, -(1/2), 1/2}, {1, -(1/2), 1/2, 
    0, -(1/2), 0, 0, -(1/2)}}; 
hadamard = KroneckerProduct @@ Table[{{1, 1}, {1, -1}}, {3}]; 
invHad = Inverse[hadamard]; 
vs = Range[8]; 
m = mm /@ vs; 
section = 
    Thread[m -> 
    p0 + {x, y, z}.Orthogonalize[{p1 - p0, p2 - p0, p3 - p0}]]; 
cons = And @@ Thread[invHad.m >= 0 /. section]; 
eqs = Apply[Equal, List @@@ Subsets[cons, {3}], {2}]; 
sols = Flatten[{x, y, z} /. Table[Solve[eq, {x, y, z}], {eq, eqs}], 
    1]; // Quiet 
pts = Select[sols, And @@ (NumericQ /@ #) &]; 
ptPic = Graphics3D[{PointSize[Large], Point[pts]}]; 
regionPic = 
    RegionPlot3D[cons, {x, -2, 2}, {y, -2, 2}, {z, -2, 2}, 
    PlotPoints -> 40]; 
Show[{regionPic, ptPic}] 

Ainsi, il y a des points qui sont finalement coupés par le plan défini par une autre contrainte. Voici un moyen (je suis absolument inefficace) de trouver ceux que vous voulez.

regionPts = regionPic[[1, 1]]; 
nf = Nearest[regionPts]; 
trimmedPts = Select[pts, Norm[# - nf[#][[1]]] < 0.2 &]; 
trimmedPtPic = Graphics3D[{PointSize[Large], Point[trimmedPts]}]; 
Show[{regionPic, trimmedPtPic}] 

Ainsi, vous pouvez utiliser la coque convexe de trimmedPts. Cela dépend finalement du résultat de RegionPlot et vous pourriez avoir besoin de rampe de la valeur de PlotPoints pour le rendre plus fiable. Googler sur un peu révèle le concept d'une région de faisabilité en programmation linéaire. Cela semble être exactement ce que vous cherchez.

Mark

Answer 3

Voir toutes les réponses précédentes,! ce qui est mal à l'aide de la construction en fonction RegionPlot3D, par exemple

RegionPlot3D[ 2*y+3*z <= 5 && x+y+2*z <= 4 && x+2*y+3*z <= 7 && 
       x >= 0 && y >= 0 && z >= 0, 
      {x, 0, 4}, {y, 0, 5/2}, {z, 0, 5/3} ]