如何正确求多项式根？

Question

考虑一个多项式，例如：

p = [1 -9 27 -27];

显然真正的根是 3:

polyval(p,3)

0

使用 roots 函数时

q = roots([1 -9 27 -27]);

与 format short:

q =

   3.0000 + 0.0000i
   3.0000 + 0.0000i
   3.0000 - 0.0000i

并检查根是否真实：

bsxfun(@eq,ones(size(q)),isreal(q))

0
0
0

更糟糕的是 format long 我得到：

roots([1 -9 27 -27])

ans =

  3.000019414068325 + 0.000000000000000i
  2.999990292965843 + 0.000016813349886i
  2.999990292965843 - 0.000016813349886i

如何正确计算多项式的根？

Answer 1

这是由于浮点数不准确造成的。看看这个 post 了解详情： Is floating point math broken?

您可以做的一件事是将 answer/s 四舍五入到一些小数位，如下所示：

q = round(roots([1 -9 27 -27]), 4) % rounding off to 4 decimal places

Answer 2

您可能需要象征性地工作。为此，您需要符号数学工具箱。

将多项式定义为符号函数。您可以 (a) 使用 poly2sym 从其系数生成符号多项式。或者 (b) 更好的是，直接使用字符串定义符号函数。这样您就可以避免将系数表示为 double.
使用solve，符号求解代数方程。

带有选项 (a) 的代码：

p = [1 -9 27 -27];
ps = poly2sym(p);
rs = solve(ps);

带有选项 (b) 的代码：

ps = sym('x^3-9*x^2+27*x-27');
rs = solve(ps);

无论哪种情况，结果都是象征性的：

>> rs
rs =
 3
 3
 3

您可能希望使用

转换为数值

r = double(rs);

在您的示例中，这给出了

>> format long
>> r
r =
     3
     3
     3

Answer 3

这对您的多项式来说非常具体。通常，您必须期望多重根 m 具有大小为 mu^(1/m) 的相对浮点误差，其中 mu=1e-15 是机器精度。在这种情况下，多重性为 m=3，因此误差范围为 10^(-5)。这正是您结果中错误的规模。

这里出现的这种情况具有明确的整数系数是matlab使用的方法的结果，它计算伴随矩阵的特征值，特征值算法会将整数矩阵转换为具有相应舍入误差的适当浮点矩阵在算法的第一步。

其他算法对多重性和相关的近似根簇进行了实证检验，因此能够纠正此错误。在这种情况下，您可以通过将每个根替换为 3 个根的平均值来实现。

在数学上，你有一些多项式

p(x)=(x-a)^m*q(x)

在多重性 m 的 x=a 处有一个根。由于浮点运算，求解器有效地 "sees" 一个多项式

p(x)+e(x)

其中 e(x) 的系数的大小是 p 乘以 mu 的系数的大小。靠近根a，这个扰动多项式可以有效地替换为

(x-a)^m*q(a)+e(a) = 0  <==>  (x-a)^m = -e(a)/q(a)

使得解形成一个以 a 为中心、半径为 |e(a)/q(a)|^(1/m) 的 m 点正多边形或星形，它应该在 |a|*mu^(1/m).

的区域内

如何正确求多项式根？

How to find polynomial roots correctly?

matlab

polynomials