Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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

After the comment/recommendation of belisarius... I edit 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 for consider only a,b,c,d in [0,1] with a+b+c+d=1 ?

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

share|improve this question
    
At least for me, it isn't clear what you're trying to do –  belisarius Dec 20 '13 at 16:11
    
@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
add comment

5 Answers

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

share|improve this answer
    
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 1:09
add comment

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

share|improve this answer
    
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 3:10
add comment

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

share|improve this answer
    
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 3:09
add comment

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

share|improve this answer
add comment

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[
                          (* I didn't find a good way to obtain the boundingbox: *)
                          1.1 Outer[
                                    RegionNearest[parareg,
                                       {Identity, Reverse}[[#1]]@{(-1)^#2 10^3, 0}
                                       ][[#1]] &,
                                    {1, 2}, {1, 2}],
                          RegionProperty[parareg, {x, y}, "ImplicitDescription"]
                         ]
                }]

GeometryPlot[meshreg, PlotRange -> {{-3, 3}, {-3, 3}}]

region plot

share|improve this answer
    
Why doesn't just GeometryPlot@RegionConvert@parareg show the result? –  belisarius Dec 22 '13 at 1:03
    
@belisarius It looks like RegionConvert will truncate the object at the boundary of a $2\times 2$ square. –  Silvia Dec 22 '13 at 1:16
    
Yep, GeometryPlot@ RegionConvert@ TransformedRegion[parareg, ScalingTransform[1/10 {1, 1}]] ... but there must be a way to tell RegionConvert to expand the domain –  belisarius Dec 22 '13 at 1: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 1:32
    
@Silvia I just added the bounding box computation to my response (above). –  A.G. Dec 22 '13 at 3:04
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.