numpy 多项式根与 Matlab 相比不正确？

Question

我想算出多项式

x^2 - 800x + 160000

的根，它的倍数根为400。

在 Matlab 中，如果我输入

roots([1 -800 160000)

，我会正确得到

ans = 400.0000, 400.0000

，而 numpy

np.roots([1, -800, 160000])

则给出

array([400.+2.69940682e-06j, 400.-2.69940682e-06j])

。

为什么 numpy 显示某种浮点精度问题，而 matlab 却没有？他们不都使用基于伴随矩阵的算法吗？我怎样才能在Python中可靠地确定“实际”根是400，就像Matlab似乎正在做的那样？

编辑：这是一个“玩具”示例，我在一个更大、更复杂的算法中面临，该算法不一定只面临实根和多个根，也就是说，在我更复杂的情况下，我实际上事先并不知道根是多重且实数的。

Answer 1

数值精度是有一定限制的，尤其是涉及到很多运算的时候，比如求根、除法等。在您的示例中，可以使用不同的方法计算根，这反过来会给出略有不同的结果。为了说明这一点，我使用以下示例，其中我不仅更改浮点精度，而且还对同一任务使用两种不同的 numpy 算法。

重要 另请注意，有时浮点数不会以完整的小数精度打印，例如你的 Matlab 示例有 4 位小数。

import numpy as np
print(f'numpy version: {np.__version__}\n')

def print_roots(roots):
    print(f'    roots (round to 5 decimals): {np.round(roots, 5)}')
    print(f'    roots (round to 6 decimals): {np.round(roots, 6)}')
    print(f'    roots (full decimals)      : {roots}')
    print(f'    roots (full decimals real) : {roots.real}')
    print(f'    roots (full decimals imag) : {roots.imag}')
    print()
    

coefficients = np.array([1, -800, 160000])
for i in [np.float32, np.float64]:
    print(f'type: {i}')
    poly = np.polynomial.Polynomial(coefficients.astype(i)[::-1])
    print(f'  np.polynomial.Polynomial: f(x) = {poly}')
    print_roots(poly.roots())

    print('  np.roots:')
    print_roots(np.roots(coefficients.astype(i)))

我的输出

numpy version: 1.26.4

type: <class 'numpy.float32'>
  np.polynomial.Polynomial: f(x) = 160000.0 - 800.0·x + 1.0·x²
    roots (round to 5 decimals): [400. 400.]
    roots (round to 6 decimals): [400. 400.]
    roots (full decimals)      : [400. 400.]
    roots (full decimals real) : [400. 400.]
    roots (full decimals imag) : [0. 0.]

  np.roots
    roots (round to 5 decimals): [400. 400.]
    roots (round to 6 decimals): [400. 400.]
    roots (full decimals)      : [400. 400.]
    roots (full decimals real) : [400. 400.]
    roots (full decimals imag) : [0. 0.]

type: <class 'numpy.float64'>
  np.polynomial.Polynomial: f(x) = 160000.0 - 800.0·x + 1.0·x²
    roots (round to 5 decimals): [400. 400.]
    roots (round to 6 decimals): [399.999995 400.000005]
    roots (full decimals)      : [399.9999954 400.0000046]
    roots (full decimals real) : [399.9999954 400.0000046]
    roots (full decimals imag) : [0. 0.]

  np.roots
    roots (round to 5 decimals): [400. 400.]
    roots (round to 6 decimals): [400.000002 399.999998]
    roots (full decimals)      : [400.00000242 399.99999758]
    roots (full decimals real) : [400.00000242 399.99999758]
    roots (full decimals imag) : [0. 0.]