Root Finding, Again

March 21, 2014

We give our solution in Python, starting with the secant method. This method is not very useful in practice, because going outside the original bounds means that some steps actually retreat from the root instead of advancing toward it:

def secant(fun, x1, x2, args=(), xtol=1e-12, maxiter=1000):
    """the fastest (superlinear), but will often fail"""
    a, b = float(x1), float(x2)
    fa, fb = fun(a, *args), fun(b, *args)
    for i in xrange(maxiter):
        if fa == fb:
            raise ZeroDivisionError("fa == fb in secant")
        dx = -(a - b) / (1 - fb / fa)
        x = a + dx
        y = fun(x, *args)
        if abs(dx) < xtol or y == 0.0:
            return x
        a, fa, b, fb = b, fb, x, y
    raise ValueError("too many iterations")

This is essentially the same as the regula falsi method, but after every iteration the last two points are kept for the next iteration, instead of trying to keep a bracket on the root.

Now the Dekker method. Points x1 and x2 should bracket a root. If disp=True the iterates are printed and one can see if a bisection or a secant step is taking place.

def dekker(fun, x1, x2, args=(), xtol=1e-8, disp=False):
    """see wiki Brent
        c is last iterate of b
        a, b bracket solution abs(fb)
        if f(s) has same sign as fa, then old b beomes new contrapoint
        the solution is within 6*EPS + 2 * xtol
    """
    a, b = float(x1), float(x2)
    fa, fb = fun(a, *args), fun(b, *args)
    if fa * fb > 0.0: raise RootNotBracketedException(a, b, fa, fb)
    c, fc = a, fa
    feval = 2
    while 1:
        if abs(fa) < abs(fb):
            a, fa, b, fb, c, fc = b, fb, a, fa, b, fb
        xm = (a - b) / 2 # bisection step
        tol1 = 2 * EPS * abs(b) + xtol
        if abs(xm) 0:
            q = -q
        p = abs(p)
        d = p / q if p tol1 else __sign(tol1, xm))
        fb = fun(b, *args)
        feval += 1
        if disp:
            msg = "b" if abs(d) == abs(xm) else "s"
            fmt = "{:3d} a={:7.4g} b={:12.9g} c={:7.4g} fb={:10.3g} xm={:8.2g} {:s}"
            print fmt.format(feval, a, b, c, fb, xm, msg)
        if fa * fb > 0:     # if new iterate has same sign as a - old b take s place of a
            a, fa = c, fc

The __sign function returns x if y > 0, −x if y < 0, and 0 otherwise:

def __sign(x, y):
    if y > 0.0 :
        return x
    elif y < 0.0:
        return -x
    else:
        return 0.0

As a test, we compare the secant, regula falsi and dekker methods finding roots of the function cos(x) - x in the interval 0 to 1; the left column is the result and the right column is the time taken:

Secant       0.739085133215 0.000127772968474
Regula Falsi 0.739085133215 0.00897438331583
Dekker       0.739085133215 0.00184064914826

The secant method is very good when it works. Here are the intermediate steps of the Dekker algorithm:

3 a=      0 b= 0.685073357 c=      1 fb=    0.0893 xm=    -0.5 s
4 a=      1 b= 0.736298998 c= 0.6851 fb=   0.00466 xm=    0.16 s
5 a=      1 b= 0.739119362 c= 0.7363 fb= -5.73e-05 xm=    0.13 s
6 a= 0.7363 b= 0.739085112 c= 0.7391 fb=  3.53e-08 xm= -0.0014 s
7 a= 0.7391 b= 0.739085133 c= 0.7391 fb=  2.67e-13 xm= 1.7e-05 s
8 a= 0.7391 b= 0.739085133 c= 0.7391 fb=         0 xm= 1.7e-05 s

By comparison, here are the results when the three methods are run on a “dirty function” that has the value -1e-10 in (0,0.99) and then rises steeply to 1 at x=1:

Secant       fa == fb in secant
Regula Falsi too many iterations
Dekker       0.990000100000 0.000124694101764

The secant method fails because the two points have the same value, and regula falsi takes too long, but the Dekker method finds the root without difficulty.

You can run the program at http://programmingpraxis.codepad.org/Xe17Hcuf.

About these ads

Pages: 1 2

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 632 other followers

%d bloggers like this: