1.3.3 Raízes de Equações — Raízes de Polinômios — O Método de Bairstow

1. Introdução

O Método de Bairstow consiste em dividir um polinômio em fatores isolados que permitam aproximar suas raízes.

Possui a grande vantagem de conseguir extrair todas as raízes de um polinômio arbitrário sem a necessidade de percorrer todo o espaço de busca,como fariam o métodos de Muller ou Newton Raphson, que sempre convergem à aproximação mais próxima.

Além disso, ainda possui o benefício de fornecer as raízes complexas à uma taxa de convergência quadrática utilizando-se apenas de aritmética real.

2. Desenvolvimento do Método

Seja a função um polinômio geral:
(eq1)


\( f_{n} = a_{0} + a_{1}x + a_{2}x^2 + \cdots + a_{n}x^n \)

Se dividirmos o polinômio 1 pelo fator \(x-t \) , obtemos outro polinômio de um grau a menos:
(eq2)

\( f_{n-1} = b_{1} + b_{2}x + b_{3}x^2 + \cdots + b_{n}x^{n-1} \)

Caso esta divisão não gere resto, então \(t \) é uma raiz do polinômio \(f \). Mas, tomando o pior caso, onde há um resto \(R \) . Nesse caso \(R = b_0 \) e os coeficientes poderão ser calculados pela relação:
(eq3)

\( b_{n} = a_{n} \)

e
(eq4)

\( b_{i} = a_{i} + b_{i+1}t \)

onde \(i=n-1, n-2, \ldots, 0 \) .

O Método de Bairstow segue um procedimento semelhante ao descrito acima, contudo, ao invés de dividir \(f \) por um fator \(x-t \) , ele divide o polinômio pelo fator \(x*x – r*x – s \) . A introdução deste fator quadrático serve para permitir a determinação de raízes complexas, enquanto que o fator linear permitiria apenas aproximações reais.

Ao dividirmos o polinômio pelo fator de segunda ordem, obtemos:
(eq5)


\(f_{n-2} = b_{2} + b_{3} x + b_{4} x^2 + \cdots + b_{n-1} x^{n-3} + b_{n} x^{n-2}\)

isso gera uma expressão para o resto que é:
(eq6)

\( R = b_{1} \left(x -r \right) + b_{0} \)

Aplicando o algoritmo de Briot-Ruffine com \(x*x – r*x – s \) como divisor, obtemos a seguinte relação de recorrência:
(eq7)


\( b_{n} = a_{n} \)


\(\Downarrow \)

(eq8)

\( b_{n-1} = a_{n-1} + r b_{n} \)


\(\Downarrow \)

(eq9)


\( b_j = a_{j} + r b_{j+1} + s b_{j+2} \)

onde \(j = n-2, n-1, \ldots, 0 \) .

Com isso, o Método de Bairstow se resume a determinar os valores de \(r \) e \(s \) que tornam o fator quadrático \((x*x – r*x – s) \) exato. Ou seja, ele utiliza uma estratégia para encontrar os valores de \(r \) e \(s \) para qual \(R \) é zero. Quando \(R=0 \) , \(r \) será uma raiz do polinômio.

Da equação do resto (Equação 6) é possível notar que quando \(R \) é nulo, \(b_{0} \) e \(b_{1} \) são nulos. A técnica usada por Bairstow segue o desenvolvimento do Método de Newton, expandindo \(b_{0} \) e \(b_{1} \) em expressões que os aproximem de zero. Como \(b_{0} \) e \(b_{1} \) são funções de \(r \) e \(s \) , é possível expandir \(R \) por Série de Taylor para duas variáveis. Bairstow utiliza uma expansão de até a primeira ordem, gerando:

(eq10)


\( b_{1}(r+\Delta r, s+\Delta s) = b_{1} + \left( \dfrac{\partial}{\partial r} b_{1} \right) \Delta r + \left( \dfrac{\partial}{\partial s} b_{1} \right) \Delta s \)

(eq11)


\( b_{0}(r+\Delta r, s+\Delta s) = b_{0} + \left( \dfrac{\partial}{\partial r} b_{0} \right) \Delta r + \left( \dfrac{\partial}{\partial s} b_{0} \right) \Delta s \)

Como \(b_{0} \) e \(b_{1} \) serão nulos quando \(R=0 \) , então é possível reescrever as equações acimas como:
(eq12)


\( -b_{1} = \left( \dfrac{\partial}{\partial r} b_{1} \right) \Delta r + \left( \dfrac{\partial}{\partial s} b_{1} \right) \Delta s \)

(eq13)


\( -b_{0} = \left( \dfrac{\partial}{\partial r} b_{0} \right) \Delta r + \left( \dfrac{\partial}{\partial s} b_{0} \right) \Delta s \)

A partir das duas últimas equações, podemos obter os valores de \(b_{0} \) e \(b_{1} \) . O problema agora está nas derivadas parciais de \(b \) em relação a \(r \) e \(s \) . Se elas forem determinadas o problema é reduzido à um sistema de duas equações lineares com apenas duas incógnitas: \(\Delta r\) e \(\Delta s \).

Tomando as derivadas como constantes, renomeamos elas em termos de um vetor \(c \) , tal que:
(eq14)


\( \dfrac{\partial}{\partial r} b_{0} = c_{1} \)

(eq15)


\( \dfrac{\partial}{\partial s} b_{0} = c_{2} \)

(eq16)


\( \dfrac{\partial}{\partial r} b_{1} = c_{2} \)

(eq17)


\( \dfrac{\partial}{\partial s} b_{1} = c_{3} \)

A partir das suposições acimas, o sistema de equações pode ser reescrito com suas derivadas parciais substituídas pela variável \(c \) , tal que
(eq18)


\( c_{2} \Delta r + c_{3} \Delta s = -b_{1} \)

(eq19)

\( c_{1} \Delta r + c_{2} \Delta s = -b_{0} \)

O Método de Bairstow utiliza este sistema de equações lineares para ajustar variáveis complexas da mesma forma que as reais. As derivadas parciais, obtidas pela divisão sintética de \(b \), são:
(eq20)

\( c_{n} = b_{n} \)

(eq21)

\( c_{n-1} = b_{n-1} + r c_{n} \)

portanto, generalizando para todos termos,
(eq22)
\( c_{i} = b_{i} + r~c_{i+1} + s~c_{i+2} \)

onde \(i = n-1, n-2, \ldots, 1 \) .

Com este conjunto de equações é possível determinar \(\Delta r \) e \(\Delta s \), necessários para as aproximações de \(r \) e \(s \) , sendo o processo de cálculo do erro relativo de cada iteração é obtido pelas expressões:
(eq23)

\( \left| err(r) \right| = \left| \dfrac{\Delta r}{r} \right| \)

e
(eq24)
\( \left| err(s) \right| = \left| \dfrac{\Delta s}{s} \right| \)

Quando ambos erros estiverem abaixo do máximo tolerado, \(r \) e \(s \) estarão aproximados de forma que as raízes \(x \) serão determinadas por

  1. Pela expressão abaixo, caso o polinômio que deseja-se encontrar a raiz seja quadrático:
    (eq25)

    \( x = \dfrac{1}{2}r + \dfrac{1}{2} \sqrt{r^2 + 4s} \)

  2. Pela expressão abaixo, caso o polinômio seja de grau um:
    (eq26)

    \( x = – \dfrac{s}{r} \)

  3. Caso o polinômio seja de grau maior ou igual a três, o Método de Bairstow será aplicado ao quociente para o cálculo de novos valores \(r \) e \(s \) usando a função abaixo. Os valores de \(r \) e \(s \) calculados anteriormente servem como parâmetro inicial para a próxima iteração.

3. Implementação

Implementação em Python:

def bairstow(a, rr=1, ss=1, errto=0.001, imax=500):
    """
    Calculate the roots of a function (using the Bairstow Method)
 
    bairstow(a, polDegree, errto, rr, ss, imax)
 
    * a[]: list with the coeficients
    * rr: initial aprox. of "r" param
    * ss: initial aprox. of "s" param
    * errto: tolerated error
    * imax: max iterCountation allowed (not inf-loop)
 
    return: a list with all roots of polynom
 
    Author: Pedro Garcia [sawp@sawp.com.br]
    see: http://www.sawp.com.br
 
    License: Creative Commons
         http://creativecommons.org/licenses/by-nc-nd/2.5/br/
 
    Jan 2010
    """
 
    def filt(l):
        r = l.real
        i = l.imag
        if abs(r) < errto:
            r = 0.0
        if abs(i) < errto:
            i = 0.0
        return complex(r, i)
 
    r = rr
    s = ss
    n = len(a) - 1
 
    ea1 = errto + 1
    ea2 = errto + 1
    iterCount = 0
    roots = [ complex(0,0) for i in xrange(len(a)) ]
    b = [ complex(0,0) for i in xrange(len(a)) ]
    c = [ complex(0,0) for i in xrange(len(a)) ]
 
    while n >= 3 or iterCount < imax:
        iterCount = 0
 
        while ea1 >= errto and ea2 >= errto or iterCount < imax:
            iterCount += 1
            b[n] = a[n]
            b[n-1] = a[n-1] + r * b[n]
            c[n] = b[n]
            c[n-1] = b[n-1] + r * c[n]
 
            for i in xrange(n-2, -1, -1):
                 b[i] = a[i] + r * b[i+1] + s * b[i+2]
                 c[i] = b[i] + r * c[i+1] + s * c[i+2]
 
            det = c[2]*c[2] - c[3]*c[1]
 
            if det != 0.0:
                dr = (-b[1]*c[2] + b[0]*c[3])/det
                ds = (-b[0]*c[2] + b[1]*c[1])/det
                r = r + dr
                s = s + ds
 
                if r != 0.0:
                    ea1 = abs(dr/r)
 
                if s != 0.0:
                    ea2 = abs(ds/s)
            else:
                r += errto
                s += errto
                iterCount = 0
 
        (r1, i1, r2, i2) = quadratic_solve(r, s)
        roots[n] = complex(r1, i1)
        roots[n-1] = complex(r2, i2)
        n = n - 2
 
        for i in xrange(n):
            a[i] = b[i+2]
 
    if n == 2:
        r = -a[1]/a[2]
        s = -a[0]/a[2]
        (r1, i1, r2, i2) = quadratic_solve(r, s)
        roots[n] = complex(r1, i1)
        roots[n-1] = complex(r2, i2)
    else:
        roots[n] = complex(-a[0]/a[1], 0.0)
    print iterCount
    roots[1:len(roots)].sort(key=lambda r: r.real)
    return map(filt, roots[1:len(roots)])

Note que a função quadratic\_solve foi incluída para deixar explícita a resolução da fórmula quadrática. A implementação desta função é:

def quadratic_solve(r, s):
    """
    Calculate the roots of a parabolic equation: a*x^2 + b*x + c = 0 (*)
 
    quadratic_solve(r, s)
 
    * r: b coeficient of (*)
    * s: a*c, coeficients of (*)
 
    return: a tuple with roots of (*)
 
    Author: Pedro Garcia [sawp@sawp.com.br]
    see: http:!www.sawp.com.br
 
    License: Creative Commons
         http://creativecommons.org/licenses/by-nc-nd/2.5/br/
 
    Jan 2009
    """
    delta = r*r + 4.0*s
 
    if delta > 0:
        r1 = (r + sqrt(delta))/2.0
        r2 = (r - sqrt(delta))/2.0
        i1 = 0.0
        i2 = 0.0
    else:
        r1 = r/2.0
        r2 = r1
        i1 = sqrt(abs(delta))/2.0
        i2 = -i1
 
    return (r1, i1, r2, i2)

Observe que o programa altera os valores de \(r \) e \(s \) com algum incremento aleatório e repete o procedimento da resolução das equações quando a determinante é nula, forçando o programa a convergir à uma aproximação melhor.

Um exemplo da utilização desta sub-rotina pode ser encontrado em: http://www.sawp.com.br/code/rootfind/bairstow.py

A utilização desta função não é tão simples, pois exige como parâmetro de entrada um vetor com os coeficientes do polinômio. Embora a utilização de listas para representar um polinômio seja algo comum para muitos programadores, certamente é uma abordagem menos intuitiva que passar uma função como referência.

Também há disponível uma outra implementação deste algoritmo em Fortran 90 no endereço http://www.sawp.com.br/code/rootfind/bairstow.f90

4. Copyright

Este documento é disponível sob a licença Creative Commons. As regras dos direitos de cópia deste conteúdo estão acessíveis em http://creativecommons.org/licenses/by-nc-nd/2.5/br/.

References


[1] http://mathworld.wolfram.com/BairstowsMethod.html
[2] http://en.wikipedia.org/wiki/Bairstow's_method
[3] Ralston, Anthony and Rabinowitz, Philip, A First Course in Numerical Analysis.,(2nd ed.) McGraw-Hill and Dover, (2001), 374.
[4] N.B Franco, Cálculo Numérico, Pearson Prentice Hall, (2006), 99.