1.3.4 Raízes de Equações – Raízes de Polinômios – O Método de Aberth

1. Introdução

Apresentado por Oliver Aberth em seu artigo “Iteration Methods for finding all zeros of a polynomial simultaneously”, este método permite encontrar todas as raízes de um polinômio, desde que seja conhecido seu grau.

Possui a vantagem do método de Bairstow de convergir à todas as raízes, sendo elas reais ou complexas, diferentemente do método de Laguerre.

Sendo derivado do método de Newton, possui convergência quadrática e pode ocorrer das raízes não convergirem para polinômios aproximados (interpoladores de funções não-lineares).

2. Desenvolvimento do Método

Seja o polinômio de grau \(n \) formado por uma composição de coeficientes reais:
(eq1)


\( p \left( x \right) = p_{{n}}{x}^{n}+p_{{n-1}}{x}^{n-1} + \cdots + p_{{1}}x+p_{{0}} \)

Supondo que este polinômio possui todas raízes complexas, ele pode ser descrito na forma fatorizada:
(eq2)


\( p \left( x \right) =p_{{n}}\prod _{i=1}^{n}(x-z_{{i}}) \)

onde \(z_{i} \) são as raízes complexas de p(x).

Podemos adicionar um termo neste produtório tal que \({C}{(x – z_{i})} = 1 \) , que implica um termo a mais no produtório da Equação 2. Contudo, ao invés de notarmos este produtório pelos índices \(i=1 \ldots n \) , vamos tomar \(i=0 \ldots n \).
(eq3)


\( p \left( x \right) =C\prod _{i=0}^{n}(z-x_{{i}}) \)

Comparando as Fórmulas 2 e 3, observamos que a segunda é um polinômio que recebe um termo a mais em \(i=0 \) . Para que \(p(x) \) continue sendo o mesmo, o novo termo deve ser nulo ou então a seguinte relação válida:
(eq4)


\( C=p_{{n}} \left( z_{{k}} \right) \)

Isolando \(p_{n} \) na Equação 2 e levando na Equação 4, temos que:
(eq5)


\( C \left( x_{{i}} \right) ={\dfrac {p \left( x_{{i}} \right) }{\prod _{k=0}^{n}(x_{{i}}-z_{{k}})}} \)

Pelo método de Newton-Raphson, sabemos que uma raiz pode ser encontrada através da seguinte função de iteração:
(eq6)


\( x_{{i+1}}=x_{{i}}-{\dfrac {f \left( x_{{i}} \right) }{{\dfrac {d}{dx}}f \left( x_{{i}} \right) }} \)

Usando a relação 5 como sendo a fórmula a ser aplicada Equação 6:
(eq7)

\( f \left( x \right) =C \left( x \right) \)

Resultando em:
(eq8)

\( z_{{i+1}}=z_{{i}}-{\dfrac {C \left( z_{{i}} \right) }{{\dfrac {d}{dx}}C \left( z_{{i}} \right) }} \)

Uma propriedade interessante da função de iteração usada no Método de Newton é que ela equivale ao inverso da derivada do logaritmo da função. Ou seja, uma fórmula facilmente obtida usando-se a Regra da Cadeia[1] é:
(eq9)


\( {\dfrac {{\dfrac {d}{dx}}C \left( x \right) }{C \left( x \right) }}={\dfrac {d}{dx}}\ln \left( C \left( x \right) \right)\)

Tomando o logaritmo da Equação 5,
(eq10)


\( \ln \left( C \left( x \right) \right) =\ln \left( {\dfrac {p \left( x_{{i}} \right) }{\prod _{k=0}^{n}(x_{{i}}-z_{{k}})}} \right)\)

onde
(eq11)

\( \ln \left( C \left( x \right) \right) =\ln \left( p \left( x_{{i}} \right) \right) -\ln \left( \prod _{k=0}^{n}(x_{{i}}-z_{{k}})\right)\)

como
(eq12)

\( \ln \left( \prod _{k=0}^{n}x_{{i}}-z_{{k}} \right) =\sum _{k=0}^{n} \ln \left( \left| x-z_{{k}} \right| \right) \)

Substituindo a equação 12 em 11,
(eq13)

\( \ln \left( C \left( x \right) \right) =\ln \left( p \left( x_{{i}} \right) \right) -\sum _{k=0}^{n}\ln \left( \left| x-z_{{k}}\right| \right)\)

Levando a Equação 13 na Equação 9, temos que:
(eq14)

\( {\dfrac {{\dfrac {d}{dx}}C \left( x \right) }{C \left( x \right) }}={\dfrac {\partial }{\partial x}} \left( \ln \left( p \left( x \right) \right) -\sum _{j=0}^{n}\ln \left( \left| x-z_{{j}} \right| \right) \right) \)

Aplicando novamente a relação mencionada para a expressão 9, conseguimos eliminar o logaritmo do polinômio, onde deixamos a relação apenas com a própria função polinomial e sua derivada:
(eq15)


\( {\dfrac {{\dfrac {d}{dx}}C \left( x \right) }{C \left( x \right) }}={\dfrac {{\dfrac {d}{dx}}p \left( x \right) }{p \left( x \right) }}-\sum _{j=0}^{n} \left( x-z_{{j}} \right) ^{-1}\)

como o inverso da Equação 15 é exatamente a função de iteração do Método de Newton, temos que:
(eq16)

\( {\dfrac {C \left( z_{{i}} \right) }{{\dfrac {d}{dx}}C \left( z_{{i}} \right) }}= \left( {\dfrac {{\dfrac {d}{dx}}p \left( x \right) }{p \left( x \right) }}-\sum _{j=0}^{n} \left( x-z_{{j}} \right) ^{-1} \right) ^{-1} \)

ou seja, do Método de Newton
(eq17)

\( z_{{i+1}}=z_{{i}}-{\dfrac {C \left( z_{{i}} \right) }{{\dfrac {d}{dx}}C \left( z_{{i}} \right) }} \)

obtemos a expressão do Método de Aberth:
(eq18)

\( z_{{i+1}}=z_{{i}}- \left( {\dfrac {{\dfrac {d}{dx}}p \left( x \right) }{ p \left( x \right) }}-\sum _{{j=0}{j \neq i}}^{n} \left( z_{{i}}-z_{{j}} \right) ^{-1} \right) ^{-1} \)

3. Implementação

def aberth(p, d1p, polDegree=1, x0=complex(1,1), errto=0.001, imax=500):
    """
    Return all roots of Polynom (using the Aberth Method)
 
    aberth(p, d1p, polDegree=1, x0=1.0, errto=0.001, imax=100)
 
    * p: polynom where find the roots
    * d1p: analytic derivative of f
    * polDegree: polynom degree
    * x0: initial estimative (complex)
    * errto: tolerated error
    * imax: max of iterations allowed
 
    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/
 
    Dec 2009
    """
 
    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)
    n = polDegree
    errno = errto + 1.0
    iterCount = 0
    if type(x0) == list:
        zold = x0[:] + abs(n-len(x0))*[ complex(0,0) ]
    elif type(x0) == complex:
        zold = [ x0+random() for i in xrange(n) ]
    else:
        zold = [ complex(random()*x0,random()*x0) for i in xrange(n) ]
    znew = n * [ complex(0, 0) ]
    while errno > errto and iterCount < imax:
        for i in xrange(len(zold)):
            pz = p(zold[i])
            while pz == 0:
                zold[i] += random()
                pz = p(zold[i])  
            dp1z = d1p(zold[i])              
            d = dp1z/pz
            s = 0
            for w in (range(i) + range(i+1, len(zold))):
                den = zold[i] - zold[w]
                while den == 0:
                    den += random()
                s += 1/(den)
            if (d - s) != 0:
                znew[i] = zold[i] - 1.0/(d - s)
            else:
                znew[i] = zold[i] - 1.0/(d - s + errto*random())
        errno = 0
        for i in xrange(len(znew)):
            if znew[i] != 0:
                t = abs(znew[i] - zold[i])/abs(znew[i])
                if errno < t:
                    errno = t
        zold = znew[:]
        iterCount += 1
    znew.sort(key=lambda r: r.real)
    return map(filt, znew)

Como todas raízes são desconhecidas, lista cujas raízes serão armazenadas é iniciada com valores aleatórios distribuídos sobre o plano complexo. Esta construção toma como base o valor passado no parâmetro \(x_{0} \) . Caso a estimativa inicial \(x_{0} \) seja uma lista, o algoritmo as utiliza e completa os outros elementos com um valor nulo. Contudo, é recomendável que todas as \(n \) raízes recebam uma estimativa inicial. Além disso, a função filt filtra os resultados cujo valor deu abaixo do erro estimado, servindo apenas para melhorar a visualização do resultado. Além disso, sob algumas restrições, esta implementação também permite que funções não-polinomiais sejam aproximadas com sucesso, como pode ser visto neste exemplo de utilização: http://www.sawp.com.br/code/rootfind/aberth.py

4. Conclusões

Assim como o método de Bairstow, o método de Aberth converge a uma taxa quadrática à todas raízes de um polinômio qualquer, sendo uma das soluções mais gerais e elegantes existentes.

Possui a vantagem sobre o método de Bairstow de não exigir que sejam passados os coeficientes do polinômio como parâmetro para o algoritmo, sendo mais conveniente quando trabalhamos com uma função conhecida, com alto grau e grande quantidade de coeficientes.

Todavia, possui a desvantagem de necessitar de uma função que compute a derivada analítica do polinômio, o que geralmente é fácil de se obter, mas que pode consumir considerável tempo de processamento.

Além disso, como a função de iteração de Aberth é obtida diretamente do Método de Newton, ela carrega a desvantagem da estimativa inicial para convergir ao resultado correto.

5. 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] B.G. Thomas, Cálculo, volume 1, Addison Wesley, V.1,(2002), 179.
[2] Ralston, Anthony and Rabinowitz, Philip, A First Course in Numerical Analysis.,(2nd ed.) McGraw-Hill and Dover, (2001), 358–359.
[3] D.A Bini, Numerical computation of polynomial zeros by means of Aberth’s method, Numerical Algorithms, V.13,(1996), 179–200.
[5] http://en.wikipedia.org/wiki/Horner_scheme
[4] http://en.wikipedia.org/wiki/Aberth_method