Programming Praxis

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