# 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.
~'''