Group Abstract Group Abstract

Message Boards Message Boards

0
|
4.1K Views
|
11 Replies
|
4 Total Likes
View groups...
Share
Share this post:

Existence of Real Roots to Coupled Polynomial Equations

Can Mathematica be used to derive the symbolic conditions for the existence of real roots of a pair of second order polynomial equations, with arbitrary coefficients, in two variables?

POSTED BY: Frank Kampas
11 Replies

Fortunately, the tangent point is not at the origin, as that would be original sin.

POSTED BY: Frank Kampas

This is what happens when misdirected comics discuss intersected conics.

POSTED BY: Daniel Lichtblau

FindInstance works better if you don't pull the trigger.

In[38]:= fi = 
 FindInstance[
   Rest[gb] == 0 && cth^2 + sth^2 == 1, {xc, yc, cth, sth}] // N

Out[38]= {{xc -> -2.70687, yc -> 0., cth -> -1., sth -> 0.}}

In[46]:= ContourPlot[{q1 == 
   0, (q2 /. th -> ArcCos[cth] /. fi[[1]] ) == 0}, {x, -5, 5}, {y, -5,
   5}]

enter image description here

POSTED BY: Frank Kampas

That's useful to know. I guess it's 'cos FindInstance is not able to handle an overabundance of sin.

POSTED BY: Daniel Lichtblau

Yeah, having explicit Sin/Cos terms makes it triggy. I generally work around that by replacements and adding appropriate polynomials to cover the trig identities needed.

q1 = eleq[{5/4, 3/2, 0, 0, \[Pi]/4}][{x, y}];
q2 = eleq[{4/3, 2, xc, yc, th}][{x, y}];

quadratics = {q1, q2};
grad1 = Grad[quadratics[[1]], {x, y}];
grad2 = Grad[quadratics[[2]], {x, y}];
polys = Flatten[{quadratics, grad1 - lambda*grad2, 
    cth^2 + sth^2 - 1}] /. {Cos[th] -> cth, Sin[th] -> sth}
elim = {x, y, lambda};

(* {-1 + 4/9 (-(x/Sqrt[2]) + y/Sqrt[2])^2 + 
  16/25 (x/Sqrt[2] + y/Sqrt[2])^2, -1 + 
  1/4 (-sth (x - xc) + cth (y - yc))^2 + 
  9/16 (cth (x - xc) + sth (y - yc))^2, -(4/9) Sqrt[
   2] (-(x/Sqrt[2]) + y/Sqrt[2]) + 
  16/25 Sqrt[2] (x/Sqrt[2] + y/Sqrt[2]) - 
  lambda (-(1/2) sth (-sth (x - xc) + cth (y - yc)) + 
     9/8 cth (cth (x - xc) + sth (y - yc))), 
 4/9 Sqrt[2] (-(x/Sqrt[2]) + y/Sqrt[2]) + 
  16/25 Sqrt[2] (x/Sqrt[2] + y/Sqrt[2]) - 
  lambda (1/2 cth (-sth (x - xc) + cth (y - yc)) + 
     9/8 sth (cth (x - xc) + sth (y - yc))), -1 + cth^2 + sth^2} *)

Timing[gb = 
   GroebnerBasis[polys, Complement[Variables[polys], elim], elim, 
    MonomialOrder -> EliminationOrder];]

(* Out[111]= {0.277867, Null} *)

Now restore the trigs.

elimpoly = Rest[gb /. {cth -> Cos[th], sth -> Sin[th]}];

This might be usable for e.g. FindInstance. Of course with the trigs present the task is now probably more difficult. (More variables and a larger solution set does not actually make the job easier.)

POSTED BY: Daniel Lichtblau

Daniel, I tried out your approach with two specific ellipses with the center of the second one being the only unknowns.

In[7]:= eleq[{a_, b_, xc_, yc_, \[Theta]_}][{x_, 
    y_}] = (((x/a)^2 + (y/b)^2 - 1 ) /. {x -> 
       Cos[\[Theta]]*x + Sin[\[Theta]] y, 
      y -> -Sin[\[Theta]]*x + Cos[\[Theta]]*y}) /. {x -> x - xc, 
    y -> y - yc};

In[8]:= q1 = eleq[{5/4, 3/2, 0, 0, \[Pi]/4}][{x, y}]

Out[8]= -1 + 4/9 (-(x/Sqrt[2]) + y/Sqrt[2])^2 + 
 16/25 (x/Sqrt[2] + y/Sqrt[2])^2

In[16]:= q2 = eleq[{4/3, 2, xc, yc, \[Pi]/3}][{x, y}]

Out[16]= -1 + 1/4 (-(1/2) Sqrt[3] (x - xc) + (y - yc)/2)^2 + 
 9/16 ((x - xc)/2 + 1/2 Sqrt[3] (y - yc))^2

In[17]:= quadratics = {q1, q2};
grad1 = Grad[quadratics[[1]], {x, y}];
grad2 = Grad[quadratics[[2]], {x, y}];
polys = Flatten[{quadratics, grad1 - lambda*grad2}];
elim = {x, y, lambda};

Timing[gb = 
   GroebnerBasis[polys, Complement[Variables[polys], elim], elim, 
    MonomialOrder -> EliminationOrder];]

Out[22]= {1.70313, Null}

In[27]:= fi = FindInstance[gb == 0, {xc, yc}] // N

Out[27]= {{xc -> -3.12216, yc -> 0.}}

In[28]:= ContourPlot[{q1 == 0, (q2 /. fi[[1]]) == 0}, {x, -5, 
  5}, {y, -5, 5}]

enter image description here

However, when I tried to make the orientation of the second ellipse a variable, calculating the GroebnerBasis took too long. All those Sin[theta] and Cos[theta] terms I guess. I don't think that's tangential to the problem.

POSTED BY: Frank Kampas

Speaking parabolically? That's fine. I had our lines crossed (wouldn't want this to degenerate...)

As for determining when two ellipses collide in say a packing problem, I think the best thing is to use NSolve and check numerically. That gives a fair mix of speed and reliability, to the point where it is probably usable in a numeric optimization as part of the constraint formulation. This assumes there are not too many interior ellipses, admittedly.

POSTED BY: Daniel Lichtblau

There are various ways to get information about real root constraints. I will show one that is more or less feasible, although the result might be too large to be of use and even so the result is incomplete.

The idea is to determine where the zeros hit multiple roots (this is a part of what is called,perhaps not surprisingly, the "discriminant variety"). It is at such values that pairs cross to/from two real values to complex conjugates. This happens precisely at points where the ellipses touch tangentially (found with a Lagrange multiplier). In order to make this viable I switched from your formulation as two explicit ellipses to a pair of "generic" guadratics.

quadratics = {a*x^2 + b*x*y + c*y^2 + d*x + e*y + f, 
   p*x^2 + q*x*y + r*y^2 + s*x + t*y + u};
grad1 = Grad[quadratics[[1]], {x, y}];
grad2 = Grad[quadratics[[2]], {x, y}];
polys = Flatten[{quadratics, grad1 - lambda*grad2}];
elim = {x, y, lambda};

Timing[gb = 
   GroebnerBasis[polys, Complement[Variables[polys], elim], elim, 
    MonomialOrder -> EliminationOrder];]

(* {2.377909, Null} *)

It contains one polynomial, and that polynomial (in the various parameters) vanishes precisely where we have the transition points. So the vanishing gives the condition for where these crossings happen. I will note that the polynomial has a leafcount in the tens of thousands.

Another way to approach this might be to mess with resultants and discriminants. A resultant of the two quadratics, say with respect to x, eliminates x and gives a (parametrized) univariate in y. The discriminant of this new polynomial, with respect to y of course, should then also provide the vanishing conditions. But this step is slow, most likely because it also spawns what are called spurious factors.

If one really wants inequality conditions in the parameter space that divide the regions of zero, two, or four real solutions (counting multiplicity), that requires CylindricalDecomposition. I doubt any formulation will run in reasonable time or produce a result of less than too-stupendous-to-use size.

Did I really chase you to ellipses, or was that just, err, hyperbole?

POSTED BY: Daniel Lichtblau

Interesting. I've been using embedded Lagrange multipliers to pack ellipses into circles, but the method breaks down when I tried to pack ellipses into ellipses due to multiple solutions for the Lagrange multiplier equations. I've been playing with GroebnerBasis but your knowledge of their applications is much more extensive than mine.

As you suspect, the problem is too large for CylindricalDecomposition.

I was speaking parabolically. Some time ago I sent you the symbolic solution for the circumscribing radius for packing circles with radii 1, sqrt(1/2) and sqrt(1/3) and you commented on my obsession with inverse square root radii circle packing.

POSTED BY: Frank Kampas

My example is (what a surprise) finding out if two ellipses overlap. As before, my equation for an ellipse with semi-major axes, a and b, centered at {xc ,yc} and rotated by an angle theta is eleq[{a,b,xc,yc,theta}][{x,y}] == 0 where

eleq[{a_, b_, xc_, yc_, \[Theta]_}][{x_, 
   y_}] = (((x/a)^2 + (y/b)^2 - 1 ) /. {x -> 
      Cos[\[Theta]]*x + Sin[\[Theta]] y, 
     y -> -Sin[\[Theta]]*x + Cos[\[Theta]]*y}) /. {x -> x - xc, 
   y -> y - yc}

-1 + ((y - yc) Cos[\[Theta]] - (x - xc) Sin[\[Theta]])^2/b^2 + ((x - 
      xc) Cos[\[Theta]] + (y - yc) Sin[\[Theta]])^2/a^2

If the two ellipses do not overlap, all 4 solutions to { eleq[{a1,b1,xc1, yc1, theta1}][{x,y}] == 0, eleq[{a2,b2,xc2,yc2,theta2}]{x,y}]==0} will be complex. What I've done so far is to ComplexExpand eleq to get two equations in the real and imaginary parts of x and y and test that they give the same solutions for pairs of ellipses as the original equations. Not sure how to proceed from there. Daniel, your concern about my obsession with packing circles has caused me to move on to ellipses, so it would be good if you can help out.

Attachments:
POSTED BY: Frank Kampas
POSTED BY: Daniel Lichtblau
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard