有没有一种方法可以比使用 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

此代码的最终目标是至少求解此多项式 100,000 次，每次使用不同的系数，因此我希望代码快速高效。我一直在使用 sympy.nroots() 来获取根，但根据 %timeit，每个循环大约需要 9.6 毫秒，与每个循环需要 60 微秒的 numpy.roots() 相比，这是非常慢的。但是我不能使用 numpy.roots() 因为它不能很好地处理复数系数并且一直错误地求解这个多项式的根。使用 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

再次证明除了最后两个尾数位之外的任何变化都会大大增加残差。

为了消除"eigenvalues of the companion matrix"在实现np.roots时可能存在的不精确性，应用"root polishing"一步牛顿法重新计算残差，

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位是正确的，实部和虚部计算位置的方式相同。

有没有一种方法可以比使用 Sympy 的 nroot() 更快地找到四阶多项式的根？

Is there a way to find the roots of a 4th order polynomial that is faster than using Sympy's nroot()?

python

numpy

sympy

numerical-methods