1. Introdução

O Método de Brent é uma combinação simples de diversos métodos computacionais para busca de raízes de função. Os métodos usados no algoritmo de Brent são: Método da Bisseção, da Secante e Interpolação Quadrática Inversa.

A ideia por trás deste algoritmo está em aproveitar as vantagens de cada uma das técnicas combinadas, procurando minimizar suas desvantagens. Desta forma, utiliza a Interpolação Quadrática Inversa sempre que as aproximações intermediárias não se sobrepõem (o que acarretaria em uma divisão por zero). Quando há a sobreposição, tenta fazer a próxima aproximação utilizando o Método da Secante.

Todavia, estas duas últimas técnicas são suscetíveis à divergir ao longo da execução. Caso este comportamento seja detectado, o algoritmo “ativa” o Método da Bisseção a fim de reconverter a aproximação ao zero da função.

Sumarizando, suas maiores vantagens são: tende a convergir rapidamente, não necessita de derivadas e sempre converge, caso a raiz esteja delimitada pelo intervalo necessário ao Método da Bisseção.

Por outro lado, possui a desvantagem de exigir que um intervalo contendo a raiz seja fornecido, o que reduz a generalidade da aplicação do método em alguns casos.

Este método, embora não apresente nenhuma novidade em termos matematicamente formais, é um algoritmo que foi vastamente usado e é presente em diversas bibliotecas e softwares como função-chave na aproximação de raízes.

2. Implementação

def brent(f, xl, xr, errto=0.001, imax=500):
    """
    Return the root of function (Using the Brent's Method)

    brent(f, xl, xr, errto=0.001, imax=500)

    * f: Function where find roots
    * xl: left bound of interval
    * xr: right bound of interval
    * errto: tolerated error
    * imax: max of iterations allowed

    return: root next xl

    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
    """

    errno = errto + 1
    iterCount = 0
    x = xl
    y = xr
    xrold = 0
    mflag = True
    delta = xr

    if f(xl) * f(xr) >= 0.0:
        return Float("nan")

    if abs(f(xl)) < abs(f(xr)):
        (xr, xl) = (xl, xr)

    while errno > errto and iterCount < imax:
        fl = f(xl)
        fr = f(xr)
        fx = f(x)

        # try Inverse Quadratic Interpolation
        if fl != fx and fr != fx:
            delta = 7 * [0] 
            delta[1] = fl - fx
            delta[2] = fl - fr
            delta[3] = fx - fl
            delta[4] = fx - fr
            delta[5] = fr - fl
            delta[6] = fr - fx

            def disturb(ecs):            
                if(ecs == 0):
                    return errto
                else:
                    return ecs

            delta = map(disturb, delta)

            s = fx * fr * xl / (delta[1] * delta[2]) + \
                fl * fr * x / (delta[3] * delta[4]) + \
                fl * fx * xr / (delta[5] * delta[6])
        # try Secant's Rule
        else:
            s = xr - fr*(xr - xl)/(fr - fl)

        if (not (s > (3 * xl + xr)/4 and s < xr)) or \
            (mflag and abs(s - xr) >= abs(xr - x)/2) or \
            (not mflag and abs(s - xr) >= abs(x - y)/2) or \
            (mflag and abs(xr - x) < abs(delta)) or \
            (not mflag and abs(x - y) < abs(delta)):
                # bissection method
                s = (xr + xl)/2
                mflag = True
        else:
            mflag = False

        fs = f(s)
        y = x
        x = xr

        if fl * fs < 0.0:
            xr = s
        else:
            xl = s

        if abs(f(xl)) < abs(f(xr)):
            (xr, xl) = (xl, xr)

        if fs == 0 or f(xr) == 0:
            break

        if abs(xr - xrold) < delta:
            delta = abs(xr - xrold)            

        xrold = xr
        iterCount += 1
        errno = fabs(xr - xl)

    return xr

Disponível para download em http://www.sawp.com.br/code/rootfind/brent.py </p> </p>

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/BrentsMethod.html
[2] http://math.fullerton.edu/mathews/n2003/BrentMethodMod.html