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.
