Brent’s Method with MATLAB | MATLAB Tutorial
by Gunjan Gupta
Brent’s method is a “hybrid” method that combines aspects of the bisection and secant methods with some additional features that make it completely robust and usually very efficient. Brent’s root-location method is a clever algorithm that does just that by applying a speedy open method wherever possible, but reverting to a reliable bracketing method if necessary. The approach was developed by Richard Brent (1973) based on an earlier algorithm of Theodorus Dekker (1969).
The bracketing technique is the trusty bisection method, whereas two different open methods are employed. The first is the secant method and the second is inverse quadratic interpolation.
Brent’s method combines bisection and quadratic interpolation into an efficient root-finding algorithm. In most problems the method is much faster than bisection alone, but it can become sluggish if the function is not smooth. It is the recommended method of root finding if the derivative of the function is difficult or impossible to compute.
Figure 2: Inverse quadratic iteration
Brent’s method assumes that a root of f (x) = 0 has been initially bracketed in the interval (x1, x2). The root-finding process starts with a bisection step that halves the interval to either (x1, x3) or (x3, x2), where x3 = (x1 + x2)/2, as shown in Fig. 2. In the course of bisection we had to compute f1 = f (x1), f2 = f (x2) and f3 = f (x3), so that we now know three points on the f (x) curve (the open circles in the figure). These points allow us to carry out the next iteration of the root by inverse quadratic interpolation (viewing x as a quadratic function of f ). If the result x of the interpolation falls inside the latest bracket (as is the case in Fig. 2), we accept the result. Otherwise, another round of bisection is applied.
The next step is to relabel x as x3 and rename the limits of the new interval x1 and x2 (x1 < x3 < x2), as indicated in Fig. 3. We have now recovered the original sequencing of points in Fig. 2, but the interval (x1, x2) containing the root has been reduced. This completes the first iteration cycle. In the next cycle another inverse quadratic interpolation is attempted and the process is repeated until the convergence criterion |x − x3| < ε is satisfied, where ε is a prescribed error tolerance.
Figure 3: Relabeling points after an iteration
The inverse quadratic interpolation is carried out with Lagrange’s three-point interpolant. Interchanging the roles of x and f , we have
Setting f = 0 and simplifying, we obtain for the estimate of the root
The change in the root is
Problem4: Solve problem1 using Brent’s method.
Calling the Brent’s function from command window,
>> y = @(m) sqrt(9.81*m/0.25)*tanh(sqrt(9.81*0.25/m)*4)-36;
M-file to implement Brent’s method.
function root = brent(func,a,b,tol) %% % Finds a root of f(x) = 0 by combining quadratic % interpolation with bisection (Brent’s method). %% input: % func = name of function % a,b = limits of the interval containing the root % tol = error tolerance (default is 1.0e6*eps) %% output: % root = real root (root = NaN if failed to converge). %% if nargin < 4; tol = 1.0e6*eps; end % First step is bisection x1 = a; f1 = feval(func,x1); if f1 == 0; root = x1; return; end x2 = b; f2 = feval(func,x2); if f2 == 0; root = x2; return; end if f1*f2 > 0.0 error('Root is not bracketed in (a,b)') end x3 = 0.5*(a + b); % Beginning of iterative loop. for i = 1:30 f3 = feval(func,x3); if abs(f3) < tol root = x3; return end % Tighten brackets (a,b) on the root. if f1*f3 < 0.0; b = x3; else; a = x3; end if (b - a) < tol*max(abs(b),1.0) root = 0.5*(a + b); return end % Try quadratic interpolation. denom = (f2 - f1)*(f3 - f1)*(f2 - f3); numer = x3*(f1 - f2)*(f2 - f3 + f1)... + f2*x1*(f2 - f3) + f1*x2*(f3 - f1); % If division by zero, push x out of bracket % to force bisection. if denom == 0; dx = b - a; else; dx = f3*numer/denom; end x = x3 + dx; % If interpolation goes out of bracket, use bisection. if (b - x)*(x - a) < 0.0 dx = 0.5*(b - a); x = a + dx; end % Let x3 <-- x & choose new x1, x2 so that x1 < x3 < x2. if x < x3 x2 = x3; f2 = f3; else x1 = x3; f1 = f3; end x3 = x; end root = NaN;