Integration in Python

An integral is most often thought of as the area under a curve, where the curve is a function (Figure 11.1a). Numerical integration amounts calculating that area using an algorithm. The area of a rectangle is easy to calculate, so if the re­gion under a curve could be reasonably approximated by a bunch of rectangles, then the problem would be solved. This is the idea behind the trapezoidal rule. The integral from x0 to x1 of a function f(x) can be approximated by the width (x1-x0) multiplied by the height (the average value of the function in that range), which is just a rectangle that approximates the area under the curve (Figure 11.1b). In mathematical notation, this is

This would generally be a pretty poor approximation of a curve, and would yield correspondingly bad approximations of the integral. However, the smaller the width x₁-x₀, the more accurate the approximation can be, and so using a great many small trapezoids would be much better than using only one (Figure 11.1c). How many? That is not known at the outset, but could be increased from an ini­tial guess until a desired accuracy was achieved.

A function that performs integration using this method would accept a func­tion, the starting x₀ and the ending x_t for the integral. The function would break the interval between x₀ and x_t into n parts, when n is an initial guess. The func­tion is evaluated for all n parts, the area of each trapezoid is computed, and they are summed to get the final result. Now increase n and do it again. If the two values are close enough (delta), then the process is complete.

This is done in two steps, first using a function trap0() that computes and returns the sum of N trapezoids. The obvious but slow way to do this is shown in the following code.

def trap0 (f, x0, x1, n):   # Slow method

The integration function trapezoid() calls this function with increasing val­ues of n until two consecutive results show a small enough difference (i.e., it is smaller than a provided delta value):

def trapezoid (f, x0, x1, delta=0.0001, niter=1024):

n = 4

resold = trap0(f, x0, x1, 2)

resnew = trap0(f, x0, x1, 4)

while abs(resnew-resold) > delta:

if n>niter: break

resold = resnew

n = n * 2

resnew = trap0 (f, x0, x1, n)

return resnew

The function trap0() can be sped up significantly by not re-computing the function twice each time through the loop, but remembering the previous value instead (Exercise 2). A more popular algorithm for integration is Simpson’s Rule, which tries to minimize the error even more by using a quadratic approximation to the curve at the top of the trapezoid, instead of a straight line.

Source: Parker James R. (2021), Python: An Introduction to Programming, Mercury Learning and Information; Second edition.