Numerical Integration in Python
Recently, I've started looking into TensorFlow and what's being called "deep learning" (i.e. neural networks), and I've discovered as I try to read through the algorithms that my calculus has gotten a bit rusty. Of course, it's entirely possible to read through the descriptions of the algorithms involved without fully understanding the theory behind them, but I understood the background just well enough that I kept being irritated that I couldn't put the whole story together. I did take three semesters of undergraduate calculus, after all — I knew this stuff at one time. As it turns out, since I hadn't been able to sell my calculus textbook back to the campus bookstore because it was an old edition, I actually still have it. So I dug it out, dusted it off, and started reading over the chapters on partial differentiation and multiple integrals. I was pleasantly surprised how quickly and easily it came back to me. (Which bodes well, because I'm going to be helping my son with his calculus homework when he starts his junior year of high school in a few years).
Of course, being a math textbook, my calculus book is full of back-references to previous chapters. I came across one back reference that I remember despising as an undergraduate: Simpson's rule. See, if you can't find a closed- form integral for a particular function (or if you don't even actually have a continuous function to integrate), you can estimate the definite integral by evaluating the function at specific points. So, for example, to compute:
(to approximate the ln function, for example), you could evaluate:
This can be geometrically interpreted as the shaded area in the graph shown in figure 1, below:
More generally, the definite integral of any function can be estimated by:
For some function f and some n — the larger the better. There's a whole theory of why this works based on approximations of parabolas in the spaces between each 1/n increment of the function evaluation that you can look up if you're curious.
As I look at some of the example homework problems, I'm reminded of exactly why I hated this topic so much: estimating definite integrals by way of Simpson's rule is tedious with a hand calculator, which is what we had to work with back in the 80's. (Well, actually we did have computers back then, too, but we weren't allowed to use them in exams). I remember hand-tabulating long tables of values like this one:
I was struck immediately by how naturally this fits the functional programming paradigm, so I put together a quick Python program to work some of the examples:
from __future__ import division def simpson(a, b, n, f): sum = 0 inc = (b - a) / n for k in range(n + 1): x = a + (k * inc) summand = f(x) if (k != 0) and (k != n): summand *= (2 + (2 * (k % 2))) sum += summand return ((b - a) / (3 * n)) * sum # Examples of use print(simpson(1, 2, 10, lambda x: 1 / x)) print(simpson(1, 4, 6, lambda x: 1 / x)) import math print(simpson(0, 1, 6, lambda x: 1 / math.sqrt(1 + x * x)))
simpson function keeps track of the running sum and the
multiplicands associated with each evaluation of the function, but the
function to be evaluated is, itself, an actual python function that takes in
a single argument x and returns its value.
A related rule called the trapezoidal rule yields similar results:
from __future__ import division def trapezoid(a, b, n, f): sum = 0 inc = (b - a ) / n for k in range(n + 1): x = a + (k * inc) summand = f(x) if (k != 0) and (k != n): summand *= 2 sum += summand return sum * ((b - a) / (2 * n)) # Example of use print(trapezoid(1, 2, 10, lambda x: 1 / x))
As well as the "midpoint" rule that evaluates the function in between each point in the intervals:
from __future__ import division def midpoint(a, b, n, f): sum = 0 x_int = ((2 * n + 1) * a - b) / (2 * n) inc = (b - a) / n for k in range(1, n + 1): x = x_int + (k * inc) sum += f(x) return sum * inc # Example of use print midpoint(1,2,10,lambda x: 1/x)
Of course, all of this is already implemented in
much faster than this implementation, but I thought it was interesting to put