solving simultaneous quadratics gives inconsistent results

[Russell Young]

I'm a newbie with sympy, so this may be something I am doing wrong.

I've written code that finds the intersection of 2 hyperbolas. The equations are correct - when I plot them they are as expected. Using test data where I know there is a solution, sometimes it finds it, sometimes it doesn't. I've found and fixed a few issues so far: no division of python numbers, normalizing the graphs so no numbers are too small, adding a .simplify() call to the equations being solved. Still, it ends up missing solutions. The order the equations are passed to the solve function also affects the solutions found.

Here is some output that describes the hyperbolas, shows the equations, and gives the results:

$ python hyperbola.py -1 0 10.033445 1 1 4.729811 0 -1 7.478488
h12
  a = 0.792893, b = 0.788239, c = 1.118034
  foci are Point(-1, 0) and Point(1, 1)
  center is Point(0, 1/2)
  equation is -1.60947597191687*(-sqrt(5)*x/5 + 2*sqrt(5)*(y - 1/2)/5)**2 + 1.59063495669236*(2*sqrt(5)*x/5 + sqrt(5)*(y - 1/2)/5)**2 == 1
h13
  a = 0.381966, b = 0.595065, c = 0.707107
  foci are Point(-1, 0) and Point(0, -1)
  center is Point(-1/2, -1/2)
  equation is 6.8540998039692*(sqrt(2)*(x + 1/2)/2 - sqrt(2)*(y + 1/2)/2)**2 - 2.82404568540787*(sqrt(2)*(x + 1/2)/2 + sqrt(2)*(y + 1/2)/2)**2 == 1
h23
  a = 0.410927, b = 1.039778, c = 1.118034
  foci are Point(1, 1) and Point(0, -1)
  center is Point(1/2, 0)
  equation is -0.924950593916484*(-sqrt(5)*y/5 - 2*sqrt(5)*(x - 1/2)/5)**2 + 5.92202447335121*(2*sqrt(5)*y/5 - sqrt(5)*(x - 1/2)/5)**2 == 1
h12 and h13: [(2.00000069353477, -1.05724449909221e-7)]
h13 and h13: []

h12 and h23: [(-0.872442811259504, -1.67266454797664), (-0.671975994196070, 0.0582722523367054), (0.819499463572107, 0.688825355380073), (2.00000032674867, 4.62032762769592e-8)]
h23 and h12: [(2.00000032674867, 4.62032762919501e-8)]

h13 and h23: []
h23 and h13: []
$

Parameters for these 3 equations were chosen so that the point (2, 0) is nearly on all of them. It is found in one of the orderings of the first pair, in both of the second pair, and in neither of the third pair.

Following is an excerpt of the code that makes and solves the equation:

_x, _y, _a, _b = symbols('x, y, a, b')
_standard = Eq(_x**2/_a**2 - _y**2/_b**2, 1)
_degenerate = Eq(_x, 0)


class Hyperbola(object):
    def __init__(self, f0, f1, diff):
        self.foci = [Point(f0), Point(f1)]
        self.a = S(diff)/2
        self.c = self.foci[0].distance(self.foci[1])/2
        try:
            self.b = sqrt(self.c**2 - self.a**2)
        except ValueError:
            print 'Computed minor axis for %s %s %f is imaginary' % (str(f0), str(f1), diff)
            exit(1)
        theta, x_temp, y_temp = symbols('theta, x_temp, y_temp')
        # make equation in standard form
        self.equation = _standard if self.a else _degenerate
        self.equation = self.equation.subs({_a: self.a, _b: self.b})
        # rotate it to the proper orientation
        self.equation = self.equation.subs({_x: x_temp*cos(theta) + y_temp*sin(theta),
                                            _y: -x_temp*sin(theta) + y_temp*cos(theta)})\
                                            .subs({x_temp: _x, y_temp: _y}).\
                                            subs({theta: self.angle(*self.foci)})
        # translate to the right position
        center = self.center()
        self.equation = self.equation.subs({_x: _x - center.x, _y: _y - center.y})

    def solve(self, other):
        # def is_real(x):
        #     return not x[0].is_imaginary and not x[1].is_imaginary
        # return [x for x in solve_poly_system([self.equation.simplify(), other.equation.simplify()], [_x, _y]) or []
        #         if is_real(x)]
        return [x for x in solve_poly_system([self.equation.simplify(),
                                              other.equation.simplify()], [_x, _y]) or []]

Am I doing something wrong not to get all the roots? Is there a way to get reliably the intersections of these curves, or is sympy not the right tool for this job? If it isn't, are there any suggestions for what might be better?

Thanks,

russell

[Chris Smith]

Did you just try using solve? Sympifying the equations you give for h13 and h23 gives the following result for me:

>>> solve((h23,h13),x,y)

[(-0.557666440391888 + 4.91385686350788e-31*I, -1.35411398328542 - 2.27087209183

879e-32*I), (-0.555228172386711 + 2.35128780795143e-29*I, 0.0820809601337735 - 1

.08661602193259e-30*I), (0.871420006337371 + 4.47416209082826e-30*I, -0.28263412

4887546 - 2.06767454992846e-31*I), (2.00000355907055 - 2.84784242544466e-29*I, 5

.70736401553117e-7 + 6.42623143478428e-34*I)]

/c

...