Programa Policeros
Si lo que se quiere es realizar un programa computacional que utilice el Teorema de Sturm 4 para aproximar los ceros de un polinomio se deben seguir una serie de pasos que trataremos de enunciar a continuación, nosotros realizamos el programa en Java, sin embargo se puede programar en cualquier otro lenguaje, aquí la idea no es realizar un programa completo donde se muestre paso a paso como realizarlo, sino mostrar las ideas principales sobre cómo se debe realizar, es decir, que no se va a mostrar el código del programa completo, sino solo algunas de las funciones que son necesarias, el programa lo puede bajar al final de este artículo para analizarlo y tenerlo completo.
Un primer problema que se nos puede presentar es cómo leer el polinomio, leer un polinomio como una tira de caracteres y transformarlo a algo que la computadora entienda no es algo trivial, a esto se le llama parsear una expresión y ésto será algo de lo que escribiré en otra ocasión, ahora se debe buscar una forma sencilla, se puede poner una caja de texto para cada coeficiente que se quiera leer y hacer el polinomio de un tamaño fijo (otra forma es preguntar el tamaño primero y después construir las cajas de texto que se requieran).
La mejor forma de guardar el polinomio es en un arreglo de números, en la posición cero se guarda el coeficiente que está solo (es decir que va multiplicado por x0), en la posición 1 se guarda el coeficiente que acompaña a x1 y así consecutivamente hasta guardar todo el polinomio.
Una vez que se ha guardado el polinomio deseado, se necesitan varias funciones para trabajar con él. La primera función será la que evalúe un valor de x en el polinomio, esta función debe recibir el arreglo del polinomio y devolver la evaluación, para esta simplemente hay que sumar los término multiplicando cada uno por xn, donde n es la posición en el arreglo (hay otras formas mejores para hacerlo, pero la idea aquí no es mostrar la mejor forma, sino la más simple). Veamos el código:
/*La funcion que sigue es la que evalua el polinomio*/
public double f(Vector fun, double x) {
double total = (double)((Double)fun.elementAt(0)).doubleValue();
for(int i = 1; i < fun.size(); i++){
total += ((double)((Double)fun.elementAt(i)).doubleValue())*Math.pow(x, i*1.0);
}
return total;
}//Fin de la funcion f
Una segunda función calcula la derivada de un polinomio, esta también recibe un arreglo con el polinomio y devuelve otro arreglo con su derivada, es muy sencillo pues cada valor se debe multiplicar por su posición en el arreglo y se elimina el primer término (el que está multiplicado por x0).
//Funcion que deriva un polinomio
public Vector derivePoli(Vector fun) {
Vector derivado = new Vector();
for(int i = 1; i < fun.size(); i++)
derivado.addElement(new Double((double)((Double)fun.elementAt(i)).doubleValue()*i));
derivado.trimToSize();
return derivado;
}//derivePoli
La tercera función es la que calcula la división de dos polinomios, esta no es tan sencilla, se deben recibir dos arreglos que contengan el dividendo y el divisor (son polinomios) y se puede hacer que solo devuelva el cociente, el residuo o ambos, en este caso sólo se necesita el residuo, así que se declara un arreglo para guardarlo, se realiza la división entre los últimos coeficientes del cociente y el divisor (recuerde que los guardamos empezando con el que acompaña a x0, luego x1,etc), este valor se multiplica por cada elemento del divisor y se le suma a su correspondiente elemento del dividendo, estos valores se van guardando en el arreglo declarado.
//Funcion que divide dos polinomios y devuelve el residuo
//Si no se puede dividir devuelve el vector vacío
public Vector dividaPolis(Vector dividendo, Vector divisor) {
Vector dividendoTemp = new Vector();
Vector cocienteTemp = new Vector();
Vector residuoTemp = new Vector();
dividendoTemp = (Vector)dividendo.clone();
//el grado del dividendo debe ser mayor o igual que el grado del divisor
if(dividendo.size() >= divisor.size()){
//Cantidad de veces que se debe hacer el algoritmo
for(int i = 0; i <= (dividendo.size()-divisor.size()); i++){
if(residuoTemp.size() > 0)
residuoTemp.removeAllElements();
cocienteTemp.addElement(new Double((double)((Double)dividendoTemp.elementAt(
dividendoTemp.size()-1)).doubleValue()/(double)((Double)divisor.elementAt(
divisor.size()-1)).doubleValue()));
//Hace la división de los mas grandes para ir llenando el cociente
for(int j = 2; j <= dividendoTemp.size(); j++){
if (j <= divisor.size())
residuoTemp.addElement(new Double((double)((Double)dividendoTemp.elementAt(
dividendoTemp.size()-j)).doubleValue()-((double)((Double)divisor.elementAt(
divisor.size()-j)).doubleValue()*(double)((Double)cocienteTemp.elementAt(i)
).doubleValue())));
else
residuoTemp.addElement(new Double((double)((Double)dividendoTemp.elementAt(
dividendoTemp.size()-j)).doubleValue()));
}//for
residuoTemp = volverVector(residuoTemp);
dividendoTemp = (Vector)residuoTemp.clone();
}//for
cocienteTemp = volverVector(cocienteTemp);
residuoTemp = quitaCeros(residuoTemp);
if(residuoTemp.size() == 0)
residuoTemp.addElement(new Double(0.0));
}//if
residuoTemp.trimToSize();
cocienteTemp.trimToSize();
return residuoTemp; //Aquí se puede retornar el cociente o el residuo cambiando la palabra
}//dividaPolis
Se podría hacer otra función para simplificar algunos cálculos, esta recibe un polinomio y le quita los ceros que se encuentran al final, es decir, que el coeficiente que acompaña a la variable de mayor grado sea diferente de cero.
//Procedimiento que elimina el número de grado mayor si es 0
public Vector quitaCeros(Vector fun){
Vector sinCeros = new Vector();
boolean salir = false; //para ver si encontró un valor de cierto grado o si no poner 0
int i = fun.size()-1;
sinCeros = (Vector)fun.clone();
while(i >= 0 && salir == false){
if((double)((Double)sinCeros.elementAt(i)).doubleValue() == 0.0)
sinCeros.removeElementAt(i);
else
salir = true;
i--;
}
sinCeros.trimToSize();
return sinCeros;
}//quitaCeros
Otra función que puede ser necesaria es la que calcula - f (x), es decir, una función que le cambie el signo a todos los coeficientes del polinomio. También es necesaria la que calcula f (- x), esta función debe cambiarle el signo a los campos impares (1, 3, 5, 7, ...).
//Función que le cambia el signo a todas las entradas del polinomio
public Vector cambieSignoPoli(Vector fun) {
Vector signoCambiado = new Vector();
for(int i = 0; i < fun.size(); i++)
signoCambiado.addElement(new Double(-1*(double)((Double)fun.elementAt(i)).doubleValue()));
signoCambiado.trimToSize();
return signoCambiado;
}//cambieSignoPoli
//Funcion que le cambia el signo a las potencias impares del polinomio
public Vector cambiaSignoImpares(Vector fun){
Vector nuevafun = new Vector();
for(int i = 0; i < fun.size(); i++){
if(i%2 == 0)
nuevafun.addElement((Double)fun.elementAt(i));
else
nuevafun.addElement(new Double(-1*((double)((Double)fun.elementAt(i)).doubleValue())));
}
return nuevafun;
}//cambiaSignoImpares
Para el teorema de Sturm puede ser útil realizar una función que devuelva el signo de algún valor que recibe, si es negativo puede devolver -1, si es positivo devuelve 1 y si es cero entonces devuelve ese valor, es decir 0.
//Busca la cota superior de lo ceros positivos de una funcion
public double cota(Vector fun){
double valorCota = 0, mayorCoef;
mayorCoef = Math.abs((double)(((Double)fun.elementAt(0)).doubleValue()));
for(int i = 1; i < fun.size()-1; i++)
mayorCoef = Math.max(mayorCoef, Math.abs((double)(((Double)fun.elementAt(i)).doubleValue())));
if(Math.abs((double)(((Double)fun.elementAt(fun.size()-1)).doubleValue())) == 0)
valorCota = 0;
else
valorCota = 1+(mayorCoef/Math.abs((double)(((Double)fun.elementAt(fun.size()-1)).doubleValue())));
return valorCota;
}//cota
Falta también la función (esta ya usa las anteriores) que calcula el sistema de Sturm para un polinomio dado, esta es bastante simple si ya tenemos la funciones anteriores. Otra función es la que evalúa un valor en el sistema de Sturm y devuelve la cantidad de cambios de signo que hubieron.
/*La función que sigue es la que hace el teorema de sturm
primero deriva el polinomio, divide el polinomio entre la derivada
y recoge el residuo, luego divide la derivada entre el residuo
y nuevamente recoge el residuo y sigue dividiendo los dos últimos
residuos para obtener otro residuo hasta que ya no se pueda.
Luego evalúa las funciones en -infinito y +infinito y con la resta
de los cambios de signo se sabe cuantos ceros hay...*/
public Vector sturm(Vector fun) {
boolean salir = false;
int cerosPositivos = 0, cerosNegativos = 0;
Vector funciones = new Vector();
Vector raices = new Vector();
Vector cotas = new Vector();
if(fun.size() > 1){
fun = quitaFallo(fun);
if(fun.size() > 1){
//Se buscan las funciones del sistema de sturm
funciones.addElement(fun);
funciones.addElement(derivePoli(fun));
int i = 0;
while(!salir && i < (fun.size()-2)){
funciones.addElement(cambieSignoPoli(dividaPolis((Vector)funciones.elementAt(i),
(Vector)funciones.elementAt(i+1))));
if(((Vector)funciones.elementAt(funciones.size()-1)).size() == 1){
if(((double)(((Double)(((Vector)funciones.elementAt(funciones.size()-1))).
elementAt(0))).doubleValue()) == 0){
salir = true;
}
}
i++;
}
//Se buscan los intervalos con el sistema de sturm
cerosPositivos = numeroCerosPositivos(funciones);
cerosNegativos = numeroCerosNegativos(funciones);
cotas = busqueCotas((Vector)funciones.elementAt(0));
raices = unaVectores(busqueRaices(funciones, cerosNegativos, (double)(((Double)cotas.
elementAt(0)).doubleValue()), (double)(((Double)cotas.elementAt(1)).doubleValue())),
busqueRaices(funciones, cerosPositivos, (double)(((Double)cotas.elementAt(2)).doubleValue()),
(double)(((Double)cotas.elementAt(3)).doubleValue())));
if(falloSturm)
raices.addElement(new Double(0.0));
}
else
raices.addElement(new Double(0.0));
}
raices.trimToSize();
return raices;
}//sturm
En este programa se presentan algunas funciones de las que no hemos hablado, por ejemplo, el número de ceros positivos y negativos, estas funciones simplemente se evalúa el sistema de Sturm de cero a infinito para saber los ceros positivos y de menos infinito a cero para los negativos. busqueRaices es la función que se encarga de ir acotando las raices y buscarlas. La función quitaFallo lo que hace es eliminar la multiplicidad del cero
//Función que busca los intervalos donde hay ceros
public Vector busqueRaices(Vector fun, int numCeros, double minimo, double maximo){
Vector raiz = new Vector();
int cerosEncontrados = 0;
double ci = minimo, cs = maximo, cotaux;
cotaux = (ci+cs)/2;
if(numCeros == 1)
raiz.addElement(new Double(miRedondeo(buscarCero((Vector)fun.elementAt(0), ci, cs))));
else{
while(cerosEncontrados < numCeros){
if(Math.abs(cambiosDeSigno(fun, cotaux)-cambiosDeSigno(fun, ci)) == 0){
ci = cotaux;
cotaux = (cs+cotaux)/2;
}
else if(Math.abs(cambiosDeSigno(fun, ci)-cambiosDeSigno(fun, cotaux)) > 1)
cotaux = (ci+cotaux)/2;
else if(Math.abs(cambiosDeSigno(fun, ci)-cambiosDeSigno(fun, cotaux)) == 1){
raiz.addElement(new Double(miRedondeo(buscarCero((Vector)fun.elementAt(0), ci, cotaux))));
ci = cotaux;
cotaux = (cs+ci)/2;
cerosEncontrados++;
}
}//while
}//if
raiz.trimToSize();
return raiz;
}//busqueraices
//Funcion que quita el fallo de Sturm, especificamente, la multiplicidad del cero
public Vector quitaFallo(Vector fun){
Vector nuevaFun = new Vector();
falloSturm = false;
int i = 0;
while((double)(((Double)fun.elementAt(i)).doubleValue()) == 0.0){
i++;
}
if(i > 0){
falloSturm = true;
while(i < fun.size()){
nuevaFun.addElement((Double)fun.elementAt(i));
i++;
}
}
else
nuevaFun = fun;
return nuevaFun;
}//quitaFallo
Por último, se debe hacer el programa principal que maneje todas las funciones anteriores, es decir que debe leer el polinomio, acomodarlo, calcular el sistema de sturm e imprimir los ceros.
public boolean action(Event evt, Object arg){
if (evt.target instanceof Button){
if (evt.target == panel.b1){
panel.lista.removeAll();
leaParametros();
funcion = leePoli(formulaPoli);
panel.info.setText("");
if(funcion.size() > 0){
funcion = acomodaPoli(funcion);
raicesPol = sturm(funcion);
imprimaCeros(raicesPol);
}
lienzo.repaint();
}
}//if del botón
return true;
}//action
Hay que hacer la funcion que reciba el polinomio y los valores iniciales y aproxime el cero, es decir, la funcion que realice el método de bisección, bisección acelerada o el de Newton.
//Buscar cero es la función que aproxima el cero, utiliza el método de bisección
public double buscarCero(Vector fun, double a, double b){
double xmax = b, xmin = a;
double xn = a;
boolean cero = false;
if(f(fun, xmin) == 0)
xn = xmin;
else if(f(fun, xmax) == 0)
xn = xmax;
else{
int i = 0;
while(!cero && i < 1000){
xn = (xmin+xmax)/2;
if(Math.abs(f(fun, xn)) == 0)
cero = true;
else if(f(fun, xn)*f(fun, xmin) < 0)
xmax = xn;
else
xmin = xn;
i++;
}//while
}
//Si este método no funciona se hace la llamada al método newton
if(Math.abs(f(fun, xn)) > 0.000001)
xn = newton(fun, a, b);
return xn;
}//buscarCero
//Buscar cero es la función que aproxima el cero, utiliza el método de Newton
public double newton(Vector fun, double a, double b){
double xn = (a+b)/2;
Vector derivada = derivePoli(fun);
boolean cero = false;
int i = 0;
while(!cero && i < 500){
xn = xn-(f(fun, xn)/f(derivada, xn));
if(Math.abs(f(fun, xn)) == 0)
cero = true;
i++;
}//while
return xn;
}//buscarCero
Ya para finalizar, se deben mostrar en pantalla los ceros obtenidos.
//Imprime los ceros en una lista
public void imprimaCeros(Vector raices){
if(raices.size() == 0)
panel.lista.addItem("No hay ceros reales.");
else
for(int i = 0; i < raices.size(); i++)
panel.lista.addItem("x = "+(double)((Double)raices.elementAt(i)).doubleValue());
}//imprimaCeros
Espero que estas reflexiones hayan servido para que puedan realizar un programa computacional eficaz y eficiente, cualquier duda se pueden comunicar a mi correo y con mucho gusto les puedo contestar y ayudar con algún problema que tengan. Más adelante voy a escribir como hacer un parseador simple para leer un polinomio y otro más complicado sobre como leer cualquier expresión metemática y evaluar un valor en ella. El programa que se puede bajar tiene mucho más que solo la búsqueda de ceros, pues grafica el polinomio y marca los ceros en la gráfica.
Revista Matemáticas, Educación
e Internetâ
Derechos Reservados