Método de Brent e implementación en Python

El método de Brent es un método numérico para encontrar las raíces de una función que combina los métodos de la interpolación cuadrática inversa, la secante y la bisección. Convirtiéndolo en un método más robusto y eficiente que los anteriores. Veamos más en detalle los fundamentos de este método y como se puede hacer una implementación en Python.

El método de Brent

Cuando se combinan varios métodos numéricos en uno se puede obtener un algoritmo más eficiente, como sucede en el método de Brent. Un algoritmo que combina los métodos de la interpolación cuadrática inversa, la secante y la bisección de la siguiente manera:

• Interpolación cuadrática inversa: es el método que se utiliza para obtener una mejor aproximación de la raíz cuando se puede aplicar al intervalo. Para ello, este método utiliza información sobre la función en tres puntos distintos para construir un polinomio cuadrático que se ajusta a la función. Siendo el punto donde el polinomio se hace cero la que se usa como aproximación de la raíz.
• Secante: cuando no es posible utilizar el método de la interpolación cuadrática inversa para aproximar la raíz en una iteración se usa el método de la secante. Calculando para ello la secante y usado el punto donde se cruza con el origen para mejorar la aproximación de la raíz.
• Bisección: finalmente, para garantizar la convergencia del algoritmo se usa el método de bisección. En la primera iteración, el método de Brent utiliza el de la bisección para reducir el tamaño del intervalo que contiene la raíz. Después de eso, el método de bisección se utiliza solo cuando el método de interpolación cuadrática inversa o el método de la secante no logran mejorar la aproximación de la raíz.

Implementación en Python del método de Brent

Una posible implementación en Python del método de Brent es la que se muestra en el siguiente código.

En Analytics Lane

Nuevo test de normalidad interactivo en el laboratorio de Analytics Lane

def brent(f, a=0.0, b=1.0, tol=1e-12, max_iter=100):
    """
    Encuentra la raíz de la función f(x) en el intervalo [a, b] utilizando el método de Brent.

    Parámetros
    ----------
    f : callable
        Función a la que se desea encontrar su raíz.
    a : float, opcional (valor por defecto = 0.0)
        Extremo inferior del intervalo.
    b : float, opcional (valor por defecto = 1.0)
        Extremo superior del intervalo.
    tol : float, opcional (valor por defecto = 1e-6)
        Tolerancia para el error en la aproximación de la raíz.
    max_iter : int, opcional (valor por defecto = 100)
        Número máximo de iteraciones permitidas.

    Devuelve
    --------
    x : float
        Aproximación de la raíz de f(x).
    """
    fa = f(a)
    fb = f(b)
    if fa*fb >= 0:
        raise ValueError("La función debe tener signos opuestos en los extremos del intervalo.")

    if abs(fa) < abs(fb):
        a, b = b, a
        fa, fb = fb, fa

    c = a
    fc = fa
    d = None
    mflag = True
    i = 0
    while i < max_iter and abs(b-a) > tol:
        if fa != fc and fb != fc:
            # Interpolación cuadrática inversa
            s = a * fb * fc / ((fa - fb) * (fa - fc)) + \
                b * fa * fc / ((fb - fa) * (fb - fc)) + \
                c * fa * fb / ((fc - fa) * (fc - fb))
        elif fa != fb:
            # Método de la secante
            s = b - fb * (b - a)/(fb - fa)

        if (s < (3 * a + b) / 4 or s > b) or \
                (mflag and abs(s - b) >= abs(b - c) / 2) or \
                (not mflag and abs(s - b) >= abs(c - d) / 2) or \
                (mflag and abs(b - c) < tol) or (not mflag and abs(c - d) < tol):
            # Método de la bisección
            s = (a + b) / 2
            mflag = True
        else:
            mflag = False

        fs = f(s)
        d = c
        c = b
        fc = fb
        if fa * fs < 0:
            b = s
            fb = fs
        else:
            a = s
            fa = fs

        if abs(fa) < abs(fb):
            a, b = b, a
            fa, fb = fb, fa

        i += 1

    if i == max_iter:
        raise RuntimeError("El método de Brent no converge")
    else:
        return b

Siendo los pasos que se siguen en el código los siguientes:

1. Verificación de que los signos de f(a) y f(b) sean diferentes. Si no es así, se lanzará una excepción.
2. Comprobación de que |f(a)| < |f(b)|, para intercambiar los valores si es necesario.
3. Inicialización de los valores usados en el cálculo.
4. Si los valores de f(a), f(b) y f(c) son distintos, se utiliza la interpolación cuadrática inversa para obtener una nueva aproximación s de la raíz de f(). En caso contrario, se utiliza el método de la secante.
5. Si la nueva aproximación s está fuera del intervalo de búsqueda o si la diferencia entre s y c es mayor que la diferencia entre las dos últimas aproximaciones, se utiliza la bisección para obtener una nueva aproximación de la raíz.
6. Se evalúa la aproximación s en la función y se actualizan los límites del intervalo.
7. El proceso se repite hasta que se alcance el límite de iteraciones o el intervalo sea inferior a la tolerancia indicada.

Evaluación de la implementación

Al igual que en entradas anteriores se puede usar diferentes funciones para comprobar que se obtienen buenas aproximaciones de las raíces con esta implementación.

f = lambda x: x - 4
g = lambda x: x**2 + 2*x - 8
h = lambda x: np.cos(2 * x)

brent(f, 0, 5) # 4.0
brent(g, 0, 5) # 1.9999999999997755
brent(h, 0, 2) # 0.7853981633974483

Conclusiones

El método de Brent combina en uno diferentes algoritmos para obtener la raíz de funciones: la interpolación cuadrática inversa, la secante y la bisección. Lo que le permite conseguir una rápida convergencia hacia la solución, aprovechad lo mejor de cada uno.

Imagen de Marcos Garzo en Pixabay