Root Finding
March 14, 2014
Finding the roots of an equation is an important operation in mathematics, and there are many algorithms for finding roots numerically. We will study two of those algorithms in today’s exercise, and perhaps we will look at other algorithms in future exercises.
Let’s begin with the rules of the game. We are given a univariate function f(x) and want to find one or more values of x such that f(x) = 0. The function could be a polynomial, or could use trigonometric or exponential functions, or integrals, or any other mathematical operation.
The bisection algorithm starts with two points lo and hi that bound the root; thus, one of f(lo) and f(hi) is positive and the other is negative. The algorithm is iterative; at each step, the midpoint mid = (lo + hi) / 2 is found, the function f(mid) is evaluated at the midpoint, and then it replaces either lo or hi, whichever has the same sign. Iteration stops if f(mid) = 0 or if the difference between lo and hi is sufficiently small.
The regula falsi method is similar, but instead of calculating the center midpoint it calculates the midpoint at the point where a line connecting the current lo and hi crosses the axis. The method dates to the ancient Egyptians and Babylonians.
Both methods work only when the function f is continuous, with no discontinuities between lo and hi.
Your task is to write functions that find roots by the bisection and regula falsi methods. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.
Programming Praxis 
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