# Eliminate several variables between five Conic_section equations

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) *)

-
Maybe you have too many equations ? Try this for instance : Solve[{x1^2/a^2 + y1^2/b^2 == 1, x2^2/a^2 + y2^2/b^2 == 1(*,x1 x2+ y1 y2\[Equal]0*), y2/(x2 - c) == k, y1/(x1 - c) == k}, {x1, x2, y1, y2}]. – b.gatessucks Feb 7 '13 at 8:23
Alternatively you might want to solve for more variables, say Solve[{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) == k, y1/(x1 - c) == k}, {x1, x2, y1, y2, k}]. – b.gatessucks Feb 7 '13 at 8:25
Thank you. But I expected the result is k^2 == (a^2 b^2)/(-a^2 b^2 + (a^2 + b^2) c^2). – chyaong Feb 7 '13 at 8:50
There are other solutions as well as k^2==... – Daniel Lichtblau Feb 7 '13 at 19:35
The solutions find configurations of ellipses and lines where the two intersection points and the center of the ellipse form a right angle at the center. There's a much nicer way to parameterize the ellipses for this purpose: use $(x/s)^2+(y/t)^2-2\cos(2\theta)xy/(st)=1$. These are centered at the origin and contain the points $(s,0)$ and $(0,t)$, which obviously form a solution--and all possible solutions are in this form. If desired, rotate this configuration back to make the ellipse parallel with the coordinate axes: now you have all solutions to the original equations. – whuber Feb 7 '13 at 20:50

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 ---

-

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]}}

-