含有浮点数的变系数多项式的Sympy判别式



我想知道当两个椭圆的"半径"增加时,它们第一次相交的时间。在那一刻,它们应该是外部的余切。我们有两个函数f和g,它们的形式都是a x ^2+b y ^2+c xy+d x+f y+(e-t(。注意常数项中的变量t。f、 g(x,y(=0描述了"半径"t的第一个、第二个椭圆。Contangent曲线的意思是复根。我的方法是:取消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,但还是没什么成功。

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

我最终按照建议做了。通过在四次方中插入一个特定的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,有时我会因为使用空列表而出错。那么,二进制搜索方法似乎是可行的。

最新更新