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 region 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 x1-x0, 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 initial guess until a desired accuracy was achieved.
A function that performs integration using this method would accept a function, the starting x0 and the ending xt for the integral. The function would break the interval between x0 and xt into n parts, when n is an initial guess. The function 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 values 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)
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.