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:
\( 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 \)
\( 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 \)
\( -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} \)
\( \dfrac{\partial}{\partial s} b_{0} = c_{2} \)
\( \dfrac{\partial}{\partial r} b_{1} = c_{2} \)
\( \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} \)
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)
portanto, generalizando para todos termos,
(eq22)
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)
e
(eq24)
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
- 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} \)
- Pela expressão abaixo, caso o polinômio seja de grau um:
(eq26)
\( x = – \dfrac{s}{r} \)
- 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/.