[Edu-sig] Intro to calculus through Python

[Kirby Urner]

We learn in intro the calculus that the derivative of a single-variable
function at some point, is the ratio of change in f(x) to change in x,
where this change approaches zero as its limiting value:

    D(f(x)) = lim  (f(x+h)-f(x))/h
              h->0

What's maybe not so clear to some students, because the book isn't clear,
is that D(f) -- derivative of f -- is itself a function, with a function
as input.  The variable x stands for "any x" in the domain, not a specific
value.  It's a parameter.

Given Python's top-level functions, and the ability to take functions as
arguments, and return them as results, we have the option to illustrate the
above concept with some definite value of h.  This is discrete mathematics,
using floating point numbers, so not exactly formal calculus -- but such
approximations serve a valid pedagogical purpose.

def D(f,h=1e-3):
     """
     Return derivative of function f (also a function)
     """
     def df(x):
         deriv=(f(x+h)-f(x))/h
         return round(deriv,3)
     return df

Let's pause to admire Python's flexibility here.  D( ) is taking an
arbitrary function f as its input (with a default value of h), and
returning df.  Python's scoping rules (as of 1.6) will build df
using whatever function we pass in (f).

For example, let's define f(x) = x**2 and pass that to D( ) above:

   >>> def g(x): return x*x

   >>> g          # g is a function
   <function g at 0x00B50030>

   >>> dg = D(g)  # we pass g, a function, to another function

   >>> dg         # ... and we get back a function
   <function df at 0x0147C4B0>

Now compare g(x), for x = 0,1,2,3...10 and dg(x) for the same values:

   >>> [g(x) for x in range(11)]
   [0, 1, 4, 9, 16, 25, 36, 49, 64, 81, 100]
   >>> [dg(x) for x in range(11)]
   [0.001, 2.0009999999999999, 4.0010000000000003, 6.0010000000000003,
   8.0009999999999994, 10.000999999999999, 12.000999999999999,
   14.000999999999999, 16.001000000000001, 18.001000000000001,
   20.001000000000001]

Again, we're dealing in approximations here, but clearly dg(x) closely
approximates 2*x, which is the actual derivative of x**2.  We can highlight
this in a couple of ways.  By further rounding:

   >>> [round(dg(x)) for x in range(11)]
   [0.0, 2.0, 4.0, 6.0, 8.0, 10.0, 12.0, 14.0, 16.0, 18.0, 20.0]

...or by redefining dg using a smaller value for h (h yet closer to zero):

   >>> dg = D(g,h=1e-8)
   >>> [dg(x) for x in range(11)]
   [0.0, 2.0, 4.0, 6.0, 8.0, 10.0, 12.0, 14.0, 16.0, 18.0, 20.0]

How about integration?  The definite integral is the sum of the little 
rectangles
fitted to the curve, as their widths become vanishingly small.  The rectangles
stand next to each other starting from a, running up to b.  And because f(x+h)
and f(x) will generally differ, because f(x) is a curve of some kind, it makes
sense to compute the area of the rectangles using the *average* of these two
values.  Again, an approximation.  And again, we use functions for inputs and
outputs:

def I(f,h=1e-3):
     """
     Return definite integral function
     """
     def intf(b,a=0):
         sum=0
         x = a
         while x<=b:                    # keep accumulating
             sum += h*(f(x+h)+f(x))/2.0 # area of rectangle
             x += h
         return round(sum,3)
     return intf

Our definite integral function is inherently less precise, and has to do more
work, than our derivative function.  Use of the average is precise only insofar
as our straight line segments from (x,f(x)) to (x+h, f(x+h)) are good 
approximations
of the curve.  Clearly this gets better as h gets smaller, but since we have to
sum a lot of these rectangles, we can't afford to make h too small -- we're
depending on successive additions of h to get us all the way from a to b, 
after
all.

But let's admire Python some more.  The fundamental theorem of the calculus
states that the definite integral and derivative are inverse functions, i.e.
if we take the derivative of the integral of f, or the integral of the
derivative of f, over the interval [a,b], we should get back our original
function: f.  And Python allows us to show this rather directly (in an
approximate form):

   >>> inv1 = I(D(g))    # a function of a function of a function!
   >>> [inv1(x) for x in range(11)]
   [0.0, 1.0009999999999999, 4.0060000000000002, 9.0090000000000003,
   16.012, 25.004999999999999, 36.006, 49.006999999999998, 64.007999999999996,
   81.009, 100.03]

   >>> inv2 = D(I(g))     # reversing I and D from the above case
   >>> [inv2(x) for x in range(11)]
   [0.0, 1.0, 4.0, 9.0, 16.0, 25.0, 36.0, 49.0, 64.0, 81.0, 100.0]

In both of the above examples, we get back a function reasonably close
to g(x) = x**2, as demonstrated by feeding it the values 0...10.

inv1 was defined by applying the integral function to the function returned
by the derivative function, which got function g as its argument.

inv2 was defined by applying the derivative function to the definite
integral of g.  The interval [a,b] is arbitrary, with a=0 by default.

Python's ability to support high level functional programming helps us
appreciate that differentiation and integration apply to functions to
create other functions, and that these resulting functions have an
inverse relationship.

Kirby