Roots of a Quartic Function

Question

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?

I don't know why you didn't just implement the beautiful equations from this PlanetMath site: http://planetmath.org/encyclopedia/QuarticFormula.html :) — Justin Peel, Jun 04 '10 at 17:52
Are you interested in the complex solutions too? If you are trying to find the Real solutions any standard numerical root finding method will do with adequate seeds. — Dr. belisarius, Jun 04 '10 at 21:27
Don't use this method! It's beautiful mathematics but a bad algorithm. You probably only want real roots. Count up how many real operations are in your code, including expanding out the complex arithmetic. It's expensive compared to iterative algorithms. Additionally, can you say how accurate the final result is after all those operations? In fact, the iterative methods will probably be more accurate. Other people have listed those methods. — sigfpe, Jun 04 '10 at 22:49

Olivier Verdier · Answer 1 · 2010-06-04T15:54:13.557

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.

Can you suggest a resource for learning how to use this method? — Scott Langendyk, Jun 04 '10 at 15:46

score 4 · Accepted Answer · edited Jan 04 '16 at 07:16

4

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.

edited Jan 04 '16 at 07:16

kristianp

5,496
37
56

answered Jun 04 '10 at 16:00

user168715

5,469
1
31
42

Interesting, I assume I can then use one function to compute the roots for polynomials of varying degrees? – Scott Langendyk Jun 04 '10 at 16:14
Just implemented this method, and it seems to be working great! Thanks a lot! – Scott Langendyk Jun 04 '10 at 23:07
@Scott Langendyk, you should accept the answer then by clicking the check mark off to the side. – Justin Peel Jun 10 '10 at 04:23
FYI there is an open source c++ implementation (BSD license) of the Jenkins-Traub algorithm called [RPOLY++](https://github.com/sweeneychris/RpolyPlusPlus) – kip622 Sep 20 '15 at 15:01
Or see http://www.akiti.ca/PolyRootRe.html for a javascript translation of RPOLY. – Lutz Lehmann Jan 31 '16 at 18:59

score 2 · Answer 3 · edited May 23 '17 at 11:54

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.

Notice that my original answer was for a cubic polynomial, but the approach is identical for degree `n`; in your case `n=4`. — Escualo, Jun 04 '10 at 16:55
If you need the subroutine for computing QR using a header-only C++ library, download the "TNT - Template Numerical Toolkit" - it defines a `Matrix` object, on which you can call the `QR` method, or, even better, you can form the companion Matrix and call the `eigenvalue` method. — Escualo, Jun 04 '10 at 16:59

marcel hendrix · Answer 4 · 2012-02-26T09:54:48.850

The other answers are good and sound advice. However, recalling my experience with the implementation of Ferrari's method in Forth, I think your wrong results are probably caused by 1. wrong implementation of the necessary and rather tricky sign combinations, 2. not realizing yet that ".. == beta" in floating-point should become "abs(.. - beta) < eps, 3. not yet having found out that there are other square roots in the code that may return complex solutions.

For this particular problem my Forth code in diagnostic mode returns:

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

The text after "-->" follows from backsubstituting the root into the original equation.

For reference, here are Mathematica/Alpha's results to the highest possible precision I managed to set it:

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

user5108_Dan · Answer 5 · 2016-02-11T14:55:17.017

A good alternative to the methods already mentioned is the TOMS Algorithm 326, which is based on the paper "Roots of Low Order Polynomials" by Terence R.F.Nonweiler CACM (Apr 1968).

This is an algebraic solution to 3rd and 4th order polynomials that is reasonably compact, fast, and quite accurate. It is much simpler than Jenkins Traub.

Be warned however that the TOMS code doesn't work all that well.

This Iowa Hills Root Solver page has code for a Quartic / Cubic root finder that is a bit more refined. It also has a Jenkins Traub type root finder.

Roots of a Quartic Function

5 Answers5

Linked