四次函数的根

Question

我在进行一些高级碰撞检测时遇到了一种情况，需要计算四次函数的根。

我使用法拉利的通用解决方案编写了一个似乎工作正常的函数，如下所示： http:// /en.wikipedia.org/wiki/Quartic_function#Ferrari.27s_solution。

这是我的功能：

    private function solveQuartic(A:Number, B:Number, C:Number, D:Number, E:Number):Array{          
        // For paramters: Ax^4 + Bx^3 + Cx^2 + Dx + E
        var solution:Array = new Array(4);

        // Using Ferrari's formula: http://en.wikipedia.org/wiki/Quartic_function#Ferrari.27s_solution
        var Alpha:Number = ((-3 * (B * B)) / (8 * (A * A))) + (C / A);
        var Beta:Number = ((B * B * B) / (8 * A * A * A)) - ((B * C) / (2 * A * A)) + (D / A);          
        var Gamma:Number = ((-3 * B * B * B * B) / (256 * A * A * A * A)) + ((C * B * B) / (16 * A * A * A)) - ((B * D) / (4 * A * A)) + (E / A);

        var P:Number = ((-1 * Alpha * Alpha) / 12) - Gamma; 
        var Q:Number = ((-1 * Alpha * Alpha * Alpha) / 108) + ((Alpha * Gamma) / 3) - ((Beta * Beta) / 8);

        var PreRoot1:Number = ((Q * Q) / 4) + ((P * P * P) / 27);
        var R:ComplexNumber = ComplexNumber.add(new ComplexNumber((-1 * Q) / 2), ComplexNumber.sqrt(new ComplexNumber(PreRoot1)));

        var U:ComplexNumber = ComplexNumber.pow(R, 1/3);

        var preY1:Number = (-5 / 6) * Alpha;
        var RedundantY:ComplexNumber = ComplexNumber.add(new ComplexNumber(preY1), U);

        var Y:ComplexNumber;

        if(U.isZero()){
            var preY2:ComplexNumber = ComplexNumber.pow(new ComplexNumber(Q), 1/3);

            Y = ComplexNumber.subtract(RedundantY, preY2);
        } else{
            var preY3:ComplexNumber = ComplexNumber.multiply(new ComplexNumber(3), U);
            var preY4:ComplexNumber = ComplexNumber.divide(new ComplexNumber(P), preY3);

            Y = ComplexNumber.subtract(RedundantY, preY4);
        }

        var W:ComplexNumber = ComplexNumber.sqrt(ComplexNumber.add(new ComplexNumber(Alpha), ComplexNumber.multiply(new ComplexNumber(2), Y)));

        var Two:ComplexNumber = new ComplexNumber(2);
        var NegativeOne:ComplexNumber = new ComplexNumber(-1);

        var NegativeBOverFourA:ComplexNumber = new ComplexNumber((-1 * B) / (4 * A));
        var NegativeW:ComplexNumber = ComplexNumber.multiply(W, NegativeOne);

        var ThreeAlphaPlusTwoY:ComplexNumber = ComplexNumber.add(new ComplexNumber(3 * Alpha), ComplexNumber.multiply(new ComplexNumber(2), Y));
        var TwoBetaOverW:ComplexNumber = ComplexNumber.divide(new ComplexNumber(2 * Beta), W);

        solution["root1"] = ComplexNumber.add(NegativeBOverFourA, ComplexNumber.divide(ComplexNumber.add(W, ComplexNumber.sqrt(ComplexNumber.multiply(NegativeOne, ComplexNumber.add(ThreeAlphaPlusTwoY, TwoBetaOverW)))), Two));
        solution["root2"] = ComplexNumber.add(NegativeBOverFourA, ComplexNumber.divide(ComplexNumber.subtract(NegativeW, ComplexNumber.sqrt(ComplexNumber.multiply(NegativeOne, ComplexNumber.subtract(ThreeAlphaPlusTwoY, TwoBetaOverW)))), Two));
        solution["root3"] = ComplexNumber.add(NegativeBOverFourA, ComplexNumber.divide(ComplexNumber.subtract(W, ComplexNumber.sqrt(ComplexNumber.multiply(NegativeOne, ComplexNumber.add(ThreeAlphaPlusTwoY, TwoBetaOverW)))), Two));
        solution["root4"] = ComplexNumber.add(NegativeBOverFourA, ComplexNumber.divide(ComplexNumber.add(NegativeW, ComplexNumber.sqrt(ComplexNumber.multiply(NegativeOne, ComplexNumber.subtract(ThreeAlphaPlusTwoY, TwoBetaOverW)))), Two));

        return solution;
    }

唯一的问题是我似乎遇到了一些例外。最值得注意的是，当我有两个实根和两个虚根时。

例如，这个方程： y = 0.9604000000000001x^4 - 5.997600000000001x^3 + 13.951750054511718x^2 - 14.326264455924333x + 5.474214401412618

返回根： 1.7820304835380467 + 0i 1.34041662585388 + 0i 1.3404185025061823 + 0i 1.7820323472855648 + 0i

如果我绘制该特定方程，我可以看到实际根更接近 1.2 和 2.9（大约）。我不能将四个不正确的根视为随机的，因为它们实际上是方程一阶导数的两个根：

y = 3.8416x^3 - 17.9928x^2 + 27.9035001x - 14.326264455924333

请记住，我是实际上并没有寻找我发布的方程的具体根。我的问题是是否有某种特殊情况我没有考虑到。

有什么想法吗？

Answer 1

为了求次数 >= 3 的多项式的根，我总是使用 Jenkins-Traub 得到更好的结果（ http://en.wikipedia.org/wiki/Jenkins-Traub_algorithm）而不是显式公式。

Answer 2

我不知道为什么法拉利的解决方案不起作用，但我尝试使用标准数值方法（创建伴随矩阵并计算其特征值），并且我获得了正确的解决方案，即在1.2和1.9处的两个实根。

这种方法不适合胆小的人。构造多项式的伴随矩阵后，运行QR 算法 查找该矩阵的特征值。这些是多项式的零点。

我建议您使用 QR 算法的现有实现，因为它的大部分内容比算法更接近厨房食谱。但我相信，它是计算特征值以及多项式根时使用最广泛的算法。

Answer 3

您可以查看我对相关问题的回答。我支持奥利维尔的观点：要走的路可能就是伴随矩阵/特征值方法（非常稳定、简单、可靠和快速）。

编辑

我想，为了方便起见，如果我在这里重现答案并没有什么坏处：

以可靠、稳定的方式多次执行此操作的数值解决方案，涉及：（1）形成伴随矩阵，（ 2) 求伴随矩阵的特征值。

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

对于多项式：

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

伴随矩阵是：

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

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

为了快速完成此操作，请从 LAPACK 下载奇异值子例程，编译它们，并将它们链接到您的代码。如果您有太多（例如，大约一百万）组系数，请并行执行此操作。您可以使用 QR 分解或任何其他稳定的方法来计算特征值（请参阅有关“矩阵特征值”的维基百科条目）。

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

祝你好运。

编辑：Numpy 和 Octave 也依赖于这种方法来计算多项式的根。例如，请参阅此链接。

Answer 4

其他答案都是好的、合理的建议。然而，回顾我在 Forth 中实现 Ferrari 方法的经验，我认为您的错误结果可能是由于 1. 错误地实现了必要且相当棘手的符号组合，2. 尚未意识到“.. == beta”浮点应该变成“abs(.. - beta) < eps, 3.尚未发现代码中还有其他平方根可能返回复杂的解决方案。

对于这个特定问题，我在诊断模式下的第四代码返回：

x1 =  1.5612244897959360787072371026316680470492e+0000 -1.6542769593216835969789894020584464029664e-0001 i
 --> -4.2123274051525879873007970023884313331788e-0054  3.4544674220377778501545407451201598284464e-0077 i
x2 =  1.5612244897959360787072371026316680470492e+0000  1.6542769593216835969789894020584464029664e-0001 i
 --> -4.2123274051525879873007970023884313331788e-0054 -3.4544674220377778501545407451201598284464e-0077 i
x3 =  1.2078440724224197532447709413299479764843e+0000  0.0000000000000000000000000000000000000000e-0001 i
 --> -4.2123274051525879873010733597821943554068e-0054  0.0000000000000000000000000000000000000000e-0001 i
x4 =  1.9146049071693819497220585618954851525216e+0000 -0.0000000000000000000000000000000000000000e-0001 i
 --> -4.2123274051525879873013497171759573776348e-0054  0.0000000000000000000000000000000000000000e-0001 i

“-->”之后的文本是通过将根代入原始方程得出的，

作为参考，以下是我设法设置的尽可能高的精度的 Mathematica/Alpha 结果：

Mathematica:
x1 = 1.20784407242
x2 = 1.91460490717
x3 = 1.56122449 - 0.16542770 i
x4 = 1.56122449 + 0.16542770 i 

Answer 5

已经提到的方法的一个很好的替代方法是 TOMS 算法 326，它基于论文“Roots低阶多项式”作者：Terence RFNonweiler CACM（1968 年 4 月）。

这是三阶和四阶多项式的代数解，相当紧凑、快速且相当准确。它比 Jenkins Traub 简单得多。

但请注意，TOMS 代码的工作效果并不那么好。

这个爱荷华山根求解器页面包含更精致的四次/三次根查找器代码。它还具有 Jenkins Traub 型寻根器。