涉及浮点数的可变系数多项式的 Sympy 判别式答案

【问题标题】：Sympy discriminant of polynomial with variable coefficients involving floats涉及浮点数的可变系数多项式的 Sympy 判别式
【发布时间】：2021-08-06 13:02:29
【问题描述】：

我想弄清楚两个椭圆何时第一次相交，因为它们的“半径”增加了。在那一刻，它们应该是外部余切的。我们有两个函数 f 和 g，它们的形式都是 a x^2+b y^2+ c xy +d x+f y+(e-t)。注意常数项中的变量 t。 f, g(x,y) = 0 描述了“半径”t 的第一个、第二个椭圆。正切曲线意味着双根。我的方法：取消 y^2 项，然后将 y 作为 x 的有理函数求解。将其代入 f，并清除分母：这会产生 x 中的四次多项式，其系数涉及 t（其中 x 为零）。现在我应该能够使用判别式计算出该多项式具有双根的 t 值，从而恢复椭圆余切的最小 t。

我在 sympy live shell 中尝试了一个简单的示例，其中 f 和 g 中涉及的 a、b 等是整数。经过一些试验和错误，我得到了以下工作：

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)

丢弃复杂的根源后，我得到了我想要的。如何让它与非整数输入一起使用？

我的尝试：第一个问题是乘以 h.as_expr().coeff(y,1)**2 无法取消 f.subs(y, solve(h,y)[0]) 中的分母，大概是由于舍入错误。我通过分别提取表达式的分子和分母制作了一个解决方法，用系数运算替换替换：

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)

但是，我只收到一堆 sympy.polys.polyerrors.PolynomialDivisionFailed 错误。如果我打印出 f2 的表达式，我会得到 Poly(10108.448*x**3 + 155858.976*x**2 + (707853.024000001 - 949.184000000001*t)*x + 918397.152000001 - 6870.336*t, x, domain='RR[t]') .000001尖叫浮点错误...我尝试切换到QQ，但仍然收效甚微。

我也对其他方法持开放态度：也许有一种比我正在做的更容易从这对多项式中提取我想要的信息的方法。

【问题讨论】：

在开始时将系数替换为有理近似值或表示。或者使用数值求解器（来自 scipy.optimize）而不使用 sympy。 Afaik 在标准 python 包中没有实现“混合”求解器方法。
有没有办法让这个工作与浮动？
关于近似结果，请参阅 Kalthofen（和其他人），本质上，这对 Sylvester 矩阵的 SVD 进行了维度分析。这不是在 sympy 中实现的东西。
另外：我看不出数值求解器有什么帮助。我正在求解 t 值，使得 (x,y) 系统具有双根。直到我有一些表达“这些是具有双根的 t 值”（判别式 0）我不明白数值求解器是如何适用的
应该首先用变量作为系数来解决这个问题，然后将结果中的变量替换为它们的浮点值。

标签： python sympy

【解决方案1】：

我最终按照建议做了。通过将特定的 t 插入四次并检查正实根，您可以查看椭圆是否以特定大小相交。然后我进行了二分搜索以找到椭圆第一次相交的最小 t 值。

后来，我意识到了一个不优雅的解决方法：在https://mathworld.wolfram.com/PolynomialDiscriminant.html 处有一个四次判别式的显式表达式。所以我花了 10 分钟把它输入到我的程序中。然后我可以代入 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])

但是，扩展这个庞大的表达式似乎会产生舍入误差。使用这种方法，如果椭圆相距约 0.1，有时我会因为占用一个空列表而出错。那么，二分查找法似乎是可行的方法。

【讨论】：