# integration.py Phys 305 David MacQuigg 1/28/10 ''' Illustrating the use of Python to teach math and science. The main purpose of Python is clarity. Later we will re-write these functions in C for 100X boost in speed. For now, we are just trying to understand some algorithms used in numerical integration, including their accuracy and relative speed. There is a lot of "overhead" in most languages if you want nice printouts, timings, etc. Python minimizes that overhead. ~''' def f(x): return 0.5 + x*x # This is the function we will integrate. ''' A function can be as simple as the one-liner above. By the way, we will be mixing text and executable code freely in this tutorial. The triple quotes make the text a "multi-line string". At the beginning of a module or function, these are the "docstrings" we will use later in automated documentation and testing. Elsewhere, they are just ignored by the Python iterpreter. You should be reading this file in IDLE (Python's Integrated Development Environment. Many text editors won't display all the colors that make readability so nice. ~''' def integLeft(a, b, f, nbins=10): '''Return the integral from a to b of function f using the left hand rule >>> integLeft(0, 1, f, 4) 0.71875 ''' h = float(b-a)/nbins # allow input of integer a and b assert h > 0 # warn if b not > a assert type(nbins) == int sum = 0.0 for n in range(nbins): # [0, 1, ... nbins-1] sum = sum + h * f(a + n*h) # left hand rule return sum ''' The first thing you may notice about Python, is that there are no brackets and semicolons or other punctuation to define statements and code blocks. Thus "program structure" is established by indentation. While this may seem like a bad idea to anyone who has not tried Python, it actually is better for teaching (and for experienced programmers who might make simple blunders). The interpreter catches indentation errors, and won't let the student continue with a subtle eror because a curly bracket is out of place. Teachers no longer have to nag students about proper style. Notes on the code above: 1) Python has a very flexible argument-passing syntax. If you include a "default value" with a parameter (nbins=10), that makes it a "keyword" parameter. Keyword parameters are optional in a call to the function. 2) Python variables have no type. Parameters a, b and f above are expected to be two floats and a function. While type declarations would catch some errors, Python expects the programmer write explicit and complete checks in those situations where defensive programming is important (like user input). The checks above do a lot more than just type checking. They also serve as documentation. Students will learn early about good programming practices, if teachers include some bad inputs in their automated test cases. On the other hand, if the emphasis is on learning algorithms, not programming, these tests can be left out. Python allows that flexibility. ~''' def integMid(a, b, f, nbins=10): '''Return the integral from a to b of function f using the midpoint rule >>> integMid(0, 1, f, 4) 0.828125 ''' h = float(b-a)/nbins assert h > 0 assert type(nbins) == int sum = 0.0 x = a + h/2 # first midpoint while (x < b): sum += h * f(x) x += h return sum def integTrap(a, b, f, nbins=10): '''Return the integral from a to b of function f using trapezoidal rule >>> integTrap(0, 1, f, 4) 0.84375 ''' h = float(b-a)/nbins assert h > 0 assert type(nbins) == int sum = (h/2) * (f(a) + f(b)) # endpoints are special for n in range(1, nbins): # [1, 2, ... nbins-1] sum += h * f(a + n*h) return sum if __name__ == '__main__': ''' Everything from here down runs only if this "integration.py" module is run directly from the interpreter. Here is the place for all the messy "overhead" needed to run tests and generate nicely-formatted output. The mess is ignored if this module is imported into another, for example, to make this module part of a library you might call to do integrations. ''' def table_of_errors(a, b, f, exact): ''' Print timings and a table of errors comparing these integration methods to an exact result. Show how the error decreases as o(1/nbins) or o(1/nbins**2), depending on the method. ~''' print print "nbins Left Mid Trap" for nbins in [4, 40, 400]: print '%4s %13.9f %13.9f %13.9f' % (nbins, integLeft(a,b,f,nbins), integMid(a,b,f,nbins), integTrap(a,b,f,nbins)) print for nbins in [4, 40, 400]: print '%4s %13.9f %13.9f %13.9f' % (nbins, integLeft(a,b,f,nbins) - exact, integMid(a,b,f,nbins) - exact, integTrap(a,b,f,nbins) - exact ) print from time import time # from time module import time function t0 = time() # seconds integLeft (0, 1, f, 1000) t1 = time() integMid (0, 1, f, 1000) t2 = time() integTrap (0, 1, f, 1000) t3 = time() print ("Time for 1000 bins:\n %13.6f %13.6f %13.6f" % ( (t1-t0), (t2-t1), (t3-t2) ) ) table_of_errors(0, 1, f, 5/6.) ''' >>> table_of_errors(0, 1, f, 5/6.) nbins Left Mid Trap 4 0.718750000 0.828125000 0.843750000 40 0.820937500 0.833281250 0.833437500 400 0.832084375 0.833332812 0.833334375 4 -0.114583333 -0.005208333 0.010416667 40 -0.012395833 -0.000052083 0.000104167 400 -0.001248958 -0.000000521 0.000001042 Time for 1000 bins: 0.000865 0.000763 0.000814 ''' from doctest import testmod testmod(verbose=True) # run the doctests every time this module is run. ''' Doctests serve two purposes: 1) They explain by example how your functions are called. 2) They serve as "unit tests". If you ever make a change that breaks an existing test "testmod" will holler at you :>) Don't let subtle errors go by. If a program is to fail, let it fail immediately and loudly. Notes: "testmod" requires the output to be exactly as in the example. We had to move the table_of_errors test outside of the function (where testmod doesn't see it) to allow for run-to-run variations in the timings. ~'''