如何正确地找到多项式的根？

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

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

- NKN

---

小提示：您检查根是否为实数的方法不正确。如果数组q是复数，则isreal(q)会返回false。但是，某些条目可能具有零虚部。实际上，isreal(q)返回false，而for x = q(:).', isreal(x), end返回true，false，false。 q的第一个条目是实数，其他条目不是，整个q也不是实数。 - Luis Mendo

Answer 1

您可能需要使用符号计算工具箱来进行符号计算。

以下是翻译的内容：

您可能需要以符号形式工作。这需要使用符号数学工具箱。

1. 将多项式定义为符号函数。您可以 (a) 使用 poly2sym 将其系数转换为符号多项式。或者 (b) 更好的方法是，直接使用字符串定义符号函数。这样可以避免将系数表示为 double 时可能导致的精度损失。

2. 使用 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

- Luis Mendo

---

谢谢您的回答。我正在对1000个具有不同系数和次数的多项式执行此操作。就效率而言，是使用您答案中的选项（a）还是只使用@sardar_usama提到的round函数？ - NKN

我认为 round 函数更快，但是精度和效率之间存在权衡。 - NKN

不需要使用 poly2sym。这样可以正常工作：roots(sym([1 -9 27 -27]))。roots函数似乎没有对符号类型进行重载，但它调用的基础函数是重载的。在符号数学工具箱中还有 root 函数，它可以代替更通用的 solve 函数。 - horchler

@horchler 谢谢你的信息！ - Luis Mendo

@NKN 是的，round 更快（但可能不够准确）是有道理的。 - Luis Mendo

Answer 2

这是由于浮点数不精确性导致的。请参阅此帖子了解详细信息：浮点数运算有问题吗？

你可以做的一件事情是将答案四舍五入到某些小数位，像这样：

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

- Sardar Usama

---

谢谢您的回答，请查看我下面的评论。 - NKN

两个答案都很好，但很遗憾我不能都接受。所以我决定接受你的答案，因为我会在我的代码中使用它。同时给予@LuisMendo的回答一个赞(+1)，因为我学到了新的东西。谢谢你们俩。 - NKN

Answer 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)

这样，解决方案将形成一个以a为中心的m角形规则多边形或星形，半径为|e（a）/q（a）|^(1/m)，应该在|a|*mu^(1/m)区域内。

- Lutz Lehmann