r - 实三次多项式的最快数值解？

Question

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

score 6 · Accepted Answer

以可靠、稳定的方式多次执行此操作的数值解包括：（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 也依赖于这种计算多项式根的方法。例如，请参阅此链接。

score 4 · Accepted Answer

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

也许这个问题更适合mathoverflow？

score 2 · Accepted Answer

在上面充实 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 所建议的）是在您的系统上对其进行基准测试的问题。

score 1 · Accepted Answer

1) Solve for the derivative polynomial P' to locate your three roots. See there to know how to do it properly. Call those roots a and b (with a < b)

2) For the middle root, use a few steps of bisection between a and b, and when you're close enough, finish with Newton's method.

3) For the min and max root, "hunt" the solution. For the max root:

Start with x0 = b, x1 = b + (b - a) * lambda, where lambda is a moderate number (say 1.6)
do x_n = b + (x_{n - 1} - a) * lambda until P(x_n) and P(b) have different signs
Perform bisection + newton between x_{n - 1} and x_n

score 1 · Accepted Answer

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

score 0 · Accepted Answer

0

于 2010-01-09T19:37:06.127 回答

score 0 · Accepted Answer

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

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

7 回答 7