3

I'm attempting to create a function that calculates the 4 roots of a 4th degree polynomial with included complex numbers. In my search for a formula, I came across a rather simple one contained in this discussion, described by Tito Piezas III towards the bottom of the page.

Now, I believe that the real error isn't really fault in my code (as I'm sure will be annoying to proofread) but in my understanding of the method involved. My problem is, that the quadratic roots are complex, and I don't know how to use the complex numbers when calculating the quartic roots, programatically.

He suggests deriving the quartic roots using the roots of two quadratic equations. I attempted to mimic the formula the best I can with my code below. The idea is that I calculate the two quadratic roots (under the premise they are positive only --- I don't know how otherwise), then, using those results, I can calculate the qurtic roots which I then save the real and complex values into x1,x2,x3,x4 into r1,r2,r3,r4,c1,c2,c3,c4 respectively. But, when calculating the quadratic roots, u, a value used later to calculate quartic roots: is complex!

Here is an image of his formula and steps. Blow is my code with captions on most steps.

double a, b, c, d;
double c1, c2, c3, c4;      //complex values
double r1, r2, r3, r4;      //real values

//      x^4+ax^3+bx^2+cx+d=0
a = 3;     
b = 4; 
c = 5;     //<--- example coefficients
d = 6;
if (a != 0) {

    double u,v1,v2;       
    double x,y,z;        //essentially a,b,c that he uses

    x=1;
    y= -2*b*b*b+9*a*b*c-27*c*c-27*a*a*d+72*b*d;
    z= Math.pow((b*b-3*a*c+12*d),3);

      //calculation of the v roots
    v1 = -y+(Math.sqrt(y*y-4*x*z))/(2*x);  // < negative root
    v2 = -y-(Math.sqrt(y*y-4*x*z))/(2*x);  // < negative root

//---calculations after this are invalid since v1 and v2 are NaN---       

    u = (a*a)/4 + ((-2*b+Math.pow(v1,1/3)+Math.pow(v2,1/3))/3);

    double x12sub,x34sub;

    x12sub= 3*a*a-8*b-4*u+((-a*a*a+4*a*b-8*c)/(Math.sqrt(u)));
    x34sub= 3*a*a-8*b-4*u-((-a*a*a+4*a*b-8*c)/(Math.sqrt(u)));

    r1 = -(1/4)*a +(1/2)*(Math.sqrt(u));
    r2 = -(1/4)*a +(1/2)*(Math.sqrt(u));
    r3 = -(1/4)*a -(1/2)*(Math.sqrt(u));
    r4 = -(1/4)*a -(1/2)*(Math.sqrt(u));

//--casting results into their orderly variables--

    if(x12sub<0){
        x12sub= x12sub*-1;
        x12sub = Math.sqrt(x12sub);
        x12sub = x12sub*(1/4);
        c1=x12sub;
        c2=x12sub;
    }
    else{
        r1=r1+x12sub;
        r2=r2-x12sub;
    }
    if(x34sub<0){
        x34sub= x34sub*-1;
        x34sub = Math.sqrt(x34sub);
        x34sub = x34sub*(1/4);
        c3=x34sub;
        c4=x34sub;  
    }
    else{
        r3=r3+x34sub;
        r4=r4+x34sub;
    }
}

I'm open for ANY solution. Even ones which involve the use of libraries that could help me. Thanks for the help.

Nikolai Shevchenko
  • 7,083
  • 8
  • 33
  • 42
daedsidog
  • 1,732
  • 2
  • 17
  • 36
  • [This](http://math.stackexchange.com/questions/44406/how-do-i-get-the-square-root-of-a-complex-number) might be useful for part of your problem. – President James K. Polk Jan 19 '16 at 00:48
  • 1
    Possible duplicate of [Finding real roots of quartic equation using ferrari's method](http://stackoverflow.com/questions/27217113/finding-real-roots-of-quartic-equation-using-ferraris-method) – Teepeemm Jan 19 '16 at 01:19
  • You're probably better off using a non-explicit method such as the Jenkins-Traub method. Here is a [pretty good article](http://math-blog.com/2008/03/06/polynomial-root-finding-with-the-jenkins-traub-algorithm/) about the Jenkins-Traub method. – Robert Dodier Jan 19 '16 at 06:15
  • I see that there is an implementation of Laguerre's method in [Apache Commons-Math](https://commons.apache.org/proper/commons-math/). Here's the javadoc for [LaguerreSolver](https://commons.apache.org/proper/commons-math/apidocs/org/apache/commons/math3/analysis/solvers/LaguerreSolver.html). – Robert Dodier Jan 19 '16 at 06:19
  • @RobertDodier: This article tells you almost nothing about Jenkins-Traub. Compare with the wikipedia article. Synthetic division of linear factors is the fundamental building block, not just the method for deflation. -- However, the companion javascript translation http://www.akiti.ca/PolyRootRe.html of RPoly is often useful for a quick root computation. – Lutz Lehmann Jan 19 '16 at 08:09

2 Answers2

1

Try using the Efficient Java Matrix Library. You can download the jars here: https://sourceforge.net/projects/ejml/files/v0.28/

You need to have this method in your class:

public static Complex64F[] findRoots(double... coefficients) {
    int N = coefficients.length-1;

    // Construct the companion matrix
    DenseMatrix64F c = new DenseMatrix64F(N,N);

    double a = coefficients[N];
    for( int i = 0; i < N; i++ ) {
        c.set(i,N-1,-coefficients[i]/a);
    }
    for( int i = 1; i < N; i++ ) {
        c.set(i,i-1,1);
    }

    // use generalized eigenvalue decomposition to find the roots
    EigenDecomposition<DenseMatrix64F> evd =  DecompositionFactory.eig(N,false);

    evd.decompose(c);

    Complex64F[] roots = new Complex64F[N];

    for( int i = 0; i < N; i++ ) {
        roots[i] = evd.getEigenvalue(i);
    }

    return roots;
}

Then you can use this to find, for example, the roots of x^2 + 4x + 4:

Complex64F[] c = findRoots(4, 4, 1);
    for(Complex64F f : c)
        System.out.println(f.toString());

This will print out:

-2
-2

Which is the desired result.

sguan
  • 1,039
  • 10
  • 26
1

You can deplete the polynomial az^4+bz^3+cz^2+dz+e to the form t^4+pt^2+qt+r by a linear change of variable. Then you factor as (t^2 + ux + v)(t^2 - ux + w). This leads to the identification p = v + w - u^2, q = u (w - v), r = vw. From this, 2w = p + u^2 + q / u, 2v = p + u^2 - q / u and finally, (p + u^2 + q / u)(p + u^2 - q / u) = 4r. After rearranging, this is a cubic equation in u^2.

To solve a cubic az^3+bz^2+cz+d, you can again deplete to t^3+pt+q, and write t = u - p / 3u. This gives u^3 + q - p^3 / 27u^3 = 0, a quadratic equation in u^3.

To solve a quadratic polynomial ax^2 + bx + c, you can deplete to the form t^2 + p = 0. The roots are obviously ±√(-p).

Hence you solve the quadratic for u^3 and take the cubic root of one of the roots. From it, you apply the inverse depletion transformation and obtain a solution of the cubic. Now you can compute the coefficients of the factored quartic, and solve the two factors separately, obtaining four roots. Finally, undo the initial depletion.

Easy, no ?

To deplete a polynomial, divide all coefficients by the leading one, then set x = t + s. Then substituting in the polynomial and expanding, you cancel the coefficient of the second term. This gives s = - b / ka where k is the degree of the polynomial. The other coefficients follow from s.

For a polynomial of real coefficients, note that the cubic equation always have a real solution, giving a factorization in two quadratic polynomials.

Anyway, the resolution of a cubic of real coefficients may involve an unescapable resort to complex numbers.