Roots of a Quartic Function

Question

Hey guys,

I came across a situation doing some advanced collision detection, where I needed to calculate the roots of a quartic function.

I wrote a function that seems to work fine using Ferrari's general solution as seen here: http://en.wikipedia.org/wiki/Quartic_function#Ferrari.27s_solution.

Here's my function:

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

The only issue is that I seem to get a few exceptions. Most notably when I have two real roots, and two imaginary roots.

For example, this equation: y = 0.9604000000000001x^4 - 5.997600000000001x^3 + 13.951750054511718x^2 - 14.326264455924333x + 5.474214401412618

Returns the roots: 1.7820304835380467 + 0i 1.34041662585388 + 0i 1.3404185025061823 + 0i 1.7820323472855648 + 0i

If I graph that particular equation, I can see that the actual roots are closer to 1.2 and 2.9 (approximately). I can't dismiss the four incorrect roots as random, because they're actually two of the roots for the equation's first derivative:

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

Keep in mind that I'm not actually looking for the specific roots to the equation I posted. My question is whether there's some sort of special case that I'm not taking into consideration.

Any ideas?

Answer 1

Shouldn't you be posting this on MathOverflow?

Answer 2

I do not know why Ferrari's solution does not work, but I tried to use the standard numerical method (create a companion matrix and compute its eigenvalues), and I obtain the correct solution, i.e., two real roots at 1.2 and 1.9.

This method is not for the faint of heart. After constructing the companion matrix of the polynomial, you run the QR algorithm to find the eigenvalues of that matrix. Those are the zeroes of the polynomial.

I suggest you to use an existing implementation of the QR algorithm since a good deal of it is closer to kitchen recipe than algorithmics. But it is, I believe, the most widely used algorithm to compute eigenvalues, and thereby, roots of polynomials.

Answer 3

For finding roots of polynomials of degree >= 3, I've always had better results using Jenkins-Traub (http://en.wikipedia.org/wiki/Jenkins-Traub_algorithm) than explicit formulas.

Answer 4

You can see my answer to a related question. I support the view of Olivier: the way to go may just be the companion matrix / eigenvalue approach (very stable, simple, reliable, and fast).

Edit

I guess it does not hurt if I reproduce the answer here, for convenience:

The numerical solution for doing this many times in a reliable, stable manner, involve: (1) Form the companion matrix, (2) find the eigenvalues of the companion matrix.

You may think this is a harder problem to solve than the original one, but this is how the solution is implemented in most production code (say, Matlab).

For the polynomial:

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

the companion matrix is:

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

Find the eigenvalues of such matrix; they correspond to the roots of the original polynomial.

For doing this very fast, download the singular value subroutines from LAPACK, compile them, and link them to your code. Do this in parallel if you have too many (say, about a million) sets of coefficients. You could use QR decomposition, or any other stable methodology for computing eigenvalues (see the Wikipedia entry on "matrix eigenvalues").

Notice that the coefficient of t^3 is one, if this is not the case in your polynomials, you will have to divide the whole thing by the coefficient and then proceed.

Good luck.

Edit: Numpy and octave also depend on this methodology for computing the roots of polynomials. See, for instance, this link.