%pylab inline
Populating the interactive namespace from numpy and matplotlib
from __future__ import division
import ad
import math
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.
def bisection(fun, x1, x2, args=(), xtol=1e-12, aftol=1e-12, rftol=1e-6, maxiter=1000):
steps = []
fx1 = fun(x1, *args)
y = fun(x2, *args)
neval = 2
if fx1*y > 0.0:
msg = "No bracket: f({:f})={:f}, f({:f})={:f}".format(x1, fx1, x2, y)
raise ValueError(msg)
for i in xrange(maxiter):
x = (x1 + x2)/2.0
steps.append(x)
y = fun(x, *args)
neval += 1
if abs(x1 - x2) < xtol or allclose((y,), (0.0,), atol=aftol, rtol=rftol):
return x, neval, steps
else:
(x1, x2) = (x1, x) if fx1*y < 0 else (x, x2)
raise ValueError("Exceeded MAXIT")
def regfalsi(fun, x1, x2, args=(), xtol=1e-12, aftol=1e-12, rftol=1e-6, maxiter=1000):
steps = []
fx1 = fun(x1, *args)
y = fun(x2, *args)
neval = 2
if fx1*y > 0.0:
msg = "No bracket: f({:f})={:f}, f({:f})={:f}".format(x1, fx1, x2, y)
raise ValueError(msg)
for i in xrange(maxiter):
x = (x1 + x2)/2.0
steps.append(x)
y = fun(x, *args)
neval += 1
if abs(x1 - x2) < xtol or allclose((y,), (0.0,), atol=aftol, rtol=rftol):
return x, neval, steps
else:
(x1, x2) = (x1, x) if fx1*y < 0 else (x, x2)
raise ValueError("Exceeded MAXIT")
from ad import *
def newton(fun, x1, args=(), xtol=1e-12, aftol=1e-12, rftol=1e-6, maxiter=1000):
steps = []
fx1 = fun(x1, *args)
g, h = gh(fun)
neval = 1
for i in xrange(maxiter):
m = g(x1, *args)
x1 += -fx1/m[0]
steps.append(x1)
fx1 = fun(x1, *args)
neval += 1
if allclose(fx1, 0.0, atol=aftol, rtol=rftol):
return x1, neval, steps
raise ValueError("Exceeded MAXIT")
print "Solving an example function:\nf(x) = 1/20 * (x + 4)*(x + 2)*(x + 1)*(x - 1)*(x - 3) - 2"
f = lambda x: 1/20 * (x + 4)*(x + 2)*(x + 1)*(x - 1)*(x - 3) - 2
print "\n Solution iterations"
x1, x2 = 2, 4
x_b, neval_b, steps_b = bisection(f, x1, x2)
print x_b, neval_b
x_s, neval_s, steps_s = regfalsi(f, x1, x2)
print x_s, neval_s
x_n, neval_n, steps_n = newton(f, x1-(x1 - x2)/2)
print x_n, neval_n
Solving an example function: f(x) = 1/20 * (x + 4)*(x + 2)*(x + 1)*(x - 1)*(x - 3) - 2 Solution iterations 3.12497159176 44 3.12497159176 44 3.12497159176 5
print "Timing test:"
print " Bisection method: ",
%timeit x, neval, steps_b = bisection(f, x1, x2)
print "Regfalsi method: ",
%timeit x, neval, steps_s = regfalsi(f, x1, x2)
print "Newton's method: ",
%timeit x, neval, steps_n = newton(f, x1-(x1 - x2)/2)
Timing test: Bisection method: 1000 loops, best of 3: 1.65 ms per loop Regfalsi method: 100 loops, best of 3: 2.03 ms per loop Newton's method: 1000 loops, best of 3: 1.23 ms per loop
figure(figsize=(20,10))
f = np.vectorize(f)
xdd = arange(x1, x2, .0100)
subplot(311)
title("REGFALSI")
plot(xdd, np.zeros_like(xdd), 'g-')
plot(xdd, f(xdd), lw=1.5)
plot([x1, x2], f([x1, x2]), 'ko')
steps_s = array(steps_s)
step(steps_s, f(steps_s), 'yo')
plot([x_s], f([x_s]), 'ro')
xlim([x1, x2])
subplot(312)
title("BISECTION")
plot(xdd, f(xdd), lw=1.5)
plot(xdd, np.zeros_like(xdd), 'g-')
plot([x1, x2], f([x1, x2]), 'ko')
steps_b = array(steps_b)
step(steps_b, f(steps_b), 'yo')
plot([x_b], f([x_b]), 'ro')
xlim([x1, x2])
subplot(313)
title("AD-NEWTON")
plot(xdd, f(xdd), lw=1.5)
plot(xdd, np.zeros_like(xdd), 'g-')
plot([x1-(x1 - x2)/2], f(x1-(x1 - x2)/2), 'ko')
step(steps_n, f(steps_n), 'yo')
plot([x_n], f([x_n]), 'ro')
xlim([x1, x2])
show()
loglog(abs(f(steps_b) - f(x_b)), 'o')
plot(100/exp(.7*arange(1, 50)), '-.')
hold(True)
plot(abs(f(steps_s) - f(x_s)), 'o')
plot(abs(f(steps_n) - f(x_n)), 'o')
grid(True)