## 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.

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
```