In
numerical analysis
Numerical analysis is the study of algorithms that use numerical approximation (as opposed to symbolic computation, symbolic manipulations) for the problems of mathematical analysis (as distinguished from discrete mathematics). It is the study of ...
, Brent's method is a hybrid
root-finding algorithm
In mathematics and computing, a root-finding algorithm is an algorithm for finding zeros, also called "roots", of continuous functions. A zero of a function , from the real numbers to real numbers or from the complex numbers to the complex numbers ...
combining the
bisection method
In mathematics, the bisection method is a root-finding method that applies to any continuous function for which one knows two values with opposite signs. The method consists of repeatedly bisecting the interval defined by these values and the ...
, the
secant method
In numerical analysis, the secant method is a root-finding algorithm that uses a succession of roots of secant lines to better approximate a root of a function ''f''. The secant method can be thought of as a finite-difference approximation of ...
and
inverse quadratic interpolation In numerical analysis, inverse quadratic interpolation is a root-finding algorithm, meaning that it is an algorithm for solving equations of the form ''f''(''x'') = 0. The idea is to use quadratic interpolation to approximate the inverse of ''f''. ...
. It has the reliability of bisection but it can be as quick as some of the less-reliable methods. The algorithm tries to use the potentially fast-converging secant method or inverse quadratic interpolation if possible, but it falls back to the more robust bisection method if necessary. Brent's method is due to
Richard Brent and builds on an earlier algorithm by
Theodorus Dekker
Theodorus Jozef Dekker (Dirk Dekker, 1 March 1927 - 25 November 2021) was a Dutch mathematician.
Dekker completed his Ph.D. degree from the University of Amsterdam in 1958. His thesis was titled "Paradoxical Decompositions of Sets and Spaces".
...
. Consequently, the method is also known as the Brent–Dekker method.
Modern improvements on Brent's method include
Chandrupatla's method, which is simpler and faster for functions that are flat around their roots;
Ridders' method In numerical analysis, Ridders' method is a root-finding algorithm based on the false position method and the use of an exponential function to successively approximate a root of a continuous function f(x). The method is due to C. Ridders.
Ridders' ...
, which performs exponential interpolations instead of quadratic providing a simpler closed formula for the iterations; and the
ITP method which is a hybrid between regula-falsi and bisection that achieves optimal worst-case and asymptotic guarantees.
Dekker's method
The idea to combine the bisection method with the secant method goes back to .
Suppose that we want to solve the equation ''f''(''x'') = 0. As with the bisection method, we need to initialize Dekker's method with two points, say ''a''
0 and ''b''
0, such that ''f''(''a''
0) and ''f''(''b''
0) have opposite signs. If ''f'' is continuous on
0, ''b''0">'a''0, ''b''0 the
intermediate value theorem
In mathematical analysis, the intermediate value theorem states that if f is a continuous function whose domain contains the interval , then it takes on any given value between f(a) and f(b) at some point within the interval.
This has two import ...
guarantees the existence of a solution between ''a''
0 and ''b''
0.
Three points are involved in every iteration:
* ''b''
''k'' is the current iterate, i.e., the current guess for the root of ''f''.
* ''a''
''k'' is the "contrapoint," i.e., a point such that ''f''(''a''
''k'') and ''f''(''b''
''k'') have opposite signs, so the interval
''k'', ''b''''k''">'a''''k'', ''b''''k''contains the solution. Furthermore, , ''f''(''b''
''k''), should be less than or equal to , ''f''(''a''
''k''), , so that ''b''
''k'' is a better guess for the unknown solution than ''a''
''k''.
* ''b''
''k''−1 is the previous iterate (for the first iteration, we set ''b''
''k''−1 = ''a''
0).
Two provisional values for the next iterate are computed. The first one is given by linear interpolation, also known as the secant method:
::
and the second one is given by the bisection method
::
If the result of the secant method, ''s'', lies strictly between ''b''
''k'' and ''m'', then it becomes the next iterate (''b''
''k''+1 = ''s''), otherwise the midpoint is used (''b''
''k''+1 = ''m'').
Then, the value of the new contrapoint is chosen such that ''f''(''a''
''k''+1) and ''f''(''b''
''k''+1) have opposite signs. If ''f''(''a''
''k'') and ''f''(''b''
''k''+1) have opposite signs, then the contrapoint remains the same: ''a''
''k''+1 = ''a''
''k''. Otherwise, ''f''(''b''
''k''+1) and ''f''(''b''
''k'') have opposite signs, so the new contrapoint becomes ''a''
''k''+1 = ''b''
''k''.
Finally, if , ''f''(''a''
''k''+1), < , ''f''(''b''
''k''+1), , then ''a''
''k''+1 is probably a better guess for the solution than ''b''
''k''+1, and hence the values of ''a''
''k''+1 and ''b''
''k''+1 are exchanged.
This ends the description of a single iteration of Dekker's method.
Dekker's method performs well if the function ''f'' is reasonably well-behaved. However, there are circumstances in which every iteration employs the secant method, but the iterates ''b''
''k'' converge very slowly (in particular, , ''b''
''k'' − ''b''
''k''−1, may be arbitrarily small). Dekker's method requires far more iterations than the bisection method in this case.
Brent's method
proposed a small modification to avoid the problem with Dekker's method. He inserts an additional test which must be satisfied before the result of the secant method is accepted as the next iterate. Two inequalities must be simultaneously satisfied:
Given a specific numerical tolerance
, if the previous step used the bisection method, the inequality
must hold to perform interpolation, otherwise the bisection method is performed and its result used for the next iteration.
If the previous step performed interpolation, then the inequality
is used instead to perform the next action (to choose) interpolation (when inequality is true) or bisection method (when inequality is not true).
Also, if the previous step used the bisection method, the inequality
must hold, otherwise the bisection method is performed and its result used for the next iteration. If the previous step performed interpolation, then the inequality
is used instead.
This modification ensures that at the kth iteration, a bisection step will be performed in at most
additional iterations, because the above conditions force consecutive interpolation step sizes to halve every two iterations, and after at most
iterations, the step size will be smaller than
, which invokes a bisection step. Brent proved that his method requires at most ''N''
2 iterations, where ''N'' denotes the number of iterations for the bisection method. If the function ''f'' is well-behaved, then Brent's method will usually proceed by either inverse quadratic or linear interpolation, in which case it will converge
superlinearly.
Furthermore, Brent's method uses
inverse quadratic interpolation In numerical analysis, inverse quadratic interpolation is a root-finding algorithm, meaning that it is an algorithm for solving equations of the form ''f''(''x'') = 0. The idea is to use quadratic interpolation to approximate the inverse of ''f''. ...
instead of
linear interpolation
In mathematics, linear interpolation is a method of curve fitting using linear polynomials to construct new data points within the range of a discrete set of known data points.
Linear interpolation between two known points
If the two known poin ...
(as used by the secant method). If ''f''(''b''
''k''), ''f''(''a''
''k'') and ''f''(''b''
''k''−1) are distinct, it slightly increases the efficiency. As a consequence, the condition for accepting ''s'' (the value proposed by either linear interpolation or inverse quadratic interpolation) has to be changed: ''s'' has to lie between (3''a''
''k'' + ''b''
''k'') / 4 and ''b''
''k''.
Algorithm
input ''a'', ''b'', and (a pointer to) a function for ''f''
calculate ''f''(''a'')
calculate ''f''(''b'')
if ''f''(''a'')''f''(''b'') ≥ 0 then
exit function because the root is not bracketed.
end if
if , ''f''(''a''), < , ''f''(''b''), then
swap (''a'',''b'')
end if
''c'' := ''a''
set mflag
repeat until ''f''(''b'' or ''s'') = 0 or , ''b'' − ''a'', is small enough ''(convergence)''
if ''f''(''a'') ≠ ''f''(''c'') and ''f''(''b'') ≠ ''f''(''c'') then
''(
inverse quadratic interpolation In numerical analysis, inverse quadratic interpolation is a root-finding algorithm, meaning that it is an algorithm for solving equations of the form ''f''(''x'') = 0. The idea is to use quadratic interpolation to approximate the inverse of ''f''. ...
)''
else
''(
secant method
In numerical analysis, the secant method is a root-finding algorithm that uses a succession of roots of secant lines to better approximate a root of a function ''f''. The secant method can be thought of as a finite-difference approximation of ...
)''
end if
if ''(condition 1)'' ''s'' is not or
''(condition 2)'' (mflag is set and
, ''s''−''b'', ≥ , ''b''−''c'', /2) or
''(condition 3)'' (mflag is cleared and
, ''s''−''b'', ≥ , ''c''−''d'', /2) or
''(condition 4)'' (mflag is set and
, ''b''−''c'', < , , ) or
''(condition 5)'' (mflag is cleared and
, ''c''−''d'', < , , ) then
''(
bisection method
In mathematics, the bisection method is a root-finding method that applies to any continuous function for which one knows two values with opposite signs. The method consists of repeatedly bisecting the interval defined by these values and the ...
)''
set mflag
else
clear mflag
end if
calculate ''f''(''s'')
''d'' := ''c'' ''(d is assigned for the first time here; it won't be used above on the first iteration because mflag is set)''
''c'' := ''b''
if ''f''(''a'')''f''(''s'') < 0 then
''b'' := ''s''
else
''a'' := ''s''
end if
if , ''f''(''a''), < , ''f''(''b''), then
swap (''a'',''b'')
end if
end repeat
output ''b'' ''or s (return the root)''
Example
Suppose that we are seeking a zero of the function defined by ''f''(''x'') = (''x'' + 3)(''x'' − 1)
2.
We take
0, ''b''0">'a''0, ''b''0=
minus;4, 4/3'' as our initial interval.
We have ''f''(''a''
0) = −25 and ''f''(''b''
0) = 0.48148 (all numbers in this section are rounded), so the conditions ''f''(''a''
0) ''f''(''b''
0) < 0 and , ''f''(''b''
0), ≤ , ''f''(''a''
0), are satisfied.
# In the first iteration, we use linear interpolation between (''b''
−1, ''f''(''b''
−1)) = (''a''
0, ''f''(''a''
0)) = (−4, −25) and (''b''
0, ''f''(''b''
0)) = (1.33333, 0.48148), which yields ''s'' = 1.23256. This lies between (3''a''
0 + ''b''
0) / 4 and ''b''
0, so this value is accepted. Furthermore, ''f''(1.23256) = 0.22891, so we set ''a''
1 = ''a''
0 and ''b''
1 = ''s'' = 1.23256.
# In the second iteration, we use inverse quadratic interpolation between (''a''
1, ''f''(''a''
1)) = (−4, −25) and (''b''
0, ''f''(''b''
0)) = (1.33333, 0.48148) and (''b''
1, ''f''(''b''
1)) = (1.23256, 0.22891). This yields 1.14205, which lies between (3''a''
1 + ''b''
1) / 4 and ''b''
1. Furthermore, the inequality , 1.14205 − ''b''
1, ≤ , ''b''
0 − ''b''
−1, / 2 is satisfied, so this value is accepted. Furthermore, ''f''(1.14205) = 0.083582, so we set ''a''
2 = ''a''
1 and ''b''
2 = 1.14205.
# In the third iteration, we use inverse quadratic interpolation between (''a''
2, ''f''(''a''
2)) = (−4, −25) and (''b''
1, ''f''(''b''
1)) = (1.23256, 0.22891) and (''b''
2, ''f''(''b''
2)) = (1.14205, 0.083582). This yields 1.09032, which lies between (3''a''
2 + ''b''
2) / 4 and ''b''
2. But here Brent's additional condition kicks in: the inequality , 1.09032 − ''b''
2, ≤ , ''b''
1 − ''b''
0, / 2 is not satisfied, so this value is rejected. Instead, the midpoint ''m'' = −1.42897 of the interval
2, ''b''2">'a''2, ''b''2is computed. We have ''f''(''m'') = 9.26891, so we set ''a''
3 = ''a''
2 and ''b''
3 = −1.42897.
# In the fourth iteration, we use inverse quadratic interpolation between (''a''
3, ''f''(''a''
3)) = (−4, −25) and (''b''
2, ''f''(''b''
2)) = (1.14205, 0.083582) and (''b''
3, ''f''(''b''
3)) = (−1.42897, 9.26891). This yields 1.15448, which is not in the interval between (3''a''
3 + ''b''
3) / 4 and ''b''
3). Hence, it is replaced by the midpoint ''m'' = −2.71449. We have ''f''(''m'') = 3.93934, so we set ''a''
4 = ''a''
3 and ''b''
4 = −2.71449.
# In the fifth iteration, inverse quadratic interpolation yields −3.45500, which lies in the required interval. However, the previous iteration was a bisection step, so the inequality , −3.45500 − ''b''
4, ≤ , ''b''
4 − ''b''
3, / 2 need to be satisfied. This inequality is false, so we use the midpoint ''m'' = −3.35724. We have ''f''(''m'') = −6.78239, so ''m'' becomes the new contrapoint (''a''
5 = −3.35724) and the iterate remains the same (''b''
5 = ''b''
4).
# In the sixth iteration, we cannot use inverse quadratic interpolation because ''b''
5 = ''b''
4. Hence, we use linear interpolation between (''a''
5, ''f''(''a''
5)) = (−3.35724, −6.78239) and (''b''
5, ''f''(''b''
5)) = (−2.71449, 3.93934). The result is ''s'' = −2.95064, which satisfies all the conditions. But since the iterate did not change in the previous step, we reject this result and fall back to bisection. We update ''s'' = -3.03587, and ''f''(''s'') = -0.58418.
# In the seventh iteration, we can again use inverse quadratic interpolation. The result is ''s'' = −3.00219, which satisfies all the conditions. Now, ''f''(''s'') = −0.03515, so we set ''a''
7 = ''b''
6 and ''b''
7 = −3.00219 (''a''
7 and ''b''
7 are exchanged so that the condition , ''f''(''b''
7), ≤ , ''f''(''a''
7), is satisfied). (''Correct'' : linear interpolation )
# In the eighth iteration, we cannot use inverse quadratic interpolation because ''a''
7 = ''b''
6. Linear interpolation yields ''s'' = −2.99994, which is accepted. (''Correct'' : )
# In the following iterations, the root ''x'' = −3 is approached rapidly: ''b''
9 = −3 + 6·10
−8 and ''b''
10 = −3 − 3·10
−15. (''Correct'' : Iter 9 : ''f''(''s'') = −1.4 × 10
−7, Iter 10 : ''f''(''s'') = 6.96 × 10
−12)
Implementations
* published an
Algol 60
ALGOL 60 (short for ''Algorithmic Language 1960'') is a member of the ALGOL family of computer programming languages. It followed on from ALGOL 58 which had introduced code blocks and the begin and end pairs for delimiting them, representing a k ...
implementation.
*
Netlib Netlib is a repository of software for scientific computing maintained by AT&T, Bell Laboratories, the University of Tennessee and Oak Ridge National Laboratory. Netlib comprises many separate programs and libraries. Most of the code is written in ...
contains a Fortran translation of this implementation with slight modifications.
* The
PARI/GP
PARI/GP is a computer algebra system with the main aim of facilitating number theory computations. Versions 2.1.0 and higher are distributed under the GNU General Public License. It runs on most common operating systems.
System overview
The P ...
method implements the method.
* Other implementations of the algorithm (in C++, C, and Fortran) can be found in the
Numerical Recipes
''Numerical Recipes'' is the generic title of a series of books on algorithms and numerical analysis by William H. Press, Saul A. Teukolsky, William T. Vetterling and Brian P. Flannery. In various editions, the books have been in print since 1 ...
books.
* The
Apache Commons The Apache Commons is a project of the Apache Software Foundation, formerly under the Jakarta Project. The purpose of the Commons is to provide reusable, open source Java software. The Commons is composed of three parts: proper, sandbox, and dorm ...
Math library implements the algorithm in
Java
Java (; id, Jawa, ; jv, ꦗꦮ; su, ) is one of the Greater Sunda Islands in Indonesia. It is bordered by the Indian Ocean to the south and the Java Sea to the north. With a population of 151.6 million people, Java is the world's List ...
.
* The
SciPy
SciPy (pronounced "sigh pie") is a free and open-source Python library used for scientific computing and technical computing.
SciPy contains modules for optimization, linear algebra, integration, interpolation, special functions, FFT, signal ...
optimize module implements the algorithm in
Python (programming language)
Python is a high-level, general-purpose programming language. Its design philosophy emphasizes code readability with the use of significant indentation.
Python is dynamically-typed and garbage-collected. It supports multiple programming para ...
* The Modelica Standard Library implements the algorithm in
Modelica
Modelica is an object-oriented, declarative, multi-domain modeling language for component-oriented modeling of complex systems, e.g., systems containing mechanical, electrical, electronic, hydraulic, thermal, control, electric power or process-o ...
.
* The function implements the algorithm in
R (software).
* The function implements the algorithm in
MATLAB
MATLAB (an abbreviation of "MATrix LABoratory") is a proprietary multi-paradigm programming language and numeric computing environment developed by MathWorks. MATLAB allows matrix manipulations, plotting of functions and data, implementation ...
.
* The
Boost (C++ libraries)
Boost, boosted or boosting may refer to:
Science, technology and mathematics
* Boost, positive manifold pressure in turbocharged engines
* Boost (C++ libraries), a set of free peer-reviewed portable C++ libraries
* Boost (material), a material ...
implements two algorithms based on Brent's method in
C++
C++ (pronounced "C plus plus") is a high-level general-purpose programming language created by Danish computer scientist Bjarne Stroustrup as an extension of the C programming language, or "C with Classes". The language has expanded significan ...
in the Math toolkit:
# Function minimization a
minima.hppwith an exampl
# Root finding implements the newer TOMS748, a more modern and efficient algorithm than Brent's original, a
TOMS748 an
tha
wit
examples
* Th
Optim.jlpackage implements the algorithm in
Julia (programming language)
Julia is a high-level, dynamic programming language. Its features are well suited for numerical analysis and computational science.
Distinctive aspects of Julia's design include a type system with parametric polymorphism in a dynamic programmi ...
* Th
SICMUtilscomputer algebra system (written in
Clojure (programming language)
Clojure (, like ''closure'') is a dynamic and functional dialect of the Lisp programming language on the Java platform. Like other Lisp dialects, Clojure treats code as data and has a Lisp macro system. The current development process is commu ...
) implements a variant of the algorithm designed for univariate function minimization.
Root-Finding in C#library hosted in Code Project.
References
*
*
Further reading
*
*
*
External links
zeroin.fat
Netlib Netlib is a repository of software for scientific computing maintained by AT&T, Bell Laboratories, the University of Tennessee and Oak Ridge National Laboratory. Netlib comprises many separate programs and libraries. Most of the code is written in ...
.
module brent in C++ (also C, Fortran, Matlab)by John Burkardt
implementation.
implementation.
implementation
{{Root-finding algorithms
Root-finding algorithms