Eliminate several variables between five Conic_section equations

Question

I want to eliminate x1, x2, y1, y2 between these 5 equations. I tried

Eliminate[{x1^2/a^2 + y1^2/b^2 == 1, x2^2/a^2 + y2^2/b^2 == 1, 
  x1 x2 + y1 y2 == 0, y2/(x2 - c) == y1/(x1 - c) == k}, {x1, x2, y1, y2}]

Mathematica takes a long time, I couldn't even wait for it to finish.Can you recommend an efficient method?

Updated:

http://en.wikipedia.org/wiki/Proportionality_(mathematics)

ForAll[{a, b}, b d != 0 && b^2 != d^2, 
 Equivalent[a/b == c/d, 
  a/b == c/d == (a - c)/(b - d) == (a + c)/(b + d)]]
Resolve[%]
(*True*)

  k^2 /. Solve[
     Eliminate[{x1^2/a^2 + y1^2/b^2 == 1, x2^2/a^2 + y2^2/b^2 == 1, 
     x1 x2 + y1 y2 == 0, y2/(x2 - c) == (y1 - y2)/(x1 - x2) == k}, {x1, x2, y1, y2}], 
   k] // DeleteDuplicates//Tr

(* (a^2 b^2)/(-a^2 b^2 + a^2 c^2 + b^2 c^2) *)

Answer 1

Can use GroebnerBasis for this.

eqs = {x1^2/a^2 + y1^2/b^2 == 1, x2^2/a^2 + y2^2/b^2 == 1, 
   x1 x2 + y1 y2 == 0, y2/(x2 - c) == y1/(x1 - c) == k};

First take differences to make expressions from the equations.

exprs = Flatten[Map[#[[1]] - Rest[#] &, Apply[List, eqs, 1]]]

(* Out[81]= {-1 + x1^2/a^2 + y1^2/b^2, -1 + x2^2/a^2 + y2^2/b^2, 
 x1 x2 + y1 y2, -(y1/(-c + x1)) + y2/(-c + x2), -k + y2/(-c + x2)} *)

Now get common denominators and extract numerators, to get polynomials.

polys = Numerator[Together[exprs]]

(* Out[82]= {-a^2 b^2 + b^2 x1^2 + a^2 y1^2, -a^2 b^2 + b^2 x2^2 + 
  a^2 y2^2, x1 x2 + y1 y2, 
 c y1 - x2 y1 - c y2 + x1 y2, -c k + k x2 - y2} *)

We now find the variables to keep.

elims = {x1, x2, y1, y2};
vars = Complement[Variables[polys], elims]

(* Out[84]= {a, b, c, k} *)

Now find the polynomial resulting from eliminating the x and y variables.

Timing[
 gb = GroebnerBasis[polys, vars, elims, 
   MonomialOrder -> EliminationOrder]]

(* Out[85]= {0.100000, {a^8 b^8 - a^6 b^8 c^2 + 3 a^8 b^8 k^2 - 
   3 a^8 b^6 c^2 k^2 - 2 a^6 b^8 c^2 k^2 + 3 a^6 b^6 c^4 k^2 - 
   a^4 b^8 c^4 k^2 + 3 a^8 b^8 k^4 - 2 a^8 b^6 c^2 k^4 - 
   5 a^6 b^8 c^2 k^4 + 3 a^8 b^4 c^4 k^4 + a^4 b^8 c^4 k^4 - 
   3 a^6 b^4 c^6 k^4 + 2 a^4 b^6 c^6 k^4 + a^2 b^8 c^6 k^4 + 
   a^8 b^8 k^6 + a^8 b^6 c^2 k^6 - 4 a^6 b^8 c^2 k^6 - 
   a^8 b^4 c^4 k^6 - 3 a^6 b^6 c^4 k^6 + 6 a^4 b^8 c^4 k^6 - 
   a^8 b^2 c^6 k^6 + 2 a^6 b^4 c^6 k^6 + 3 a^4 b^6 c^6 k^6 - 
   4 a^2 b^8 c^6 k^6 + a^6 b^2 c^8 k^6 - a^4 b^4 c^8 k^6 - 
   a^2 b^6 c^8 k^6 + b^8 c^8 k^6}} *)

Thus far I did not try to force denominators not to vanish. One could do that by creating a new polynomial of the form product(denoms)*newvariable-1, adding that to the polynomial list, and adding newvariable to the list of elimination variables. We'll do this below.

denoms = 
 Apply[Times, 
  FactorSquareFreeList[Times @@ Denominator[Together[exprs]]][[All, 
    1]]]

(* Out[93]= a b (c - x1) (c - x2) *)

newpoly = denoms*newvar - 1

(* Out[94]= -1 + a b newvar (c - x1) (c - x2) *)

This will be slower...

Timing[
 gb2 = GroebnerBasis[Append[polys, newpoly], vars, 
   Append[elims, newvar], MonomialOrder -> EliminationOrder]]

(* Out[96]= {27.900000, {a^6 b^6 + 3 a^6 b^6 k^2 - 3 a^6 b^4 c^2 k^2 + 
   a^4 b^6 c^2 k^2 + 3 a^6 b^6 k^4 - 2 a^6 b^4 c^2 k^4 - 
   2 a^4 b^6 c^2 k^4 + 3 a^6 b^2 c^4 k^4 - 2 a^4 b^4 c^4 k^4 - 
   a^2 b^6 c^4 k^4 + a^6 b^6 k^6 + a^6 b^4 c^2 k^6 - 
   3 a^4 b^6 c^2 k^6 - a^6 b^2 c^4 k^6 - 2 a^4 b^4 c^4 k^6 + 
   3 a^2 b^6 c^4 k^6 - a^6 c^6 k^6 + a^4 b^2 c^6 k^6 + 
   a^2 b^4 c^6 k^6 - b^6 c^6 k^6}} *)

So the "correct" result is actually of lower degree. Indeed, the new one is a proper factor of the old.

Factor[gb[[1]]/gb2[[1]]]

(* Out[98]= b^2 (a - c) (a + c) *)

(I guess the pros call it a "quotient ideal" for a reason.)

--- edit ---

Alternatively, since it is a principle ideal, I guess I should have just factored the generating polynomial for gb and removed all factors that divide any denominator. That would also have given gb2 and at much lower computational cost.

--- end edit ---

Answer 2

This is a matter of how you define your assumptions probably. As @b.gatessucks noted, you have a carbon copy of the same set of two equations:

$ \frac{x^2}{a^2} + \frac{y^2}{b^2} = 1 $ and $ \frac{y}{x-c}=k $

which yield two solutions each, call them $ (x_1, y_1), (x_2, y_2). $

Solve[{x^2/a^2 + y^2/b^2 == 1, y/(x - c) == k}, {x, y}];
{{x1, y1}, {x2, y2}}/.%

Mind you, overall, you have four solutions now - a pair for each carbon copy of equations. Then you want to couple these two sets using the third equation that @b.gatessucks commented out $ x_{FIRST} x_{SECOND} + y_{FIRST} y_{SECOND} = 0. $

So there are three possible combinations of these:

$ x_{i}^2 + y_{i}^2 = 0,\;\;i=1,2 $

or

$ x_1 x_2 + y_1 y_2 = 0 $

So, overall, solving the former for $ k $ you have:

(Solve[x1^2 + y1^2 == 0, k] // FullSimplify)

which gives:

{{k -> -Sqrt[((a^4 b^4)/(-a^4 b^4 + a^2 (a - b) b^2 (a + b) c^2 + 2 Sqrt[a^6 b^6 (-a^2 + b^2) c^2]))]}, 
 {k -> Sqrt[(a^4 b^4)/(-a^4 b^4 + a^2 (a - b) b^2 (a + b) c^2 + 2 Sqrt[a^6 b^6 (-a^2 + b^2) c^2])]}, 
 {k -> -Sqrt[-((a^4 b^4)/(a^2 b^4 c^2 + 2 Sqrt[a^6 b^6 (-a^2 + b^2) c^2] + a^4 b^2 (b - c) (b + c)))]}, 
 {k -> Sqrt[-((a^4 b^4)/(a^2 b^4 c^2 + 2 Sqrt[a^6 b^6 (-a^2 + b^2) c^2] + a^4 b^2 (b - c) (b + c)))]}}

this is the same as for $ i=2 $ (which you can check):

(Solve[x1^2 + y1^2 == 0, k] // FullSimplify) == (Solve[x2^2 + y2^2 == 0, k] // FullSimplify)

and also you get the third combination (that you quote):

Solve[x1 x2 + y1 y2 == 0, k] // FullSimplify

which gives

{{k -> -((a b)/Sqrt[-a^2 b^2 + (a^2 + b^2) c^2])}, {k -> (a b)/Sqrt[-a^2 b^2 + (a^2 + b^2) c^2]}}