5

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

我使用法拉利的通用解决方案编写了一个似乎可以正常工作的函数，如下所示：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

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

有任何想法吗？

4

5 回答 5

5

为了找到度数 >= 3 的多项式的根，我总是使用 Jenkins-Traub ( http://en.wikipedia.org/wiki/Jenkins-Traub_algorithm ) 比显式公式获得更好的结果。

于 2010-06-04T16:00:28.913 回答
4

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

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

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

于 2010-06-04T15:31:36.280 回答
2

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

编辑

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

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

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

对于多项式：

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

伴随矩阵是：

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

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

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

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

祝你好运。

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

于 2010-06-04T16:50:24.880 回答
2

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

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

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
于 2012-02-26T08:29:02.667 回答
0

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

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

但是请注意，TOMS 代码不能很好地工作。

这个Iowa Hills Root Solver 页面有一个更精细的四次/三次根查找器的代码。它还有一个 Jenkins Traub 类型的根查找器。

于 2016-01-29T14:40:06.560 回答