绘制Mathematica中的线性不等式 [英] Plotting linear inequalities in Mathematica

查看:161
本文介绍了绘制Mathematica中的线性不等式的处理方法,对大家解决问题具有一定的参考价值,需要的朋友们下面随着小编来一起学习吧!

问题描述

我有3个变量的不等式的线性系统,我想绘制这些区域。
理想情况下,我希望看起来像PolyhedronData中的对象。我尝试了RegionPlot3D,但结果是视觉上很差,太多边形以至于无法实时旋转。下面是我的意思,下面的代码生成10组线性约束,并且绘制它们

 
randomCons:= Module [{},
hadamard = KroneckerProduct @@ Table [{{1,1}, {1,-1}},{3}];
invHad =反[hadamard];
vs =范围[8];
m = mm / @ vs;
sectionAnchors =子集[vs,{1,7}];
randomSection:=
平均值[hadamard [[#]]&/ @#]&/ @
前置[RandomChoice [sectionAnchors,3],vs]; {p0,p1,p2,p3} =
randomSection;
section =
Thread [m - >
p0 + {x,y,z} .Orthogonalize [{p1 - p0,p2 - p0,p3 - p0}]];
和@@线程[invHad.m> = 0 /。部分]
];
表[RegionPlot3D @@ {randomCons,{x,-3,3},{y,-3,3},{z,-3,
3}},{10}]

任何建议? 更新 ,这里是我最终用来绘制一个线性不等式系统的可行区域的版本

 
(*绘制线性程序的可行区域in 3个变量,\
指定为cons [[1]]> = 0,cons [[2]]> = 0,...
每个元素必须是\
是* b
$ b plotFeasible3D [cons_]:=
模块[{maxVerts = 20,vcons,vertCons,polyCons},
(*找到飞机所有三元组的交点,并去除不是点的\
交点*)

vcons =线程[#== 0]&/ @ Subsets [cons,{ 3}];
vcons =选择[vcons,长度[Reduce [#]] == 3&];
(*将顶点约束与不等式约束结合起来,并找到\
直到maxVerts可行点*)
vertCons =或@@(和@@@ vcons);
polyCons =和@@线程[cons> = 0];
verts = {x,y,z} /。
FindInstance [polyCons && vertCons,{x,y,z},maxVerts];
ComputationalGeometry`Methods`ConvexHull3D [verts,
Graphics`Mesh`FlatFaces - > False]
]

测试代码

 
randomCons:= Module [{},
hadamard = KroneckerProduct @@ Table [{{1, 1},{1,-1}},{3}];
invHad =反[hadamard];
vs =范围[8];
m = mm / @ vs;
sectionAnchors =子集[vs,{1,7}];
randomSection:=
平均值[hadamard [[#]]&/ @#]&/ @
前置[RandomChoice [sectionAnchors,3],vs]; {p0,p1,p2,p3} =
randomSection;
section =
Thread [m - >
p0 + {x,y,z} .Orthogonalize [{p1 - p0,p2 - p0,p3 - p0}]];
和@@线程[invHad.m> = 0 /。部分]
];
Table [plotFeasible3D [List @@ randomCons [[All,1]]],{50}];


解决方案

这是一个小程序,它可以做你想做的事情: b
$ b

  rstatic = randomCons; (*调用你的函数*)
randeq = rstatic /。 x_> = y_ - > x == y; (*做一组平面方程
用== *替换不等式)

eqset =子集[randeq,{3}]; (*使所有可能的3个平面的子集*)

(*现在找到顶点候选
求解三个方程的集合*)
vertexcandidates =
Flatten [表[求解[eqset [[i]]],{i,长度[eqset]}],1];

(*现在选择满足所有原始方程*的候选

vertex = Union [选择[vertexcandidates,rstatic /。 #&]];

(*现在使用一个未记录的Mathematica
函数绘制曲面*)

gr1 = ComputationalGeometry`Methods`ConvexHull3D [{x,y,z} / 。顶点];
(*您的绘图如下*)
gr2 = RegionPlot3D [rstatic,
{x,-3,3},{y,-3,3},{z,-3,3 },
PerformanceGoal - > 质量,PlotPoints - > 50]

显示[gr1,gr2](*显示叠加的图形*)



结果是:





缺点:无证函数并不完美。当脸部不是三角形时,它会显示一个三角形:



编辑



摆脱犯规的三角剖分

  ComputationalGeometry`Methods`ConvexHull3D [{x,y,z} /。顶点,
Graphics`Mesh`FlatFaces - > False]

可以做到这一点。示例:



编辑2



正如我在评论中提到的,这里有两组退化顶点,您的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平方[2/5],x - > Sqrt [3/5]}
}

  {{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)}
}

Still试图看清楚如何轻松应对这些......

编辑3



此代码对于完整问题来说不够全面,但可以消除样本数据生成器的圆柱体退化问题。我使用了这样的事实:致病病例总是与其轴线平行于其中一个坐标轴的圆柱体,然后使用RegionPlot3D绘制它们。
我不确定这是否对您的一般情况有用:(。

  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]];

如果[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 - >质量,PlotPoints - > 50,
PlotLabel - >我,PlotStyle - >指令[黄色,不透明度[0.5]],
网格 - >无]
];


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

Here 可以找到生成的输出图像,退化的案例(所有圆柱体)都是透明的黄色


I have linear systems of inequalities in 3 variables and I'd like to plot these regions. Ideally, I'd like something that looks like objects in PolyhedronData. I tried RegionPlot3D, but the results are visually poor and too polygon-heavy to rotate in real time

Here's what I mean, the code below generates 10 sets of linear constraints and plots them

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}]

Any suggestions?

Update: Incorporating suggestions below, here's the version I ended up using to plot feasible region of a system of linear inequalities

(* 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 for testing

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}];

解决方案

Here is a small program that seems to do what you want:

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 *)

The result is:

Downside: the undocumented function is not perfect. When the face is not a triangle, it will show a triangulation:

Edit

There is an option to get rid of the foul triangulation

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

does the magic. Sample:

Edit 2

As I mentioned in a comment, here are two sets of degenerate vertex generated by your 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]}
 }

and

{{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)}
}

Still trying to see how to cope gently with those ...

Edit 3

This code is not general enough for the full problem, but eliminates the cylinder degenerancy problem for your sample data generator. I used the fact that the pathogenic cases are always cylinders with their axis paralell to one of the coordinate axis, and then used RegionPlot3D to plot them. I'm not sure if this will be useful for your general case :(.

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 you can find an image of the generated output, the degenerated cases (all cylinders) are in transparent yellow

HTH!

这篇关于绘制Mathematica中的线性不等式的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!

查看全文
登录 关闭
扫码关注1秒登录
发送“验证码”获取 | 15天全站免登陆