El método de Brent

En análisis numérico , el método de Brent es un algoritmo híbrido para la búsqueda de raíces que combina el método de bisección , el método de la secante y la interpolación cuadrática inversa . Posee la fiabilidad de la bisección, pero puede ser tan rápido como algunos de los métodos menos fiables. El algoritmo intenta utilizar el método de la secante, que converge rápidamente, o la interpolación cuadrática inversa si es posible, pero recurre al método de bisección, más robusto, si es necesario. El método de Brent se debe a Richard Brent ^{[ 1 ]} y se basa en un algoritmo anterior de Theodorus Dekker ^{[ 2 ]} . Por consiguiente, también se le conoce como el método de Brent-Dekker .

Las mejoras modernas del método de Brent incluyen el método de Chandrupatla, que es más simple y rápido para funciones que son planas alrededor de sus raíces; ^{[ 3 ] [ 4 ]} el método de Ridders , que realiza interpolaciones exponenciales en lugar de cuadráticas, proporcionando una fórmula cerrada más simple para las iteraciones; y el método ITP , que es un híbrido entre regula-falsi y bisección que logra garantías óptimas en el peor de los casos y asintóticas.

La idea de combinar el método de bisección con el método de la secante se remonta a Dekker (1969) .

Supongamos que queremos resolver la ecuación f ( x ) = 0. Al igual que con el método de bisección, es necesario inicializar el método de Dekker con dos puntos, digamos a ₀ y b ₀ , tales que f ( a ₀ ) y f ( b ₀ ) tengan signos opuestos. Si f es continua en [ a ₀ , b ₀ ], el teorema del valor intermedio garantiza la existencia de una solución entre a ₀ y b ₀ .

En cada iteración intervienen tres puntos:

• b _k es la iteración actual, es decir, la estimación actual para la raíz de f .
• a _k es el "contrapunto", es decir, un punto tal que f ( a _k ) y f ( b _k ) tienen signos opuestos, por lo que el intervalo [ a _k , b _k ] contiene la solución. Además, | f ( b _k )| debe ser menor o igual que | f ( a _k )|, de modo que b _k es una mejor aproximación para la solución desconocida que a _k .
• b _{k − 1} es la iteración anterior (para la primera iteración, se establece b _{k − 1} = a ₀ ).

Se calculan dos valores provisionales para la siguiente iteración. El primero se obtiene mediante interpolación lineal, también conocida como método de la secante:

${\displaystyle s={\begin{cases}b_{k}-{\frac {b_{k}-b_{k-1}}{f(b_{k})-f(b_{k-1})}}f(b_{k}),&{\mbox{si }}f(b_{k})\neq f(b_{k-1})\\m&{\mbox{en otro caso }}\end{cases}}}$

y el segundo se obtiene mediante el método de bisección.

${\displaystyle m={\frac {a_{k}+b_{k}}{2}}.}$

Si el resultado del método de la secante, s , se encuentra estrictamente entre b _k y m , entonces se convierte en la siguiente iteración ( b _{k +1} = s ), de lo contrario se utiliza el punto medio ( b _{k +1} = m ).

Luego, se elige el valor del nuevo contrapunto de manera que f ( a _{k +1} ) y f ( b _{k +1} ) tengan signos opuestos. Si f ( a _k ) y f ( b _{k +1} ) tienen signos opuestos, entonces el contrapunto permanece igual: a _{k +1} = a _k . De lo contrario, f ( b _{k +1} ) y f ( b _k ) tienen signos opuestos, por lo que el nuevo contrapunto se convierte en a _{k +1} = b _k .

Finalmente, si | f ( a _{k +1} )| < | f ( b _{k +1} )|, entonces a _{k +1} es probablemente una mejor estimación para la solución que b _{k +1} , y por lo tanto los valores de a _{k +1} y b _{k +1} se intercambian.

Con esto finaliza la descripción de una sola iteración del método de Dekker.

El método de Dekker funciona bien si la función f se comporta razonablemente bien. Sin embargo, existen circunstancias en las que cada iteración emplea el método de la secante, pero las iteraciones b _k convergen muy lentamente (en particular, | b _k − b _{k − 1} | pueden ser arbitrariamente pequeñas). En este caso, el método de Dekker requiere muchas más iteraciones que el método de bisección.

Brent (1973) propuso una pequeña modificación para evitar el problema del método de Dekker. Inserta una prueba adicional que debe cumplirse antes de que el resultado del método de la secante se acepte como la siguiente iteración. Deben cumplirse simultáneamente dos desigualdades:

Dada una tolerancia numérica específica${\displaystyle \delta }$, si el paso anterior utilizó el método de bisección, la desigualdad${\estilo de texto |\delta |<|b_{k}-b_{k-1}|}$debe cumplirse para realizar la interpolación; de lo contrario, se realiza el método de bisección y su resultado se utiliza para la siguiente iteración.

Si el paso anterior realizó una interpolación, entonces la desigualdad${\textstyle |\delta |<|b_{k-1}-b_{k-2}|}$se utiliza en su lugar para realizar la siguiente acción (elegir) interpolación (cuando la desigualdad es verdadera) o método de bisección (cuando la desigualdad no es verdadera).

Además, si el paso anterior utilizó el método de bisección, la desigualdad${\textstyle |s-b_{k}|<{\begin{matrix}{\frac {1}{2}}\end{matrix}}|b_{k}-b_{k-1}|}$ debe cumplirse, de lo contrario se realiza el método de bisección y su resultado se utiliza para la siguiente iteración. Si el paso anterior realizó una interpolación, entonces la desigualdad${\textstyle |s-b_{k}|<{\begin{matrix}{\frac {1}{2}}\end{matrix}}|b_{k-1}-b_{k-2}|}$ se utiliza en su lugar.

Esta modificación garantiza que en la k -ésima iteración se realizará un paso de bisección en como máximo${\displaystyle 2\log _{2}(|b_{k-1}-b_{k-2}|/\delta )}$iteraciones adicionales, porque las condiciones anteriores obligan a que los tamaños de paso de interpolación consecutivos se reduzcan a la mitad cada dos iteraciones, y después de como máximo${\displaystyle 2\log _{2}(|b_{k-1}-b_{k-2}|/\delta )}$iteraciones, el tamaño del paso será menor que${\displaystyle \delta }$, que invoca un paso de bisección. Brent demostró que su método requiere como máximo N ² iteraciones, donde N denota el número de iteraciones para el método de bisección. Si la función f se comporta bien, entonces el método de Brent generalmente procederá mediante interpolación cuadrática inversa o lineal, en cuyo caso convergerá superlinealmente .

Además, el método de Brent utiliza interpolación cuadrática inversa en lugar de interpolación lineal (como la que utiliza el método de la secante). Si f ( b _k ), f ( a _k ) y f ( b _{k − 1} ) son distintos, aumenta ligeramente la eficiencia. En consecuencia, la condición para aceptar s (el valor propuesto por la interpolación lineal o la interpolación cuadrática inversa) debe cambiarse: s debe estar entre (3 a _k + b _k ) / 4 y b _k .

Ingrese a , b y (un puntero a) una función para f. Calcule f ( a ). Calcula f ( b ) si f ( a ) f ( b ) ≥ 0 entonces Función de salida porque la raíz no está entre corchetes. fin si si | f ( a )| < | f ( b )| entonces intercambiar ( a , b ) fin si c := a establecer mflag repetir hasta que f ( b o s ) = 0 o | b − a | sea suficientemente pequeño (convergencia) si f ( a ) ≠ f ( c ) y f ( b ) ≠ f ( c ) entonces${\textstyle s:={\frac {af(b)f(c)}{(f(a)-f(b))(f(a)-f(c))}}+{\frac {bf(a)f(c)}{(f(b)-f(a))(f(b)-f(c))}}+{\frac {cf(a)f(b)}{(f(c)-f(a))(f(c)-f(b))}}}$( interpolación cuadrática inversa ) de lo contrario${\textstyle s:=b-f(b){\frac {b-a}{f(b)-f(a)}}}$( método de la secante ) fin si si (condición 1) s no está entre${\displaystyle (3a+b)/4}$y b o (condición 2) (mflag está activado y | s − b | ≥ | b − c |/2) o (condición 3) (mflag está desactivado y | s − b | ≥ | c − d |/2) o (condición 4) (mflag está activado y | b − c | < | δ |) o (condición 5) (mflag está desactivado y | c − d | < | δ |) entonces${\textstyle s:={\frac {a+b}{2}}}$( método de bisección ) establecer mflag sino borrar mflag fin si calcular f ( s ) d := c (d se asigna por primera vez aquí; no se usará arriba en la primera iteración porque mflag está establecido) c := b si f ( a ) f ( s ) < 0 entonces b := s sino a := s fin si si | f ( a )| < | f ( b )| entonces intercambiar ( a , b ) fin si fin repetir imprimir b o s (devolver la raíz)

Supongamos que buscamos un cero de la función definida por f ( x ) = ( x + 3)( x − 1) ² .

Tomamos [ a ₀ , b ₀ ] = [ − 4, 4/3] como nuestro intervalo inicial.

Tenemos f ( a ₀ ) = − 25 y f ( b ₀ ) = 0.48148 (todos los números en esta sección están redondeados), por lo que se cumplen las condiciones f ( a ₀ ) f ( b ₀ ) < 0 y | f ( b ₀ )| ≤ | f ( a _{0 )|.}

1. En la primera iteración, usamos interpolación lineal entre ( b _{− 1} , f ( b _{− 1} )) = ( a ₀ , f ( a ₀ )) = ( − 4, − 25) y ( b ₀ , f ( b ₀ )) = (1.33333, 0.48148), lo que produce s = 1.23256. Esto se encuentra entre (3 a ₀ + b ₀ ) / 4 y b ₀ , por lo que este valor es aceptado. Además, f (1.23256) = 0.22891, por lo que establecemos a ₁ = a ₀ y b ₁ = s = 1.23256.
2. En la segunda iteración, utilizamos la interpolación cuadrática inversa entre ( a ₁ , f ( a ₁ )) = ( − 4, − 25) y ( b ₀ , f ( b ₀ )) = (1.33333, 0.48148) y ( b ₁ , f ( b ₁ )) = (1.23256, 0.22891). Esto da como resultado 1.14205, que se encuentra entre (3 a ₁ + b ₁ ) / 4 y b ₁ . Además, se cumple la desigualdad |1.14205 − b ₁ | ≤ | b ₀ − b _{− 1} | / 2, por lo que se acepta este valor. Además, f (1.14205) = 0.083582, por lo que establecemos a ₂ = a ₁ y b ₂ = 1.14205.
3. En la tercera iteración, usamos interpolación cuadrática inversa entre ( a ₂ , f ( a ₂ )) = ( − 4, − 25) y ( b ₁ , f ( b ₁ )) = (1.23256, 0.22891) y ( b ₂ , f ( b ₂ )) = (1.14205, 0.083582). Esto produce 1.09032, que se encuentra entre (3 a ₂ + b ₂ ) / 4 y b ₂ . Pero aquí entra en juego la condición adicional de Brent: la desigualdad |1.09032 − b ₂ | ≤ | b ₁ − b ₀ | / 2 no se satisface, por lo que este valor se rechaza. En cambio, se calcula el punto medio m = − 1,42897 del intervalo [ a ₂ , b ₂ ]. Tenemos f ( m ) = 9,26891, así que establecemos a ₃ = a ₂ y b ₃ = − 1,42897.
4. En la cuarta iteración, usamos interpolación cuadrática inversa entre ( a ₃ , f ( a ₃ )) = ( − 4, − 25) y ( b ₂ , f ( b ₂ )) = (1.14205, 0.083582) y ( b ₃ , f ( b ₃ )) = ( − 1.42897, 9.26891). Esto produce 1.15448, que no está en el intervalo entre (3 a ₃ + b ₃ ) / 4 y b ₃ ). Por lo tanto, se reemplaza por el punto medio m = − 2.71449. Tenemos f ( m ) = 3.93934, así que establecemos a ₄ = a ₃ y b ₄ = − 2.71449.
5. En la quinta iteración, la interpolación cuadrática inversa produce − 3.45500, que se encuentra en el intervalo requerido. Sin embargo, la iteración anterior fue un paso de bisección, por lo que se debe satisfacer la desigualdad | − 3.45500 − b ₄ | ≤ | b ₄ − b ₃ | / 2. Esta desigualdad es falsa, por lo que usamos el punto medio m = − 3.35724. Tenemos f ( m ) = − 6.78239, por lo que m se convierte en el nuevo contrapunto ( a ₅ = − 3.35724) y la iteración permanece igual ( b ₅ = b ₄ ).
6. En la sexta iteración, no podemos usar la interpolación cuadrática inversa porque b ₅ = b ₄ . Por lo tanto, usamos la interpolación lineal entre ( a ₅ , f ( a ₅ )) = ( − 3.35724, − 6.78239) y ( b ₅ , f ( b ₅ )) = ( − 2.71449, 3.93934). El resultado es s = − 2.95064, que satisface todas las condiciones. Pero como la iteración no cambió en el paso anterior, rechazamos este resultado y volvemos a la bisección. Actualizamos s = − 3.03587, y f ( s ) = − 0.58418.
7. En la séptima iteración, podemos usar nuevamente la interpolación cuadrática inversa. El resultado es s = − 3.00219, que satisface todas las condiciones. Ahora, f ( s ) = − 0.03515, por lo que establecemos a ₇ = b ₆ y b ₇ = − 3.00219 ( a ₇ y b ₇ se intercambian de modo que se cumpla la condición | f ( b ₇ )| ≤ | f ( a ₇ )|). ( Correcto  : interpolación lineal )${\displaystyle s=-2.99436,f(s)=0.089961}$)​
8. En la octava iteración, no podemos usar la interpolación cuadrática inversa porque a ₇ = b ₆ . La interpolación lineal produce s = − 2.99994, que es aceptado. ( Correcto  :${\displaystyle s=-2.9999,f(s)=0.0016}$)​
9. En las siguientes iteraciones, la raíz x = − 3 se aproxima rápidamente: b ₉ = − 3 + 6·10 ^{− 8} y b ₁₀ = − 3 − 3·10 ^{− 15} . ( Correcto  : Iteración 9  : f ( s ) = −1,4 × 10 ⁻⁷ , Iteración 10  : f ( s ) = 6,96 × 10 ⁻¹² )

• Brent (1973) publicó una implementación de Algol 60 .
• Netlib contiene una traducción al Fortran de esta implementación con ligeras modificaciones.
• El método PARI/GPsolve implementa el método.
• Otras implementaciones del algoritmo (en C++, C y Fortran) se pueden encontrar en los libros de Numerical Recipes .
• La biblioteca Apache Commons Math implementa el algoritmo en Java .
• El módulo de optimización de SciPy implementa el algoritmo en Python (lenguaje de programación).
• La biblioteca estándar de Modelica implementa el algoritmo en Modelica .
• La unirootfunción implementa el algoritmo en R (software) .
• La fzerofunción implementa el algoritmo en MATLAB .
• La biblioteca Boost (bibliotecas de C++) implementa dos algoritmos basados ​​en el método de Brent en C++ en el conjunto de herramientas matemáticas:
  1. Minimización de funciones en minima.hpp con un ejemplo que localiza los mínimos de la función .
  2. La búsqueda de raíces implementa el nuevo TOMS748, un algoritmo más moderno y eficiente que el original de Brent, en TOMS748 , y la búsqueda de raíces de Boost.Math que utiliza TOMS748 internamente con ejemplos .
• El paquete Optim.jl implementa el algoritmo en Julia (lenguaje de programación).
• El sistema de álgebra computacional Emmy (escrito en Clojure (lenguaje de programación) ) implementa una variante del algoritmo diseñado para la minimización de funciones univariadas.
• Búsqueda de raíz en la biblioteca C# alojada en Code Project .

• Brent, RP (1973), "Capítulo 4: Un algoritmo con convergencia garantizada para encontrar un cero de una función", Algoritmos para la minimización sin derivadas , Englewood Cliffs, NJ: Prentice-Hall, ISBN 0-13-022335-2
• Dekker, TJ (1969), "Encontrar un cero mediante interpolación lineal sucesiva", en Dejon, B.; Henrici, P. (eds.), Aspectos constructivos del teorema fundamental del álgebra , Londres: Wiley-Interscience, ISBN 978-0-471-20300-1

• Atkinson, Kendall E. (1989). «Sección 2.8.». Introducción al análisis numérico (2.ª  ed.). John Wiley and Sons. ISBN 0-471-50023-2.
• Press, WH; Teukolsky, SA; Vetterling, WT; Flannery, BP (2007). «Sección 9.3. Método de Van Wijngaarden–Dekker–Brent» . Numerical Recipes: The Art of Scientific Computing (3.ª  ed.). Nueva York: Cambridge University Press. ISBN 978-0-521-88068-8Archivado del original el 11 de agosto de 2011. Consultado el 28 de febrero de 2012 .
• Alefeld, GE; Potra, FA; Shi, Yixun (septiembre de 1995). "Algoritmo 748: encerrar ceros de funciones continuas" . Transacciones ACM sobre software matemático . 21 (3): 327– 344. doi : 10.1145/210089.210111 . S2CID 207192624 . 

• zeroin.f en Netlib .
• Módulo brent en C++ (también C, Fortran, Matlab). Archivado el 5 de abril de 2018 en Wayback Machine por John Burkardt.
• Implementación de GSL .
• Implementación de Boost en C++ .
• Implementación en Python (Scipy)