Department of Mathematics

Math 300: Mathematical Computing

A Tale of Two Trapezoids

In our discussion of flow control we wrote a python function to compute a trapezoidal rule approximation to the integral function we specified over some interval we also specified. That approximation was given by ab fct x dx = h2 fcta+ fctb+ 2 i=1 n-1 fctxi = h2 -fcta -fctb+ 2 i=0 n fctxi , where a=x0 and b=xn. The function we used to evaluate this approximation was

def trapezoid(fct,a,b,n): if n<0 or n==0: print "Error: n must be positive" return False h = (double(b)-double(a))/double(n) total = 0.0 for i in arange(n+1): total += fct(a+i*h) total *= 2.0 total -= fct(a)+fct(b) return total*h/2.0

It seems interesting to find out how fast this function is. We can do that using the clock function from the time module. This function gives the current time to some level of precision. On our computer it gave the time to the nearest hundredth of a second. We should be able to find the time before we start our trapezoid function, then run the trapezoid function, then find the time again. The difference between the two times is roughly how long our function took to run. We can write a little function to do this for us.

>>> import myfunctions as mf >>> import time >>> def ttrap(fct): ... t1=time.clock() ... fct(mf.f,0,1,1000000) ... t2=time.clock()-t1 ... print t2 ... >>> ttrap(mf.trapezoid) 11.58

Our function is called ttrap, and it takes one argument: the name of the function it is going to time. It finds the current time and puts it in a variable called t1. Then it runs the function you supplied as an argument. After that it finds the time again, subtracts t1 to find the elapsed time, and prints that. We find that running our trapezoid function with a million subintervals (n=1000000) takes eleven and a half seconds.

You might wonder why we were so cagey about the name of the function we wanted to time - why not just type trapezoid instead of fct and then supplying trapezoid as an argument. You might also wonder whether we could write a faster function. That, of course is the point. We can write a function that runs faster, and we'll give it a different name so we can use this ttrap function to time it.

The thing that takes all the time in the trapezoid function is the "for" construction. A million times it has to evaluate f(x) for some scalar x. We could write this instead to do the calculations on a numpy array. Maybe that would go faster. Our new function looks like this.

def arraytrapezoid(fct,a,b,n): if n<0 or n==0: print "Error: n must be positive" return False h = (double(b)-double(a))/double(n) x = arange(a,b+h,h) y = 2.0*fct(x) y[0] = y[0]/2.0 y[n] = y[n]/2.0 return sum(y)*h/2.0

Everything is the same up through the definition of h. At that point, instead of starting a "for" construction, we define a vector x of points of the partition, in such a way that xi= a+ih for each i=0, 1, ..., n. After that, we have only to evaluate the function for each point of that array, which we do simply by feeding the entire array to fct as an argument. The properties of arrays take care of the elementwise multiplication that is required. After that we just subtract the extra function value at a and b, and return the result.

Evidently this version of the code needs a great deal more memory, since it has to store the entire vector x as well as the function evaluation over that: y. However, this version of the code is spectacularly faster than our earlier version. Remember that our other version took over eleven seconds to run.

>>> ttrap(mf.arraytrapezoid) 0.05

Why is this more than twenty times faster? It is the difference between interpreted code and compiled code. These are essentially the two ways to run computer programs. Interpreters look at the text commands you type, parse their meaning, translate them into machine code, and then execute them. Compilers translate an entire file of code into machine language, which then only needs to be loaded and run full speed. Interpreted code thus runs much more slowly than that compiled, as it were pretranslated code.

In our case, the trapezoid function relies on that "for" construction. Thus each time it executes the code, the python interpreter must translate that command to add to the total, then execute it. By contrast, in the arraytrapezoid function, everything happens through the numpy array construction. This, like many things in numpy, is compiled. Thus in this case, instead of one million interpreted lines, there is only one, and then the computation happens in compiled code, so it is much faster. In general, because the numpy package contains much compiled code, the more we rely on it and not native python, the faster our programs will run.

Assignment C is posted.

The test solution is available.

Department of Mathematics, PO Box 643113, Neill Hall 103, Washington State University, Pullman WA 99164-3113, 509-335-3926, Contact Us
Copyright © 1996-2015 Kevin Cooper