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.

Four-sided convex land plots are usually denoted with edge lengths and area.

How do I create a Polygon Object with these parameters in Mathematica ?

share|improve this question
    
What have you tried so far? –  s0rce Jun 12 '13 at 21:33
    
I tried tinkering on Wolfram Alpha to get the inner angles but so far I couldn't think of any way to accomplish it :/ –  a3f Jun 12 '13 at 22:17
    
You can write the area of a convex quadrilateral in terms of the angle between the diagonals and edge lengths or angle between different sides and the side lengths, etc. Try formulating the problem using those... –  rm -rf Jun 12 '13 at 23:11
    
Writing some code to bruteforce it at the moment. Are there any 'more elegant' solutions ? –  a3f Jun 12 '13 at 23:13
    
As noted by rm, you might be able to use Bretschneider's formula to reckon out admissible ranges of diagonal lengths. From this, you can build up the two triangles that comprise your quadrilateral. –  J. M. Jun 13 '13 at 6:52

2 Answers 2

up vote 7 down vote accepted

I'll illustrate one way with an example. The idea is to fix the first two points at origin and on positive x axis, set up equations for the lengths and the area, solve them, then get rid of some symmetric solutions by insisting the third point have positive y value.

I cribbed the area formula from some random Internet site so I'm certain it's correct.

lens = {3, 1, 4, 2};
area = 3;

p1 = {0, 0};
p2 = {First[lens], 0};
p3 = {x3, y3};
p4 = {x4, y4};
{a, b, c, d} = lens;
psq = (p3 - p1).(p3 - p1);
qsq = (p4 - p2).(p4 - p2);
polys = {(p3 - p2).(p3 - p2) - lens[[2]]^2,
   (p4 - p3).(p4 - p3) - lens[[3]]^2,
   p4.p4 - lens[[4]]^2, 
   area^2 - 1/16 (4*psq*qsq - (a^2 + c^2 - b^2 - d^2)^2)};

sol = Solve[polys == 0];
Select[sol, (y3 /. #) >= 0 &] // N

(* Out[198]= {{y4 -> -1.88729833462, x4 -> 0.661895003862, 
  y3 -> 0.536554646919, x3 -> 3.8438655763}, {y4 -> 1.11270166538, 
  x4 -> -1.66189500386, y3 -> 0.733275958395, x3 -> 2.32006884993}} *)

--- edit --- Motivated by the response from @Michael E2, I packaged this as a Module and used NSolve (which, under the hood, also computes a Groebner basis over InexactNumbers).

quad[lens_, area_, x_, y_] := Module[
  {p1, p2, p3, p4, a, b, c, d, psq, qsq, polys, sol},
  p1 = {0, 0};
  p2 = {First[lens], 0};
  p3 = {x[3], y[3]};
  p4 = {x[4], y[4]};
  {a, b, c, d} = lens;
  psq = (p3 - p1).(p3 - p1);
  qsq = (p4 - p2).(p4 - p2);
  polys = {(p3 - p2).(p3 - p2) - lens[[2]]^2,
    (p4 - p3).(p4 - p3) - lens[[3]]^2, p4.p4 - lens[[4]]^2,
    area^2 - 1/16 (4*psq*qsq - (a^2 + c^2 - b^2 - d^2)^2)};
  sol = NSolve[polys];
  Select[sol, With[{y3 = y[3] /. #}, Im[y3] == 0 && y3 >= 0] &]
  ]

I cannot say much about the relative speeds since I have not succeeded in getting both methods to work. My choice of an area formula might lead to suboptimal performance. I show the example from the @Michael E2 response; one of my answers is the same as the one indicated therein.

qsols = quad[{2, 3, 1, 2}, Pi, x, y]

Out[8]= {{y[4] -> 1.72915590151, x[4] -> -1.0049974469, 
  y[3] -> 2.14622307243, 
  x[3] -> -0.0961217816172}, {y[4] -> 1.95571758855, 
  x[4] -> 0.418531616287, y[3] -> 2.80744298296, 
  x[3] -> 0.942520025049}}

--- end edit ---

share|improve this answer
    
I added some timings to my answer. Please verify their relative correctness. Thanks. –  Michael E2 Jun 14 '13 at 1:58

I had roughly the same idea as Daniel Lichtblau, but with a playful twist that turned out to improve the speed. (I was surprised.) So on that basis I feel it's worth sharing.

We'll set up some polynomials that will define the desired quadrilateral. I made the origin an explicit argument, even though later the quadrilateral function will use it only as {0, 0}. The Total[..] line creates equations that specify the distance of three sides to be s2, s3,s4. (The first side from the origin topwill be preassigned top = {s1, 0}`. One could add the equation here, but we have the freedom to determine the position of one side of the quadrilateral.) To these are added the equation for the area.

area[orig_, p_, q_, r_] := Det[{p - orig, r - orig}]/2 + Det[{r - q, p - q}]/2;

quadIdeal[s1_, s2_, s3_, s4_, A_, orig_, p_, q_, r_] := 
  Join[
   Total[Subtract @@@ Partition[{p, q, r, orig}, 2, 1]^2, {2}] - {s2, s3, s4}^2,
   {area[orig, p, q, r] - A}];

We use GroebnerBasis to reduce these polynomials. It puts the polynomials in an order such that their roots can be solved sequentially -- a nice setup for Fold. We apply Reduce and convert the equations it outputs to Set, to assign the value to the coordinate variable in the equation. In GroebnerBasis we set CoefficientDomain -> InexactNumbers, so that any numeric arguments may be processed. One caveat: if the area is too great, the figure is impossible and complex roots will be returned.

quadrilateral[s1_?NumericQ, s2_?NumericQ, s3_?NumericQ, s4_?NumericQ, A_?NumericQ] := 
  Polygon @
   Block[{x2, y2, x3, y3},
    Fold[
      Reduce[#2 == 0, Variables[#2]] /. Equal -> Set &,
      0,
      GroebnerBasis[#, Reverse@Variables[#], CoefficientDomain -> InexactNumbers] &@
        quadIdeal[s1, s2, s3, s4, A, {0, 0}, {s1, 0}, {x2, y2}, {x3, y3}]];
    {{0, 0}, {s1, 0}, {x2, y2}, {x3, y3}}
    ];

Example -- "circle" (π) quadrature :-)

N @ quadrilateral[2, 3, 1, 2, Pi]
Polygon[{{0., 0.}, {2., 0.}, {0.94252, 2.80744}, {0.418532, 1.95572}}]
Graphics[quadrilateral[2, 3, 1, 2, Pi], Frame -> True]

Polygon

Surprisingly using Fold etc. this way is faster than Solve or NSolve. Edit - updated timings. On a fresh kernel everything ran a little faster today, so I'll report them with the other updates. Also GroebnerBasis with InexactNumbers uses very high precision arithmetic apparently (perhaps a hundred digits). One can save a little if at least one machine precision number is passed to quadrilateral:

Do[quadrilateral[2, 3, 1, 3, N @ Pi], {100}] // Timing // First

Do[quadrilateral[2, 3, 1, 3, Pi], {100}] // Timing // First

Do[Solve[quadIdeal[2, 3, 1, 3, Pi, {0, 0}, {2, 0}, {x2, y2}, {x3, y3}] == 0,
         {x2, y2, x3, y3}], {100}] // Timing // First

Do[NSolve[quadIdeal[2, 3, 1, 3, Pi, {0, 0}, {2, 0}, {x2, y2}, {x3, y3}] == 0,
          {x2, y2, x3, y3}], {100}] // Timing // First
0.202163
0.220306
1.774986
1.049178

Using @DanielLichtblau's quad (also a little faster with machine precision input):

Do[quad[{2, 3, 1, 2}, N @ Pi, x, y], {100}] // Timing // First
Do[quad[{2, 3, 1, 2}, Pi, x, y], {100}] // Timing // First
1.548769
1.581755

To see how it works, let's examine the steps. Here are the initial polynomials.

quadIdeal[1, 4, 2, 5, 6, {0, 0}, {x1, 0}, {x2, y2}, {x3, y3}]
{-16 + (1 - x2)^2 + y2^2,  -4 + (x2 - x3)^2 + (y2 - y3)^2,
 -25 + x3^2 + y3^2,  -6 + y3/2 + 1/2 (y2 - x3 y2 - y3 + x2 y3)}

GroebnerBasis reduces the set to a simpler one that has the same roots. Because of the nature of this problem and the way GroebnerBasis works, it turns out to be easy to solve for the roots.

GroebnerBasis[%, Reverse@Variables[%]]
{21 - 26 x2 + 5 x2^2,  33 - 20 x2 + 17 x3, -9 + x2 + 2 y2,  -89 + 5 x2 + 17 y3}

The root of the first polynomial yields x2; here Reduce returns two equations, and when they are converted to Set, the second assignment overrides the first. The second polynomial, with x2 being found, yields x3. Finally y2 and y3 are set, again x2 having been set earlier.

share|improve this answer
    
(Not to display too many sour grapes, but...) is it possible you have a typo somewhere? In a fresh session I get: In[4]:= N@quadrilateral[2, 3, 1, 2, Pi] During evaluation of In[4]:= Reduce::ratnz: Reduce was unable to solve the system with inexact coefficients. The answer was obtained by solving a corresponding exact system and numericizing the result. >> Out[4]= Polygon[{{0., 0.}, {2.14331906606 + 2.93152121243 I, 0.}, {-0.355115569688 + 2.70725977699 I, 1.7075673539 - 0.328129099269 I}, {-1.14934382499 - 0.693069787178 I, -1.82998174396 + 0.435291490098 I}}] –  Daniel Lichtblau Jun 13 '13 at 18:22
    
@DanielLichtblau Yes there was a typo (an x1 that should have been an s1). Fixed now. Sorry for the trouble. –  Michael E2 Jun 14 '13 at 1:42
    
Very nice now. Offhand I am not sure what accounts for the relative speed gain. It could be any of a few things. (1) A lex basis might be faster to find than degree-reverse-lex, for this particular type of system. (2) The one-at-a-time solving might be faster, since it discards all but one root as it goes. (3) This form of back-solving might be generally faster than what NSolve does, for this type of system but for other reasons unknown. Interesting. –  Daniel Lichtblau Jun 14 '13 at 3:49

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.