如何正确地找到多项式根

Question

考虑如下多项式:

p = [1 -9 27 -27];

显然实根是3:

polyval(p,3)
0

使用roots函数时

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

with 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

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

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

Answer 2

这是由于浮点数不准确造成的。看看这篇文章的细节:浮点数学坏了吗?

你可以做的一件事是将答案四舍五入到小数点后几位，像这样:

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

Answer 3

这是您的多项式所特有的。一般来说，您必须期望多重根m具有mu^(1/m)量级的相对浮点误差，其中mu=1e-15是机器精度。在本例中，多重度为m=3，因此误差在10^(-5)的范围内。这正是你的结果中误差的大小。

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

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

---

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

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为圆心，半径为|e(a)/q(a)|^(1/m)的m点正多边形或星形，该正多边形或星形应位于|a|*mu^(1/m)区域内。