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:

Figure 1: The ln function

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:

kxkf(xk)mmf(xk)
01.01.011.0
11.10.9090909143.636363
21.20.8333333321.666667
31.30.7692307643.076923
41.40.7142857121.428571
51.50.6666666742.666667
61.60.62521.25
71.70.5882352942.352941
81.80.5555555621.111111
91.90.5263157842.105263
102.00.510.5

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)))

Listing 1: Simpson's rule in Python

The 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))

Listing 2: Trapezoidal rule

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)

Listing 3: Midpoint rule in Python

Of course, all of this is already implemented in scipy, surely much faster than this implementation, but I thought it was interesting to put together.

Add a comment:

Completely off-topic or spam comments will be removed at the discretion of the moderator.

You may preserve formatting (e.g. a code sample) by indenting with four spaces preceding the formatted line(s)

Name: Name is required
Email (will not be displayed publicly):
Comment:
Comment is required
My Book

I'm the author of the book "Implementing SSL/TLS Using Cryptography and PKI". Like the title says, this is a from-the-ground-up examination of the SSL protocol that provides security, integrity and privacy to most application-level internet protocols, most notably HTTP. I include the source code to a complete working SSL implementation, including the most popular cryptographic algorithms (DES, 3DES, RC4, AES, RSA, DSA, Diffie-Hellman, HMAC, MD5, SHA-1, SHA-256, and ECC), and show how they all fit together to provide transport-layer security.

My Picture

Joshua Davies

Past Posts