Quadrature

What is it?

Quadrature is a form of numerical integration, where the area under a curve is approximated by forming quadrilaterals using the function values at certain points (in the following, all points are evenly spaced or dictated by interpolation points).

Quadrature, combined with polynomial interpolation, allow us to quickly estimate the integral of complex functions.  While cos(\pi x) + x is straightforward to integrate, we will use it as our example function in the following methods.

What are some methods we use?

Midpoint

The midpoint rule approximates the area under the curve f on the interval [x_0, x_1] as a rectangle of length x_1-x_0 and height f(\frac{x_0+x_1}{2}).  In the image below, the midpoint rule is applied to cos(\pi x) + x over the interval [0, 1], approximated using 4 squares.

The midpoint rule can be described as follows, courtesy Wikipedia:

Trapezoidal

The trapezoidal rule creates a quadrilateral whose points on the interval [x_0, x_1] are: (x_0, 0), (x_0, f(x_0)), (x_1, f(x_1)), (x_1, 0).  In other words, it takes the function’s values at the endpoints (whereas the midpoint takes the function’s value only at the midpoint of the interval).  The picture below again uses the function cos(\pi x) +x over the interval [0, 1], approximated using the trapezoidal rule.

And again courtesy Wikipedia, the trapezoidal rule:

Simpson’s Rule

Simpson’s rule can be derived in multiple ways.  One such was is as a composite of the midpoint and trapezoidal rules.  Simpson’s rule is exactly \frac{2}{3}Midpoint Rule + \frac{1}{3} Trapezoidal Rule.  The below picture shows Simpson’s Rule using 4 points to approximate the integral of the function cos(\pi x) + x.

Once more, Wikipedia provides the description of Simpson’s Rule:

Code

NOTE: the most recent version of JUtil (available here) is only required to generate the sample output.  It is not required to use the following methods.  A version of the Quadrature.py that requires only numpy can be found here.

This python code has three methods:

midpoint(f, x0, xn, n)

trapezoid(f, x0, xn, n)

simpson(f, x0, xn, n)

For each method, the following input are taken:

f is a function that takes one number and returns one number (the function value at that point). This function necessarily returns values for input between x0 and xn, at equally spaced points x_0, x_1, ..., xn where x_{i}-x_{i-1} = \frac{x_n-x_0}{n} for i \in [1, n].

x0, xn are the left and right endpoints, respectively, of the interval over which we will approximate the integral of the function f.

n is the number of sections to divide the interval [x0,xn] into to approximate the integral of f on [x0,xn].

Each function returns three vectors: X,Y,VLines

X is the set of x coordinates of the quadrilaterals estimating the integral of f.  For an estimate using 3 rectangles, there would therefore be 3*4 = 12 values in X.

Y is the set of y coordinates of the quadrilaterals estimating the integral of f.  For an estimate using 3 rectangles, there would therefore be 3*4 = 12 values in Y.

(Note: in Python, simply do zip(X,Y) to return a vector of coordinates in the form (x,y))

VLines is a vector of x-coordinates representing the evenly-spaced sampling over the interval [x0,xn] which can be used to easily create vertical lines (for plotting).

Sample Code Output:

The code linked above contains a long section starting with:

if “__main__” in __name__:

The code in that section generates a plot which overlays the Midpoint, Trapezoidal, and Simpson’s Rule approximations of the function cos(\pi x) + x using 4 quadrilaterals.  It then plots dashed vertical lines at each left and right endpoint of the quadrilaterals to show each section.  Finally, the actual function is plotted in black.

The following image is produced:

Where the red line is the midpoint rule, green line is the trapezoidal rule, and blue line is Simpson’s rule.

The following image is a closer look at a similar function, to show the differences between methods.  The coloring remains unchanged.

Leave a comment