The code below is my implementation of finding roots using Durand-Kerner's method more information can be obtained at wikipedia.
The source code is available at the bottom of the page.
/**A class that defines a complex number
*@author Tan Hong Cheong
*@version 20101202
*/
class ComplexNumber
{
/**constructor
*creates an imaginary number with real and imaginary component 0
*/
public ComplexNumber()
{
real = 0;
imaginary = 0;
}
/**copy constructor
*@param rhs the complex number to be copied
*/
public ComplexNumber(ComplexNumber rhs)
{
this.real = rhs.real;
this.imaginary = rhs.imaginary;
}
/**constructor
*@param real the real component
*@param imgainary the imaginary component
*/
public ComplexNumber(double real, double imaginary)
{
this.real = real;
this.imaginary = imaginary;
}
/**set the complex number to be the same value as rhs
*@param rhs the complex number to be copied
*/
public void set(ComplexNumber rhs)
{
this.real = rhs.real;
this.imaginary = rhs.imaginary;
}
/**@return the conjugate of this complex number
*/
public ComplexNumber congugate()
{
return new ComplexNumber(real,-imaginary);
}
/**return the complex number this+rhs
*@param rhs the complex number
*/
public ComplexNumber add(ComplexNumber rhs)
{
return new ComplexNumber(this.real+rhs.real,this.imaginary+rhs.imaginary);
}
/**return the complex number this-rhs
*@param rhs the complex number
*/
public ComplexNumber minus(ComplexNumber rhs)
{
return new ComplexNumber(this.real-rhs.real,this.imaginary-rhs.imaginary);
}
/**return the complex number this*a
*@param a the real value
*/
public ComplexNumber times(double a)
{
return new ComplexNumber(this.real*a,this.imaginary*a);
}
/**return the complex number this*rhs
*@param rhs the complex number
*/
public ComplexNumber times(ComplexNumber rhs)
{
double x = this.real*rhs.real-this.imaginary*rhs.imaginary;
double y = this.real*rhs.imaginary+rhs.real*this.imaginary;
return new ComplexNumber(x,y);
}
/**return the complex number this/rhs
*@param rhs the complex number
*/
public ComplexNumber divide(ComplexNumber rhs)
{
double l = rhs.real*rhs.real+rhs.imaginary*rhs.imaginary;
double x = (this.real*rhs.real+this.imaginary*rhs.imaginary)/l;
double y = (rhs.real*this.imaginary-this.real*rhs.imaginary)/l;
return new ComplexNumber(x,y);
}
/**@return true if this complex number equals to rhs
*/
public boolean equals(Object o)
{
try
{
ComplexNumber rhs = (ComplexNumber) o;
return (this.real==rhs.real) && (this.imaginary==rhs.imaginary);
}
catch(ClassCastException e)
{
return false;
}
}
/**@return true if this.real==real and this.imaginary==imaginary
*/
public boolean equals(double real, double imaginary)
{
return this.real==real && this.imaginary==imaginary;
}
/**return the length of the complex number
*/
public double length()
{
return Math.sqrt(real*real+imaginary*imaginary);
}
/**@return the argument of the complex number
*/
public double argument()
{
return Math.atan2(imaginary,real);
}
/**
*the real component
*@param read the read value
*/
public void setReal(double real)
{
this.real = real;
}
/**
*@return the read value
*/
public double getReal()
{
return real;
}
/**
*the imaginary component
*@param imaginary the imaginary value
*/
public void setImaginary(double imaginary)
{
this.imaginary = imaginary;
}
/**
*@return the imaginary value
*/
public double getImaginary()
{
return imaginary;
}
/**print a string represntation of it
*/
public String toString()
{
String s = ""+real;
if (imaginary<0)
{
s = s+"-i"+(-imaginary);
}
else
{
s = s+"+i"+imaginary;
}
return s;
}
/**
*the imaginary component
*/
private double imaginary;
/**
*the real component
*/
private double real;
public static void main(String[] args)
{
ComplexNumber numerator = new ComplexNumber(-4,0);
ComplexNumber denominator = new ComplexNumber(0.34199999999999975,1.9169999999999998);
System.out.println("numerator="+numerator);
System.out.println("denominator="+denominator);
System.out.println("fraction="+numerator.divide(denominator));
}
}
/**a class that defines a polynomial with complex coefficients
*modified frm http://sites.google.com/site/drjohnbmatthews/polyrots
*@author Tan Hong Cheong
*@version 20101202
*/
public class ComplexPolynomial
{
/**Constructor
*@param coefficients the array of coefficients with the highest power at coefficientss[0]
*/
public ComplexPolynomial(ComplexNumber[] coefficients)
{
a = new ComplexNumber[coefficients.length];
for(int i=0;i<a.length;i++)
{
a[i] = new ComplexNumber(coefficients[i]);
}
}
/**constructor, create a polynomial with real coefficients
*@param coefficients the array of coefficients with the highest power at coefficientss[0]
*/
public ComplexPolynomial(double... coefficients)
{
a = new ComplexNumber[coefficients.length];
for(int i=0;i<a.length;i++)
{
a[i] = new ComplexNumber(coefficients[i],0);
}
}
/**Evaluate the polynomial at x using Horner's scheme.
*@return the Complex value of the function at x.
*/
public ComplexNumber evaluate(ComplexNumber x)
{
ComplexNumber result = a[0];
for(int i=1;i<a.length;i++)
{
result = result.times(x).add(a[i]);
}
return result;
}
/**
*This implementation uses the Durand-Kerner-Weierstrass method
*to find the rots of a polynomial with complex coefficients.
*The method requires a monic polynomial; some errr may occur
*when dividing by the leading coefficient.
*The array should have the highest order coefficient first.
*
*@param max the maximum iteration
*@param errr the errr that is acceptable
*@return an array of the Complex rots of the polynomial.
*/
public ComplexNumber[] roots(int max,double errr)
{
ComplexNumber one = new ComplexNumber(1,0);
if (!a[0].equals(1,0))
{
//copy the coefficients so that it will not be changed
ComplexNumber[] ca = new ComplexNumber[a.length];
for(int i=0;i<a.length;i++)
{
ca[i] = new ComplexNumber(a[i]);
}
for (int i=1;i<ca.length;i++)
{
ca[i] = ca[i].divide(ca[0]);
}
ca[0].setReal(1);
ca[0].setImaginary(0);
ComplexPolynomial poly = new ComplexPolynomial(ca);
return poly.roots(max,errr);
}
//the polynomial is monix
//i.e coefficient of x^n is 1
ComplexNumber[] r = new ComplexNumber[a.length - 1];//the array of rots
//Initialize r
//the n-1 rots are
//(0.4, 0.9)^0
//(0.4, 0.9)^1
//(0.4, 0.9)^2
//(0.4, 0.9)^3
//(0.4, 0.9)^4
//etc
//no particular reason to choose (0.4, 0.9)
ComplexNumber multiplier = new ComplexNumber(0.4, 0.9);
r[0] = new ComplexNumber(1,0);
for (int i=1;i<r.length;i++)
{
r[i] = r[i-1].times(multiplier);
}
//Iteration
int count = 0;
boolean done = false;
while(!done)
{
boolean unchanged = true;
/*uncomment to see the roots for each iterations
System.out.println("Iteration = "+count);
for(int i=0;i<r.length;i++)
{
System.out.println(r[i]);
}
*/
for(int i=0;i<r.length;i++)
{
ComplexNumber numerator = evaluate(r[i]);
ComplexNumber denominator = new ComplexNumber(1,0);
for (int j=0;j<r.length;j++)
{
if (i!=j)
{
denominator = r[i].minus(r[j]).times(denominator);
}
}
ComplexNumber newRi = r[i].minus(numerator.divide(denominator));
double realDiff = r[i].getReal()-newRi.getReal();
double imgDiff = r[i].getImaginary()-newRi.getImaginary();
if (realDiff<0)
{
realDiff = -realDiff;
}
if (imgDiff<0)
{
imgDiff= -imgDiff;
}
if ((realDiff>errr)||(imgDiff>errr))
{
unchanged = false;
}
//replace r[i] with new value
r[i].set(newRi);
}
count++;
done = (count>max||unchanged);
if (count>max)
{
System.out.println("Exceeded");
}
}
return r;
}
/**the coefficients the array of coefficients with the highest power at coefficientss[0]
*/
private ComplexNumber[] a;
public static void main(String[] args)
{
{
ComplexPolynomial poly = new ComplexPolynomial(1,-3,3,-5);
ComplexNumber[] roots = poly.roots(999,1E-15);
System.out.println("roots are:");
for(ComplexNumber root:roots)
{
System.out.println(root);
}
}
{
ComplexNumber[] coeff = new ComplexNumber[4];
coeff[0] = new ComplexNumber(1,1);
coeff[1] = new ComplexNumber(2,2);
coeff[2] = new ComplexNumber(3,3);
coeff[3] = new ComplexNumber(4,4);
ComplexPolynomial poly = new ComplexPolynomial(coeff);
ComplexNumber[] roots = poly.roots(999,1E-15);
System.out.println("roots are:");
for(ComplexNumber root:roots)
{
System.out.println(root);
}
}
{
ComplexNumber[] coeff = new ComplexNumber[6];
coeff[0] = new ComplexNumber(1,0);
coeff[1] = new ComplexNumber(-13.999,-5);
coeff[2] = new ComplexNumber(74.99,55.998);
coeff[3] = new ComplexNumber(-159.959,-260.982);
coeff[4] = new ComplexNumber(1.95,463.934);
coeff[5] = new ComplexNumber(150,-199.95);
ComplexPolynomial poly = new ComplexPolynomial(coeff);
ComplexNumber[] roots = poly.roots(999,1E-15);
System.out.println("roots are:");
for(ComplexNumber root:roots)
{
System.out.println(root);
}
}
}
}