如何正确求多项式根？答案

【问题标题】：How to find polynomial roots correctly?如何正确求多项式根？
【发布时间】：2017-01-08 18:30:20
【问题描述】：

考虑一个多项式，例如：

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

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

【问题讨论】：

次要注意：您检查根是否真实是不正确的。如果 array q 很复杂，isreal(q) 会给出 false。但有些条目的虚部可能为零。事实上，isreal(q) 给出了false，而for x = q(:).', isreal(x), end 给出了true、false、false。 q的第一个条目是真实的，其他的不是，q整体不是真实的

标签： matlab polynomials

【解决方案1】：

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

将多项式定义为符号函数。您可以 (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

【讨论】：

感谢您的回答。我正在为具有不同系数和度数的 1000 个多项式执行此操作。关于效率，在您的答案中使用 option(a) 或仅使用 @sardar_usama 提到的圆形函数是否有意义？
我认为round函数更快，但是在准确性和效率之间存在权衡。
我认为不需要使用poly2sym。这很好用：roots(sym([1 -9 27 -27]))。 roots 函数似乎没有重载符号类型，但它调用的底层函数是。符号数学工具箱中还有root，可以用来代替更通用的solve。
@horchler 感谢您的信息！
@NKN 是的，round 更快是有道理的（但可能不太准确）

【解决方案2】：

这是由于浮点数不准确造成的。看看这篇文章的详细信息： Is floating point math broken?

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

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

【讨论】：

感谢您的回答，请查看我下面的评论。
两个答案都很好，不幸的是我不能同时接受。所以我决定接受你的回答，因为我在我的代码中使用了这个。并给@LuisMendo 答案（+1），因为我学到了新东西。谢谢你们。

【解决方案3】：

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

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

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

数学上，你有一些多项式

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

以x=a 为根，多重m。由于浮点运算，求解器有效地“看到”多项式

p(x)+e(x)

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

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

这样解形成一个 m 点正多边形或星形，以 a 为中心，半径为 |e(a)/q(a)|^(1/m)，应该在 |a|*mu^(1/m) 的区域内。

【讨论】：