SAWP > 1.3.3 Raízes de Equações — Raízes de Polinômios

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

March 9th, 2010 by SAWP | Filed under Métodos Computacionais

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

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/.

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.

Tags:

[ Métodos Computacionais ]

Você pode ler as eventuais respostas desta entrada através do RSS 2.0 feed.
As respostas estão encerradas, mas você pode trackback a partir do seu próprio site.

SAWP

Software Advice Working Party