Root Finding

March 14, 2014

We begin with the bisection method; the midpoint is just the average of lo and hi, and iteration stops when the difference between lo and hi is less than epsilon:

(define (bisection f lo hi epsilon)
  (let loop ((lo lo) (f-lo (f lo)) (hi hi) (f-hi (f hi)))
    (let* ((mid (/ (+ lo hi) 2)) (f-mid (f mid)))
      (cond ((or (zero? f-mid) (< (abs (- hi lo)) epsilon)) mid)
            ((= (signum f-lo) (signum f-mid))
              (loop mid f-mid hi f-hi))
            (else (loop lo f-lo mid f-mid))))))

There are two changes for the regula falsi method. First, the midpoint is linear interpolation between lo and hi, instead of just their simple average. Second, since one of the endpoints never moves during regula falsi, the stopping criterion occurs when the difference between the new midpoint and one of the old endpoints is less than epsilon:

(define (regula-falsi f lo hi epsilon)
  (let loop ((lo lo) (f-lo (f lo)) (hi hi) (f-hi (f hi)))
    (let* ((mid (/ (- (* f-lo hi) (* f-hi lo)) (- f-lo f-hi))) (f-mid (f mid)))
      (cond ((or (zero? f-mid) (< (abs (- mid lo)) epsilon) (< (abs (- mid hi)) epsilon)) mid)
            ((= (signum f-lo) (signum f-mid))
              (loop mid f-mid hi f-hi))
            (else (loop lo f-lo mid f-mid))))))

For example, we can find the root between 1 and 10 of the function sin(x) * x3 in any of four ways:

> (bisection f 1. 10. 1e-12)
3.141592653589697
> (bisection f 10. 1. 1e-12)
3.141592653589697
> (regula-falsi f 1. 10. 1e-12)
3.1415926535883956
> (regula-falsi f 10. 1. 1e-12)
3.1415926535883956

The signum function is +1 if its argument is positive, −1 if its argument is negative, and 0 if its argument is zero; it’s implementation is given at http://programmingpraxis.codepad.org/jtw5uNay, where you can run the program.

About these ads

Pages: 1 2

One Response to “Root Finding”

  1. Paul said

    I have used these methods quite a lot in my work. Good methods are IMO bisection and secant (only if close to the root). Regula Falsi is not so good as bisection, as it is in general not faster and is not so robust as bisection. Try Regula Falsi on the function called dirty and you have to wait a long time (with this code it fails, as too many iterations are taken). The root finding functions return the solution and the number of function evaluations used. Combining the bisection and secant methods it is possible to make a method, that works fast for all functions, that do not underflow.

    def bisection(fun, x1, x2, args=(), xtol=1e-12, maxiter=1000):
        """ very robust - always works - not the fastest (linear)"""
        a, b   = float(x1), float(x2)
        fa, y = fun(a, *args), fun(b, *args)
        ncalls = 2
        if fa * y > 0.0:
            msg = "No bracket: f({:f})={:f}, f({:f})={:f}".format(a, fa, b, y)
            raise ValueError(msg)
        for i in xrange(maxiter):
            x = (a + b) / 2.0
            y = fun(x, *args)
            ncalls += 1
            if abs(a - b) < xtol or y == 0.0:
                return x, ncalls
            a, b = (x, b) if y * fa > 0 else (a, x)
        raise ValueError("too many iterations")
    
    def regfalsi(fun, x1, x2, args=(), xtol=1e-12, maxiter=1000):
        """ not so robust - also not fast (linear)"""
        a, b   = float(x1), float(x2)
        fa, fb = fun(a, *args), fun(b, *args)
        ncalls = 2
        if fa * fb > 0.0:
            msg = "No bracket: f({:f})={:f}, f({:f})={:f}".format(a, fa, b, fb)
            raise ValueError(msg)
        for i in xrange(maxiter):
            dx = -(a - b) / (fa - fb) * fa
            x = a + dx
            y = fun(x, *args)
            ncalls += 1
            if abs(dx) < xtol or y == 0.0:
                return x, ncalls
            a, fa, b, fb = (x, y, b, fb) if y * fa > 0 else (a, fa, x, y)
        raise ValueError("too many iterations")
    
    def secant(fun, x1, x2, args=(), xtol=1e-12, maxiter=1000):
        """the fastest (superlinear), but will often fail - not a bracketing method"""
        a, b   = float(x1), float(x2)
        fa, fb = fun(a, *args), fun(b, *args)
        ncalls = 2
        for i in xrange(maxiter):
            dx = -(a - b) / (1 - fb / fa)
            x = a + dx
            y = fun(x, *args)
            ncalls += 1
            if abs(dx) < xtol or y == 0.0:
                return x, ncalls
            a, fa, b, fb = b, fb, x, y
        raise ValueError("too many iterations")
    
    def f(x):
        return x ** 2 - 9
        
    from math import sin
    
    def g(x):
        return sin(x) * x ** 3
        
    def dirty(x, eps=1e-10):
        return ((x * 1000.0) ** 2 if x < 0 else 0) - eps
        
    fun, a, b = g, 1, 4
    print bisection(fun, a, b), secant(fun, a, b), regfalsi(fun, a, b) 
    #(3.1415926535897825, 45) (1.898432588672972e-12, 139) (3.141592653589149, 48)
    fun, a, b = f, 0, 10
    print bisection(fun, a, b), secant(fun, a, b), regfalsi(fun, a, b) 
    #(2.9999999999998295, 47) (3.0, 12) (2.9999999999992526, 50)
    fun, a, b = dirty, -0.001, 1
    print bisection(fun, a, b), secant(fun, a, b), regfalsi(fun, a, b) 
    #(-9.999935628022207e-09, 43) FAILS FAILS    
    

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 621 other followers

%d bloggers like this: