Rewriting the Newton Method example to use functions
The last example from the previous lecture is an obvious candidate to introduce functions. The example makes use of two functions, the polynomial f(x) whose root we are trying to find, and its derivative.
Here is that example rewritten to make use of those two functions.
# Program to find the roots of f(x) = x^4+2 x^3-14 x^2+2 x+1
# by Newton's method
# Define the function and its derivative
def f(x):
return (((x+2)*x-14)*x+2)*x + 1
def fp(x):
return ((4*x+6)*x-28)*x+2
x = input("Enter a starting guess for the root: ")
x = float(x)
y = f(x)
tolerance = 10e-6
while y < -tolerance or y > tolerance:
# Compute f'(x)
yp = fp(x)
# Do one round of Newton's method
x = x - y/yp
y = f(x)
print("The root estimate is " + str(x))
Constructing a useful function: computing exp(x)
Most programming systems include a library of useful functions, and Python is no exception. For example, to make use of the standard mathematical functions such as exp(x) or sin(x), we can import the math module.
>>> import math >>> print(str(math.exp(3.3))) 27.112638920657883
As an exercise in function writing, let us see if we can construct our own function to compute exp(x). The key to doing this is the power series expression for exp(x):
Here is a first draft of a function that can compute this power series.
def exp1(x):
"""A first attempt to compute exp(x)."""
sum = 1
for n in range(1,16):
sum += x**n/math.factorial(n)
return sum
Since the sum in the power series is an infinite sum, we will have to settle for a finite approximation. I determined the cutoff of n = 15 by trial and error, since this produced the best possible results with the fewest terms in the sum. Here is some output from a simple test program that confirms that this works well enough.
math.exp(0.3)= 1.3498588075760032
exp1(0.3)= 1.349858807576003
math.exp(0.5)= 2.0137527074704766
exp1(0.5)= 2.0137527074704766
Next, we introduce a simple optimization. Since x**n and math.factorial(n) are both expensive to compute, we take advantage of a simple relationship between successive terms in the sum.
This means that we can easily compute each new term in the sum from the term that came before it.
Here is a version of the exp function that uses this optimization.
def exp2(x):
"""A more refined version of the exponential function."""
sum = 1 + x
term = x
for n in range(2,16):
term = term*x/n
sum += term
return sum
Here is a quick test run to confirm that this of course produces the same numerical results as the earlier, less efficient version.
math.exp(0.3)= 1.3498588075760032
exp1(0.3)= 1.349858807576003
exp2(0.3)= 1.349858807576003
math.exp(0.5)= 2.0137527074704766
exp1(0.5)= 2.0137527074704766
exp2(0.5)= 2.0137527074704766
Our exponential computing function seems to be working OK, but it still has one last, serious weakness. The problem with power series approximations is that they don't converge equally well for all values of x. Here is another test run that shows that when we plug in a larger x the approximation gets much worse.
math.exp(3.3)= 27.112638920657883
exp2(3.3)= 27.11262722548732
We could work around this problem by increasing the number of terms in the sum as the input x gets larger. Instead, I am going to pursue an alternative strategy. Since exp2 seems to do a pretty good job with inputs x in the range from 0 to 1, the trick is to map any x we are given into that range. We can do that by the following strategy:
- Compute the smallest power n such that x/2n < 1.
- Compute exp(x/2n).
- Note that
- We can compute the latter expression by repeated squaring:
exp(x) = (exp(x/2n))2n
(exp(x/2n))2n = ((exp(x/2n))2)2n-1
Here now is the final version of the exponential computing function that makes use of this strategy.
def exp3(x):
"""Final version of exp - works for a broader range of x."""
# Initially we work with the absolute value of x
z = math.fabs(x)
# Compute the smallest n such that z/2**n <= 1
n = 1
while z > 1:
n += 1
z /= 2.0
# Since z is now in the range 0 <= z <= 1, exp2 will do a good job with it
result = exp2(z)
# Since now result = exp(z/2**n), to compute exp(z) we now need to compute
# result**(2**n). We do this by repeated squaring of result, since
# (result**2)**(2**(n-1)) = result**(2**n)
while n > 1:
result = result*result
n -= 1
# Finally, we handle the case of a negative x.
if(x < 0):
return 1/result
else:
return result
and a test run to confirm that this works much better than exp2:
math.exp(3.3)= 27.112638920657883
exp2(3.3)= 27.11262722548732
exp3(3.3)= 27.112638920657805
Creating and using a module
As programs get larger and more complex we will want to start organizing useful functions into modules that we can import into our programs. In the next example I am going to write a set of useful functions for working with prime numbers and then import those functions into a simple test program.
First, here is the module with the useful function definitions:
def isPrime(n,known_primes):
"""
Test n to determine whether or not it is a prime number.
The second parameter must be a list of prime numbers that
includes all of the prime numbers smaller than the square
root of n.
"""
for p in known_primes:
if p*p > n:
return True
if n % p == 0:
return False
return True
def primeFactors(n,primes):
"""Compute and return a list of the prime factors of n."""
factors = []
x = n
for p in primes:
while x % p == 0:
factors.append(p)
x = x//p
if x > 1:
factors.append(x)
return factors
def findPrimes(m,n,primes):
"""
Construct and return a list of all the prime numbers less than n.
primes must be a seed list of prime numbers less than or equal to m.
"""
myPrimes = primes[:]
for num in range(m,n):
if isPrime(num,myPrimes) and not num in myPrimes:
myPrimes.append(num)
return myPrimes
Here is the test program that uses the functions in this module.
from prime_stuff import *
p = [2,3,5,7]
allPrimes = findPrimes(10,100,p)
x = input("Enter an integer between 2 and 10000: ")
if isPrime(int(x),allPrimes):
print(x + " is a prime number.")
else:
factors = primeFactors(int(x),allPrimes)
print(x + " is not a prime number.")
print("Its factors are " + str(factors))
Extended example - numerical integration
Integration as a geometry problem
At the end of calculus I you were introduced to the definite integral of a function f(x).
This quantity can be defined geometrically as an area: the area between the function f(x) and the x-axis over the range of x values a ≤ x ≤ b. As you saw back in calculus 1, you can use geometric techniques to estimate this area. The simplest geometric technique involves dividing the interval from a to b into a mesh of N + 1 equally spaced points xi = a + h i, where the step size h is equal to (b-a)/N.
This mesh of points divides the interval [a,b] into subintervals [xi , xi+1] for i = 0 to N. Over each of those subintervals we can place a rectangle that approximates the area underneath f(x) over the subinterval [xi , xi+1]. One way to get a halfway decent approximating rectangle is to give the rectangle that sits over the interval [xi , xi+1] a height f(xi). The drawing below illustrates this set-up.
By summing the areas of all of these approximating rectangles we can get a crude approximation for the area we want.
A program to compute this area
Writing a Python program to estimate areas using the method I described above is relatively easy. Below you will see the code for a simple program that uses this method to estimate the area
This is a nice example to work with, because the area underneath the curve in this example is a quarter circle with radius 2, which has an area of π. Since we can predict what the exact area should be, we can construct an error estimate each time we do the rectangle sum. The program recomputes the sum of the rectangle areas for a variety of different values for N, so we can see how quickly the error shrinks as we increase N.
import math
# This is the function we will work with for this example
def f(x):
return math.sqrt(4.0 - (2-x)*(2-x))
# Compute an approximation to the integral
# of f(x) from a to b by using rectangles
# with steps of size (b-a)/N
def rectangleArea(a,b,N):
h = (b - a)/N;
sum = 0.0;
for k in range(N):
term = h*(f(a + k*h))
sum += term
return sum
print(" N estimate error");
N = 10
power = 1
while power <= 6:
estimate = rectangleArea(0.0,2.0,N)
error = math.fabs(math.pi - estimate)
print("{:>10d} {:f} {:.16f}".format(N,estimate,error))
N = N*10
power += 1
Here are the results that this program produces.
N estimate error
10 2.904518 0.2370743273414746
100 3.120417 0.0211756218107477
1000 3.139555 0.0020371866787654
10000 3.141391 0.0002011759784653
100000 3.141573 0.0000200371878285
1000000 3.141591 0.0000020011759911
As you can clearly see, the size of the error decreases by a factor of 10 each time we increase the number of subdivisions by 10. A somewhat disappointing aspect of these results is that the error does not decrease quickly enough. To get the error to decrease to, say, 10-12 we would need something like 1012 subintervals. The time needed to compute the total area of all of those subintervals would be unreasonably large.
A more sophisticated method
To improve our error estimate we need to do a better job of matching our estimated areas to the actual area under the curve. One way to do this is to replace rectangles with trapezoids. Over each subinterval [xi , xi+1] we place a trapezoid whose left height is given by f(xi) and whose right height is given by f(xi+1). The picture below illustrates that the trapezoid does a better job of matching the area under the curve than the rectangle did.
The area of one such trapezoid is h (f(xi) + f(xi+1))/2.
It is a relatively straightforward matter to rewrite the area computation program shown above to use this formula for the area estimate. Here are the results that produces.
N estimate error
10 3.104518 0.0370743273414744
100 3.140417 0.0011756218107473
1000 3.141555 0.0000371866787638
10000 3.141591 0.0000011759784582
100000 3.141593 0.0000000371877436
1000000 3.141593 0.0000000011760504
This is much better, but not yet perfect. A Python float can represent numbers with a precision of 16 digits, but we only managed to get our result accurate to about 10 digits. Squeezing out those last digits of the result is going to be harder to do. One way to proceed would be to keep increasing the number of subdivisions. This is not practical, because computing the last row in the table took more than a second, and it looks like we would need to use 100 times as many subdivisions to get a result accurate to 16 digits.
The Richardson extrapolation
The situation above illustrates a classic situation in computational science. A relatively crude method can often be significantly improved by making a few simple adjustments. However, to reach a high standard of performance one has to deploy more exotic, specialized techniques. The rest of these notes are going to take you through the development of one such technique, called the Richardson extrapolation. Through the use of this technique we will be able to get improved performance with less computation.
The Richardson extrapolation begins with a consideration of the error term in the trapezoid rule. When we deploy the trapezoid rule to estimate an integral with N subintervals and a step size of h, the trapezoid rule typically produces results that are accurate to a factor that is proportional to h2.
The notation O(h2) I used here means that the error is of order h2. Another way to write this fact is to write the error term as a combination of various powers of h:
(This is somewhat beyond the scope of the discussion here, but it does turn out that the error term from the trapezoid rule contains only even powers of h.)
The next step is deploy the Richardson extrapolation, which is a strange, but effective algebraic trick. In this trick we write down the error relationship twice, once for N trapezoids as we saw above, and a second time for 2 N trapezoids. Increasing the number of trapezoids by a factor of 2 reduces the step size from h to h/2, so we get these two expressions for the error:
Multiplying the second equation by 4 and subtracting the first equation from it produces
We can rewrite this as
or
By using this simple algebraic trick, we have switched from a method that has an error term with magnitude O(h2) to a method with an error term with magnitude O(h4).
Better yet, this trick can be iterated. I will spare you the details, but iterating this trick is the basis for an integration technique called Romberg integration. Here is a summary of the Romberg technique. The function R(k,j) represents the jth iteration of this technique starting from a trapezoid rule with N = 2k-1 steps.
R(a,b,k,1) = trapezoidArea(a,b,2k-1)
Here is a Python program that uses this recursive formula to compute integral estimates.
import math
# This is the function we will work with for this example
def f(x):
return math.sqrt(4.0 - (2-x)*(2-x))
# Compute an approximation to the integral
# of f(x) from a to b by using trapezoids
# with steps of size (b-a)/N
def trapezoidArea(a,b,N):
h = (b - a)/N;
sum = 0.0;
for k in range(N):
term = h*(f(a + k*h)+f(a + (k+1)*h))/2
sum += term
return sum
# Recursively defined Romberg formula
def romberg(a,b,k,j):
if j == 1:
return trapezoidArea(a,b,2**(k-1))
elif k < j:
return 0
else:
return romberg(a,b,k,j-1) + (romberg(a,b,k,j-1)-romberg(a,b,k-1,j-1))/(4**(j-1)-1)
print(" j estimate error");
for j in range(2,6):
estimate = romberg(0.0,2.0,4*j,j)
error = math.fabs(math.pi - estimate)
print("{:>3d} {:f} {:.16f}".format(j,estimate,error))
It takes a little bit of effort to unpack what this program is doing. The loop at the bottom is computing R(a,b,4 j,j) for j in the range from 2 to 5. For a given value of j, this method starts with a trapezoid rule estimate using 2k-1 = 24 j - 1 steps. When j = 5 this means using 220 - 1 = 524288 steps. This looks like a lot, but remember that in the last example program we pushed the number of steps up to a million.
Here are the results this program returns
j estimate error 2 3.141275 0.0003172346542017 3 3.141588 0.0000043511837551 4 3.141593 0.0000000660130737 5 3.141593 0.0000000010239907
This program has one last flaw in it, and this flaw seriously bogs down the computation starting at around j= 4. The problem here is one of excessive recursion. As you can see from the code above, a typical call to the romberg function triggers three additional recursive calls. If each of those in turn triggers more calls to the romberg function, you can end up making a huge number of calls this function.
To counteract this problem, we can use a simple caching strategy. Each time we compute a value for romberg(a,b,k,j) for a particular combination of k and j, we should make a note of that result. The next time we end up needing that value we can simply look up the result without having to recompute anything.
The final version of the program makes use of this caching strategy. We set up a two dimensional list R that can store values for many combinations of k and j. The romberg function then hands most of the work off to a recursive internal function, r_int. r_int does an optimized calculation by first checking the cache to see if a result has already been recorded for this combination of k and j. If so, it immediately stops and returns the cached value. If not, it falls back to the recursive formula to compute a result, and then caches that result before returning an answer.
import math
# This is the function we will work with for this example
def f(x):
return math.sqrt(4.0 - (2-x)*(2-x))
# Compute an approximation to the integral
# of f(x) from a to b by using trapezoids
# with steps of size (b-a)/N
def trapezoidArea(a,b,N):
h = (b - a)/N;
sum = 0.0;
for k in range(N):
term = h*(f(a + k*h)+f(a + (k+1)*h))/2
sum += term
return sum
# Recursively defined Romberg formula with caching
# r_int will compute R[k,j] but then remember that
# result in the R list to speed up subsequent calculations.
def r_int(a,b,k,j,R):
if R[k][j] != -1: # Have we seen this k,j before?
return R[k][j] # If so, just return the cached value
# If not, compute, cache, and return the result
if j == 1:
R[k][j] = trapezoidArea(a,b,2**(k-1))
elif k < j:
R[k][j] = 0
else:
R[k][j] = r_int(a,b,k,j-1,R) + (r_int(a,b,k,j-1,R)-r_int(a,b,k-1,j-1,R))/(4**(j-1)-1)
return R[k][j]
def romberg(a,b,k,j):
R = [[-1 for c in range(j+1)] for r in range(k+1)]
return r_int(a,b,k,j,R)
print(" j estimate error");
for j in range(3,22,2):
estimate = romberg(0.0,2.0,j,j)
error = math.fabs(math.pi - estimate)
print("{:>3d} {:f} {:.16f}".format(j,estimate,error))
Note that since we have a more efficient method now we can safely expand the range of j that we compute with. Here are the results that this program returns.
j estimate error 3 3.090764 0.0508290045413791 5 3.135506 0.0060864702273435 7 3.140835 0.0007579750725255 9 3.141498 0.0000946982204093 11 3.141581 0.0000118360657910 13 3.141591 0.0000014794727754 15 3.141592 0.0000001849330067 17 3.141593 0.0000000231165633 19 3.141593 0.0000000028895282 21 3.141593 0.0000000003611502