How to plot a 2D region when there are 3 free parameters?

Question

How could I obtain the region of payments, with the parameters a, b, c, d between 0 and 1 and the sum of all equal to 1?

After the comment/recommendation of belisarius, I edited the original question, and ask you a question in the last line of code.

P = {{a, b}, {c, d}}

A = {{2, -1}, {0, 1}}    
B = {{-1, 1}, {2, -1}}

A[[1]].P[[1]] + A[[2]].P[[2]]

f1[a_, b_, c_, d_] = A[[1]].P[[1]] + A[[2]].P[[2]]
f2[a_, b_, c_, d_] = B[[1]].P[[1]] + B[[2]].P[[2]]

ParametricPlot[{f1[a, b, c, d], f2[a, b, c, d]}, {  a, 0, 1}, {b, 0, 
  1}, {c, 0, 1}, {d, 0, 1} ]

In the last line, how can I code to consider only a, b, c, d in $[0,1]$ with $a+b+c+d=1$?

It´s similar to obatining all the convex combinations of the points $(a(i,j),b(i,j))$

@belisarius Yes It´s not very clear! :-) I want to draw the pairs (f1[a,b,c,d],f2[a,b,c,d]), with a,b,c,d between 0 and 1, and verifyng a+b+c+d=1 — Mika Ike, Dec 20 '13 at 16:17
It´s possible that it´s difficult to use Manipulate in this case, because there are 4 parameters (a,b,c,d) interconnected and depending each other. — Mika Ike, Dec 20 '13 at 16:18

score 10 · Answer 1 · edited Apr 13 '17 at 12:55

10

Being this the first answer at the site using the Region context, I think it deserves an entry on its own:

   AppendTo[$ContextPath, "Region`"]
   z = 0; 
   RegionPlot[
         RegionProperty[
                  TransformedRegion[
                        ParametricRegion[{{p1, p2, p4}, 
                                          0 < p1 < 1 && 0 < p2 < 1 - p1 && 0 < p4 < 1 - p1 - p2}], 
                         {2 #1 - #2 + #3, 2 - 3 #1 - #2 - 3 #3} &], 
                   {x, y, z}, ImplicitDescription], 
         {x, -3, 3}, {y, -3, 3}]

Mathematica graphics

edited Apr 13 '17 at 12:55

Community

1

answered Dec 21 '13 at 06:47

Dr. belisarius

115,881
13
203
453

Beautiful, many thanks! Btw the ParametricRegion can be defined in a more straightforward manner using just the problem statement : ParametricRegion[{{p1,p2,p3,p4},0<p1<1 && 0<p2<1 && 0<p3<1&& 0<p4<1 && p1+p2+p3+p4==1}],{2 #1 - #2 + #4, - #1 + #2 +2 #3 - #4}(...). – A.G. Dec 22 '13 at 01:09

score 5 · Answer 2 · edited Dec 22 '13 at 03:09

To start with, I will work in 4 dimensions to accommodate a,b,c,d.

P = {a, b, c, d};
M = {{2, -1, 0, 1}, {-1, 1, 2, -1}};

The points of your 2D set are the linear transforms, of the form

M.P  (* {2 a - b + d, -a + b + 2 c - d} *)

where P lies in some 4D set which is defined by $a+b+c+d=1$ and all positive, that is some polytope in 4D. Since you get your 2D points by a linear transform on the 4D points (M.P), and since the linear transform of a polytope is the convex hull of the transform of the extreme points of the polytope, your set will be the convex hull of the points M.P for all P being an extreme point of the 4D polytope.

Not, this 4D polytope is a simple one, you can get its extreme points by all combinations {a,b,c,d} of 0's and 1's that have a sum <= 1. You can generate all the combinations and then Select the right ones :

Flatten[Table[{a, b, c, d},
  {a, 0, 1}, {b, 0, 1}, {c, 0, 1}, {d, 0, 1}], 3]
extremePoints4D = Select[%, Total[#] <= 1 &]

(* {{0, 0, 0, 0}, {0, 0, 0, 1}, {0, 0, 1, 0}, {0, 1, 0, 0}, {1, 0, 0, 0}} *)

(yes, you could have guessed that just fine!)

Now, transform these to their 2D transforms :

Points2D = M.# & /@ extremePoints4D

(* {{0, 0}, {1, -1}, {0, 2}, {-1, 1}, {2, -1}} *)

(not all of these will end up extreme points of the 2D region but that's OK, we will take the convex hull and interior points do not matter) Now, plot the convex hull :

<< ComputationalGeometry`
Graphics[Polygon[Points2D[[ConvexHull[Points2D]]]], Frame -> True]

enter image description here

Note: This approach needs to generate all extreme points of a polytope; in higher dimensions (4D -> nD) and with constraints not so simple as $a+b+c+d=1$, getting those extreme points gets trickier.

About the bounding box

@Silvia got me thinking about the bounding box of the 2D region. With the above method it is quite obvious (min and max x and y over extreme points) but I thought I would suggest another way to do this that has the advantage of working in higher dimensions.

The idea is to get the min and max x and y over all possible (x,y). Linear programming will do just that :

(*
Min x (or -x, y, -y)
subject to linear constraints
a + b + c + d = 1
2a - b - d -x = 0
-a + b + 2c - d - y = 0
0 <= a, b, c, d <= 1
-infinity < x, y < infinity

with a 6D vector {a,b,c,d,x,y}.
*)

obj = {{0, 0, 0, 0, 1, 0}, {0, 0, 0, 0, -1, 0}, {0, 0, 0, 0, 0, 1}, {0, 0, 0, 0, 0, -1}};
(* obj will be used to min/max y/y *)
M = {{1, 1, 1, 1, 0, 0}, {2, -1, 0, 1, -1, 0}, {-1, 1, 2, -1, 0, -1}};
(* M.{a,b,c,d,x,y}=0 generates all 3 constraints *)

bounds = {{0, 1}, {0, 1}, {0, 1}, {0, 1}, {-Infinity, Infinity}, {-Infinity, Infinity}}

(* Now get the 4 solutions to min/max x/y: *)
sol = LinearProgramming[#, M, {{1, 0}, {0, 0}, {0, 0}}, bounds] & /@ obj
(* and keep the {x,y}'s
sol[[1 ;; 4, 5 ;; 6]]

(* all downhill from here *)

xmin = Min[points[[1 ;; 4, 1]]]
xmax = Max[points[[1 ;; 4, 1]]]
ymin = Min[points[[1 ;; 4, 2]]]
ymax = Max[points[[1 ;; 4, 2]]]
ListPlot[points, AxesOrigin -> {xmin - 1, ymin - 1}, 
 Epilog -> {Orange, Line[{{xmin, ymin}, {xmax, ymin}, {xmax, ymax}, {xmin, ymax}, {xmin, ymin}}]}]

The bounding box

Hi A.G., I guess you missed an image by carelessness. I added it back, hope you don't mind. — Silvia, Dec 22 '13 at 03:10

Dr. belisarius · Answer 3 · 2013-12-20T21:24:57.640

You have 3 free parameters, since the imposed condition over their sum takes one away.
So, let's plot in 3D the region for those three free parms:

RegionPlot3D[(0 < p1 < 1 && 0 < p2 < 1 && 0 < p4 < 1 &&  p1 + p2 + p4 < 1), 
             {p1, 0, 1}, {p2, 0, 1}, {p4, 0, 1}, PlotPoints -> 30]

Mathematica graphics

Which is equivalent to:

Graphics3D[gc = GraphicsComplex[{{0, 0, 0}, {1, 0, 0}, {0, 1, 0}, {0, 0, 1}}, 
           Polygon@{{1, 2, 3}, {2, 3, 4}, {1, 2, 4}, {1, 3, 4}}], 
           Axes -> True]

Mathematica graphics

Now your functions f1 and f2 can bee written as:

f[p_, x_] := Total[p x, 2]
p = {{p1, p2}, {p3, p4}};
a = {{2, -1}, {0, 1}};
b = {{-1, 1}, {2, -1}};
{f[p, a], f[p, b]}
(*
  {2 p1 - p2 + p4, -p1 + p2 + 2 p3 - p4}
*)

Which is your target transform.

Now we find a transformation to map any {p1, p2, p3, p4} (minus one parm set by the sum restriction) into your target:

r = ToRules@
    Reduce[ForAll[{p1, p2, p4}, 
                  TransformationFunction[{{a1, a2, a3, a4}, {b1, b2, b3, b4}, 
                                          {c1, c2, c3, c4}, {d1, d2, d3, d4}}][{p1, p2, p4}] == 
                  {f[p, a], f[p, b], 0} /. p3 -> 1 - p1 - p2 - p4], 
          {a1, a2, a3, a4, b1, b2, b3, b4, c1, c2, c3, c4, d1, d2, d3, d4}];

t = TransformationFunction[{{a1, a2, a3, a4}, {b1, b2, b3, b4}, 
                            {c1, c2, c3, c4}, {d1, d2, d3, d4}}] /. r /. a1 -> 1

And now we apply that transformation to our graphics complex:

Graphics3D[{FaceForm[Black], EdgeForm[None], 
           GeometricTransformation[gc, t]}, Axes -> True]

Mathematica graphics

And we can show it in 2D as:

Graphics[Normal@ GeometricTransformation[gc, t] /. {x_, y_, 0} :> {x, y} /. 
         (VertexNormals -> x_) :> Sequence[], Frame -> True]

Mathematica graphics

Your first RegionPlot3D could be simpler: RegionPlot3D[a + b + c <= 1, {a, 0, 1}, {b, 0, 1}, {c, 0, 1}] — A.G., Dec 21 '13 at 03:09

score 3 · Answer 4 · edited Apr 13 '17 at 12:55

3

Inspired by belisarius and Simon's answers, I came up with a similar method, but without burdensome RegionPlot.

Graphics`Region`RegionInit[];

parareg = TransformedRegion[
                ParametricRegion[{
                         {p1, p2, p4},
                         0 < p1 < 1 && 0 < p2 < 1 - p1 && 0 < p4 < 1 - p1 - p2
                         }],
                {2 #1 - #2 + #3, 2 - 3 #1 - #2 - 3 #3} &]

meshreg = MeshRegion@RegionToMeshObject[Evaluate@{
                  {x, y},
                  Prepend[
                   (* Thanks to belisarius, we have a good way to obtain the boundingbox: *)
                          1.1 RegionProperty[parareg, "BoundingBox"],
                          RegionProperty[parareg, {x, y}, "ImplicitDescription"]
                         ]
                }]

GeometryPlot[meshreg, PlotRange -> 1.1 RegionProperty[parareg, "BoundingBox"]]

region plot

edited Apr 13 '17 at 12:55

Community

1

answered Dec 22 '13 at 00:23

Silvia

27,556
3
84
164

Why doesn't just GeometryPlot@RegionConvert@parareg show the result? – Dr. belisarius Dec 22 '13 at 01:03
@belisarius It looks like RegionConvert will truncate the object at the boundary of a $2\times 2$ square. – Silvia Dec 22 '13 at 01:16
Yep, GeometryPlot@ RegionConvert@ TransformedRegion[parareg, ScalingTransform[1/10 {1, 1}]] ... but there must be a way to tell RegionConvert to expand the domain – Dr. belisarius Dec 22 '13 at 01:23
1

@belisarius I Trace-ed the RegionConvert and eventually arrived RegionToMeshObject, where the {{-1,1},{-1,1}} boundary is just inserted "suddenly". I think more internal trace will be needed for revealing the key point. – Silvia Dec 22 '13 at 01:32
@Silvia I just added the bounding box computation to my response (above). – A.G. Dec 22 '13 at 03:04
1

@Silvia Try RegionProperty[parareg, "BoundingBox"] – Dr. belisarius Jun 10 '14 at 13:03
@belisarius That's Nice! I've included your method in my answer :) So are you digging the version in mma 9 or the up coming one in 10? – Silvia Jun 10 '14 at 16:18
@Silvia I'm crawling on V9. In fact I have a need to solve this very problem http://mathematica.stackexchange.com/q/6357/193 using Voronoi tessellations. I'm close, but not yet done. – Dr. belisarius Jun 10 '14 at 17:58
@belisarius Aha~ That question is in my TODO list too :D But I don't have time for it yet ... Looking forward to your solution :) – Silvia Jun 11 '14 at 03:35
@Silvia I'm delayed because the v9 implementation of the Mesh and Region contexts are buggy. I'm still trying to decide if I'll implement my own or wait till v10. Oh laziness!, mother of us all- – Dr. belisarius Jun 12 '14 at 16:10
@belisarius It looks like many of them are documented function in v10 with the same name. I guess even if you write your code in v9, it won't be hard to upgrade to v10. – Silvia Jun 12 '14 at 17:12

Michael E2 · Answer 5 · 2013-12-22T03:43:50.447

Reduce can be used to turn the system

{x == 2 a - b + d, y == -a + b + 2 c - d,
 0 <= a <= 1, 0 <= b <= 1, 0 <= c <= 1, 0 <= d <= 1,
 a + b + c + d == 1}

into a system of inequalities of the form (via a CylindricalDecomposition)

(x1 < x <= x2 && f1[x] < y < f2[x] && ..) ||
..

where the omitted inequalities in the first line constrain the parameters a, b, c, and d in terms of x and y. Only the first pair of each && expression is needed to plot the {x, y} region. From there, the functions f1 and f2 form the boundary of the region, so it can be plotted.

The command that will convert the system is basically

Reduce[system, {x, y, a, b, c, d}, Reals, Backsubstitution -> True]

A key is putting x and y first. Then we can extract the x and y parts with

Reduce[system, {x, y, a, b, c, d}, Reals, Backsubstitution -> True][[All, {1, 2}]]

Some of the entries will be strict equalities. These represent boundary points or lines that bound the inequalities and may be removed. Finally, some of the cases may be combined. Reduce initially divides the region according to the x coordinates of points of intersection and then it subdivides these. If we gather the ones with the same x inequality (with Split below), Reduce can then be used to combine them.

Here is this first step:

Clear[x, y, a, b, c, d];

P = {{a, b}, {c, d}};
A = {{2, -1}, {0, 1}};
B = {{-1, 1}, {2, -1}};

redSystem = Reduce[{
     {x, y} == Inner[Dot, {Hold @@@ A, Hold @@@ B}, Hold @@@ P, Plus] /. Hold -> List,
     0 <= a <= 1, 0 <= b <= 1, 0 <= c <= 1, 0 <= d <= 1,
     a + b + c + d == 1},
   {x, y, a, b, c, d},
   Reals, 
   Backsubstitution -> True];

removeBoundaries[ineqs_] := Select[ineqs, ! MemberQ[#, _Equal] &];

pieces = removeBoundaries @ redSystem[[All, {1, 2}]];
combinedPieces = Reduce[#, {x, y}, Reals] & /@ Or @@@ SplitBy[List @@ pieces, #[[1, 1]] &]

(* {-1   < x <= 0   && -x < y < 2 + x, 
     0   < x <= 5/7 && -x < y < 1/2 (4 - 3 x), 
     5/7 < x <= 1   && -x < y < 1/2 (4 - 3 x), 
     1    < x < 2   && -1 < y < 1/2 (4 - 3 x)} *)

There are two ways to process these inequalities. If the image is a polygon, as in this case, then we can get its convex hull. If the region is bounded by curves, then we can use Plot.

Convex Hull method

vertices[ineqx_ && ineqy_] := 
  With[{ineq1 = List @@ ineqx, ineq2 = List @@ ineqy}, 
   Thread[{x, ineq2[[{1, -1}]]}] /. ({x -> #} & /@ #) &@ ineq1[[{1, -1}]]];

points = DeleteDuplicates @ Flatten[vertices /@ combinedPieces, 2]

(* {{-1, 1}, {0, 0}, {0, 2}, {5/7, -(5/7)}, {5/7, 13/14}, {1, -1}, {1, 1/2}, {2, -1}} *)

Needs["ComputationalGeometry`"];
Graphics[{
  EdgeForm[ColorData[1][1]], RGBColor[0.368417, 0.506779, 0.709798], Opacity[0.2], 
  Polygon[points[[ComputationalGeometry`ConvexHull @ points]]]},
 Frame -> True]

Mathematica graphics

Plot method

Clear[rplot];
rplot[_[x1_, ___, x2_] && _[y1_, ___, y2_], opts : OptionsPattern[Plot]] := 
 Plot[{y1, y2}, {x, x1, x2}, Filling -> {1 -> {2}}, opts]

Show[
 rplot[#,
    PlotStyle -> ColorData[1][1], 
    FillingStyle -> 
     Directive[Opacity[1], Lighter[RGBColor[0.368417, 0.506779, 0.709798], 0.8]]] & /@ 
  combinedPieces,
 PlotRange -> All, Frame -> True, Axes -> False, AspectRatio -> Automatic
 ]

Mathematica graphics

How to plot a 2D region when there are 3 free parameters?

5 Answers5

Linked