Brent's method

From Wikipedia, the free encyclopedia

In numerical analysis, Brent's method is a complicated but popular root-finding algorithm combining the bisection method, the secant method and inverse quadratic interpolation. It has the reliability of bisection but it can be as quick as some of the less reliable methods. The idea is to use the secant method or inverse quadratic interpolation if possible, because they converge faster, but to fall back to the more robust bisection method if necessary. Brent's method is due to Richard Brent (1973) and builds on an earlier algorithm of Theodorus Dekker (1969).

Contents 1 Dekker's method 2 Brent's method 3 Example 4 Algorithm 5 Implementations 6 References 7 External link

[edit] Dekker's method

The idea to combine the bisection method with the secant method goes back to Dekker.

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₀ and b₀, such that f(a₀) and f(b₀) have opposite signs. If f is continuous on [a₀, b₀], the intermediate value theorem guarantees the existence of a solution between a₀ and b₀.

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 [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₋₁ = a₀).

Two provisional values for the next iterate are computed. The first one is given by the secant method:

and the second one is given by the bisection method

If the result of the secant method, s, lies 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.

[edit] Brent'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 proposed a small modification to avoid this problem. He inserted an additional test which must be satisfied before the result of the secant method is accepted as the next iterate. Specifically, if the previous step used the bisection method, the inequality

must hold, otherwise the midpoint m is used for the next iterate. If the previous step used interpolation, then the inequality

is used instead. This modification ensures that a bisection step is done if the secant method does not progress fast enough. Brent proved that his method requires at most N² 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) have all different values. This 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 (3a_k + b_k) / 4 and b_k.

[edit] Example

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

1. In the first iteration, we use linear interpolation between (b₋₁, f(b₋₁)) = (a₀, f(a₀)) = (−4, −25) and (b₀, f(b₀)) = (1.33333, 0.48148), which yields s = 1.23256. This lies between (3a₀ + b₀) / 4 and b₀, so this value is accepted. Furthermore, f(1.23256) = 0.22891, so we set a₁ = a₀ and b₁ = s = 1.23256.
2. In the second iteration, we use inverse quadratic interpolation between (a₁, f(a₁)) = (−4, −25) and (b₀, f(b₀)) = (1.33333, 0.48148) and (b₁, f(b₁)) = (1.23256, 0.22891). This yields 1.14205, which lies between (3a₁ + b₁) / 4 and b₁. Furthermore, the inequality |1.14205 − b₁| ≤ |b₀ − b₋₁| / 2 is satisfied, so this value is accepted. Furthermore, f(1.14205) = 0.083582, so we set a₂ = a₁ and b₂ = 1.14205.
3. In the third iteration, we use inverse quadratic interpolation between (a₂, f(a₂)) = (−4, −25) and (b₁, f(b₁)) = (1.23256, 0.22891) and (b₂, f(b₂)) = (1.14205, 0.083582). This yields 1.09032, which lies between (3a₂ + b₂) / 4 and b₂. But here Brent's additional condition kicks in: the inequality |1.09032 − b₂| ≤ |b₁ − b₀| / 2 is not satisfied, so this value is rejected. Instead, the midpoint m = −1.42897 of the interval [a₂, b₂] is computed. We have f(m) = 9.26891, so we set a₃ = a₂ and b₃ = −1.42897.
4. In the fourth iteration, we use inverse quadratic interpolation between (a₃, f(a₃)) = (−4, −25) and (b₂, f(b₂)) = (1.14205, 0.083582) and (b₃, f(b₃)) = (−1.42897, 9.26891). This yields −1.15448, which is not in the interval between (3a₃ + b₃) / 4 and b₃). Hence, it is replaced by the midpoint m = −2.71449. We have f(m) = 3.93934, so we set a₄ = a₃ and b₄ = −2.71449.
5. 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₄| ≤ |b₄ − b₃| / 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₅ = −3.35724) and the iterate remains the same (b₅ = b₄).
6. In the sixth iteration, we cannot use inverse quadratic interpolation because b₅ = b₄. Hence, we use linear interpolation between (a₅, f(a₅)) = (−3.35724, −6.78239) and (b₅, f(b₅)) = (−2.71449, 3.93934). The result is s = −2.95064, which satisfies all the conditions. We calculate f(s) = 0.77032, so we set a₆ = a₅ and b₆ = −2.95064.
7. 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₇ = b₆ and b₇ = −3.00219 (a₇ and b₇ are exchanged so that the condition |f(b₇)| ≤ |f(a₇)| is satisfied).
8. In the eighth iteration, we cannot use inverse quadratic interpolation because a₇ = b₆. Linear interpolation yields s = −2.99994, which is accepted.
9. In the following iterations, the root x = −3 is approached rapidly: b₉ = −3 + 6·10⁻⁸ and b₁₀ = −3 − 3·10⁻¹⁵.

[edit] Algorithm

• input a, b, and a pointer to a subroutine for f
• calculate f(a)
• calculate f(b)
• if f(a) f(b) >= 0 then error-exit end if
• if |f(a)| < |f(b)| then swap (a,b) end if
• c := a
• set mflag
• repeat until f(b) = 0 or |b − a| is small enough (convergence)
  • if f(a) ≠ f(c) and f(b) ≠ f(c) then
    • (inverse quadratic interpolation)
  • else
    • (secant rule)
  • end if
  • if s is not between (3a + b)/4 and b or (mflag is set and |s−b| ≥ |b−c| / 2) or (mflag is cleared and |s−b| ≥ |c−d| / 2) then
    • 
    • set mflag
  • else
    • clear mflag
  • end if
  • calculate f(s)
  • d := c
  • 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 (return the root)

[edit] Implementations

Brent (1973) published an Algol 60 implementation. Netlib contains a Fortran translation of this implementation with slight modifications. The MATLAB function fzero also implements Brent's method. Other implementations of the algorithm (in C++, C, and Fortran flavours) can be found in the Numerical Recipes books.

[edit] References

• Atkinson (1989). An Introduction to Numerical Analysis (2nd ed.), Section 2.8. John Wiley and Sons. ISBN 0-471-50023-2.
• R.P. Brent (1973). Algorithms for Minimization without Derivatives, Chapter 4. Prentice-Hall, Englewood Cliffs, NJ. ISBN 0-13-022335-2.
• T.J. Dekker (1969). "Finding a zero by means of successive linear interpolation." In B. Dejon and P. Henrici (eds), Constructive Aspects of the Fundamental Theorem of Algebra, Wiley-Interscience, London. SBN 471-28300-9.

[edit] External link

• zeroin.f at Netlib.