Raíces de un polinomio. Método de Laguerre

El polinomio

a₀+a₁x+a₂x²+...a_nxⁿ=0

de grado n tiene n raíces que pueden ser reales o complejas y pueden no ser distintas.Cuando los coeficientes del polinomio son reales, las raices complejas aparecen por pares x₁=a+bi es una raíz y su conjugada x₂=a-bi es también raíz. Pero si los coeficentes del polinomio son complejos las raíces complejas no están necesariamente relacionadas.

En esta página se proporciona el código del procedimiento de Laguerre.

Habitualmente, el polinomio tiene coeficientes reales y se escribe de mayor a menor grado, por ejemplo la ecuación

x⁴+0x³+4x²-3x+3=0

Se crea el array de coeficientes

 double[] coef={1, 0, 4, -3, 3};   
	calculaRaices(coef);

y se llama a la función calculaRaices que se encarga de crear el array de coeficientes complejos a[0],a[1], a[2]... a[n] en orden inverso.

Para el caso general, el lector puede crear directamente el array de coeficientes complejos, en la forma en la que se señala en la función main entre delimitadores de comentarios /* ... */

public class Complejo {
   public double r;
   public double i;
  public Complejo(double r, double i) {
    this.r=r;
    this.i=i;
  }
  public Complejo() {
    this.r=0.0;
    this.i=0.0;
  }

public static Complejo add(Complejo a, Complejo b)   {
  double real=a.r+b.r;
  double imag=a.i+b.i;
	return new Complejo(real, imag);
}

public static Complejo sub(Complejo a, Complejo b) {
  double real=a.r-b.r;
  double imag=a.i-b.i;
	return new Complejo(real, imag);
}

public static Complejo mul(Complejo a, Complejo b)  {
	double real=a.r*b.r-a.i*b.i;
	double imag=a.i*b.r+a.r*b.i;
	return new Complejo(real, imag);
}

public static Complejo conjg(Complejo z){
	return new Complejo(z.r, -z.i);
}

public static Complejo div(Complejo a, Complejo b) {
	Complejo c=new Complejo();
  double real, imag;
	double r,den;
	if (Math.abs(b.r) >= Math.abs(b.i)) {
		r=b.i/b.r;
		den=b.r+r*b.i;
		real=(a.r+r*a.i)/den;
		imag=(a.i-r*a.r)/den;
	} else {
		r=b.r/b.i;
		den=b.i+r*b.r;
		real=(a.r*r+a.i)/den;
		imag=(a.i*r-a.r)/den;
	}
	return new Complejo(real, imag);
}

public static double abs(Complejo z){
	double x,y,ans,temp;
	x=Math.abs(z.r);
	y=Math.abs(z.i);
	if (x == 0.0)
		ans=y;
	else if (y == 0.0)
		ans=x;
	else if (x > y) {
		temp=y/x;
		ans=x*Math.sqrt(1.0+temp*temp);
	} else {
		temp=x/y;
		ans=y*Math.sqrt(1.0+temp*temp);
	}
	return ans;
}

public static Complejo sqrt(Complejo z){
	double real, imag;
	double x,y,w,r;
	if ((z.r == 0.0) && (z.i == 0.0)) {
		return new Complejo();
	} else {
		x=Math.abs(z.r);
		y=Math.abs(z.i);
		if (x >= y) {
			r=y/x;
			w=Math.sqrt(x)*Math.sqrt(0.5*(1.0+Math.sqrt(1.0+r*r)));
		} else {
			r=x/y;
			w=Math.sqrt(y)*Math.sqrt(0.5*(r+Math.sqrt(1.0+r*r)));
		}
		if (z.r >= 0.0) {
			real=w;
			imag=z.i/(2.0*w);
		} else {
			imag=(z.i >= 0) ? w : -w;
			real=z.i/(2.0*imag);
		}
		return new Complejo(real, imag);
	}
}

public static Complejo mul(double x, Complejo a){
	return new Complejo(x*a.r, x*a.i);
}

public String toString(){
     if(i>0)     return new String((double)Math.round(1000*r)/1000+" 
     + "+(double)Math.round(1000*i)/1000+"*i");
     return new String((double)Math.round(1000*r)/100+" - "
     +(double)Math.round(-1000*i)/100+"*i");
}

} 
public class Laguerre {
  static final double EPSS=1.0e-7;
  static final int MR= 8;
  static final int MT= 10;
  static final int MAXIT=(MT*MR);
  static  final double EPS=2.0e-6;
  static final int MAXM= 100;

  public Polinomio() {
  }
//ojo al pasar el objeto x, queremos que 
//modifique sus miembros dato r e i en el curso de la llamada a la función

private static void laguer(Complejo[] a, int m, Complejo x){
	int iter,j;
	double abx,abp,abm,err;
	Complejo dx,x1,b,d,f,g,h,sq,gp,gm,g2;
	final double[] frac = {0.0,0.5,0.25,0.75,0.13,0.38,0.62,0.88,1.0};

	for (iter=1;iter<=MAXIT;iter++) {
		b=a[m];
		err=Complejo.abs(b);
		d=f=new Complejo(0.0,0.0);
		abx=Complejo.abs(x);
		for (j=m-1;j>=0;j--) {
			f=Complejo.add(Complejo.mul(x,f),d);
			d=Complejo.add(Complejo.mul(x,d),b);
			b=Complejo.add(Complejo.mul(x,b),a[j]);
			err=Complejo.abs(b)+abx*err;
		}
		err *= EPSS;
		if (Complejo.abs(b) <= err){
     return ;
    }
		g=Complejo.div(d,b);
		g2=Complejo.mul(g,g);
		h=Complejo.sub(g2,Complejo.mul(2.0,Complejo.div(f,b)));
		sq=Complejo.sqrt(Complejo.mul((double) (m-1),
            Complejo.sub(Complejo.mul((double) m,h),g2)));
		gp=Complejo.add(g,sq);
		gm=Complejo.sub(g,sq);
		abp=Complejo.abs(gp);
		abm=Complejo.abs(gm);
		if (abp < abm) gp=gm;
		dx=((FMAX(abp,abm) > 0.0 ? Complejo.div(new Complejo((double) m,0.0),gp)
			: Complejo.mul(1+abx,new Complejo(Math.cos((double)iter),
                 Math.sin((double)iter)))));
		x1=Complejo.sub(x,dx);
		if (x.r == x1.r && x.i == x1.i){
     return ;
    }
		if (iter % MT!=0){
        x.r=x1.r;
        x.i=x1.i;
    }	else{
     Complejo temp=Complejo.sub(x,Complejo.mul(frac[iter/MT],dx));
     x.r=temp.r;
     x.i=temp.i;
    }
 	}
	System.out.println("too many iterations in laguer");
	return;
}
private static  double FMAX(double a, double b){
    return ((a>b) ?a: b);
}


private static Complejo[] zroots(Complejo[] a, int m, boolean polish){
	int i, j,jj;
	Complejo x,b,c;
  Complejo[] roots=new Complejo[a.length];
	for (i=0;i<a.length;i++){
        roots[i]=new Complejo();
  }
  Complejo[] ad=new Complejo[MAXM];
  for(j=0; j<MAXM; j++){
    ad[j]=new Complejo();
  }

	for (j=0;j<=m;j++) ad[j]=a[j];
	for (j=m;j>=1;j--) {
		x=new Complejo();
		laguer(ad,j, x);
		if (Math.abs(x.i) <= 2.0*EPS*Math.abs(x.r)) x.i=0.0;
		roots[j]=x;
		b=ad[j];
		for (jj=j-1;jj>=0;jj--) {
			c=ad[jj];
			ad[jj]=b;
			b=Complejo.add(Complejo.mul(x,b),c);
		}
	}
	if (polish)
		for (j=1;j<=m;j++)
			laguer(a,m,roots[j]);
	for (j=2;j<=m;j++) {
		x=roots[j];
		for (i=j-1;i>=1;i--) {
			if (roots[i].r <= x.r) break;
			roots[i+1]=roots[i];
		}
		roots[i+1]=x;
	}
  return roots;
}
public static void calculaRaices(double[] coef){
        Complejo[] a=new Complejo[coef.length];
        for(int i=0; i<coef.length; i++){
            a[coef.length-1-i]=new Complejo(coef[i], 0.0);
        }
	      Complejo[] roots=zroots(a,a.length-1, true);
	      for (int i=1;i<roots.length;i++)
		    System.out.println(i+"   "+roots[i].r+"   "+roots[i].i);
 }
public static void calculaRaices(Complejo[] a){
	      Complejo[] roots=zroots(a,a.length-1, true);
	      for (int i=1;i<roots.length;i++)
		    System.out.println(i+"   "+roots[i].r+"   "+roots[i].i);
 }
}
public class MyClass1 {
  public static void main(String[] args) {
//Ecuación con coeficientes reales
//a[n]·x^n+a[n-1]·x^(n-1)+.....+a[3]·x^3+a[2]·x^2+a[1]·x+a[0]=0
           double[] coef={1, 0, 4, -3, 3};
           Laguerre.calculaRaices(coef);
//Ecuación con coeficientes complejos en orden inverso
//a[0]+a[1]·x+a[2]·x^2+a[3]·x^3+...=0
  /*Complejo[] a={new Complejo(3,0), new Complejo(-3,0), new Complejo(4,0), 
       new Complejo(0,0), new Complejo(1,0)};
  Laguerre.calculaRaices(a);
   */
}
}

Referencias

Press W. H., Teukolsky S. A., Vetterling W. T., Flannery B. P. Numerical Recipes in C, Second edition. Root Finding and Nonlinear Sets of Equations. Sección 9.5. Roots of Polynomials. Cambridge University Press. Código en C adaptado por el autor al lenguaje Java.