MTH 337: Week 8¶
%pylab inline
Finding Roots of Equations¶
We're solving the equation $\cos(x) = x$ numerically.¶
Rewrite as the function $f(x) = \cos(x) - x$. This then becomes the "root-finding problem" of finding a value $x$ for which $f(x) = 0$.
def myf(x):
return cos(x) - x
Bisection Method¶
- Start with a "bracket" [a, b] with $sign(f(a)) \neq sign(f(b))$.
- Evaluate $f$ at the midpoint c = (a + b)/2.
- "Replace" either a or b by c, to keep $sign(f(a)) \neq sign(f(b))$.
- Repeat while bracket width b - a > tolerance.
This is guaranteed to work, but is slow.
def mybisection(f, a, b, tol):
sign_a = sign(f(a))
sign_b = sign(f(b))
# First check that we don't already have a root
if sign_a == 0:
return a
elif sign_b == 0:
return b
# Then check that we do actually have a bracket.
elif sign_a == sign_b:
print "Not a bracket!"
return
# Keep track of the current step
step = 0
# Keep bisecting the interval [a, b]
while b - a > tol:
print step, a, b
step += 1
c = (a + b)/2.
sign_c = sign(f(c))
if sign_c == 0:
return c
elif sign_a == sign_c:
a = c
else:
b = c
return (a + b)/2.
mybisection(myf, 0, pi/2., 1e-5)
Functional Iteration¶
- Start with some guess $x_0$ at the solution.
- Generate a series of approximations $x_1, x_2, x_3, \ldots$ by $x_{i + 1} = f(x_i) + x_i$
def func_iter(f, x, tol):
step = 0
fval = f(x)
while abs(fval) > tol:
print step, x
fval = f(x)
x = fval + x
step += 1
return x
func_iter(myf, 1., 1e-5)
Newton's Method¶
- Assume that f is differentiable, and that values are available for both f(x) and f'(x).
- Find the tangent to f at the current estimate $x_i$.
- Let $x_{i + 1}$ be the intersection of the tangent with the x-axis.
Then: $x_{i + 1} = x_i - \dfrac{f(x_i)}{f'(x_i)}$
This is very fast when it works, and converges quickly. But, the method is a little "fragile", and needs a good initial guess.
# Return both the function and the derivative
def myfdf(x):
return cos(x) - x, -sin(x) - 1
def newton(fdf, x, tol):
while True:
f, df = fdf(x)
print x, f
s = - f/df
x += s
if abs(s) <= tol:
return x
newton(myfdf, .5, 1e-10)
Problems with Newton's Method¶
We use Newton's Method to find the root of $f(x) = tanh(x)$.
The graph of this function is shown below.
x = linspace(-2, 2)
plot(x, tanh(x))
axhline(color='k')
axvline(color='k')
Define a function which returns the values f(x), f'(x)¶
def tanhdf(x):
return tanh(x), 1./cosh(x)**2
Modify the "newton" function to stop after a maximum number of steps¶
def newton(fdf, x, tol, maxsteps=10):
for i in range(maxsteps):
f, df = fdf(x)
print x, f
s = - f/df
x += s
if abs(s) <= tol:
return x
Testing, we find that initial values |x| < 1.08 converge, |x| > 1.09 diverge¶
newton(tanhdf, 1.08, 1e-5)
newton(tanhdf, 1.09, 1e-5)
We can find the initial value of x that neither converges nor diverges by solving:¶
$x_{i + 1} = -x_i$
$\Rightarrow x_i - \tanh(x_i)\cosh^2(x_i) = -x_i$
$\Rightarrow \sinh(x_i)\cosh(x_i) - 2 x_i = 0$
This is itself a root-finding problem, with $f(x) = \sinh(x)\cosh(x) - 2 x$
So, Newton's method can be used to find the initial values of x for which Newton's method does not converge!
def myfdf(x):
return sinh(x)*cosh(x) - 2*x, cosh(2*x) - 2
newton(myfdf, 1, 1e-6)
Colors¶
The RGB model¶
In the RGB model, colors are represented in terms of their three components - red, green, and blue.
Each component is represented as a number between fully off and fully on.
These components can be integers from 0 to 255, or floats in the interval [0, 1]
Images are 3D arrays of size (rows, columns, 3)¶
- Image height is the number of rows.
- Image width is the number of columns
- The '3' indicates the three color components
The "dtype" for the array needs to be either a float in the interval [0, 1], or an 8-bit integer of type 'uint8'. We'll use floats for report 7.
A pure red square¶
a = zeros((10, 10, 3))
a[:,:,0] = 1
imshow(a)
A pure green square¶
a = zeros((10, 10, 3))
a[:,:,1] = 1
imshow(a)
A pure blue square¶
a = zeros((10, 10, 3))
a[:,:,2] = 1
imshow(a)
Slicing can be used to set the color components in different parts of the image¶
Modern art?¶
n = 100
a = zeros((n, n, 3))
for i in range(100):
w, x, y, z = random.randint(0, n, size=4)
dim = random.randint(0, 3)
a[min(w,x):max(w,x), min(y,z):max(y,z), dim] = random.rand()
imshow(a)
Boolean Indexing¶
An array can be indexed using another array.
If a boolean array is used as an index, only the values with a matching True in the boolean array are returned.
a = arange(5)
print a
b = array([True, True, False, False, True])
print b
print a[b]
Boolean indexing can be used to select just those array elements which pass some test¶
a = arange(5)
print a
print a > 2 # This is the test
print a[a > 2] # Only the elements passing the test are returned
Boolean indexing can also be used to modify just those array elements which pass some test¶
a = arange(5)
print a
b = a % 2 == 0 # This is the test, which is only passed by even numbers
print b
a[b] = 100 # Set all even numbers to 100
print a