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 theBrent'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 instead of linear interpolation (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)'' else ''( secant method)'' 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)'' 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 implementation. * Netlib contains a Fortran translation of this implementation with slight modifications. * TheReferences
* *Further reading
* * *External links