Programming Praxis

One Response to “Root Finding”

1. Paul said

   March 17, 2014 at 9:59 AM

   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

Enter your comment here...

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

Email (required) (Address never made public)

Name (required)

Website

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

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

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

Cancel

Connecting to %s

Notify me of new comments via email.

Notify me of new posts via email.

Δ