# 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(abfnbins=10):
    '''Return the integral from a to b of function f using the left hand rule

    >>> integLeft(0, 1, f, 4)
    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(abfnbins=10):
    '''Return the integral from a to b of function f using the midpoint rule

    >>> integMid(0, 1, f, 4)
    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(abfnbins=10):
    '''Return the integral from a to b of function f using trapezoidal rule

    >>> integTrap(0, 1, f, 4)
    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(1nbins):    # [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(abfexact):
    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 "nbins       Left         Mid           Trap"

        for nbins in [440400]:
            print '%4s %13.9f %13.9f %13.9f' % (nbins,
               integLeft(a,b,f,nbins), integMid(a,b,f,nbins), integTrap(a,b,f,nbins))
        for nbins in [440400]:
            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 )

        from time import time  # from time module import time function
        t0 = time()  # seconds
        integLeft (01f1000)
        t1 = time()
        integMid  (01f1000)
        t2 = time()
        integTrap (01f1000)
        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.)

    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.

"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.