16. Numerical integration
[status: enough written to run a course, but incomplete]
16.1. The integral
The integral of a function is the “area under the curve” of that function between two given points a and b. We write it in this way:
A picture that illustrates this is in the Wikipedia article on integrals which I reprodue here in Figure 16.1.1.
Figure 16.1.1 The integral is the area under the curve of a function, between two points a and b.
Other written tutorials explain integrals quite well, but the figure above allows an instructor to explain it at a blackboard or whiteboard.
For now remember a couple of things: the integral is an area, and notice the beautiful notation in (16.1.1).
Integrals come up a lot in science. In Physics, to example, work is defined as the integral of force over a distance. This way of describing of work then leads to defining a body’s energy. Electric flux is the integral of the electric field over a surface (a 2-dimensional integral, but the same idea). A simple application in pure math is to calculate the areas of certain figures and the volumes of some solids.
Integrals are related to derivatives by the fundamental theorem of calculus, developed in the 17th century and based on work by Archimedes in ancient Greece. I will not touch on this too much at this time: it’s possible to explain derivatives and integrals to kids, but explaining the “antiderivative” might be out of our reach. We will just occasionally mention little bits of information about the link.
16.2. Calculating the integral numerically
It rare that we can calculate an integral exactly, but we can write computer programs that approximate the integral of any function. In Figure 16.1.1 we see one way of doing it by breaking the area down into many small rectangles. The area is then close to the sum of the areas of those small rectangles.
Figure 16.2.1 Approximating the area under the curve by making many small rectangles that try to fit the curve.
Figure 16.1.1 shows how we can break that region of space into a bunch of rectangles. The area of all those rectangles will approximate the integral of our function. If you insert more rectangles, they will snuggle up to the curve even better.
The drawing in Figure 16.2.1 lets us explain the notation in Equation (16.1.1). Each little rectangle in the figure has a base of \(dx\), and a height of \(f(x)\). The area of that little rectangle is height \(\times\) base, which is \(f(x) dx\). Now think of \(\int\) as a fancily drawn letter S and you see that the notation in Equation (16.1.1) reads as “the sum of \(f(x) dx\) for x between a and b”, which means “the sum of the areas of all the baby rectangles between a and b”.
The key insight is then to notice that in the limit (yes! it’s limits again!) of an increasing number of smaller and smaller rectangles we get the the integral.
Let us write a program which calculates the area under a curve using this approximation. Students are encouraged to think for a while before looking at Listing 16.2.1. I spend some time devising an algorithm with the students. I typically go to the whiteboard and, interactively with the students, I point out that we need:
A loop over the number of rectangles.
A variable with the area we have accumulated so far. (This is a standard paradigm in many programs: a variable that starts at zero, and you add stuff to it as you compute it.)
Identifying the lower left and right corners of each baby rectangle.
Finding the upper left and right corners of each baby rectangle.
Finding the height of the rectangle using \(f(x)\).
Taking that individual rectangle’s area.
Adding that rectangle’s area to the sum variable.
Now let us write the program in Listing 16.2.1 program and run it. It will report the area it found. If you put few rectangles you will get a worse approximation, and if you put in more rectangles you get a better approximation.
#! /usr/bin/env python3
import math
def main():
a = -1
b = 1
n_rectangles = 100
# A = integrate_function(curvy_cubic, a, b, n_rectangles)
A = integrate_function(upper_semicircle, a, b, n_rectangles)
print('with %d rectangles, the area is %g' % (n_rectangles, A))
def integrate_function(f, a, b, n_rectangles):
width = (b - a) / n_rectangles
total_area = 0
for i in range(n_rectangles):
midpoint = a + i * width + width/2
height = f(midpoint)
area = width * height
total_area += area
return total_area
def curvy_cubic(x):
return (x-1)**3 - 3*(x-1) + math.sin(x)
# Note that this function, when integrated from -1 to 1, should give
# an area of pi/2
def upper_semicircle(x):
return math.sqrt(1 - x*x)
main()
You can modify integrate-with-rectangles.py to use a curvy
function curvy_cubic() or to use the upper_semicircle()
function. The latter is interesting because the integral under the
upper semicircle should be half the area of a circle with that radius:
\(\pi r^2\). In our case the radius is 1, so we expect that
\(\int_{-1}^{1} {\rm upper\_semicircle}(x) dx = \pi/2\). This
means that we can use our integration program to calculate
\(\pi\)!
But integrate-with-rectangles.py was a rather dry program: it is
simple to write and understand, it gives the answer, but it does not
carry us on a wave of enthusiasm.
So let us see if we can make an animation of the process of this approximation. Enter the program in Listing 16.2.2. This version shows the plot of the curve as well as the subdivision of the area into rectangles, and it gives us an update on the estimate of the area.
#! /usr/bin/env python3
import math
import matplotlib.pyplot as plt
import numpy as np
def main():
a = -1.3
b = 3.1
func = curvy_cubic
# a = -1
# b = 1
# func = upper_semicircle
fig = plt.figure()
for n_rectangles in range(5, 200, 3):
do_single_integration(fig, func, a, b, n_rectangles)
plt.waitforbuttonpress()
def do_single_integration(fig, func, a, b, n_rectangles):
x = np.arange(a, b, 1.0/n_rectangles)
fig.clear()
plt.grid(True)
plt.plot(x, func(x), color='g')
A = integrate_function(func, a, b, n_rectangles)
info_para = """func: %s
n_rectangles: %d
area: %g
""" % (func.__name__, n_rectangles, A)
plt.annotate(
info_para,
xy=(0, -1), xytext=(0, -2))
print(info_para)
fig.canvas.draw_idle()
plt.pause(0.8)
def integrate_function(f, a, b, n_rectangles):
width = (b - a) / n_rectangles
total_area = 0
for i in range(n_rectangles):
left_point = a + i*width
right_point = a + (i+1)*width
midpoint = a + i * width + width/2
height = f(midpoint)
area = width * height
total_area += area
# now plot that rectangle; we must renormalize all to be in
# the [0,1] range for both x and y
if f(midpoint) < 0:
color = 'r'
else:
color = 'b'
norm_left = left_point
plt.bar([midpoint], height, width=width, edgecolor=color, fill=False)
return total_area
def curvy_cubic(x):
return (x-1)**3 - 3*(x-1) + np.sin(x)
# Note that this function, when integrated from -1 to 1, should give
# an area of pi/2
def upper_semicircle(x):
return np.sqrt(1 - x*x)
main()
Run the program and pay close attention to the ongoing calculation of the area and see if it converges as you get more and more rectangles.
16.3. Improving the numerical approximation
If you imagine using trapezoids instead of rectangles in Figure 16.2.1 you will notice that they hug the curve more closely. The area of a trapezoid is straightforward to calculate (I usually get the students to work it out while I write their ideas on the whiteboard). This method is called the trapezoidal rule.
Modify the program in Listing 16.2.2 to use trapezoids instead of rectangles. Compare the results to those from the rectangular method.
Even better is a method called Simpson’s rule. We take three points along the line and approximate the curve between them with a second degree polynomial. The curve then becomes a stitched together sequence of baby parabolas. Parabolas hug our function curve even more closely, so we should get a better approximation.
The procedure for Simpson’s rule has three main steps:
Break the interval from a to b into segments which have a left, middle and right point.
Find the parabolas using techniques similar to those in Section 7.7.1.
We know how to integrate polynomials exactly, so each tiny parabola has a known area, for example \(\int_a^b x^2 dx = (b^3-a^3)/3\).
Modify the program in Listing 16.2.2 to use Simpson’s rule.