实三次多项式的最快数值解？答案

【问题标题】：Fastest numerical solution of a real cubic polynomial?实三次多项式的最快数值解？
【发布时间】：2011-01-01 11:41:19
【问题描述】：

R 问题：寻找最快的方法来以数字方式求解一组已知具有实系数和三个实根的任意三次方。据报道，R 中的 polyroot 函数将 Jenkins-Traub 的算法 419 用于复杂多项式，但对于实数多项式，作者参考了他们早期的工作。对于实数三次或更一般的实数多项式，更快的选择是什么？

【问题讨论】：

【解决方案1】：

您需要全部 3 个根还是只需要一个？如果只有一个，我认为牛顿法可以。如果全部 3 个，那么在两个靠近的情况下可能会出现问题。

【讨论】：

只有一个。已知位于 (0,1) 中的那个。
那么我投票赞成牛顿的方法。从 0.5 开始作为猜测，将域限制在 0 和 1 之间，它应该可以正常工作。它应该具有二次收敛性；您遇到问题的唯一原因是该函数是否在该范围内具有最小值或最大值（导数 = 0）

【解决方案2】：

常用的方法有：牛顿法、二分法、正割法、定点迭代法等。Google 任何一种。

另一方面，如果您有一个非线性系统（例如，在 N 个未知数中的 N 个多项式 eqn 上的系统），则可以使用诸如 high-order Newton 之类的方法。

【讨论】：

【解决方案3】：

在 n 变量中的任意多项式系统中找到真正解决方案的最快已知方法（据我所知）是多面体同伦。详细的解释可能超出了 StackOverflow 的答案，但本质上它是一种使用复曲面几何利用每个方程结构的路径算法。谷歌会给你a number of papers。

也许这个问题更适合mathoverflow？

【讨论】：

【解决方案4】：

以可靠、稳定的方式多次执行此操作的数值解决方案包括：(1) 形成伴随矩阵，(2) 找到伴随矩阵的特征值。

您可能认为这是一个比原始问题更难解决的问题，但这是解决方案在大多数生产代码（例如 Matlab）中的实现方式。

对于多项式：

p(t) = c0 + c1 * t + c2 * t^2 + t^3

伴随矩阵是：

[[0 0 -c0],[1 0 -c1],[0 1 -c2]]

求该矩阵的特征值；它们对应于原始多项式的根。

要快速完成此操作，请从 LAPACK 下载奇异值子例程，编译它们，并将它们链接到您的代码。如果您有太多（例如，大约一百万）组系数，请并行执行此操作。

注意t^3 的系数是1，如果在你的多项式中不是这种情况，你将不得不将整个事情除以系数然后继续。

祝你好运。

编辑：Numpy 和 octave 也依赖于这种计算多项式根的方法。例如，请参阅this link。

【讨论】：

我已经注意到这是一个 R 问题。对不起。
我认为这仍然是一个相关的答案，使用 ''eigen()'' 函数很容易在 R 中找到伴随矩阵的特征值。
+1。 GSL 也这样做。 Jenkins-Traub 算法是直接求解多项式的最先进的算法，但它有时会遇到不稳定性（而且它非常复杂），并且当仔细观察时，似乎在对伴随矩阵进行幂迭代......这您指出的方法是防弹的。

【解决方案5】：

【讨论】：

【解决方案6】：

充实上面 Arietta 的回答：

> a <- c(1,3,-4)
> m <- matrix(c(0,0,-a[1],1,0,-a[2],0,1,-a[3]), byrow=T, nrow=3)
> roots <- eigen(m, symm=F, only.values=T)$values

这是否比使用 GSL 包中的三次求解器更快或更慢（如上面 knguyen 所建议的）是在您的系统上对其进行基准测试的问题。

【讨论】：

伴生矩阵法可能是所有多项式求根算法中最慢的（尽管对于三次方，它应该仍然相当快）。 Jenkins-Traub 的实系数变体当然更快。 Jenkins-Traub 算法有一个开源 c++ 实现（BSD 许可），称为RPOLY++

【解决方案7】：

1) 求解导数多项式 P' 以定位您的三个根。请参阅there 以了解如何正确执行此操作。将这些根称为 a 和 b（使用 a

2) 对于中间根，在 a 和 b 之间使用一些平分步骤，当你足够接近时，使用牛顿法完成。

3) 对于最小和最大根，“寻找”解决方案。对于最大根：

【讨论】：