有没有比使用Sympy的nroot（）更快的方法来找到四阶多项式的根？

Question

我想解一个复系数的四阶多项式

-0.678916793992528*w^4 + 9207096.65180878*i*w^3 
+ 1.47445911372677e+15*w^2 - 1.54212540689566e+21*i*w 
+ 2.70530138119032e+26

这段代码的最终目标将是解决这个多项式至少100000次，每次不同的系数，所以我希望代码是快速和有效的。我一直在用sympy.nroots公司（）获得根，但根据%时间，每个循环大约需要9.6毫秒，与其他循环相比非常缓慢根（）每个循环需要60µs。但是我不能用根（）因为它不能很好地处理复系数，并且一直不正确地求解这个多项式的根。使用sympy.solve公司（）的速度更慢，每环路122毫秒。你知道吗

为了加速这个过程，我想到的一件事是，我实际上只需要根的虚部，特别是最负的虚部，但我不确定是否可以利用它来加快代码的运行时间。你知道吗

我的问题是，有没有其他的函数可以用来查找根更快？或者有没有一种我可以编写的更快的根查找方法？最后，有没有一种方法可以只解复数根，这样会更快？你知道吗

Answer 1

在双精度浮点数中，没有比np.root更好的结果了。计算一个接近根的多项式需要很多灾难性的抵消。你知道吗

使用numpy的例程来尝试您的示例，得到如下根

def print_Carr(z):
    for zz in z: print(">>> % 22.17e %+.17ej"%(zz.real, zz.imag))

p=np.array([-0.678916793992528, 9207096.65180878j, 1.47445911372677e+15, -1.54212540689566e+21j, 2.70530138119032e+26])
z=np.roots(p); print_Carr(z)
>>>  4.60399640251209885e+07 +6.25409784852022864e+06j
>>> -4.60399640251209214e+07 +6.25409784852025378e+06j
>>>  6.97016694994478896e-13 +1.20627308238215139e+06j
>>>  5.23825344503222243e-11 -1.53018048966713541e+05j

这些值对于多项式求值来说相当大。这些根的评估值为

print_Carr(np.polyval(p,z))
>>> -3.48222204464332800e+15 +2.82412997568102400e+15j
>>>  5.73769835033395200e+15 -1.64254152287846400e+15j
>>> -4.12316860416000000e+11 +1.37984933104284096e+09j
>>>  6.87194767360000000e+10 -1.04451799855962357e+11j

对于残差来说，这看起来相当糟糕，然而根尾数的最后一位的变化引入了值的巨大绝对变化。记住，精确的根（对于给定的系数）在浮点数之间。这些变化对多项式值的影响可以通过用其绝对值替换系数和根来估计，因为mu*|p|(|z|)是对浮点计算误差的估计。你知道吗

print_Carr(np.polyval(abs(p),abs(z)) *2**-52)
>>>  1.63036010254646300e+15 +0.00000000000000000e+00j
>>>  1.63036010254645625e+15 +0.00000000000000000e+00j
>>>  9.53421868314746094e+11 +0.00000000000000000e+00j
>>>  1.20139515277909210e+11 +0.00000000000000000e+00j

残差几乎在这些界限的范围内。你知道吗

改变根近似的最后尾数位或多项式系数具有可以通过根位置处的导数来估计的影响

print_Carr(abs(np.polyval(np.polyder(p),z))*(2**-52*abs(z)))
>>>  1.38853576300226150e+15 +0.00000000000000000e+00j
>>>  1.38853576300225050e+15 +0.00000000000000000e+00j
>>>  5.30242273857438416e+11 +0.00000000000000000e+00j
>>>  6.77504690635207825e+10 +0.00000000000000000e+00j

再次证明，任何超过最后两个尾数位的变化都会大幅增加残差。你知道吗

为了消除np.roots实现中“伴随矩阵的特征值”可能存在的不精确性，采用牛顿法的一步“根抛光”并重新计算残差

z = z - np.polyval(p,z)/np.polyval(np.polyder(p),z); print_Carr(z)
>>>  4.60399640251209661e+07 +6.25409784852025565e+06j
>>> -4.60399640251209661e+07 +6.25409784852025472e+06j
>>>  1.00974195868289511e-28 +1.20627308238215116e+06j
>>>  0.00000000000000000e+00 -1.53018048966713570e+05j
print_Carr(np.polyval(p,z))
>>>  6.74825261547520000e+13 -7.41139556597760000e+13j
>>>  1.55993212190720000e+13 -1.15513145425920000e+14j
>>>  2.74877906944000000e+11 +1.99893600285358499e-07j
>>>  0.00000000000000000e+00 +0.00000000000000000e+00j

实际上，残差减少了一个或两个小数位，这表明使用这种浮点数据类型几乎可以达到最佳效果。你知道吗

因此，针对您的任务的新建议是使用numpy.roots和一个牛顿步骤进行根部抛光。你知道吗

---

最后与另一个多精度结果进行了比较

from mpmath import mp
mp.dps = 20; mp.pretty = True;
mp.polyroots(p, maxsteps=20, extraprec=30) # prec=bits, dps=digits, 10bits=3digits
>>> [(0.0 - 153018.04896671356797j),
>>>  (0.0 + 1206273.0823821511478j),
>>>  (-46039964.025120967306 + 6254097.8485202553318j),
>>>  ( 46039964.025120967306 + 6254097.8485202553318j)]

当以实部和虚部相同的方式计算位置时，根+牛顿结果在15个前导数字中是正确的。你知道吗