Método de Müller - Universidad EAFIT

MÉTODO DE MULLER PARA OBTENER TODAS LAS RAÍCES DE UN
POLINOMIO
RESUMEN
Se describe en detalle el algoritmo de Muller para obtener una raíz de una
función f (x) . Se discute la forma clásica (más eficiente) de evaluar polinomios
y de realizar la división sintética, procedimientos necesarios para hallar todas
las raíces de un polinomio de coeficientes reales y grado arbitrario. Los
resultados obtenidos mediante un programa desarrollado por los autores en
la Facultad de Ingeniería Eléctrica de la Universidad Tecnológica de Pereira
se comparan, en términos de precisión y tiempo de ejecución, , con los
encontrados con el MATLAB 5.3.
Alvaro Acosta Montoya
Profesor Titular
Ingeniería Eléctrica U. T. P
e-mail: [email protected]
Carlos Galván Ceballos
Estudiante X semestre
Ingeniería Eléctrica U. T. P
e-mail: [email protected]
Mario Acosta Acosta
Ingeniero Eléctricista U. T. P.
Colombiana Kimberly Colpapel S. A.
e-mail: [email protected]
ABSTRACT
The Muller’s Method to obtain a single root of a function f(x) is described in
detail. The most efficient (classical) approach to evaluate a polynomial and
to make synthetic division are discussed, which are necesary functions to get
all the roots of a polynomial of real coefficients and arbitrary degree. The results
obtained by a program developed by the authors in the Electrical Engineering
Department at Universidad Tecnológica de Pereira, are compared to those
obtained by using MATLAB 5.3
MÉTODO DE MULLER PARA HALLAR UNA RAÍZ DE LA
ECUACIÓN f (x) ' 0
Entre las propiedades más destacadas del método se pueden
mencionar[1]:
<
Se puede usar para encontrar cualquier número prefijado
de raíces.
<
Se aplica a funciones arbitrarias.
<
Encuentra raíces reales y/o complejas.
<
En el entorno de la raíz posee convergencia casi
cuadrática.
<
<
Se basa en aproximar en la proximidad de la raíz cualquier
función mediante una parábola. Para ello se escogen tres
puntos relativamente contiguos que pertenecen a la función
y se halla la ecuación de la única parábola que pasa por ellos.
Las soluciones de la ecuación cuadrática que resulta
determinan los “puntos” por donde dicha parábola corta el
eje x, uno de los cuales se toma como mejor aproximación
a la raíz que se busca. Los detalles de la deducción de las
fórmulas se muestran a continuación:
Sean x0 , x1 y x2 tres aproximaciones diferentes para una
raíz de la ecuación f (x) ' 0 . La ecuación de la única parábola
que pasa por los puntos [x0 , f (x0 )] [x1 , f(x1 )] [x2 , f(x2 )] se
supone de la forma:
p(x) ' A x 2 % B x % C
No requiere evaluación de la derivada de la función como
en el método de Newton
Cuando la función es un polinomio,[( f (x) ' p(x) ], y se
aproxima una raíz repetida no diverge como el método
de Bairstow o el de Newton1.
entonces se debe cumplir:
2
p (x0 ) ' f0 • f ( x0 ) ' A x0 % B x0 % C
p(x1 ) ' f1 • f (x1 ) '
p (x2 ) ' f2 • f ( x2 ) '
1
No es difícil verificar que cuando f(x) es un
polinomio p(x), el método de Newton en las
proximidades de una raíz repetida xk diverge
por cuanto p´(xk)=0.
(1)
2
A x1
2
A x2
(2a)
% B x1 % C
(2b)
% B x2 % C
(2c)
Combinando (2c) con (2b) y (2b) con (2a) para eliminar la
incógnita C se obtiene:
SCIENTIA ET TECHNICA No. 12
2
ABRIL /2000 80
2
A (x2 & x1 ) % B (x2 & x1 ) ' f2 & f1
2
2
A (x1 & x0 ) % B (x1 & x0 ) ' f1 & f0
(3)
la ecuación de la parábola también se puede expresar de la
siguiente manera:
p (x ) ' f2 % g (x2, x1 ) ( x & x2 ) %
%
que se puede expresar matricialmente de la siguiente manera:
2
2
( x2 & x1 ) ( x2 & x1 ) A
2
2
( x1 & x0 ) ( x1 & x0 ) B
'
g (x2, x1 ) & g( x1, x0 )
x2 & x0
( x & x2 )( x & x1 )
que permite verificar de manera inmediata que corresponde
a la ecuación de la única parábola que pasa por los puntos
[x0 , f (x0 )] [x1 , f(x1 )] [x2 , f(x2 )]
(f2 & f1 )
(f1 & f0 )
cuya solución es la siguiente:
A '
g( x2, x1 ) & g( x1, x0 )
x2 & x0
Un enfoque alternativo supone que la ecuación de la parábola
es de la forma:
(4a)
(4)
B ' g( x2, x1 ) & A (x2 % x1 )
(4b)
g (x1, x0 ) '
f2 & f1
x2 & x1
f1 & f0
c ' f (x0)
(5)
x1 & x0
a '
Reemplazando (4) en (2c) se obtiene:
C 'f2 & x2 g( x2, x1 ) & x1 A
(9)
y en este caso mediante un procedimiento similar al descrito
se obtienen las siguientes expresiones para las constantes
a, b y c
donde
g (x2, x1 ) '
p(x) ' a( x & x0 )2 % b (x & x0 ) % c
g (x1, x0) & g (x2, x0)
x1 & x2
(10)
b ' g(x1 , x0) & a (x1 & x0)
(6)
donde
En [1] se sugiere y se justifica que para evitar pérdida de
dígitos significativos en la solución a la ecuación cuadrática
p (x ) ' A x 2 % Bx % C ' 0 se haga mediante la fórmula
alternativa:
g ( x1 , x0 ) '
g ( x2 , x0 ) '
2C
x3 '
B ±
B 2 & 4AC
(7)
donde el signo antes del radical se debe escoger de tal
manera que el valor absoluto del denominador sea el más
grande lo que da la raíz de valor absoluto más pequeño.
x1 & x0
f2 & f0
(11)
x2 & x0
y se mantienen las definiciones dadas en (2a), (2b) y (2c) para
f0 , f1 y f2 , pero en este caso se utiliza (9) para p (x)
p(x0 ) ' f0 • f (x0 ) 'c
p(x1 ) ' f1 • f (x1 ) ' a(x1 &x0)2 % b(x1 &x0) % c
Si f ( x3 ) # å donde å es una tolerancia especificada de
antemano, se ha hallado una solución. En caso contrario se
repite el proceso con las siguientes asignaciones:
x2 ' x3
x1 ' x2
x0 ' x1
f1 & f0
p(x2 ) ' f2 • f (x2 ) ' a(x2 &x0)2 % b(x1 &x0) % c
(12)
ALGORITMO DE MULLER
(8)
Reemplazando (4a) en (4b) en (6) y las expresiones resultantes
para las constantes A , B y C en (1), se puede demostrar que
Se supone que los coeficientes del polinomio se almacenan
en el vector A , y que los subíndices más bajos identifican
la posiciones de las potencias mas altas.
1.
Suponer que los tres puntos iniciales contiguos se
encuentran separados entre sí INI ' ½ unidad..
SCIENTIA ET TECHNICA No. 12
Inicialmente se busca una raíz en el intervalo
& 0.5 # x # 0.5 Es decir,
x0 ' & INI
x1 ' 0
3.
Fijar tolerancia å ' 10 & 6 y máximo número de
iteraciones MAX ' 1000 e inicializar contador de
iteraciones j ' 0
Si el grado es igual a 1 hacer
a ' 0
b ' A[0]
c ' A[1]
x3 ' & c
b
4.
b ' A[1]
c ' A[2]
5.
Aplicar secuencialmente (12), (11), (10) para obtener las
constantes a , b y c
6.
Aplicar (7) con
A = a, B = b y C = c
para obtener
Evaluación eficiente de un polinomio
p (x) ' a0 % a1 x %a2 x 2 % þþ an x n
(13)
Es ineficiente evaluar cada uno de los n % 1 términos
separadamente y sumarlos después porque el esfuerzo
n (n % 1)
computacional es de n sumas y
multiplicaciones.
2
Aún si se toma como base x k para evaluar x k % 1 se
requerirían n sumas y (2n & 1) multiplicaciones.
Observando que con excepción del primero todos los
términos de (13) contienen el factor común x (12) se puede
escribir de la siguiente manera:
Si el grado es igual a 2
a ' A[0]
ASPECTOS COMPUTACIONALES
En los métodos iterativos para encontrar las raíces de un
polinomio se requiere con demasiada frecuencia la evaluación
del polinomio de grado entero n > 0 de la forma:
x0 ' INI
2.
ABRIL /2000 81
x3 - x0 , verificar si f (x3 ) # å. Si no se
cumple esta condición incrementar en 1 el número de
iteraciones j y si se excede el máximo permitido
MAX volver a iniciar el proceso reduciendo el valor deINI
a la mitad. Si tampoco se encuentra solución en MAX
iteraciones se reduce nuevamente el valor de INI a la
mitad y se reinicia el proceso. Finalmente si tampoco se
encuentra solución con INI ' 0.125 ' c después de MAX
iteraciones se imprime mensaje de no convergencia
después de 3(MAX iteraciones.
p(x) ' a0 % x(a1 % a2 x % þþ %
% an
&1
xn
&2
% anx n
&1
)
De nuevo cada término dentro del paréntesis excepto el
primero tienen el factor común x
p(x) ' a0 % x(a1 % a2 (x % þþ %
% an
&1
xn
&3
% anx n
&2
))
Continuando con este proceso se obtiene la expresión
anidada de p (x)
p(x) ' a0 % x(a1 % x (a2 % þþ %
% x(an & 1 % x(a n ))þþ))
para evaluar la cual se necesitann sumas yn multiplicaciones.
División Sintética
7.
Cuando se haya encontrado una raíz se debe proceder
a realizar la deflación (división sintética) para reiniciar
el proceso con un polinomio de orden menor. Debe
tenerse en cuenta que, puesto que el polinomio tiene
coeficientes reales, cuando ocurre una raíz compleja
también ocurre su conjugado y en ese caso la división
sintética se debe aplicar dos veces.
Dados los n % 1 coeficientes a 0 , þþ, a n del polinomio (12)
y el número z hacer
b n ' an
Para k ' n & 1, þþ, 0 hacer
bk ' ak % z bk % 1
(14)
Entonces aplicando (14) se obtiene mediante un proceso
inductivo
SCIENTIA ET TECHNICA No. 12
b 0 ' a0 % z b 1
' a0 % z (a1
' a0 % z (a1
!
b 0 ' p (z) ' a0
%
ABRIL /2000 82
% z b 2)
% z (a2 % zb 3))
! ! ! !
% z (a 1 % z (a2 % þþ %
z (an & 1 % z (an )) þþ ))
(15)
Los autores se vieron en la necesidad de desarrollarla porque
están trabajando en un paquete para la solución analítica de
circuitos eléctricos de parámetros concentrados. Se realizaron
pruebas y se obtuvieron resultados ligeramente mejores, tanto
en tiempo de ejecución como en precisión, con los arrojados
por el MATLAB 5.3. Los lectores interesados en hacer uso
de la librería completa pueden solicitarla a cualquiera de los
autores vía e-mail
La secuencia de números b n , b n & 1 , þþ, b 1 tienen una
propiedad muy útil que se deriva a continuación:
BIBLIOGRAFÍA
De (14)
(16)
[1] CONTE, S. D. y CARL de Boor, “Elementary
Numerical Analysis”, Mc-Graw Hill, 2ª edición, New
York, 1972
Por lo tanto si q (x) es un polinomio de grado n & 1 de la
forma
[2] FADEEV, D. K. y V. N. FADEEVA, “Computational
Methods of Linear Algebra”, W. H Freeman and
Company, San Francisco, 1963
ak '
k ' n
bk
bk & bk
% 1z
k < n
k ' 0, þþ, n
q(x) ' b 1 % b 2 x % þþ % b n x n
&1
[3] JOYANES, A. Luis, “C++ A SU ALCANCE: Un
enfoque orientado a objetos”, McGraw-Hill,
México, 1995
entonces
b 0 % q(x)(x & z) ' b 0 % (b 1 % b 2 x % þþ
& (b 1 % b 2 x % þ
Apéndice A
Codificación de Evaluación de polinomios
' (b 0 & b 1 z) % (b 1 & b 2 z)x %þþ
(b n
&1
& b n z)x n
' a0 % a1 x % þþ % an
&1
&1
% b0x n
xn
&1
(17)
%
' p (x)
Es decir, b 1 , þþ, b n & 1 , b n son los coeficientes del
polinomio resultante de dividir p (x) entre el polinomio de
primer orden x & z y b 0 es el residuo. Haciendo x ' z en
(16) se obtiene nuevamente b 0 ' p (z) [ver ecuación (15)].
Sip (z) # å x ' z esunaraízylasecuenciadecoeficientesb 1 , þþ, b n & 1 , b n
son los del polinomio base para continuar hallando todas
las demás raíces.
CONCLUSIONES
Se ha descrito en detalle el algoritmo de Muller para obtener
todas las raíces de un polinomio de coeficientes reales. Esta
resulta ser una herramienta útil en el análisis de sistemas
lineales. También sirve de base para realizar expansión en
fracciones parciales.
float eval(FPoly &A,float &x)
{
float resul;
resul=A.poly[0];//coeficiente de mayor potencia
for (register unsigned char j=1;j<=A.grado;j++) //limita a un
polinomio de grado 126
{
resul=resul*x+A.poly[j];
}
return resul;
}
La clase Fpoly ha sido creada para representar un polinomio
de coeficientes reales y grado arbitrarion y se le han definido,
además de los constructores de copia y sobrecarga de
operadores binarios y de flujo, funciones para encontrar sus
raíces, evaluar el polinomio y realizar la división sintética
Apéndice B
Codificación de la División Sintética
FPoly Div(FPoly &A,const float &x)
{
SCIENTIA ET TECHNICA No. 12
float resul;
Fpoly B(char(A.grado-1));
resul=A.poly[0];//coeficiente de mayor potencia
for (register unsigned char j=1;j<=A.grado;j++)
//limita a un orden del polinomio a 126
{
B[char(j-1)]=resul;
resul=resul*x+A.poly[j];
}
return B;
}
Apéndice C
Código para encontrar raíces de un polinomio
Roots roots(FPoly A)
{
Roots raices(A.grado);
const int MAX=1000;
const float TOL=1E-6,REAL=1E-5;
//Contador de las raices encontradas
unsigned char N;
const unsigned char NRAICES=char(A.grado-1);
complex a,b,c,x0,x1,x2,h1,h12,h2,D,f10,f20;
float INI;
float var;
bool FINDROOT;
N=0;
INI=0.5;
do
{
FINDROOT=false;
x0=-INI;
x1=0;
x2=INI;
for(short i=0;i<=MAX;i++)
{
if(A.grado==2)
{
bool r1;
r1=false;
D=A[1]*A[1]-4*A[0]*A[2];
D=sqrt(D);
ABRIL /2000 83
if(abs(A[1]+D)<abs(A[1]-D))
{
raices[N]=-2*A[2]/(A[1]-D);
r1=true;
}
else raices[N]=-2*A[2]/(A[1]+D);
if(fabs(imag(raices[N]))<=REAL)
{
if(r1) raices[char(N+1)]=-2*A[2]/(A[1]+D);
else raices[char(N+1)]=-2*A[2]/(A[1]-D);
raices[N]=real(raices[N]);
raices[char(N+1)]=real(raices[char(N+1)]);
}
else raices[char(N+1)]=conj(raices[N]);
N=char(N+2);
FINDROOT=true;
break;
}
else if(A.grado==1)
{
raices[N]=-A[1]/A[0];
N++;
FINDROOT=true;
break;
}
c=eval(A,x0);
h1=x1-x0;
h2=x2-x0;
h12=x1-x2;
f10=(eval(A,x1)-eval(A,x0))/h1;
f20=(eval(A,x2)-eval(A,x0))/h2;
a=(f10-f20)/h12;
b=f10-a*h1;
D=b*b-4*a*c;
D=sqrt(D);
if(abs(b+D)<abs(b-D))
{
D*=-1;
D+=b;
}
else D+=b;
raices[N]=x0-2*c/D;
if(abs(eval(A,raices[N]))<TOL)
SCIENTIA ET TECHNICA No. 12
{
//una vez se obtenga una raiz aplicamos la deflacion
FINDROOT=true;
if(fabs(imag(raices[N]))<=REAL)
{
//si la raiz es real
var=real(raices[N]);
if(N<NRAICES) A=Div(A,var);
raices[N]=var;
}
else {
//si la raiz es compleja
if(N+1<NRAICES) A=Div(A,raices[N]);
raices[char(N+1)]=conj(raices[N]);
N++;
}
N++;
INI=0.5;
break;
}
//si no se encuentra raiz preparamos la siguiente
iteracion
else
{
x0=x1;
x1=x2;
x2=raices[N];
}
}
// si no converge en MAX iteraciones,
cambiamos el valor de INI
if(!FINDROOT)
{
INI/=2;
if(INI<0.125) {
cout<<"NO CONVERGE EN
"
<<3*MAX<<"ITERACIONES";
getch();
exit(-1);
}
}
ABRIL /2000 84
}
while(N<=NRAICES);
return raices;
}
Nótese que un objeto de la clase Fpoly no sirve para
almacenar las raíces del polinomio, ya que algunas de
estas podrían ser complejos. De ahí surge la necesidad
de crear la clase Roots que representa un polinomio de
coeficientes complejos.

Download Report