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.

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) *)
share|improve this question
    
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
1  
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
add comment

2 Answers 2

up vote 3 down vote accepted

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

share|improve this answer
add comment

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]}}
share|improve this answer
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.