Brent's method

<h2 id="dekkers-method">Dekker's method</h2>
The idea to combine the bisection method with the secant method goes back to Dekker (1969).
Suppose that one wants to solve the equation f(x) = 0. As with the bisection method, one needs to initialize Dekker's method with two points, say a0 and b0, such that f(a0) and f(b0) have opposite signs. If f is continuous on [a0, b0], the <a href="/facts/Intermediate_value_theorem/SaFdYE7W">intermediate value theorem</a> guarantees the existence of a solution between a0 and b0.
Three points are involved in every iteration:

<ul><li>bk is the current iterate, i.e., the current guess for the root of f.</li>
<li>ak is the "contrapoint", i.e., a point such that f(ak) and f(bk) have opposite signs, so the interval [ak, bk] contains the solution. Furthermore, |f(bk)| should be less than or equal to |f(ak)|, so that bk is a better guess for the unknown solution than ak.</li>
<li>bk−1 is the previous iterate (for the first iteration, one sets bk−1 = a0).</li></ul>
Two provisional values for the next iterate are computed. The first one is given by linear interpolation, also known as the secant method:

s
        =
        
          
            {
            
              
                
                  
                    b
                    
                      k
                    
                  
                  −
                  
                    
                      
                        
                          b
                          
                            k
                          
                        
                        −
                        
                          b
                          
                            k
                            −
                            1
                          
                        
                      
                      
                        f
                        (
                        
                          b
                          
                            k
                          
                        
                        )
                        −
                        f
                        (
                        
                          b
                          
                            k
                            −
                            1
                          
                        
                        )
                      
                    
                  
                  f
                  (
                  
                    b
                    
                      k
                    
                  
                  )
                  ,
                
                
                  
                    
                      if 
                    
                  
                  f
                  (
                  
                    b
                    
                      k
                    
                  
                  )
                  ≠
                  f
                  (
                  
                    b
                    
                      k
                      −
                      1
                    
                  
                  )
                
              
              
                
                  m
                
                
                  
                    
                      otherwise 
                    
                  
                
              
            
            
          
        
      
    
    {\displaystyle s={\begin{cases}b_{k}-{\frac {b_{k}-b_{k-1}}{f(b_{k})-f(b_{k-1})}}f(b_{k}),&{\mbox{if }}f(b_{k})\neq f(b_{k-1})\\m&{\mbox{otherwise }}\end{cases}}}

and the second one is given by the bisection method

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

If the result of the secant method, s, lies strictly between bk and m, then it becomes the next iterate (bk+1 = s), otherwise the midpoint is used (bk+1 = m).
Then, the value of the new contrapoint is chosen such that f(ak+1) and f(bk+1) have opposite signs. If f(ak) and f(bk+1) have opposite signs, then the contrapoint remains the same: ak+1 = ak. Otherwise, f(bk+1) and f(bk) have opposite signs, so the new contrapoint becomes ak+1 = bk.
Finally, if |f(ak+1)| < |f(bk+1)|, then ak+1 is probably a better guess for the solution than bk+1, and hence the values of ak+1 and bk+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 bk converge very slowly (in particular, |bk − bk−1| may be arbitrarily small). Dekker's method requires far more iterations than the bisection method in this case.

<h2 id="brents-method">Brent's method</h2>
Brent (1973) 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 
 
 
 
 δ
 
 
 {\displaystyle \delta }
 
, if the previous step used the bisection method, the inequality 
 
 
 
 
 |
 
 δ
 
 |
 
 <
 
 |
 
 
 b
 
 k
 
 
 −
 
 b
 
 k
 −
 1
 
 
 
 |
 
 
 
 {\textstyle |\delta |<|b_{k}-b_{k-1}|}
 
 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 
 
 
 
 
 |
 
 δ
 
 |
 
 <
 
 |
 
 
 b
 
 k
 −
 1
 
 
 −
 
 b
 
 k
 −
 2
 
 
 
 |
 
 
 
 {\textstyle |\delta |<|b_{k-1}-b_{k-2}|}
 
 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 
 
 
 
 
 |
 
 s
 −
 
 b
 
 k
 
 
 
 |
 
 <
 
 
 
 
 
 
 1
 2
 
 
 
 
 
 
 
 |
 
 
 b
 
 k
 
 
 −
 
 b
 
 k
 −
 1
 
 
 
 |
 
 
 
 {\textstyle |s-b_{k}|<{\begin{matrix}{\frac {1}{2}}\end{matrix}}|b_{k}-b_{k-1}|}

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 
 
 
 
 
 |
 
 s
 −
 
 b
 
 k
 
 
 
 |
 
 <
 
 
 
 
 
 
 1
 2
 
 
 
 
 
 
 
 |
 
 
 b
 
 k
 −
 1
 
 
 −
 
 b
 
 k
 −
 2
 
 
 
 |
 
 
 
 {\textstyle |s-b_{k}|<{\begin{matrix}{\frac {1}{2}}\end{matrix}}|b_{k-1}-b_{k-2}|}

is used instead.
This modification ensures that at the kth iteration, a bisection step will be performed in at most 
 
 
 
 2
 
 log
 
 2
 
 
 ⁡
 (
 
 |
 
 
 b
 
 k
 −
 1
 
 
 −
 
 b
 
 k
 −
 2
 
 
 
 |
 
 
 /
 
 δ
 )
 
 
 {\displaystyle 2\log _{2}(|b_{k-1}-b_{k-2}|/\delta )}
 
 additional iterations, because the above conditions force consecutive interpolation step sizes to halve every two iterations, and after at most 
 
 
 
 2
 
 log
 
 2
 
 
 ⁡
 (
 
 |
 
 
 b
 
 k
 −
 1
 
 
 −
 
 b
 
 k
 −
 2
 
 
 
 |
 
 
 /
 
 δ
 )
 
 
 {\displaystyle 2\log _{2}(|b_{k-1}-b_{k-2}|/\delta )}
 
 iterations, the step size will be smaller than 
 
 
 
 δ
 
 
 {\displaystyle \delta }
 
, which invokes a bisection step. Brent proved that his method requires at most N2 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 <a href="/facts/Rate_of_convergence/JVGuzPoS">superlinearly</a>.
Furthermore, Brent's method uses <a href="/facts/Inverse_quadratic_interpolation/7fQrKtYn">inverse quadratic interpolation</a> instead of <a href="/facts/Linear_interpolation/galkvRGy">linear interpolation</a> (as used by the secant method). If f(bk), f(ak) and f(bk−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 (3ak + bk) / 4 and bk.

<h2 id="algorithm">Algorithm</h2>
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
 
 
 
 
 s
 :=
 
 
 
 a
 f
 (
 b
 )
 f
 (
 c
 )
 
 
 (
 f
 (
 a
 )
 −
 f
 (
 b
 )
 )
 (
 f
 (
 a
 )
 −
 f
 (
 c
 )
 )
 
 
 
 +
 
 
 
 b
 f
 (
 a
 )
 f
 (
 c
 )
 
 
 (
 f
 (
 b
 )
 −
 f
 (
 a
 )
 )
 (
 f
 (
 b
 )
 −
 f
 (
 c
 )
 )
 
 
 
 +
 
 
 
 c
 f
 (
 a
 )
 f
 (
 b
 )
 
 
 (
 f
 (
 c
 )
 −
 f
 (
 a
 )
 )
 (
 f
 (
 c
 )
 −
 f
 (
 b
 )
 )
 
 
 
 
 
 {\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))}}}
 
 (<a href="/facts/Inverse_quadratic_interpolation/7fQrKtYn">inverse quadratic interpolation</a>)
 else
 
 
 
 
 s
 :=
 b
 −
 f
 (
 b
 )
 
 
 
 b
 −
 a
 
 
 f
 (
 b
 )
 −
 f
 (
 a
 )
 
 
 
 
 
 {\textstyle s:=b-f(b){\frac {b-a}{f(b)-f(a)}}}
 
 (<a href="/facts/Secant_method/tkNV8Rzk">secant method</a>)
 end if
 if (condition 1) s is not between 
 
 
 
 (
 3
 a
 +
 b
 )
 
 /
 
 4
 
 
 {\displaystyle (3a+b)/4}
 
 and b 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
 
 
 
 
 s
 :=
 
 
 
 a
 +
 b
 
 2
 
 
 
 
 {\textstyle s:={\frac {a+b}{2}}}
 
 (<a href="/facts/Bisection_method/3tWmtoIp">bisection method</a>)
 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)

<h2 id="example">Example</h2>
Suppose that we are seeking a zero of the function defined by f(x) = (x + 3)(x − 1)2.
We take [a0, b0] = [−4, 4/3] as our initial interval.
We have f(a0) = −25 and f(b0) = 0.48148 (all numbers in this section are rounded), so the conditions f(a0) f(b0) < 0 and |f(b0)| ≤ |f(a0)| are satisfied.

<ol><li>In the first iteration, we use linear interpolation between (b−1, f(b−1)) = (a0, f(a0)) = (−4, −25) and (b0, f(b0)) = (1.33333, 0.48148), which yields s = 1.23256. This lies between (3a0 + b0) / 4 and b0, so this value is accepted. Furthermore, f(1.23256) = 0.22891, so we set a1 = a0 and b1 = s = 1.23256.</li>
<li>In the second iteration, we use inverse quadratic interpolation between (a1, f(a1)) = (−4, −25) and (b0, f(b0)) = (1.33333, 0.48148) and (b1, f(b1)) = (1.23256, 0.22891). This yields 1.14205, which lies between (3a1 + b1) / 4 and b1. Furthermore, the inequality |1.14205 − b1| ≤ |b0 − b−1| / 2 is satisfied, so this value is accepted. Furthermore, f(1.14205) = 0.083582, so we set a2 = a1 and b2 = 1.14205.</li>
<li>In the third iteration, we use inverse quadratic interpolation between (a2, f(a2)) = (−4, −25) and (b1, f(b1)) = (1.23256, 0.22891) and (b2, f(b2)) = (1.14205, 0.083582). This yields 1.09032, which lies between (3a2 + b2) / 4 and b2. But here Brent's additional condition kicks in: the inequality |1.09032 − b2| ≤ |b1 − b0| / 2 is not satisfied, so this value is rejected. Instead, the midpoint m = −1.42897 of the interval [a2, b2] is computed. We have f(m) = 9.26891, so we set a3 = a2 and b3 = −1.42897.</li>
<li>In the fourth iteration, we use inverse quadratic interpolation between (a3, f(a3)) = (−4, −25) and (b2, f(b2)) = (1.14205, 0.083582) and (b3, f(b3)) = (−1.42897, 9.26891). This yields 1.15448, which is not in the interval between (3a3 + b3) / 4 and b3). Hence, it is replaced by the midpoint m = −2.71449. We have f(m) = 3.93934, so we set a4 = a3 and b4 = −2.71449.</li>
<li>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 − b4| ≤ |b4 − b3| / 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 (a5 = −3.35724) and the iterate remains the same (b5 = b4).</li>
<li>In the sixth iteration, we cannot use inverse quadratic interpolation because b5 = b4. Hence, we use linear interpolation between (a5, f(a5)) = (−3.35724, −6.78239) and (b5, f(b5)) = (−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.</li>
<li>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 a7 = b6 and b7 = −3.00219 (a7 and b7 are exchanged so that the condition |f(b7)| ≤ |f(a7)| is satisfied). (Correct : linear interpolation ⁠
 
 
 
 s
 =
 −
 2.99436
 ,
 f
 (
 s
 )
 =
 0.089961
 
 
 {\displaystyle s=-2.99436,f(s)=0.089961}
 
⁠)</li>
<li>In the eighth iteration, we cannot use inverse quadratic interpolation because a7 = b6. Linear interpolation yields s = −2.99994, which is accepted. (Correct : ⁠
 
 
 
 s
 =
 −
 2.9999
 ,
 f
 (
 s
 )
 =
 0.0016
 
 
 {\displaystyle s=-2.9999,f(s)=0.0016}
 
⁠)</li>
<li>In the following iterations, the root x = −3 is approached rapidly: b9 = −3 + 6·10−8 and b10 = −3 − 3·10−15. (Correct : Iter 9 : f(s) = −1.4 × 10−7, Iter 10 : f(s) = 6.96 × 10−12)</li></ol>
<h2 id="implementations">Implementations</h2>
<ul><li>Brent (1973) published an <a href="/facts/Algol_60/QnpxUYQm">Algol 60</a> implementation.</li>
<li><a href="/facts/Netlib/4QQJ2DKk">Netlib</a> contains a Fortran translation of this implementation with slight modifications.</li>
<li>The <a href="/facts/PARI%2fGP/JhTL8ysP">PARI/GP</a> method solve implements the method.</li>
<li>Other implementations of the algorithm (in C++, C, and Fortran) can be found in the <a href="/facts/Numerical_Recipes/MIJJkEau">Numerical Recipes</a> books.</li>
<li>The <a href="/facts/Apache_Commons/wTW3nDuz">Apache Commons</a> Math library implements the algorithm in <a href="/facts/Java_(programming_language)/9ScgFyAL">Java</a>.</li>
<li>The <a href="/facts/SciPy/bNcMQGc0">SciPy</a> optimize module implements the algorithm in <a href="/facts/Python_(programming_language)/YbuGqofa">Python (programming language)</a></li>
<li>The Modelica Standard Library implements the algorithm in <a href="/facts/Modelica/pjWEmtEL">Modelica</a>.</li>
<li>The uniroot function implements the algorithm in <a href="/facts/R_(software)/LSrkr8K8">R (software)</a>.</li>
<li>The fzero function implements the algorithm in <a href="/facts/MATLAB/qPjLISCk">MATLAB</a>.</li>
<li>The <a href="/facts/Boost_(C%252B%252B_libraries)/vGLGG51A">Boost (C++ libraries)</a> implements two algorithms based on Brent's method in <a href="/facts/C%252B%252B/Et0F2qn9">C++</a> in the Math toolkit:
<ol><li>Function minimization at <a href="https://www.boost.org/doc/libs/release/boost/math/tools/minima.hpp">minima.hpp</a> with an example <a href="https://www.boost.org/doc/libs/release/libs/math/doc/html/math_toolkit/brent_minima.html">locating function minima</a>.</li>
<li>Root finding implements the newer TOMS748, a more modern and efficient algorithm than Brent's original, at <a href="https://www.boost.org/doc/libs/release/boost/math/tools/toms748_solve.hpp">TOMS748</a>, and <a href="https://www.boost.org/doc/libs/release/libs/math/doc/html/root_finding.html">Boost.Math rooting finding</a> that <a href="https://www.boost.org/doc/libs/release/libs/math/doc/html/math_toolkit/roots_noderiv/bracket_solve.html">uses TOMS748 internally</a> with <a href="https://www.boost.org/doc/libs/release/libs/math/example/root_finding_example.cpp">examples</a>.</li></ol></li>
<li>The <a href="https://github.com/JuliaNLSolvers/Optim.jl">Optim.jl</a> package implements the algorithm in <a href="/facts/Julia_(programming_language)/AoB0PJ9C">Julia (programming language)</a></li>
<li>The <a href="https://github.com/mentat-collective/emmy">Emmy</a> computer algebra system (written in <a href="/facts/Clojure_(programming_language)/3JRNxKXx">Clojure (programming language)</a>) implements a variant of the algorithm designed for univariate function minimization.</li>
<li><a href="https://www.codeproject.com/Articles/79541/Three-Methods-for-Root-finding-in-C">Root-Finding in C#</a> library hosted in Code Project.</li></ul>

<ul><li><a href="/facts/Richard_Brent_(scientist)/JYLmuA0M">Brent, R. P.</a> (1973), "Chapter 4: An Algorithm with Guaranteed Convergence for Finding a Zero of a Function", Algorithms for Minimization without Derivatives, Englewood Cliffs, NJ: Prentice-Hall, <a href="/facts/ISBN_(identifier)/15AdSPa9">ISBN</a> 0-13-022335-2</li>
<li><a href="/facts/Theodorus_Dekker/L0ALFPKz">Dekker, T. J.</a> (1969), "Finding a zero by means of successive linear interpolation", in Dejon, B.; Henrici, P. (eds.), Constructive Aspects of the Fundamental Theorem of Algebra, London: Wiley-Interscience, <a href="/facts/ISBN_(identifier)/15AdSPa9">ISBN</a> 978-0-471-20300-1</li></ul>
<h2 id="further-reading">Further reading</h2>
<ul><li>Atkinson, Kendall E. (1989). "Section 2.8.". An Introduction to Numerical Analysis (2nd ed.). John Wiley and Sons. <a href="/facts/ISBN_(identifier)/15AdSPa9">ISBN</a> 0-471-50023-2.</li>
<li>Press, W. H.; Teukolsky, S. A.; Vetterling, W. T.; Flannery, B. P. (2007). <a href="https://web.archive.org/web/20110811154417/http://apps.nrbook.com/empanel/index.html#pg=454">"Section 9.3. Van Wijngaarden–Dekker–Brent Method"</a>. Numerical Recipes: The Art of Scientific Computing (3rd ed.). New York: Cambridge University Press. <a href="/facts/ISBN_(identifier)/15AdSPa9">ISBN</a> 978-0-521-88068-8. Archived from <a href="http://apps.nrbook.com/empanel/index.html#pg=454">the original</a> on 2011-08-11. Retrieved 2012-02-28.</li>
<li>Alefeld, G. E.; Potra, F. A.; Shi, Yixun (September 1995). <a href="https://doi.org/10.1145%2F210089.210111">"Algorithm 748: Enclosing Zeros of Continuous Functions"</a>. ACM Transactions on Mathematical Software. 21 (3): 327–344. <a href="/facts/Doi_(identifier)/muM9Etpq">doi</a>:<a href="https://doi.org/10.1145%2F210089.210111">10.1145/210089.210111</a>. <a href="/facts/S2CID_(identifier)/ldJsHa2Y">S2CID</a> <a href="https://api.semanticscholar.org/CorpusID:207192624">207192624</a>.</li></ul>
<h2 id="external-links">External links</h2>
<ul><li><a href="http://www.netlib.org/go/zeroin.f">zeroin.f</a> at <a href="/facts/Netlib/4QQJ2DKk">Netlib</a>.</li>
<li><a href="http://people.sc.fsu.edu/~jburkardt/cpp_src/brent/brent.html">module brent in C++ (also C, Fortran, Matlab)</a> <a href="https://web.archive.org/web/20180405205252/http://people.sc.fsu.edu/~jburkardt/cpp_src/brent/brent.html">Archived</a> 2018-04-05 at the <a href="/facts/Wayback_Machine/nmQ3a6JC">Wayback Machine</a> by John Burkardt</li>
<li><a href="https://www.gnu.org/software/gsl/doc/html/roots.html#c.gsl_root_fsolver_brent">GSL</a> implementation.</li>
<li><a href="https://www.boost.org/doc/libs/1_67_0/libs/math/doc/html/root_finding.html">Boost C++</a> implementation.</li>
<li><a href="https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.brentq.html#scipy.optimize.brentq">Python (Scipy)</a> implementation</li></ul>

<h2 id="references">References</h2>

<ol>
<li id="fn:1">Brent 1973 - Brent, R. P. (1973), "Chapter 4: An Algorithm with Guaranteed Convergence for Finding a Zero of a Function", Algorithms for Minimization without Derivatives, Englewood Cliffs, NJ: Prentice-Hall, ISBN 0-13-022335-2 <a href="#fnref:1" class="footnote-back-ref">↩</a></li>
<li id="fn:2">Dekker 1969 - Dekker, T. J. (1969), "Finding a zero by means of successive linear interpolation", in Dejon, B.; Henrici, P. (eds.), Constructive Aspects of the Fundamental Theorem of Algebra, London: Wiley-Interscience, ISBN 978-0-471-20300-1 <a href="#fnref:2" class="footnote-back-ref">↩</a></li>
<li id="fn:3">Chandrupatla, Tirupathi R. (1997). "A new hybrid quadratic/Bisection algorithm for finding the zero of a nonlinear function without using derivatives". Advances in Engineering Software. 28 (3): 145–149. doi:10.1016/S0965-9978(96)00051-8. <a href="/wiki/Doi_(identifier)" target="_blank">/wiki/Doi_(identifier)</a> <a href="#fnref:3" class="footnote-back-ref">↩</a></li>
<li id="fn:4">"Ten Little Algorithms, Part 5: Quadratic Extremum Interpolation and Chandrupatla's Method - Jason Sachs". <a href="https://www.embeddedrelated.com/showarticle/855.php" target="_blank">https://www.embeddedrelated.com/showarticle/855.php</a> <a href="#fnref:4" class="footnote-back-ref">↩</a></li>
</ol>

Brent's method open-in-new

Brent's method