math - 四次函数的根

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;
    }

唯一的问题是我似乎得到了一些异常(exception)。最值得注意的是，当我有两个实根和两个虚根时。

例如，这个等式:
y = 0.9604000000000001x^4 - 5.997600000000001x^3 + 13.951750054511718x^2 - 14.326264455924333x + 5.4618

返回根:
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

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

这种方法不适合胆小的人。构建后companion matrix多项式，你运行 QR algorithm找到该矩阵的特征值。这些是多项式的零点。

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