PDF de programación - Implementación del esquema de deflación para ecuaciones algebraicas

Implementación del esquema de deflación para ecuaciones

algebraicas

Pablo Santamaría

v0.1 (Julio 2007)

Introducción

Consideremos el problema de determinar todos los ceros de una ecuación algebraica

p(x) = a0 + a1x + a2x2 + · · · + amxm,

a0 6= 0.

De acuerdo al teorema fundamental de la aritmética, esta ecuación tiene, en el campo complejo,
exactamente m raíces, α1, α2, . . . , αn y p(x) se factoriza como

p(x) = a0(x − α1)(x − α2) · · · (x − αm)

Si los coeficientes a0, . . . , am son todos reales, entonces las posibles raices complejas aparecerán
en pares conjugados, y el polinomio p(x) se factorizará, en el campo real, como el producto de
factores lineales y polinomios cuadráticos irreducibles.

En un apunte anterior1 vimos que para ecuaciones algebraicas el método de Newton puede
ser eficientemente implementado si la evaluación del polinomio (y su derivada) es realizada por el
método iterativo de Horner. El procedimiento resultante es conocido como método de Birge–Vieta
y su implementación como una subrutina de fortran se encuentra en dicho apunte.

Ahora bien, el método de Horner nos da también un camino para factorizar el polinomio y
hallar, así, sus raíces. En efecto, recordemos que la evaluación de p(x) en un punto xn por la regla
de Horner procede computando

(bm = am

bk = ak + bk+1xn

k = m − 1, . . . , 0

siendo, entonces

Pero además, si construimos el polinomio de grado m − 1

b0 = p(xn).

q(x) = b1 + b2x2 + · · · + bm−1xm−2 + bmxm−1

entonces

p(x) = (x − xn)q(x) + b0

Ahora bien, supongamos que hemos encontrado por el método de Bierge–Vieta una raíz xn = α
de p(x), entonces b0 = 0 y las restantes raíces son también raíces del polinomio q(x). Para calcular
una de estas raíces aplicamos el método de Bierge–Vieta ahora a q(x), obteniendo la raíz y un
factor polinómico de un grado menor cuyas raíces son también raíces de p(x). El proceso se repite
tantas veces como raíces reales tenga el polinomio, obteniendo en cada paso polinomios de grado

1Ver Subrutinas para la resolución de ecuaciones de una variable.

1

cada vez menor. Si al final del proceso obtenemos un polinomio de grado mayor que 1, éste será
irreductible en el campo real y sus raíces complejas deben ser determinadas con algún método
apropiado2. El procedimiento descrito se conoce como deflación.

Nótese que en la discusión anterior no hemos tenido en cuenta el hecho de que xn no es exacta-
mente una raíz de p(x), y hemos ignorado cualquier error de redondeo introducido en el cómputo
de los coeficientes b0, b1, . . . , bm. Estos errores pueden producir que las raíces que calculamos de
los sucesivos polinomios factorizados se alejen cada vez más de las raíces del polinomio original
p(x). Puede mostrarse que los errores cometidos son despreciables si: (i) las raíces son determi-
nadas ordenándolas por su incremento en valor absoluto (es decir, primero las que tienen valores
absolutos más chicos) y (ii) cada raíz es determinada a su máxima precisión alcanzable. En todo
caso, una forma simple de solventar el problema consiste en utilizar cada raíz determinada por el
procedimiento de deflación como aproximación inicial para la ecuación original p(x) = 0 y efectuar
al menos una última iteración.

Implementación del esquema de deflación

Para implementar el esquema de deflación tenemos que modificar primero nuestra subrutina
original del método de Bierge–Vieta de manera tal que ahora devuelva también los coeficientes bk
para construir el polinomio factorizado q(x) requerido en cada paso. A continuación presentamos
la subrutina incorporando tal modificación.

C

C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C

---------------------------------------------------
SUBROUTINE BIRGE_VIETA (A,B,M,X0,TOL,NMAX,RAIZ,ICLAVE)
---------------------------------------------------
METODO DE BIRGE-VIETA para resolver ECUACIONES
ALGEBRAICAS:

P (x) = 0

donde P es un polinomio de grado n de coeficientes
reales; basado en el método de Newton-Raphson
implementando el esquema de Horner para la evalua-
ción del polinomio y su derivada.
---------------------------------------------------
Bloque de identificación de argumentos
---------------------------------------------------
A

= Arreglo unidimensional de 0 a M que

contiene los coeficientes del polinomio

B

= Arreglo unidimensional de 0 a M que

contiene los coeficientes del polinomio
factor resultante de la division sintetica
(Nota: B(0) = P(RAIZ)=~ 0)

X0
M
NMAX
TOL
RAIZ
ICLAVE = Clave de éxito

= Aproximación inicial
= Grado del polinomio
= Número máximo de iteraciones
= Tolerancia para el error relativo
= Estimación de la solución

0 = éxito
1 = iteraciones excedidas

2Si el polinomio es de grado 2 la determinación de las raíces, claramente, se reduce a utilizar la fórmula usual

para ecuaciones cuadráticas.

2

C
C
C
C

C

C

C
C
C

C
C
C

2 = división por cero

---------------------------------------------------
Bloque de declaración de tipo
---------------------------------------------------
IMPLICIT NONE
INTEGER M, NMAX
DOUBLE PRECISION A(0:M), B(0:M)
DOUBLE PRECISION X0,TOL,RAIZ
INTEGER ICLAVE
---------------------------------------------------
INTEGER I, J
DOUBLE PRECISION C, ERROR
DOUBLE PRECISION TINY
PARAMETER (TINY=1.0D-25)
---------------------------------------------------
DO I=1,NMAX

-------------------------
Esquema de Horner
-------------------------
B(M) = A(M)
C
= A(M)
DO J=M-1,1,-1

= B(J)+X0*C

B(J) = B(J+1)*X0+A(J)
C
ENDDO
B(0) = B(1)*X0+A(0)
IF (ABS(C).LE.TINY) THEN

ICLAVE = 2
RETURN

ENDIF
-------------------------
Método de Newton
-------------------------
RAIZ = X0 - B(0)/C
ERROR = ABS((RAIZ-X0)/RAIZ)
IF (ERROR.LT.TOL) RETURN
X0 = RAIZ

END DO
ICLAVE = 1
RETURN
END

Con la subrutina anterior podemos construir fácilmente un programa para implementar el
método de deflación. Inicialmente el usuario ingresa el grado del polinomio, sus coeficientes y
aproximaciones iniciales, suficientemente buenas, para cada una de las raíces reales. Entonces el
programa obtiene mejores aproximaciones de las raíces (dentro de una tolerancia prefijada en el
programa) y, paso a paso, factoriza el polinomio.

C

C

PROGRAM DEFLACION
--------------------------------------------
IMPLICIT NONE
--------------------------------------------
INTEGER MMAX
PARAMETER (MMAX=50)

! Grado máximo del polinomio

3

C

C
C
C

C
C
C

C
C
C

C
C
C

C
C
C

C
C
C

! Precisión para las raíces

! Número máximo de iteraciones

INTEGER NMAX
PARAMETER (NMAX=100)
DOUBLE PRECISION TOL
PARAMETER (TOL=1.0D-20)
--------------------------------------------
INTEGER I, J, NRAIZ, M, ICLAVE
DOUBLE PRECISION RAIZ
DOUBLE PRECISION A(0:MMAX), B(0:MMAX), X(MMAX)
--------------------------------------------
Leer coeficientes del polinomio
--------------------------------------------
WRITE(*,’(A,$)’) ’ Grado del polinomio: ’
READ(*,*) M
WRITE(*,’(A,$)’) ’ Ingrese coeficientes: ’
READ(*,*) (A(I), I=0,M)
--------------------------------------------
Leer las aproximaciones de las raices
--------------------------------------------
WRITE(*,’(A,$)’) ’ Número de raíces reales a considerar: ’
READ(*,*) NRAIZ
IF (NRAIZ.GT.M) THEN

WRITE(*,*) ’ El número de raíces no puede exceder al grado’
STOP

ENDIF
WRITE(*,’(A,$)’) ’ Ingrese aproximaciones iniciales: ’
READ(*,*) (X(I), I=1,NRAIZ)
--------------------------------------------
Imprimir polinomio original
--------------------------------------------
WRITE(*,*) ’---------------------------’
WRITE(*,*) ’Grado del polinomio = ’, M
WRITE(*,*) ’Coeficientes:’
WRITE(*,*) (A(J),J=0,M)
--------------------------------------------
Proceder con la factorización
--------------------------------------------
DO I=1,NRAIZ

----------------------------------------
Efectuar el método de Birge Vieta
----------------------------------------
RAIZ
ICLAVE = 0
CALL BIRGE_VIETA(A,B,M,X(I),TOL,NMAX,RAIZ,ICLAVE)
IF (ICLAVE.NE.0) THEN

= 0.0D0

WRITE(*,*) ’Error. No se puede continuar’
STOP

ENDIF
----------------------------------------
Imprimir resultados
----------------------------------------
WRITE(*,*) ’--------------------’
WRITE(*,*) ’Factorizacion ’, I
WRITE(*,*) ’--------------------’
WRITE(*,*) ’Raiz = ’, RAIZ

4

C
C
C

WRITE(*,*) ’Grado del factor = ’, M-1
WRITE(*,*) ’Coeficientes:’
WRITE(*,*) (B(J),J=1,M)
----------------------------------------
Tomar nuevo polinomio
----------------------------------------
M = M-1
DO J = 0,M

A(J) = B(J+1)

ENDDO

ENDDO
STOP
END

Un ejemplo

Como ejemplo, consideremos la ecuación algebraica

p(x) = −2 − x − x2 − x3 + x4 = 0

El gráfico del polinomio nos muestra que hay dos raíces reales, y que ˜x1 = −1,5 y ˜x2 = 2,5 son,
respectivamente, buenas aproximaciones para ellas3. Al ejecutar nuestro programa obtenemos lo
siguiente.

4

Grado del polinomio: 4
Ingrese coeficientes: -2.0 -1.0 -1.0 -1.0 1.0
Número de raíces reales a considerar: 2
Ingrese aproximaciones iniciales: -1.5 2.5
---------------------------
Grado del polinomio =
Coeficientes:
-2. -1. -1. -1. 1.
--------------------
Factorizacion 1
--------------------
Raiz = -1.
Grado del factor = 3
Coeficientes:
-2. 1. -2. 1.
--------------------
Factorizacion 2
--------------------
Raiz =
Grado del factor = 2
Coeficientes:

2.

1.

0. 1.

De este modo, la factorización de nuestro polinomio es

p(x) = (x + 1)(x − 2)(1 + x2)

siendo x1 = −1 y x2 = 2 las dos raíces reales. Puesto que 1 + x2 = (x + i)(x − i), resulta que las
dos restantes raíces (complejas) son x3 = i, x4 = −i.

3Los ceros exactos son (no se lo digan a nadie) x1 = −1, x2 = 2, x3 = i, x4 = −i.

5