CODEWITHSUNDEEP
pythonsympy
Luke Kiernan
Luke Kiernan

Reputation: 1

Sympy discriminant of polynomial with variable coefficients involving floats

I wish to figure out when two ellipses first intersect, as their 'radii' increase. They should be externally cotangent at that instant. We have two functions f and g, both of the form a x^2+b y^2+ c xy +d x+f y+(e-t). Note the variable t in the constant term. f, g(x,y) = 0 describes the 1st, 2nd ellipse of 'radius' t. Contangent curves means double root. My approach: cancel the y^2 terms and then solve for y as a rational function of x. Substitute this into f, and clear denominators: this yields a quartic polynomial in x with coefficients involving t (of which x is a zero). Now I should be able to use discriminants to figure out the t values for which that polynomial has a double root, recovering the smallest t for which the ellipses are cotangent.

I tried a simple example in a sympy live shell, where the a,b etc involved in f and g were integers. After some trial and error, I got the following to work:

f = 2*x**2+y**2-t
g = 3*(1/sqrt(2)*(x-1)+1/sqrt(2)*(y+1))**2+5*(1/sqrt(2)*(x-1)-1/sqrt(2)*(y+1))**2-t
h = expand(g).coeff(y,2)*f - expand(f).coeff(y,2)*g
h = Poly(h)
f2=Poly((f.subs(y, solve(h,y)[0])*h.as_expr().coeff(y,1)**2).simplify(),x)
roots = roots(discriminant(f2),t)

After discarding the complex roots, I had what I wanted. How do I get this to work with non-integer inputs?

My attempts: the first issue is that multiplying by h.as_expr().coeff(y,1)**2 fails to cancel the denominator in f.subs(y, solve(h,y)[0]) presumably due to rounding error. I crafted a workaround via extracting the numerator and denominator of the expression separately, replacing substitution with a coefficient-by-coefficient operation:

num = -1*h.as_expr().coeff(y,0)
denom = h.as_expr().coeff(y,1)
noYTerm = f.as_expr().coeff(y,0)
linearYTerm = f.as_expr().coeff(y,1)
ySquareTerm = f.as_expr().coeff(y,2)
f2 = Poly(simplify(ySquareTerm*num**2+linearYTerm*num*denom+noYTerm*denom**2),x)

However, I just get a bunch of sympy.polys.polyerrors.PolynomialDivisionFailed errors. If I print off the expression for f2, I get Poly(10108.448*x**3 + 155858.976*x**2 + (707853.024000001 - 949.184000000001*t)*x + 918397.152000001 - 6870.336*t, x, domain='RR[t]') The .000001 screams floating point error...I've tried switching to QQ, but still little success.

I'm also open to other approaches: perhaps there's an easier way to extract the information I want from this pair of polynomials than what I'm doing.

Upvotes: 0

Views: 767

Answers (1)

Luke Kiernan
Luke Kiernan

Reputation: 1

I ended up doing as suggested. By plugging a particular t into the quartic and checking for positive real roots, you can see if the ellipses intersect at a certain size. I then did a binary search to find the smallest t value at which the ellipses first intersect.

Later, I realized an inelegant workaround: there's an explicit expression for the quartic discriminant, at https://mathworld.wolfram.com/PolynomialDiscriminant.html. So I spent 10 minutes typing that into my program. Then I could plug in the coefficients of f_2:

a0 = f2.coeff(x,0) # a1 =...etc.
# formula for quartic discriminant from https://mathworld.wolfram.com/PolynomialDiscriminant.html
discr = Poly( a1**2*a2**2*a3**2 - 4*a1**3*a3**3 - 4*a1**2*a2**3*a4 + 18*a1**3*a2*a3*a4-27*a1**4*a4**2 + 256*a0**3*a4**3
       + a0*(-4*a2**3*a3**2 + 18*a1*a2*a3**3 + 16*a2**4*a4 - 80*a1*a2**2*a3*a4 - 6*a1**2*a3**2*a4 + 144*a1**2*a2*a4**2)
       + a0**2 * (-27*a3**4 + 144* a2* a3**2 * a4 - 128 * a2**2 * a4**2 - 192 * a1 * a3 * a4**2)), t)
roots = np.roots(discr.all_coeffs())
postiveAndReal  = numpy.logical_and(numpy.isreal(roots), numpy.greater(roots, numpy.zeros(roots.shape[0])))
firstIntersectAt = numpy.amin(roots[postiveAndReal])

However, it seems like expanding that massive expression creates rounding error. With this method, if the ellipses were ~0.1 apart, sometimes I'd get errors due to taking amin of an empty list. The binary search method seems to be the way to go, then.

Upvotes: 0

Related Questions