Introduction
This Java program implements the Durand-Kerner-Weierstrass method to find the roots of a univariate polynomial with complex coefficients (listing):
For better numerical stability, the implementation uses Horner's scheme to evaluate the polynomial. A similar implentation in Ada is cited here. Sample output is shown below. The implementation uses the JScience library for complex arithmetic.
Implementation
import org.jscience.mathematics.function.Polynomial;
import org.jscience.mathematics.function.Term;
import org.jscience.mathematics.function.Variable;
import org.jscience.mathematics.number.Complex;
/**
* An implementation of the Durand-Kerner-Weierstrass method to
* determine the roots of a complex univariate polynomial, as described
* <a href="http://en.wikipedia.org/wiki/Durand-Kerner_method">here</a>.
*
* @author John B. Matthews; distribution per LGPL.
*/
final class ComplexPolynomial {
private static final int MAX_COUNT = 999;
private static final boolean DEBUG = false;
private static double epsilon = 1E-15;
private ComplexPolynomial() {}
/**
* Create a complex polynomial from a number of Complex values.
* The array should have the highest order coefficient first.
*
* @param a a variable number of Complex coefficients.
* @return a Polynomial having those Complex coefficients.
*/
public static Polynomial<Complex> create(Complex... a) {
Variable<Complex> x = new Variable.Local<Complex>("x");
Polynomial<Complex> px = Polynomial.valueOf(Complex.ZERO, x);
for (int i = 0, e = a.length - 1; i < a.length; i++, e--) {
px = px.plus(Polynomial.valueOf(a[i], Term.valueOf(x, e)));
}
return px;
}
/**
* Create a complex array from an array of double.
* The array should have the highest order coefficient first.
*
* @param a a variable number of double coefficients.
* @return an array of the corresponding Complex coefficients.
*/
public static Complex[] complexArray(double... a) {
Complex[] ca = new Complex[a.length];
for (int i = 0; i < a.length; i++) {
ca[i] = Complex.valueOf(a[i], 0);
}
return ca;
}
/**
* This implementation uses the Durand-Kerner-Weierstrass method
* to find the roots of a polynomial with complex coefficients.
* The method requires a monic polynomial; some error may occur
* when dividing by the leading coefficient.
* The array should have the highest order coefficient first.
*
* @param ca a variable number of Complex polynomial coefficients.
* @return an array of the Complex roots of the polynomial.
*/
public static Complex[] roots(Complex[] ca) {
Complex[] a0 = new Complex[ca.length - 1];
Complex[] a1 = new Complex[ca.length - 1];
// Divide by leading coefficient if not monic
Complex leading = ca[0];
if (!ca[0].equals(Complex.ONE)) {
for (int i = 0; i < ca.length; i++) {
ca[i] = ca[i].divide(leading);
}
}
// Initialize a0
Complex result = Complex.valueOf(0.4, 0.9);
a0[0] = Complex.ONE;
for (int i = 1; i < a0.length; i++) {
a0[i] = a0[i - 1].times(result);
}
// Iterate
int count = 0;
while (true) {
for (int i = 0; i < a0.length; i++) {
result = Complex.ONE;
for (int j = 0; j < a0.length; j++) {
if (i != j) {
result = a0[i].minus(a0[j]).times(result);
}
}
a1[i] = a0[i].minus(ComplexPolynomial.
eval(ca, a0[i]).divide(result));
}
count++;
if (count > MAX_COUNT || done(a0, a1)) break;
System.arraycopy(a1, 0, a0, 0, a1.length); // a0 := a1
}
if (DEBUG) {
System.out.println("Iterations: " + count);
for (Complex c : a0) System.out.println(c);
}
return a1;
}
// Determine if the components of two arrays are unchanging
private static boolean done(Complex[] a, Complex[] b) {
boolean unchanged = true;
Complex delta;
for (int i = 0; i < a.length; i++) {
delta = a[i].minus(b[i]);
unchanged &= Math.abs(delta.getReal()) < epsilon
&& Math.abs(delta.getImaginary()) < epsilon;
}
return unchanged;
}
/**
* Evaluate the polynomial at x using
* <a href="http://en.wikipedia.org/wiki/Horner_scheme">Horner's scheme</a>.
* The array should have the highest order coefficient first.
*
* @param ca an array of Complex polynomial coefficients.
* @param x the function argument.
* @return the Complex value of the function at x.
*/
public static Complex eval(Complex[] ca, Complex x) {
Complex result = ca[0];
for (int i = 1; i < ca.length; i++) {
result = result.times(x).plus(ca[i]);
}
return result;
}
/** Return the nominal tolerance value. */
public static double getEpsilon() { return epsilon; }
/** Set the nominal tolerance value. */
public static void setEpsilon(double e) { epsilon = e; }
}
Sample test output
The test driver produces the following output using double complex arithmetic. On the test platform, double offers more than 15 decimal digits of precision, and a corresponding value of epsilon is used. For each polynomial P and root r, the maximal value of abs(P(r) – 0) is shown. The first few entries are easily factored polynomials with integer coefficients. These are followed by a series of all one polynomials. Finally, there are a few examples with complex coefficients. On PowerPC (G4), the largest useful order appears to be about 38; on Intel (Pentium 4) and AMD (Athlon), extended precision may give correct results up to order 48. Beyond that, round-off errors cause the roots to converge to zero.
Poly: [1.0 + 0.0i]x^2 + [-3.0 + 0.0i]x + [2.0 + 0.0i]
Real: 1.00000000000000E0
Real: 2.00000000000000E0
Error: < 1.0E-15
Poly: [1.0 + 0.0i]x^2 + [-2.0 + 0.0i]x + [1.0 + 0.0i]
Real: 1.00000000000000E0
Real: 1.00000000000000E0
Error: < 1.0E-15
Poly: [1.0 + 0.0i]x^2 + [4.0 + 0.0i]
Comp: 0.00000000000000E0 +- 2.00000000000000E0i
Error: < 1.0E-15
Poly: [4.0 + 0.0i]x^2 + [1.0 + 0.0i]
Comp: 0.00000000000000E0 +- 5.00000000000000E-1i
Error: < 1.0E-15
Poly: [1.0 + 0.0i]x^3 + [-3.0 + 0.0i]x^2 + [2.0 + 0.0i]x
Real: 0.00000000000000E0
Real: 1.00000000000000E0
Real: 2.00000000000000E0
Error: < 1.0E-15
Poly: [1.0 + 0.0i]x^4 + [-13.0 + 0.0i]x^2 + [36.0 + 0.0i]
Real: -3.00000000000000E0
Real: -2.00000000000000E0
Real: 2.00000000000000E0
Real: 3.00000000000000E0
Error: < 1.0E-15
Poly: [1.0 + 0.0i]x^5 + [-13.0 + 0.0i]x^3 + [36.0 + 0.0i]x
Real: -3.00000000000000E0
Real: -2.00000000000000E0
Real: 0.00000000000000E0
Real: 2.00000000000000E0
Real: 3.00000000000000E0
Error: < 1.0E-15
Poly: [1.0 + 0.0i]x^2 + [1.0 + 0.0i]x + [1.0 + 0.0i]
Comp: -5.00000000000000E-1 +- 8.66025403784439E-1i
Error: < 1.0E-15
Poly: [1.0 + 0.0i]x^3 + [1.0 + 0.0i]x^2 + [1.0 + 0.0i]x + [1.0 + 0.0i]
Real: -1.00000000000000E0
Comp: 0.00000000000000E0 +- 1.00000000000000E0i
Error: < 1.0E-15
Poly: [1.0 + 0.0i]x^4 + [1.0 + 0.0i]x^3 + [1.0 + 0.0i]x^2 + [1.0 + 0.0i]x + [1.0 + 0.0i]
Comp: -8.09016994374948E-1 +- 5.87785252292473E-1i
Comp: 3.09016994374947E-1 +- 9.51056516295154E-1i
Error: < 1.0E-15
Poly: [1.0 + 0.0i]x^5 + [1.0 + 0.0i]x^4 + [1.0 + 0.0i]x^3 + [1.0 + 0.0i]x^2 + [1.0 + 0.0i]x + [1.0 + 0.0i]
Real: -1.00000000000000E0
Comp: -5.00000000000000E-1 +- 8.66025403784439E-1i
Comp: 5.00000000000000E-1 +- 8.66025403784439E-1i
Error: < 1.0E-15
Poly: [1.0 + 0.0i]x^6 + [1.0 + 0.0i]x^5 + [1.0 + 0.0i]x^4 + [1.0 + 0.0i]x^3 + [1.0 + 0.0i]x^2 + [1.0 + 0.0i]x + [1.0 + 0.0i]
Comp: -9.00968867902419E-1 +- 4.33883739117558E-1i
Comp: -2.22520933956314E-1 +- 9.74927912181824E-1i
Comp: 6.23489801858734E-1 +- 7.81831482468030E-1i
Error: < 1.0E-15
Poly: [1.0 + 0.0i]x^7 + [1.0 + 0.0i]x^6 + [1.0 + 0.0i]x^5 + [1.0 + 0.0i]x^4 + [1.0 + 0.0i]x^3 + [1.0 + 0.0i]x^2 + [1.0 + 0.0i]x + [1.0 + 0.0i]
Real: -1.00000000000000E0
Comp: -7.07106781186548E-1 +- 7.07106781186548E-1i
Comp: 0.00000000000000E0 +- 1.00000000000000E0i
Comp: 7.07106781186548E-1 +- 7.07106781186548E-1i
Error: < 1.0E-15
Poly: [1.0 + 0.0i]x^8 + [1.0 + 0.0i]x^7 + [1.0 + 0.0i]x^6 + [1.0 + 0.0i]x^5 + [1.0 + 0.0i]x^4 + [1.0 + 0.0i]x^3 + [1.0 + 0.0i]x^2 + [1.0 + 0.0i]x + [1.0 + 0.0i]
Comp: -9.39692620785908E-1 +- 3.42020143325669E-1i
Comp: -5.00000000000000E-1 +- 8.66025403784439E-1i
Comp: 1.73648177666930E-1 +- 9.84807753012208E-1i
Comp: 7.66044443118978E-1 +- 6.42787609686539E-1i
Error: < 1.0E-15
Poly: [1.0 + 0.0i]x^9 + [1.0 + 0.0i]x^8 + [1.0 + 0.0i]x^7 + [1.0 + 0.0i]x^6 + [1.0 + 0.0i]x^5 + [1.0 + 0.0i]x^4 + [1.0 + 0.0i]x^3 + [1.0 + 0.0i]x^2 + [1.0 + 0.0i]x + [1.0 + 0.0i]
Real: -1.00000000000000E0
Comp: -8.09016994374948E-1 +- 5.87785252292473E-1i
Comp: -3.09016994374947E-1 +- 9.51056516295154E-1i
Comp: 3.09016994374947E-1 +- 9.51056516295154E-1i
Comp: 8.09016994374948E-1 +- 5.87785252292473E-1i
Error: < 1.0E-15
Poly: [1.0 + 0.0i]x^10 + [1.0 + 0.0i]x^9 + [1.0 + 0.0i]x^8 + [1.0 + 0.0i]x^7 + [1.0 + 0.0i]x^6 + [1.0 + 0.0i]x^5 + [1.0 + 0.0i]x^4 + [1.0 + 0.0i]x^3 + [1.0 + 0.0i]x^2 + [1.0 + 0.0i]x + [1.0 + 0.0i]
Comp: -9.59492973614497E-1 +- 2.81732556841430E-1i
Comp: -6.54860733945285E-1 +- 7.55749574354258E-1i
Comp: -1.42314838273285E-1 +- 9.89821441880933E-1i
Comp: 4.15415013001886E-1 +- 9.09631995354518E-1i
Comp: 8.41253532831181E-1 +- 5.40640817455598E-1i
Error: 1.01904927555717E-15
Poly: [1.0 + 0.0i]x^11 + [1.0 + 0.0i]x^10 + [1.0 + 0.0i]x^9 + [1.0 + 0.0i]x^8 + [1.0 + 0.0i]x^7 + [1.0 + 0.0i]x^6 + [1.0 + 0.0i]x^5 + [1.0 + 0.0i]x^4 + [1.0 + 0.0i]x^3 + [1.0 + 0.0i]x^2 + [1.0 + 0.0i]x + [1.0 + 0.0i]
Real: -1.00000000000000E0
Comp: -8.66025403784439E-1 +- 5.00000000000000E-1i
Comp: -5.00000000000000E-1 +- 8.66025403784439E-1i
Comp: 0.00000000000000E0 +- 1.00000000000000E0i
Comp: 5.00000000000000E-1 +- 8.66025403784439E-1i
Comp: 8.66025403784439E-1 +- 5.00000000000000E-1i
Error: 1.14304456355485E-15
Poly: [1.0 + 0.0i]x^12 + [1.0 + 0.0i]x^11 + [1.0 + 0.0i]x^10 + [1.0 + 0.0i]x^9 + [1.0 + 0.0i]x^8 + [1.0 + 0.0i]x^7 + [1.0 + 0.0i]x^6 + [1.0 + 0.0i]x^5 + [1.0 + 0.0i]x^4 + [1.0 + 0.0i]x^3 + [1.0 + 0.0i]x^2 + [1.0 + 0.0i]x + [1.0 + 0.0i]
Comp: -9.70941817426052E-1 +- 2.39315664287558E-1i
Comp: -7.48510748171101E-1 +- 6.63122658240795E-1i
Comp: -3.54604887042536E-1 +- 9.35016242685415E-1i
Comp: 1.20536680255323E-1 +- 9.92708874098054E-1i
Comp: 5.68064746731156E-1 +- 8.22983865893656E-1i
Comp: 8.85456025653210E-1 +- 4.64723172043769E-1i
Error: < 1.0E-15
Poly: [1.0 + 0.0i]x^13 + [1.0 + 0.0i]x^12 + [1.0 + 0.0i]x^11 + [1.0 + 0.0i]x^10 + [1.0 + 0.0i]x^9 + [1.0 + 0.0i]x^8 + [1.0 + 0.0i]x^7 + [1.0 + 0.0i]x^6 + [1.0 + 0.0i]x^5 + [1.0 + 0.0i]x^4 + [1.0 + 0.0i]x^3 + [1.0 + 0.0i]x^2 + [1.0 + 0.0i]x + [1.0 + 0.0i]
Real: -1.00000000000000E0
Comp: -9.00968867902419E-1 +- 4.33883739117558E-1i
Comp: -6.23489801858734E-1 +- 7.81831482468030E-1i
Comp: -2.22520933956314E-1 +- 9.74927912181824E-1i
Comp: 2.22520933956314E-1 +- 9.74927912181824E-1i
Comp: 6.23489801858734E-1 +- 7.81831482468030E-1i
Comp: 9.00968867902419E-1 +- 4.33883739117558E-1i
Error: < 1.0E-15
Poly: [1.0 + 0.0i]x^14 + [1.0 + 0.0i]x^13 + [1.0 + 0.0i]x^12 + [1.0 + 0.0i]x^11 + [1.0 + 0.0i]x^10 + [1.0 + 0.0i]x^9 + [1.0 + 0.0i]x^8 + [1.0 + 0.0i]x^7 + [1.0 + 0.0i]x^6 + [1.0 + 0.0i]x^5 + [1.0 + 0.0i]x^4 + [1.0 + 0.0i]x^3 + [1.0 + 0.0i]x^2 + [1.0 + 0.0i]x + [1.0 + 0.0i]
Comp: -9.78147600733806E-1 +- 2.07911690817759E-1i
Comp: -8.09016994374948E-1 +- 5.87785252292473E-1i
Comp: -5.00000000000000E-1 +- 8.66025403784439E-1i
Comp: -1.04528463267653E-1 +- 9.94521895368273E-1i
Comp: 3.09016994374947E-1 +- 9.51056516295154E-1i
Comp: 6.69130606358858E-1 +- 7.43144825477394E-1i
Comp: 9.13545457642601E-1 +- 4.06736643075800E-1i
Error: 1.25607396694702E-15
Poly: [1.0 + 0.0i]x^15 + [1.0 + 0.0i]x^14 + [1.0 + 0.0i]x^13 + [1.0 + 0.0i]x^12 + [1.0 + 0.0i]x^11 + [1.0 + 0.0i]x^10 + [1.0 + 0.0i]x^9 + [1.0 + 0.0i]x^8 + [1.0 + 0.0i]x^7 + [1.0 + 0.0i]x^6 + [1.0 + 0.0i]x^5 + [1.0 + 0.0i]x^4 + [1.0 + 0.0i]x^3 + [1.0 + 0.0i]x^2 + [1.0 + 0.0i]x + [1.0 + 0.0i]
Real: -1.00000000000000E0
Comp: -9.23879532511287E-1 +- 3.82683432365090E-1i
Comp: -7.07106781186548E-1 +- 7.07106781186548E-1i
Comp: -3.82683432365090E-1 +- 9.23879532511287E-1i
Comp: 0.00000000000000E0 +- 1.00000000000000E0i
Comp: 3.82683432365090E-1 +- 9.23879532511287E-1i
Comp: 7.07106781186548E-1 +- 7.07106781186548E-1i
Comp: 9.23879532511287E-1 +- 3.82683432365090E-1i
Error: 2.71039981919101E-15
Poly: [1.0 + 0.0i]x^37 + [1.0 + 0.0i]x^36 + [1.0 + 0.0i]x^35 + [1.0 + 0.0i]x^34 + [1.0 + 0.0i]x^33 + [1.0 + 0.0i]x^32 + [1.0 + 0.0i]x^31 + [1.0 + 0.0i]x^30 + [1.0 + 0.0i]x^29 + [1.0 + 0.0i]x^28 + [1.0 + 0.0i]x^27 + [1.0 + 0.0i]x^26 + [1.0 + 0.0i]x^25 + [1.0 + 0.0i]x^24 + [1.0 + 0.0i]x^23 + [1.0 + 0.0i]x^22 + [1.0 + 0.0i]x^21 + [1.0 + 0.0i]x^20 + [1.0 + 0.0i]x^19 + [1.0 + 0.0i]x^18 + [1.0 + 0.0i]x^17 + [1.0 + 0.0i]x^16 + [1.0 + 0.0i]x^15 + [1.0 + 0.0i]x^14 + [1.0 + 0.0i]x^13 + [1.0 + 0.0i]x^12 + [1.0 + 0.0i]x^11 + [1.0 + 0.0i]x^10 + [1.0 + 0.0i]x^9 + [1.0 + 0.0i]x^8 + [1.0 + 0.0i]x^7 + [1.0 + 0.0i]x^6 + [1.0 + 0.0i]x^5 + [1.0 + 0.0i]x^4 + [1.0 + 0.0i]x^3 + [1.0 + 0.0i]x^2 + [1.0 + 0.0i]x + [1.0 + 0.0i]
Real: -1.00000000000000E0
Comp: -9.86361303402722E-1 +- 1.64594590280734E-1i
Comp: -9.45817241700635E-1 +- 3.24699469204684E-1i
Comp: -8.79473751206489E-1 +- 4.75947393037074E-1i
Comp: -7.89140509396394E-1 +- 6.14212712689668E-1i
Comp: -6.77281571625741E-1 +- 7.35723910673132E-1i
Comp: -5.46948158122427E-1 +- 8.37166478262529E-1i
Comp: -4.01695424652969E-1 +- 9.15773326655057E-1i
Comp: -2.45485487140799E-1 +- 9.69400265939330E-1i
Comp: -8.25793454723323E-2 +- 9.96584493006670E-1i
Comp: 8.25793454723323E-2 +- 9.96584493006670E-1i
Comp: 2.45485487140799E-1 +- 9.69400265939330E-1i
Comp: 4.01695424652970E-1 +- 9.15773326655058E-1i
Comp: 5.46948158122427E-1 +- 8.37166478262528E-1i
Comp: 6.77281571625741E-1 +- 7.35723910673132E-1i
Comp: 7.89140509396394E-1 +- 6.14212712689668E-1i
Comp: 8.79473751206489E-1 +- 4.75947393037074E-1i
Comp: 9.45817241700635E-1 +- 3.24699469204683E-1i
Comp: 9.86361303402722E-1 +- 1.64594590280734E-1i
Error: 1.36549533872468E-14
Poly: [1.0 + 3.0i]x^2 + [2.0 + 2.0i]x + [3.0 + 1.0i]
Comp: -8.00000000000000E-1 - 6.00000000000000E-1i
Comp: 0.00000000000000E0 + 1.00000000000000E0i
Error: < 1.0E-15
Poly: [1.0 + 1.0i]x^3 + [2.0 + 1.0i]x^2 + [3.0 + 1.0i]x + [4.0 + 1.0i]
Comp: -1.40135938330275E0 + 2.88269653138075E-1i
Comp: -2.84985631785343E-1 - 1.30378640290474E0i
Comp: 1.86345015088091E-1 + 1.51551674976666E0i
Error: < 1.0E-15
Poly: [1.0 + 1.0i]x^3 + [2.0 + 2.0i]x^2 + [3.0 + 3.0i]x + [4.0 + 4.0i]
Real: -1.65062919143939E0
Comp: -1.74685404280306E-1 +- 1.54686888723140E0i
Error: < 1.0E-15
License
This is free software; you can redistribute it and/or modify it under terms of the GNU Lesser General Public License (LGPL) as published by the Free Software Foundation, either version 2.1, or (at your option) any later version. It is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the LGPL for more details. A copy of the LGPL, incorporated by reference, may be obtained from gnu.org.
Click here to view the source code. Copyright 2007 John B. Matthews
Distributed without warranty under the terms of the GNU Lesser General Public License.
Last updated 9-Mar-2009