## Brent’s Root Finder

### March 28, 2014

Here is our implementation of Brent’s algorithm: d is the value of the last step and e is the value of the second-last step; if the last step was a bisection then e = d. The bolded tests are those that force the bisection every three steps:

```def brent(fun, x1, x2, args=(), xtol=1e-12, disp=False, maxiter=1000):     """     Bisection: k = log2((b-a)/tol) function eval.     Dekker : upper bound: 2**k function eval.     Brent : theoretical upper bound: (k+1)**2 ? 2                practical upper bound: 3*(k+1)     Example:         tol = 1e-6         b-a = 1         Then k = 20         Upperbounds             Bisection: 20 (upperbound is also number iterations)             Dekker : 2**20             Brent : 63 (practical bound)     """     if disp: print "BRENT"     a, b = float(x1), float(x2)     fa, fb = fun(a, *args), fun(b, *args)     ncalls = 2     if fa * fb > 0.0:         raise RootNotBracketedException(a, b, fa, fb)     c, fc = b, fb     for iters in range(maxiter):         if disp: msg = ""         if fc * fb > 0.0:             if disp: msg += "c=a "             c, fc = a, fa             e = d = b - a         if abs(fc) < abs(fb):             if disp: msg += "swap "             a, b, c, fa, fb, fc = b, c, b, fb, fc, fb         tol1 = 2 * EPS * abs(b) + xtol         xm = (c - b) / 2         if abs(xm) <= tol1 or fb == 0:             return b, iters, ncalls, 0         if abs(e) >= tol1 and abs(fa) > abs(fb):             s = fb / fa             if a == c: # if a==c we cannot do quadratic interpolation, so try linear                 p = 2 * xm * s                 q = 1 - s             else:                 # try quadratic interpolation                 q = fa / fc                 r = fb / fc                 p = s * (2 * xm * q * (q - r) - (b - a) * (r - 1))                 q = (q - 1) * (r - 1) * (s - 1)         if p > 0:             q = -q         p = abs(p)         if 2 * p < abs(e * q) and 2 * p < 3 * xm * q - abs(tol1 * q):             if disp: msg += ("LINEAR" if a == c else "QUAD ")             e = d             d = p / q         else:             e = d = xm             if disp: msg += "Bisect"     else:         e = d = xm         if disp:             msg += "Bisect"     a, fa = b, fb     if abs(d) > tol1:         b += d     else:         b += __sign(tol1, xm)     fb = fun(b, *args)     ncalls += 1     if disp:         print "%15s a %8.5g %s b %8.5g %s c %8.5g %s " "fb %12.5g xm %8.5g" % (msg, a, __sss(fa), b, __sss(fb), c, __sss(fc), fb, xm)     raise TooManyIterationsException(maxiter)```

You can run the program at http://ideone.com/QXD932.

Pages: 1 2