The roots of the quartz function

Question

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

I wrote a function that works fine using the generic Ferrari solution, as shown here: http://en.wikipedia.org/wiki/Quartic_function#Ferrari.27s_solution .

Here is 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 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 problem is that I seem to get a few exceptions. First of all, 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 draw this particular equation, I can see that the actual roots are closer to 1.2 and 2.9 (approximately). I cannot discard the four irregular roots as random, because they are actually two roots for the first derivative of the equation:

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

Keep in mind that I'm not really looking for specific roots for the equation I published. My question is: is there any special case that I do not take into account.

Any ideas?

Answer 1

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

Answer 2

I do not know why the Ferrari solution does not work, but I tried to use the standard numerical method (create an accompanying matrix and calculate its eigenvalues), and I get the correct solution, i.e. two material roots in 1,2 and 1,9.

This method is not for the faint of heart. After building the companion matrix of the polynomial, you run the QR algorithm to find the eigenvalues ​​of this matrix. These are zeros of a polynomial.

I suggest you use the existing implementation of the QR algorithm, since many of them are closer to the kitchen recipe than algorithmic ones. But this, I believe, is the most widely used algorithm for calculating eigenvalues ​​and, therefore, the roots of polynomials.

Answer 3

Other answers are good and sound advice. However, recalling my experience of using the Ferrari method in Fort, I think that your incorrect results are probably caused by the incorrect implementation of the necessary and rather complex combinations of signs, 2. not yet realizing that the ".. == beta" in floating point should become "abs (.. - beta)), but didn’t find out that there are other square roots in the code that can 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 the inverse of the root equation.

For reference, here are the results of Mathematica / Alpha with the highest possible accuracy I was able to install:

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

Answer 4

You can see my answer to the corresponding question . I support Olivier's point of view: the way to go can be just an accompanying matrix / approach based on eigenvalues ​​(very stable, simple, reliable and fast).

Edit

I think this will not hurt if I reproduce the answer here, for convenience:

The numerical solution for this many times in a reliable, stable way includes: (1) Forms a matrix of companions, (2) we find the eigenvalues ​​of the matrix of related ones.

You might think that solving this problem is more difficult than the original, but this is how the solution is implemented in most production codes (for example, Matlab).

For polynomial:

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

companion matrix:

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

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

To do this, load the routines of singular values ​​from LAPACK very quickly, compile them and bind them to your code. Do this in parallel if you have too many (say, about a million) sets of odds. You can use a QR decomposition or any other stable methodology to calculate eigenvalues ​​(see Wikipedia entry “Matrix eigenvalues”).

Please note that the coefficient at t^3 is equal to one, if this is not the case in your polynomials, you will need to divide everything by a coefficient and then continue.

Good luck.

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

Answer 5

A good alternative to the methods already mentioned is TOMS Algorithm 326 , which is based on the article "Roots of Low Order Polynomial Orders" on the Terence RFNonweiler CACM (April 1968).

This is an algebraic solution of 3rd and 4th order polynomials, which is compact enough, fast and fairly accurate. This is a lot easier than Jenkins Traub.

It should be warned that TOMS code does not work so well.

This Iota Hills Root Solver page contains the Quartic / Cubic root finder code, which is a bit more sophisticated. It also has a root crawler like Jenkins Traub.

Answer 6

I also implemented this function and got a normal result, just like you. After debugging, for some time this state check "\ beta <0" should be passed, but it was not.